python3的sympy

print（“字符串”），5/2和5//2的结果是不同的5/2为2.5,5//2为2.

创新互联是少有的网站建设、做网站、营销型企业网站、小程序定制开发、手机APP,开发、制作、设计、卖友情链接、推广优化一站式服务网络公司,2013年至今,坚持透明化,价格低,无套路经营理念。让网页惊喜每一位访客多年来深受用户好评

python2需要导入from_future_import division执行普通的除法。

1/2和1//2的结果0.5和0.

%号为取模运算。

乘方运算为2**3，-2**3和-（2**3）是等价的。

from sympy import*导入库

x,y,z=symbols('x y z'),定义变量

init_printing(use_unicode=True)设置打印方式。

python的内部常量有pi，

函数simplify，simplify(sin(x)**2 + cos(x)**2)化简结果为1，

simplify((x**3 + x**2 - x - 1)/(x**2 + 2*x + 1))化简结果为x-1。化简伽马函数。simplify(gamma(x)/gamma(x - 2))得（x-2）(x-1)。

expand((x + 1)**2)展开多项式。

expand((x + 1)*(x - 2) - (x - 1)*x)

因式分解。factor(x**2*z + 4*x*y*z + 4*y**2*z)得到z*(x + 2*y)**2

from_future_import division

x,y,z,t=symbols('x y z t')定义变量，

k, m, n = symbols('k m n', integer=True)定义三个整数变量。

f, g, h = symbols('f g h', cls=Function)定义的类型为函数。

factor_list(x**2*z + 4*x*y*z + 4*y**2*z)得到一个列表，表示因式的幂，(1, [(z, 1), (x + 2*y, 2)])

expand((cos(x) + sin(x))**2)展开多项式。

expr = x*y + x - 3 + 2*x**2 - z*x**2 + x**3，collected_expr = collect(expr, x)将x合并。将x元素按阶次整合。

collected_expr.coeff(x, 2)直接取出变量collected_expr的x的二次幂的系数。

cancel()is more efficient thanfactor().

cancel((x**2 + 2*x + 1)/(x**2 + x))

，expr = (x*y**2 - 2*x*y*z + x*z**2 + y**2 - 2*y*z + z**2)/(x**2 - 1)，cancel(expr)

expr = (4*x**3 + 21*x**2 + 10*x + 12)/(x**4 + 5*x**3 + 5*x**2 + 4*x)，apart(expr)

asin(1)

trigsimp(sin(x)**2 + cos(x)**2)三角函数表达式化简，

trigsimp(sin(x)**4 - 2*cos(x)**2*sin(x)**2 + cos(x)**4)

trigsimp(sin(x)*tan(x)/sec(x))

trigsimp(cosh(x)**2 + sinh(x)**2)双曲函数。

三角函数展开，expand_trig(sin(x + y))，acos(x)，cos(acos(x))，expand_trig(tan(2*x))

x, y = symbols('x y', positive=True)正数，a, b = symbols('a b', real=True)实数，z, t, c = symbols('z t c')定义变量的方法。

sqrt(x) == x**Rational(1, 2)判断是否相等。

powsimp(x**a*x**b)幂函数的乘法，不同幂的乘法，必须先定义a和b。powsimp(x**a*y**a)相同幂的乘法。

powsimp(t**c*z**c)，注意，powsimp()refuses to do the simplification if it is not valid.

powsimp(t**c*z**c, force=True)这样的话就可以得到化简过的式子。声明强制进行化简。

(z*t)**2，sqrt(x*y)

第一个展开expand_power_exp(x**(a + b))，expand_power_base((x*y)**a)展开，

expand_power_base((z*t)**c, force=True)强制展开。

powdenest((x**a)**b)，powdenest((z**a)**b)，powdenest((z**a)**b, force=True)

ln(x)，x, y ,z= symbols('x y z', positive=True)，n = symbols('n', real=True)，

expand_log(log(x*y))展开为log(x) + log(y)，但是python3没有。这是因为需要将x定义为positive。这是必须的，否则不会被展开。expand_log(log(x/y)),expand_log(log(x**n))

As withpowsimp()andpowdenest(),expand_log()has aforceoption that can be used to ignore assumptions。

expand_log(log(z**2), force=True),强制展开。

logcombine(log(x) + log(y))，logcombine(n*log(x))，logcombine(n*log(z), force=True)。

factorial(n)阶乘，binomial(n, k)等于c（n，k），gamma(z)伽马函数。

hyper([1, 2], [3], z)，

tan(x).rewrite(sin)得到用正弦表示的正切。factorial(x).rewrite(gamma)用伽马函数重写阶乘。

expand_func(gamma(x + 3))得到，x*(x + 1)*(x + 2)*gamma(x)，

hyperexpand(hyper([1, 1], [2], z))，

combsimp(factorial(n)/factorial(n - 3))化简，combsimp(binomial(n+1, k+1)/binomial(n, k))化简。combsimp(gamma(x)*gamma(1 - x))

自定义函数

def list_to_frac(l):

expr = Integer(0)

for i in reversed(l[1:]):

expr += i

expr = 1/expr

return l[0] + expr

list_to_frac([x, y, z])结果为x + 1/z，这个结果是错误的。

syms = symbols('a0:5')，定义syms，得到的结果为(a0, a1, a2, a3, a4)。

这样也可以a0, a1, a2, a3, a4 = syms， 可能是我的操作错误 。发现python和自动缩进有关，所以一定看好自动缩进的距离。list_to_frac([1, 2, 3, 4])结果为43/30。

使用cancel可以将生成的分式化简，frac = cancel(frac)化简为一个分数线的分式。

(a0*a1*a2*a3*a4 + a0*a1*a2 + a0*a1*a4 + a0*a3*a4 + a0 + a2*a3*a4 + a2 + a4)/(a1*a2*a3*a4 + a1*a2 + a1*a4 + a3*a4 + 1)

a0, a1, a2, a3, a4 = syms定义a0到a4，frac = apart(frac, a0)可将a0提出来。frac=1/(frac-a0)将a0去掉取倒。frac = apart(frac, a1)提出a1。

help("modules"),模块的含义，help("modules yourstr")模块中包含的字符串的意思。，

help("topics")，import os.path + help("os.path")，help("list")，help("open")

# -*- coding: UTF-8 -*-声明之后就可以在ide中使用中文注释。

定义

l = list(symbols('a0:5'))定义列表得到[a0, a1, a2, a3, a4]

fromsympyimport*

x,y,z=symbols('x y z')

init_printing(use_unicode=True)

diff(cos(x),x)求导。diff(exp(x**2), x)，diff(x**4, x, x, x)和diff(x**4, x, 3)等价。

diff(expr, x, y, 2, z, 4)求出表达式的y的2阶，z的4阶，x的1阶导数。和diff(expr, x, y, y, z, 4)等价。expr.diff(x, y, y, z, 4)一步到位。deriv = Derivative(expr, x, y, y, z, 4)求偏导。但是不显示。之后用deriv.doit()即可显示

integrate(cos(x), x)积分。定积分integrate(exp(-x), (x, 0, oo))无穷大用2个oo表示。integrate(exp(-x**2-y**2),(x,-oo,oo),(y,-oo,oo))二重积分。print(expr)print的使用。

expr = Integral(log(x)**2, x)，expr.doit()积分得到x*log(x)**2 - 2*x*log(x) + 2*x。

integ.doit()和integ = Integral((x**4 + x**2*exp(x) - x**2 - 2*x*exp(x) - 2*x -

exp(x))*exp(x)/((x - 1)**2*(x + 1)**2*(exp(x) + 1)), x)连用。

limit(sin(x)/x,x,0)，not-a-number表示nan算不出来，limit(expr, x, oo)，，expr = Limit((cos(x) - 1)/x, x, 0)，expr.doit()连用。左右极限limit(1/x, x, 0, '+')，limit(1/x, x, 0, '-')。。

Series Expansion级数展开。expr = exp(sin(x))，expr.series(x, 0, 4)得到1 + x + x**2/2 + O(x**4)，，x*O(1)得到O(x)，，expr.series(x, 0, 4).removeO()将无穷小移除。exp(x-6).series(x,x0=6)，，得到

-5 + (x - 6)**2/2 + (x - 6)**3/6 + (x - 6)**4/24 + (x - 6)**5/120 + x + O((x - 6)**6, (x, 6))最高到5阶。

f=Function('f')定义函数变量和h=Symbol('h')和d2fdx2=f(x).diff(x,2)求2阶，，as_finite_diff(dfdx)函数和as_finite_diff(d2fdx2,[-3*h,-h,2*h])，，x_list=[-3,1,2]和y_list=symbols('a b c')和apply_finite_diff(1,x_list,y_list,0)。

Eq(x, y)，，solveset(Eq(x**2, 1), x)解出来x，当二式相等。和solveset(Eq(x**2 - 1, 0), x)等价。solveset(x**2 - 1, x)

solveset(x**2 - x, x)解，solveset(x - x, x, domain=S.Reals)解出来定义域。solveset(exp(x), x)    # No solution exists解出EmptySet()表示空集。

等式形式linsolve([x + y + z - 1, x + y + 2*z - 3 ], (x, y, z))和矩阵法linsolve(Matrix(([1, 1, 1, 1], [1, 1, 2, 3])), (x, y, z))得到{(-y - 1, y, 2)}

A*x = b 形式，M=Matrix(((1,1,1,1),(1,1,2,3)))，system=A,b=M[:,:-1],M[:,-1]，linsolve(system,x,y,z)，，solveset(x**3 - 6*x**2 + 9*x, x)解多项式。roots(x**3 - 6*x**2 + 9*x, x)，得出，{3: 2, 0: 1}，有2个3的重根，1个0根。solve([x*y - 1, x - 2], x, y)解出坐标。

f, g = symbols('f g', cls=Function)函数的定义，解微分方程diffeq = Eq(f(x).diff(x, x) - 2*f(x).diff(x) + f(x), sin(x))再和dsolve(diffeq,f(x))结合。得到Eq(f(x), (C1 + C2*x)*exp(x) + cos(x)/2)，dsolve(f(x).diff(x)*(1 - sin(f(x))), f(x))解出来Eq(f(x) + cos(f(x)), C1)，，

Matrix([[1,-1],[3,4],[0,2]])，，Matrix([1, 2, 3])列表示。M=Matrix([[1,2,3],[3,2,1]])

N=Matrix([0,1,1])

M*N符合矩阵的乘法。M.shape显示矩阵的行列数。

M.row(0)获取M的第0行。M.col(-1)获取倒数第一列。

M.col_del(0)删掉第1列。M.row_del(1)删除第二行，序列是从0开始的。M = M.row_insert(1, Matrix([[0, 4]]))插入第二行，，M = M.col_insert(0, Matrix([1, -2]))插入第一列。

M+N矩阵相加，M*N，3*M，M**2，M**-1，N**-1表示求逆。M.T求转置。

eye(3)单位。zeros(2, 3)，0矩阵，ones(3, 2)全1，diag(1, 2, 3)对角矩阵。diag(-1, ones(2, 2), Matrix([5, 7, 5]))生成Matrix([

[-1, 0, 0, 0],

[ 0, 1, 1, 0],

[ 0, 1, 1, 0],

[ 0, 0, 0, 5],

[ 0, 0, 0, 7],

[ 0, 0, 0, 5]])矩阵。

Matrix([[1, 0, 1], [2, -1, 3], [4, 3, 2]])

一行一行显示，，M.det()求行列式。M.rref()矩阵化简。得到结果为Matrix([

[1, 0,  1,  3],

[0, 1, 2/3, 1/3],

[0, 0,  0,  0]]), [0, 1])。

M = Matrix([[1, 2, 3, 0, 0], [4, 10, 0, 0, 1]])，M.nullspace()

Columnspace

M.columnspace()和M = Matrix([[1, 2, 3, 0, 0], [4, 10, 0, 0, 1]])

M = Matrix([[3, -2,  4, -2], [5,  3, -3, -2], [5, -2,  2, -2], [5, -2, -3,  3]])和M.eigenvals()得到{3: 1, -2: 1, 5: 2}，，This means thatMhas eigenvalues -2, 3, and 5, and that the eigenvalues -2 and 3 have algebraic multiplicity 1 and that the eigenvalue 5 has algebraic multiplicity 2.

P, D = M.diagonalize()，P得Matrix([

[0, 1, 1,  0],

[1, 1, 1, -1],

[1, 1, 1,  0],

[1, 1, 0,  1]])，，D为Matrix([

[-2, 0, 0, 0],

[ 0, 3, 0, 0],

[ 0, 0, 5, 0],

[ 0, 0, 0, 5]])

P*D*P**-1 == M返回为True。lamda = symbols('lamda')。

lamda = symbols('lamda')定义变量，p = M.charpoly(lamda)和factor(p)

expr = x**2 + x*y，srepr(expr)可以将表达式说明计算法则，"Add(Pow(Symbol('x'), Integer(2)), Mul(Symbol('x'), Symbol('y')))"。。

x = symbols('x')和x = Symbol('x')是一样的。srepr(x**2)得到"Pow(Symbol('x'), Integer(2))"。Pow(x, 2)和Mul(x, y)得到x**2。x*y

type(2)得到class 'int'，type(sympify(2))得到class 'sympy.core.numbers.Integer'..srepr(x*y)得到"Mul(Symbol('x'), Symbol('y'))"。。。

Add(Pow(x, 2), Mul(x, y))得到"Add(Mul(Integer(-1), Pow(Symbol('x'), Integer(2))), Mul(Rational(1, 2), sin(Mul(Symbol('x'), Symbol('y')))), Pow(Symbol('y'), Integer(-1)))"。。Pow函数为幂次。

expr = Add(x, x)，expr.func。。Integer(2).func，class 'sympy.core.numbers.Integer',,Integer(0).func和Integer(-1).func，，，expr = 3*y**2*x和expr.func得到class 'sympy.core.mul.Mul',,expr.args将表达式分解为得到(3, x, y**2)，，expr.func(*expr.args)合并。expr == expr.func(*expr.args)返回True。expr.args[2]得到y**2，expr.args[1]得到x，expr.args[0]得到3.。

expr.args[2].args得到(y, 2)。。y.args得到空括号。Integer(2).args得到空括号。

from sympy import *

E**(I*pi)+1，可以看出，I和E，pi已将在sympy内已定义。

x=Symbol('x')，，expand( E**(I*x) )不能展开，expand(exp(I*x),complex=True)可以展开，得到I*exp(-im(x))*sin(re(x)) + exp(-im(x))*cos(re(x))，，x=Symbol("x",real=True)将x定义为实数。再展开expand(exp(I*x),complex=True)得到。I*sin(x) + cos(x)。。

tmp = series(exp(I*x), x, 0, 10)和pprint(tmp)打印出来可读性好，print(tmp)可读性不好。。pprint将公式用更好看的格式打印出来，，pprint( series( cos(x), x, 0, 10) )

integrate(x*sin(x), x)，，定积分integrate(x*sin(x), (x, 0, 2*pi))。。

用双重积分求解球的体积。

x, y, r = symbols('x,y,r')和2 * integrate(sqrt(r*r-x**2), (x, -r, r))计算球的体积。计算不来，是因为sympy不知道r是大于0的。r = symbols('r', positive=True)这样定义r即可。circle_area=2*integrate(sqrt(r**2-x**2),(x,-r,r))得到。circle_area=circle_area.subs(r,sqrt(r**2-x**2))将r替换。

integrate(circle_area,(x,-r,r))再积分即可。

expression.sub([(x,y),(y,x)])又换到原来的状况了。

expression.subs(x, y)，，将算式中的x替换成y。。

expression.subs({x:y,u:v}) : 使用字典进行多次替换。。

expression.subs([(x,y),(u,v)]) : 使用列表进行多次替换。。

python 的scipy 里的 odeint 这个求微分方程的函数怎么用啊

scipy.integrate.odeint(func,y0,t,args=(),dfun=none,col_deriv=0,full_output=0,ml=none,mu=none,rtol=none,atol=none,tcrit=none,h0=0.0,hmax=0.0,hmin=0.0,ixpr=0,mxstep=0,mxhnil=0,mxordn=12,mxords=5,printmessg=0)

实际使用中，还是主要使用前三个参数，即微分方程的描写函数、初值和需要求解函数值对应的的时间点。接收数组形式。这个函数，要求微分方程必须化为标准形式，即dy/dt=f(y,t,)。

fromscipyimportodeint

y=odeint(dy/dt=r*y*(1-y/k),y(0)=0.1,t)

对于微分方程全还给老师了，

二级Python----Python的内置函数及标准库（DAY 8）

python的内置函数（68个）

Python考核31个内置函数，

python内置了很多内置函数、类方法属性及各种模块。当我们想要当我们想要了解某种类型有哪些属性方法以及每种方法该怎么使用时，我们可以使用dir()函数和help()函数在python idle交互式模式下获得我们想要的信息。

• dir()函数获得对象中可用属性的列表

Python中的关键词有哪些？

dir(__builtins__)：查看python内置函数

help(‘keywords‘)：查看python关键词

如微分积分方程的求解程序、访问互联网、获取日期和时间、机器学习算法等。这些程序往往被收入程序库中，构成程序库。

只有经过严格检验的程序才能放在程序库里。检验，就是对程序作充分的测试。通常进行的有正确性测试、精度测试、速度测试、边界条件和出错状态的测试。经过检验的程序不但能保证计算结果的正确性，而且对错误调用也能作出反应。程序库中的程序都是规范化的。所谓规范化有三重含义：①同一库里所有程序的格式是统一的；② 对这些程序的调用方法是相同的；③ 每个程序所需参数的数目、顺序和类型都是严格规定好的。

Python的库包含标准库和第三方库

标准库：程序语言自身拥有的库，可以直接使用。help('modules')

第三方库：第三方者使用该语言提供的程序库。

标准库： turtle 库（必选）、 random 库（必选）、 time 库（可选）。

• turtle 库：图形绘制库

原理如同控制一只海龟，以不同的方向和速度进行位移而得到其运动轨迹。

使用模块的帮助时，需要先将模块导入。

例如：在IDLE中输入import turtle

dir(turtle)

help(turtle.**)

1.画布

画布就是turtle为我们展开用于绘图区域, 我们可以设置它的大小和初始位置。

setup()方法用于初始化画布窗口大小和位置，参数包括画布窗口宽、画布窗口高、窗口在屏幕的水平起始位置和窗口在屏幕的垂直起始位置。

参数:width, height: 输入宽和高为整数时,表示 像素 ;为小数时,表示占据电脑屏幕的比例。(startx,starty):这一坐标表示

矩形窗口左上角顶点的位置,如果为空,则窗口位于屏幕中心:

例如：setup(640,480,300,300)表示在桌面屏幕(300,300)位置开始创建640×480大小的画布窗体。

2、画笔

• color() 用于设置或返回画笔颜色和填充颜色。

例如：color(‘red’)将颜色设为红色，也可用fillcolor()方法设置或返回填充颜色，或用pencolor()方法设置或返回笔触颜色。

如何用python求导数

打开python运行环境。

导入微分的模块包：from sympy import *。

定义符号变量：x = symbols('x')

定义一个函数：f = x**9

diff = diff(f,x)求导

最后输入diff，即可显示其变量值了。

众多python培训视频，尽在python学习网，欢迎在线学习！

如何使用python计算常微分方程？

常用形式

odeint(func, y0, t,args,Dfun)

一般这种形式就够用了。

下面是官方的例子，求解的是

D(D(y1))-t*y1=0

为了方便，采取D=d/dt。如果我们令初值

y1(0) = 1.0/3**(2.0/3.0)/gamma(2.0/3.0)

D(y1)(0) = -1.0/3**(1.0/3.0)/gamma(1.0/3.0)

这个微分方程的解y1=airy(t)。

令D(y1)=y0，就有这个常微分方程组。

D(y0)=t*y1

D(y1)=y0

Python求解该微分方程。

from scipy.integrate import odeint

from scipy.special import gamma, airy

y1_0 = 1.0/3**(2.0/3.0)/gamma(2.0/3.0)

y0_0 = -1.0/3**(1.0/3.0)/gamma(1.0/3.0)

y0 = [y0_0, y1_0]

def func(y, t):

... return [t*y[1],y[0]]

def gradient(y,t):

... return [[0,t],[1,0]]

x = arange(0,4.0, 0.01)

t = x

ychk = airy(x)[0]

y = odeint(func, y0, t)

y2 = odeint(func, y0, t, Dfun=gradient)

print ychk[:36:6]

[ 0.355028 0.339511 0.324068 0.308763 0.293658 0.278806]

print y[:36:6,1]

[ 0.355028 0.339511 0.324067 0.308763 0.293658 0.278806]

print y2[:36:6,1]

[ 0.355028 0.339511 0.324067 0.308763 0.293658 0.278806]

得到的解与精确值相比，误差相当小。

=======================================================================================================

args是额外的参数。

用法请参看下面的例子。这是一个洛仑兹曲线的求解，并且用matplotlib绘出空间曲线图。（来自《python科学计算》）

from scipy.integrate import odeint

import numpy as np

def lorenz(w, t, p, r, b):

# 给出位置矢量w，和三个参数p, r, b 计算出

# dx/dt, dy/dt, dz/dt 的值

x, y, z = w

# 直接与lorenz 的计算公式对应

return np.array([p*(y-x), x*(r-z)-y, x*y-b*z])

t = np.arange(0, 30, 0.01) # 创建时间点

# 调用ode 对lorenz 进行求解, 用两个不同的初始值

track1 = odeint(lorenz, (0.0, 1.00, 0.0), t, args=(10.0, 28.0, 3.0))

track2 = odeint(lorenz, (0.0, 1.01, 0.0), t, args=(10.0, 28.0, 3.0))

# 绘图

from mpl_toolkits.mplot3d import Axes3D

import matplotlib.pyplot as plt

fig = plt.figure()

ax = Axes3D(fig)

ax.plot(track1[:,0], track1[:,1], track1[:,2])

ax.plot(track2[:,0], track2[:,1], track2[:,2])

plt.show()

===========================================================================

scipy.integrate.odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0)

计算常微分方程（组）

使用 FORTRAN库odepack中的lsoda解常微分方程。这个函数一般求解初值问题。

参数：

func : callable(y, t0, ...) 计算y在t0 处的导数。

y0 : 数组 y的初值条件（可以是矢量）

t : 数组 为求出y，这是一个时间点的序列。初值点应该是这个序列的第一个元素。

args : 元组 func的额外参数

Dfun : callable(y, t0, ...) 函数的梯度(Jacobian)。即雅可比多项式。

col_deriv : boolean. True，Dfun定义列向导数（更快），否则Dfun会定义横排导数

full_output : boolean 可选输出，如果为True 则返回一个字典，作为第二输出。

printmessg : boolean 是否打印convergence 消息。

返回： y : array, shape (len(y0), len(t))

数组，包含y值，每一个对应于时间序列中的t。初值y0 在第一排。

infodict : 字典，只有full_output == True 时，才会返回。

字典包含额为的输出信息。

键值：

‘hu’ vector of step sizes successfully used for each time step.

‘tcur’ vector with the value of t reached for each time step. (will always be at least as large as the input times).

‘tolsf’ vector of tolerance scale factors, greater than 1.0, computed when a request for too much accuracy was detected.

‘tsw’ value of t at the time of the last method switch (given for each time step)

‘nst’ cumulative number of time steps

‘nfe’ cumulative number of function evaluations for each time step

‘nje’ cumulative number of jacobian evaluations for each time step

‘nqu’ a vector of method orders for each successful step.

‘imxer’index of the component of largest magnitude in the weighted local error vector (e / ewt) on an error return, -1 otherwise.

‘lenrw’ the length of the double work array required.

‘leniw’ the length of integer work array required.

‘mused’a vector of method indicators for each successful time step: 1: adams (nonstiff), 2: bdf (stiff)

其他参数，官方网站和文档都没有明确说明。相关的资料，暂时也找不到。

常微分方程的解析解(方法归纳)以及基于Python的微分方程数值解算例实现

本文归纳常见常微分方程的解析解解法以及基于Python的微分方程数值解算例实现。

考虑常微分方程的解析解法，我们一般可以将其归纳为如下几类：

这类微分方程可以变形成如下形式：

两边同时积分即可解出函数，难点主要在于不定积分，是最简单的微分方程。

某些方程看似不可分离变量，但是经过换元之后，其实还是可分离变量的，不要被这种方程迷惑。

形如

的方程叫做一阶线性微分方程，若 为0，则方程齐次，否则称为非齐次。

解法： (直接套公式)

伯努利方程

形如

的方程称为伯努利方程，这种方程可以通过以下步骤化为一阶线性微分方程：

令 , 方程两边同时乘以 ，得到

即 .

这就将伯努利方程归结为可以套公式的一阶线性微分方程。

形如

的方程称为二阶常系数微分方程，若 ，则方程称为齐次的，反之称为非齐次的。以下默认方程是非齐次的。

求解此类方程分两步：

原方程的解 = 齐次通解 + 非齐次特解

首先假设 .用特征方程法，写出对应的特征方程并且求解：

解的情况分为以下三种：

情况一：方程有两个不同的实数解

假设两个实数解分别是 , 此时方程的通解是

情况二：方程有一个二重解

假设该解等于 ，此时方程的通解是

情况三：方程有一对共轭复解

假设这对解是 , 此时方程的通解是

对于 和特征根的情况，对特解的情况做如下归纳：

形如

的方程叫做高阶常系数微分方程，若 ，则方程是齐次的，否则是非齐次的。下面默认方程是非齐次的。

求解此类方程分两步：

原方程的解 = 齐次通解 + 非齐次特解

考虑带有第三类边界条件的二阶常系数微分方程边值问题

问题一:两点边值问题的解析解

由于此方程是非齐次的,故 求解此类方程分两步：

原方程的解 = 齐次通解 + 非齐次特解

首先假设 . 用特征方程法，写出对应的特征方程

求解得到两个不同的实数特征根: .

此时方程的齐次通解是

由于 . 所以非齐次特解形式为

将上式代入控制方程有

求解得: , 即非齐次特解为 .

原方程的解 = 齐次通解 + 非齐次特解

于是,原方程的全解为

因为该问题给出的是第三类边界条件,故需要求解的导函数

且有

将以上各式代入边界条件

解此方程组可得: .

综上所述,原两点边值问题的解为

对一般的二阶微分方程边值问题

假定其解存在唯一,

为求解的近似值, 类似于前面的做法,

考虑带有第三类边界条件的二阶常系数微分方程边值问题

问题二:有限差分方法算出其数值解及误差

对于 原问题 ,取步长 h=0.2 ,用 有限差分 求其 近似解 ,并将结果与 精确解y(x)=-x-1 进行比较.

因为

先以将区间划分为5份为例，求出数值解

结果:

是不是解出数值解就完事了呢？当然不是。我们可以将问题的差分格式解与问题的真解进行比较，以得到解的可靠性。通过数学计算我们得到问题的真解为 ，现用范数来衡量误差的大小：

结果:

接下来绘图比较 时数值解与真解的差距：

结果:

将区间划分为 份, 即 时.

结果:

绘图比较 时数值解与真解的差距：

最后，我们还可以从数学的角度分析所采用的差分格式的一些性质。因为差分格式的误差为 ， 所以理论上来说网格每加密一倍，与真解的误差大致会缩小到原来的 . 下讨论网格加密时的变化：

结果：

本文标题：python求函数微分 如何用python求解微分
本文地址：http://myzitong.com/article/hhpgsg.html