[Python] 科学计算:Python VS. MATLAB(5)----常微分方程数值解

一、常微分方程的一般理论

凡含有参数,未知函数和未知函数导数 (或微分) 的方程,称为微分方程,有时简称为方程,未知函数是一元函数的微分方程称作常微分方程,未知数是多元函数的微分方程称作偏微分方程。微分方程中出现的未知函数最高阶导数的阶数,称为微分方程的阶。

二、使用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,)。而我们会发现,高阶方程的标准化工作,其实是解微分方程最主要的工作。下面是一个二阶常微分方程的例子,通过这个例子,学会变化标准形式的方法。
典型的二阶方程:y”+ay’+by=0 ,不满足标准式。变量替换,假设输入是Y=[y,y’],输出是Y’=[y’,y”],那么方程就化为了Y’=F(Y,t)的形式。我们知道Y[0]=y, Y[1]=y’,所以输出形式就是:

y’=Y[1]

y”=-ay’-by=-aY[1]-BY[0]

有了这个形式,我们很容易写出Python代码:

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用的。
hold(‘on’)
plot(time,y[:,1],label=”y’”) #y[:,1]即返回值的第二列,是y’的值
xlabel(‘t’)
ylabel(‘y’)
legend()
show()


三、使用MATLAB

同样,MATLAB解微分方程的第一步工作也是化简微分方程为一阶的标准形式,这个和Python是完全一样的。只是函数的编写不一样。MATLAB可以在m文件函数中定义一个函数和多个子函数,但是子函数只能由同一m文件中的函数调用。并且还要求m文件名必须和主函数同名;而Python无此限制,在此可见一斑,也能看出Python代码的相对更优美。另外,Python的索引从0开始,MATLAB的从1开始,为了这个习惯,编写函数的时候也要有所注意。

MATLAB在这里也有一个优势,就是提供了多个求解函数,以应对时间需求、精度需求和刚性问题等要求。这就比scipy中的函数灵活了不少。如果使用scipy,只能自己编写函数了!

这些函数几乎统一的调用方式是:[to,yo]=solver(func,tspan,y0,options)

其中,solver是ode45(最常用)、ode23、ode23t、ode23s、ode23tb等

例1:y”+4y=3sin(t) y(0)=1,y’(0)=0

定义一个函数,描写这个微分方程:

function dy=fun1(t,y)

dy(1,1)=y(2)

dy(2,1)=3sin(t)-4y(1)

将这个函数存为fun1.m, 然后写一个脚本调用它,解方程。不妨这个脚本就叫做callfun1.m

tspan=[0,5];
y0=[1,0];
[t,y]=ode45(@fun1,tspan,y0);
figure;
plot(t,y(:,1),’k-‘);
hold on
plot(t,y(:,2),’k:’);
xlabel(‘t’);

画出来的图是这样子的:

三、对比

下面用Python实现MATLAB部分的例子。

from scipy.integrate import odeint
from pylab import *

#定义函数

def deriv(y,t):
return array([ y[1], 3sin(t)-4y[0] ])

#设置初值和t区间

time = linspace(0.0,5.0,1000)
yinit = array([1,0])

#解方程并绘图

y = odeint(deriv,yinit,time)

figure()
plot(time,y[:,0],label=’y’)
hold(‘on’)
plot(time,y[:,1],label=”y’”)
xlabel(‘t’)
ylabel(‘y’)
legend()
show()
结果如图,可以和MATLAB绘制出来的对比一下:


下面用MATLAB实现Python部分的例子。

函数脚本(fun1.m):

function dy=fun1(t,y)
dy(1,1)=y(2)
dy(2,1)=-2.0y(1)-0.1y(2)

调用脚本(callfun1.m):

tspan=[0,50];
y0=[0.0005,0.2];
[t,y]=ode45(@fun1,tspan,y0);
figure;
plot(t,y(:,1),’k-‘);
hold on
plot(t,y(:,2),’k:’);
xlabel(‘t’);

绘图结果:

补充说明:

MATLAB还可以使用inline的形式调用函数。比如:

tspan=[0,5];
y0=[1,0]
fun1=inline(‘[y(2);3sin(t)-4y(1)]’,’t’,’y’); %注意这一句
[t,y]=ode45(fun1,tspan,y0);
figure;
plot(t,y(:,1),’k-‘);
hold on
plot(t,y(:,2),’k:’);
xlabel(‘t’);

这时候,不需要关心文件名的问题了,但是这种写法明显不利于推广和模块化,所以不推荐使用这样的形式。