常微分方程数值求解的命令:
求常微分方程的数值解,MATLAB的命令格式为:
[t,y]=solver('odefun',tspan,y0,options)
其中solver选择ode45等函数名,odefun为根据待解方程或方程组编写的m文件名,tspan为自变量的区间[t0,tf],即准备在那个区间上求解,y0表示初始值,options用于设定误差限制。命令格式为:
options=odeset('reltol',rt,'abstol',at)
rt输入相对误差,at输入绝对误差。
常用的函数:
函数名 | 简介 | 适用对象 |
ode45 | 单步,4/5阶龙格库塔法 | 大部分ODE |
ode23 | 单步,2/3阶龙格库塔法 | 快速、精度不高的求解 |
0de113 | 多步,Adams算法 | 误差要求严格或计算复杂 |
注:上述函数仅适用于非刚性方程(组)。
函数名 | 简介 | 适用对象 |
ode23t | 采用梯形算法 | 具有一定的刚性特点 |
ode15s | 多步,反向数值积分法 | ode45失效时可以使用 |
ode23s | 单步,2阶Rosebrock算法 | 精度设定较低时,速度快 |
ode23tb | 采用梯形算法 | 精度设定较低时,速度快 |
例1:
函数的M文件:
function dx=human(t,x)
dx=0.05*x;
使用命令:
>> [t,x]=ode45('human',[0,100],1000);
>> plot(t,x)
>> hold on
>> x=1000*exp(0.05*t);
>> plot(t,x,'*')
绘图比较真实解和数值解。
例2:
函数的M文件:
function dx=compete(t,x)
dx=zeros(2,1);
dx(1)=0.01*x(1)*(1-x(1)/50000-0.1*x(2)/60000);
dx(2)=0.02*x(2)*(1-0.2*x(1)/50000-x(2)/60000);
使用命令:
>> [t,x]=ode45('compete',[0,500],[10,10]);
>> plot(t,x(:,1),t,x(:,2))
>> plot(x(:,1),x(:,2))