Matlab解常微分方程的初值问题

Matlab解常微分方程的初值问题以下类容来源于:精通matlab5.0－张易华；清华出版社；1999年。1：问题常微分方程的初值问题的标准数学表述为：0)(),,('yaybtaytfy；我们要求解的任何高阶常微分方程都可以用替换法化为上式所示的一阶形式，其中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'1yyyyyyyy满足初值条件：1)0(31)0(20)0(1yyy已把该方程化成了标准形式。其中：y’-（y1’,y2’,y3’）,a-(0,0,0),y0-(0,1,-1),f(t,y)-(y2,y3,3y3+y2y1).(2):把微分方程组编成m函数文件。如：functiondy=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'1yyyyyy初值为：3)0(22)0(1yy○2编写函数文件ode.m,类容为：functiondy=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)的值。○4绘制解的曲线，结果如图。plot(T,Y(:,1),’-’,T,Y(:,2),’--’)title(‘SolutionofODEEquation’);xlabel(‘timeT’)ylabel(‘solutionY’);legend(‘Y1’,’Y2’)Matlab利用数值方法来求解常微分方程的解，其思路如下：把求解的时间区间划分成有限步，对应于每一步将计算出一个解，如果求得的解不满足误差限制，则减少步长，再求解。如此重复，直到满足误差限为止。○a刚性问题（stiff）:方程组的解不同分量的数量级差别较大，对于数值求解是一大困难。Matlab既能解决非刚性问题，也能解决刚性问题。○b三个解决非刚性问题的函数：ode45,ode23,ode113○c两个解刚性问题的函数：ode15s和ode23s2005-03-18