重庆分公司,新征程启航
为企业提供网站建设、域名注册、服务器等服务
有很多大学生问我,学习python有什么用呢?我说:你至少可以用来解微分方程,如下面的例子,就是解决微分方程:
大新ssl适用于网站、小程序/APP、API接口等需要进行数据传输应用场景,ssl证书未来市场广阔!成为成都创新互联公司的ssl证书销售渠道,可以享受市场价格4-6折优惠!如果有意向欢迎电话联系或者加微信:18980820575(备注:SSL证书合作)期待与您的合作!
y"+a*y'+b*y=0
代码如下:
[python] view plain copy
#y"+a*y'+b*y=0
from scipy.integrate import odeint
from pylab import *
def deriv(y,t): # 返回值是y和y的导数组成的数组
a = -2.0
b = -0.1
return array([ y[1], a*y[0]+b*y[1] ])
time = linspace(0.0,50.0,1000)
yinit = array([0.0005,0.2]) # 初值
y = odeint(deriv,yinit,time)
figure()
plot(time,y[:,0],label='y') #y[:,0]即返回值的第一列,是y的值。label是为了显示legend用的。
plot(time,y[:,1],label="y'") #y[:,1]即返回值的第二列,是y’的值
xlabel('t')
ylabel('y')
legend()
show()
输出结果如下:
本文归纳常见常微分方程的解析解解法以及基于Python的微分方程数值解算例实现。
考虑常微分方程的解析解法,我们一般可以将其归纳为如下几类:
这类微分方程可以变形成如下形式:
两边同时积分即可解出函数,难点主要在于不定积分,是最简单的微分方程。
某些方程看似不可分离变量,但是经过换元之后,其实还是可分离变量的,不要被这种方程迷惑。
形如
的方程叫做一阶线性微分方程,若 为0,则方程齐次,否则称为非齐次。
解法: (直接套公式)
伯努利方程
形如
的方程称为伯努利方程,这种方程可以通过以下步骤化为一阶线性微分方程:
令 , 方程两边同时乘以 ,得到
即 .
这就将伯努利方程归结为可以套公式的一阶线性微分方程。
形如
的方程称为二阶常系数微分方程,若 ,则方程称为齐次的,反之称为非齐次的。以下默认方程是非齐次的。
求解此类方程分两步:
原方程的解 = 齐次通解 + 非齐次特解
首先假设 .用特征方程法,写出对应的特征方程并且求解:
解的情况分为以下三种:
情况一:方程有两个不同的实数解
假设两个实数解分别是 , 此时方程的通解是
情况二:方程有一个二重解
假设该解等于 ,此时方程的通解是
情况三:方程有一对共轭复解
假设这对解是 , 此时方程的通解是
对于 和特征根的情况,对特解的情况做如下归纳:
形如
的方程叫做高阶常系数微分方程,若 ,则方程是齐次的,否则是非齐次的。下面默认方程是非齐次的。
求解此类方程分两步:
原方程的解 = 齐次通解 + 非齐次特解
考虑带有第三类边界条件的二阶常系数微分方程边值问题
问题一:两点边值问题的解析解
由于此方程是非齐次的,故 求解此类方程分两步:
原方程的解 = 齐次通解 + 非齐次特解
首先假设 . 用特征方程法,写出对应的特征方程
求解得到两个不同的实数特征根: .
此时方程的齐次通解是
由于 . 所以非齐次特解形式为
将上式代入控制方程有
求解得: , 即非齐次特解为 .
原方程的解 = 齐次通解 + 非齐次特解
于是,原方程的全解为
因为该问题给出的是第三类边界条件,故需要求解的导函数
且有
将以上各式代入边界条件
解此方程组可得: .
综上所述,原两点边值问题的解为
对一般的二阶微分方程边值问题
假定其解存在唯一,
为求解的近似值, 类似于前面的做法,
考虑带有第三类边界条件的二阶常系数微分方程边值问题
问题二:有限差分方法算出其数值解及误差
对于 原问题 ,取步长 h=0.2 ,用 有限差分 求其 近似解 ,并将结果与 精确解y(x)=-x-1 进行比较.
因为
先以将区间划分为5份为例,求出数值解
结果:
是不是解出数值解就完事了呢?当然不是。我们可以将问题的差分格式解与问题的真解进行比较,以得到解的可靠性。通过数学计算我们得到问题的真解为 ,现用范数来衡量误差的大小:
结果:
接下来绘图比较 时数值解与真解的差距:
结果:
将区间划分为 份, 即 时.
结果:
绘图比较 时数值解与真解的差距:
最后,我们还可以从数学的角度分析所采用的差分格式的一些性质。因为差分格式的误差为 , 所以理论上来说网格每加密一倍,与真解的误差大致会缩小到原来的 . 下讨论网格加密时的变化:
结果:
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,)。
from scipy import odeint
y = odeint(dy/dt=r*y*(1-y/k) ,y(0)=0.1,t)
对于微分方程全还给老师了,
这个地址有很多关于python做科学计算的文档,你可以去查查