MATLAB-6-5常微分方程数值求解

1. 常微分方程数值求解的一般概念

2. 常微分方程数值求解函数

[t,y]=solver(filename,tspan,y0,option)
t：时间向量； y：数值解；
filename：定义f(t,y)的函数名，该函数必须返回一个列向量；
tspan形式为[t0,tf]，表示求解区间； y0：初始状态向量；
option：可选参数，用于设置求解属性，常用的属性包括相对误差值RelTol（默认值为10^-3和绝对误差值AbsTol（默认值为10^-6）；

3. 常微分方程数值求解函数统一命名格式

odennxx

• ode：Ordinary Differential Equation，常微分方程
• nn：数字，代表所用方法的阶数。
  ode23采用2阶龙格-库塔（Runge-Kutta）算法，用3阶公式做误差估计来调节步长，具有低等精度。
  ode45采用4阶龙格-库塔算法，用5阶公式做误差估计来调节步长，具有中等精度。
• xx：字母，用于标注方法的专门特征。
  ode15s、ode23s中的“s”代表（Stiff），表示函数适用于刚性方程。
  
  e.g.求解微分方程初值问题，并与精确解y(t)=sqrt(t+1)+1，进行比较。
  

>> f=@(t,y)(y^2-t-2)/4/(t+1);
>> [t,y]=ode23(f,[0,10],2);
>> y1=sqrt(t+1)+1;
>> plot(t,y,'b:',t,y1,'r')

e.g.已知一个二阶线性系统的微分方程为：

>> f=@(t,x) [-2,0;0,1]*[x(2);x(1)];
>> [t,x]=ode45(f,[0,20],[1,0]);
>> subplot(2,2,1);
>> plot(t,x(:,2));
>> subplot(2,2,2);
>> plot(x(:,2),x(:,1));

4. 刚性问题

有一类常微分方程，其解的分量有的变化很快，有的变化很慢，且相差悬殊，这就是所谓的刚性问题(Stiff)。对于刚性问题，数值解算法必须取很小步长才能获得满意的结果，导致计算量会大大增加。解决刚性问题需要有专门方法。非刚性算法可以求解刚性问题，只不过需要很长的计算时间。

>> lambda=0.01;
>> f=@(t,y) y^2-y^3;
>> tic;
>> [t,y]=ode45(f,[0,2/lambda],lambda);
>> toc
时间已过 32.701709 秒。
>> disp(['ode45计算的点数' num2str(length(t))]);
ode45计算的点数157
%tic和toc函数用来记录微分方程求解命令执行的时间，使用tic函数启动计时器，使用toc函数显示从计时器启动到当前所经历的时间。最后还输出计算的点数，运行结果表明这时常微分方程不算很刚性。

>> lambda=1e-5;
>> f=@(t,y) y^2-y^3;
>> tic;
>> f=@(t,y) y^2-y^3;
>> toc
时间已过 7.096559 秒。
>> disp(['ode45计算的点数' num2str(length(t))]);
ode45计算的点数157
这时计算时间明显加长，计算的点数剧增，表明这时常微分方程表现为刚性

>> lambda=1e-5;
>> f=@(t,y) y^2-y^3;
>> tic;
>> [t,y]=ode15s(f,[0,2/lambda],lambda);
>> toc
时间已过 8.699999 秒。
>> disp(['ode15s计算的点数' num2str(length(t))]);
ode15s计算的点数98
对于刚性问题，选择以“s”结尾的函数，专门用于求解刚性方程。计算时间大大缩短，计算的点数大大减少。
相同程序，每次的计算时间都不相同