锐客网

  1. 首页 > 睿知 > it技术 > >

python计算微分函数 python常微分方程( 二 )

IT技术计算


求解得到两个不同的实数特征根:.
【python计算微分函数 python常微分方程】 此时方程的齐次通解是
由于. 所以非齐次特解形式为
将上式代入控制方程有
求解得:, 即非齐次特解为.
原方程的解 = 齐次通解 + 非齐次特解
于是,原方程的全解为
因为该问题给出的是第三类边界条件,故需要求解的导函数
且有
将以上各式代入边界条件
解此方程组可得:.
综上所述,原两点边值问题的解为
对一般的二阶微分方程边值问题
假定其解存在唯一,
为求解的近似值, 类似于前面的做法,
考虑带有第三类边界条件的二阶常系数微分方程边值问题
问题二:有限差分方法算出其数值解及误差
对于 原问题 ,取步长 h=0.2 ,用 有限差分 求其 近似解 ,并将结果与 精确解y(x)=-x-1 进行比较.
因为
先以将区间划分为5份为例，求出数值解
结果:
是不是解出数值解就完事了呢？当然不是 。我们可以将问题的差分格式解与问题的真解进行比较，以得到解的可靠性 。通过数学计算我们得到问题的真解为 ， 现用范数来衡量误差的大?。?
结果:
接下来绘图比较时数值解与真解的差距：
结果:
将区间划分为份, 即时.
结果:
绘图比较时数值解与真解的差距：
最后，我们还可以从数学的角度分析所采用的差分格式的一些性质 。因为差分格式的误差为，所以理论上来说网格每加密一倍，与真解的误差大致会缩小到原来的. 下讨论网格加密时的变化：
结果：
如何使用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.3550280.3395110.3240680.3087630.2936580.278806]
print y[:36:6,1]
[ 0.3550280.3395110.3240670.3087630.2936580.278806]
print y2[:36:6,1]
[ 0.3550280.3395110.3240670.3087630.2936580.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))

推荐阅读


上一篇：oracle怎么改表名，oracle修改表的名字

下一篇：农村金融如何推广发展，如何促进农村金融健康发展