Matlab 解常微分方程的初值问题
以下类容来源于: 精通matlab 5.0-张易华;清华出版社;1999年。
1:问题
常微分方程的初值问题的标准数学表述为:
?
??=<=<==0)(),,('y a y b t a y t f y ;我们要求解的任何高阶常微分方程都可以用替换法化为上式所示的一阶形式,其中y 为向量,yo 为初始值。
2:Matlab 中解决以上问题的步骤
(1):化方程组为标准形式。
例如:y ’’’-3y ’’-y ’y=0,y(0)=0,y ’(0)=1,y ’’(0)=-1.
把微分方程的高阶导数写为低阶导数的算式,即:
y ’’’=3y ’’+y ’y ,设:y1=y,y2=y ’,y3=y ’’,则原方程化为下列等价的方程组:
?????+===1233'33'22'1y y y y y y y y 满足初值条件:??
???-===1)0(31)0(20
)0(1y y y 已把该方程化成了标准形式。
其中:y ’->(y1’,y2’,y3’),a->(0,0,0),y0->(0,1,-1),f(t,y)->(y2,y3,3y3+y2y1).
(2):把微分方程组编成m 函数文件。
如:function dy=F(t,y)
dy=[y(2);y(3);3*y(3)+y(2)*y(1)];
注意:A:在函数文件里,虽然写微分方程时并不同时包含参数t 和y ,但第一行必须包含这两个输入变量。B:向量dy 必须为列向量。
(3):调用一个微分方程的求解函数求解。
[T,Y]=solver(‘F ’,tspan,y0);
其中:solver:求解函数名;
F:包含微分方程的m 文件;
tspan 为积分的数据范围,其格式为:[t0,tfinal];
y0为t0时刻的初值列向量。
输出参数T 和Y 为列向量
T 为时刻向量。
Y 表是不同时刻的函数值。
3:一个求解常微分方程初值问题的完整过程。
问题:求解方程y ’’-3(1-y^2)y ’+y=0在初值y ’(0)=3,y(0)=2的解。
○
1化成标准形式: 设y1=y,y2=y ’,则:???--==12)2^11(3'22
'1y y y y y y 初值为:???==3)0(22
)0(1y y
○
2编写函数文件ode.m,类容为: function dy=ode(t,y)
dy=[y(2);3*(1-y(1)^2)*y(2)-y(1)];
○
3调用函数ode45求解,时间区间为[0,20]: [T,Y]=ode45(‘ode ’,[0,20],[2;3]);
输出结果[T,Y]中T 为时间点组成的向量。Y 为对应于T 中时间点的y(1)和y(2)的值。