MATLAB中ODE的使用
2013-10-17 12:33阅读：
ode23 解非刚性微分方程，低精度，使用Runge-Kutta法的二三阶算法。
ode45
解非刚性微分方程，中等精度，使用Runge-Kutta法的四五阶算法。
ode113
解非刚性微分方程，变精度变阶次Adams-Bashforth-Moulton PECE算法。
ode23t
解中等刚性微分方程，使用自由内插法的梯形法则。
ode15s 解刚性微分方程，使用可变阶次的数值微分（NDFs）算法。
ode23s
解刚性微分方程，低阶方法，使用修正的Rosenbrock公式。
ode23tb
解刚性微分方程，低阶方法，使用TR-BDF2方法，即Runger-Kutta公式的第一级采用梯形法则，第二级采用Gear法。
[t，YY]=solver('F',tspan,Yo)
解算ODE初值问题的最简调用格式。
solver指上面的指令。
tspan=[0,30];
％时域t的范围
y0=[1;0];
％y（1）y(2)的初始值
[tt,yy]=ode45(@DyDt,tspan,y0);
plot(tt,yy(:,1)),title('x(t)')
function
ydot=DyDt(t,y)
ydot=[y(2);
2*(1-y(1)^2)*y(2)-y(1)]
刚性方程：刚性是指其Jacobian矩阵的特征值相差十分悬殊。在解的性态上表现为，其中一些解变化缓慢，另一些变化快，且相差较悬殊，这类方程常常称为刚性方程，又称为Stiff方程。
刚性方程和非刚性方程对解法中步长选择的要求不同。
刚性方程一般不适合由ode45这类函数求解，而应该采用ode15s等。
如果不能分辨是否是刚性方程，先试用ode45，再用ode15s。
[t，YY，Te，Ye,Ie]
= solver('F',tspan,Yo,options,p1,p2,…)
解算ODE初值问题的最完整调用格式。
为了能够解出方程，要用指令odeset确定求解的条件和要求。在MATLAB中，求解方程组的指令都有默认的求解的条件和要求（由结构数组options表示），但可以用odeset修改或重新建立，也可以用odeget去获取已有的“优化选项”的信息。指令odeset和odeget用法介绍如下：
语句格式如下：
options=odeset(‘name1’,value1,’name2’,value2,…)
options=odeset(oldopts,‘name1’,value1,’name2’,value2,…)
options=odeset(oldopts,newopts)
odeset
第一种调用格式是指定各个参数的取值，对不指定取值的参数，取默认值。在不引起混淆的情况下，参数名可以只键入前面的几个字母，也不必区分大小写，如用“abst”表示AbsTol.但数值的输入必须格式正确，否则仍采用默认值。
第二种格式使用了原来的优化选项，但对其中的参数1等指定了新值。
第三种格式合并了两个优化选项oldopts
newopts，重复部分取newopts的指定值）：
第四种格式可在屏幕上显示如下全部可设置的参数及其默认值。
键入help odeset可查看全部参数的说明，下面对其中几个参数举例说明。
RelTol
相对误差，默认值为1e-3
AbsTol
绝对误差，默认值为1e-6
OutputFcn
输出方式，
默认值为‘odeplot’,其它选项有：
odeplot
按时间顺序画出全部变量的解
odephase2 二维相空间中两个变量的图形
odephase2 三维相空间中三个变量的图形
odeprint
打印输出
语句格式：
val=odeget(options,’name’)
这里options是由odeset设定的优化选项。
该语句从优化选项中提取指定的参数的取值。如果该参数没有指定值，则返回空阵［］。
options
odeget(options,'Reltol')
options=odeset(options,'Reltol',1e-6)
odeget(options,'Reltol')
function ydot=lorenfcn(t,y)
ydot=[-8/3*y(1)+y(2)*y(3);
-10*y(2)+10*y(3);
-y(2)*y(1)+35*y(2)-y(3)];
axis([10 50 -50 50 -50
50])
view(3)
hold
on
title('Lorenz
Attractor')
options=odeset('OutputFcn','odephas3');
[t,y]=ode23(@lorenfcn,[0
20],[0,0,eps],options);
％ EPS
Floating point relative
accuracy.
EPS returns the
distance from 1.0 to the
next largest
floating point number.
EPS is used as a default
tolerance by
PINV and RANK, as
well as several other MATLAB
functions.
微分方程的相关知识
1、微分方程的概念
含有未知的函数及其某些阶的导数以及自变量本身的方程称为微分方程。如果未知函数是一元函数，称为常微分方程。如果未知函数是多元函数，称为偏微分方程。联系一些未知函数的一组微分方程称为微分方程组。微分方程中出现的未知函数的导数的最高阶数称为微分方程的阶。如果方程中未知函数及其各阶导数都是一次的，称为线性常微分方程。若各系数为常数，称之为常系数（或定常、自治、时不变）的。
2、初等积分法
有些方程可以直接通过积分求解。例如，一阶常系数线性常微分方程
y’=ay+b (a!=0)
可化为
dy/(ay+b)=dt
两边积分可得通解为：
y(t)=Cexp(at)-a^-1b
其中C为任意常数
3、常系数线性微分方程
例1 求x’’+0.2x’+3.92x=0的通解。
解： 特征方程为
λ²+0.2λ+3.92=0
>> roots([1 0.2 3.92])
ans =
-0.1000
+ 1.9774i
-0.1000 - 1.9774i
求得共轭复根-0.1000
±1.9774i，从而通解为：
x(t)=Aexp(-0.1t)cos(1.9774t)+Bexp(-0.1t)sin(1.9774t)
其中A,B为任意常数。
4、初值问题数值解
除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程无显式解，应用中主要依靠数值解法。
考虑一阶常微分方程组初值问题
y’=f(t,y), t0
其中y=(y1,y2,…,ym)T，f=(f1,f2,…,fm)T，y0=(y10,y20,…,ym0)T，这里T表示转置。
所谓数值解就是寻求解y(t)在一系列离散点t0上的近似值yk(k=0,1,…,n)。称hk=t(k+1)-t(k)为步长，通常取为常量h。
高阶微分方程初值问题可以化为一阶常微分方程组，已给一个n阶方程，即
![MATLAB中ODE的使用]( "MATLAB中ODE的使用")
![MATLAB中ODE的使用]( "MATLAB中ODE的使用")
，则上式可化为一阶方程组
![MATLAB中ODE的使用]( "MATLAB中ODE的使用")
解：编写如下程序：
clear
f=inline('y-2*x/y','x','y');
a=0;
b=1;
h=0.1;
n=(b-a)/h;
x=zeros(1,n+1);
y=zeros(1,n+1);
y(1)=1;
for i=1:n+1
x(i)=a+(i-1)*h;
y(i+1)=y(i)+h*f(x(i),y(i));
end
x
y=y(1:n+1)
得到结果如下：
x =
Columns
1 through 10
0 0.1000
0.2000 0.3000
0.4000
0.5000 0.6000
0.7000 0.8000
0.9000
Column 11
1.0000
y =
Columns
1 through 10
1.0000
1.1000 1.1918
1.2774 1.3582
1.4351 1.5090
1.5803
1.6498 1.7178
Column
11
1.7848
>> % 为求精确解，输入
>>
f=dsolve('Dy=y-2*x/y','y(0)=1','x')
f =
(2*x +
1)^(1/2)
>> %
在一个坐标系中画出数值解与精确解的图形：
>>
plot(x,y,'.b')
>> hold on
>>
ezplot('y-sqrt(1+2*x)',[0,1,1,1.8])
得到图1。
![MATLAB中ODE的使用]( "MATLAB中ODE的使用")
图1
由图1可知，数值解和准确解比较接近，但还是有误差。
![MATLAB中ODE的使用]( "MATLAB中ODE的使用")
解：
>> fun=inline('-2*y+2*x^2+2*x','x','y');
>>
[x,y]=ode23(fun,[0,0.5],1);
>> x';
>> y';
>>
plot(x,y,'o-');
>> x'
ans =
Columns 1 through 10
0
0.0400 0.0900
0.1400 0.1900
0.2400
0.2900 0.3400
0.3900 0.4400
Columns
11 through 12
0.4900 0.5000
>> y'
ans =
Columns
1 through 10
1.0000
0.9247 0.8434
0.7754 0.7199
0.6764 0.6440
0.6222
0.6105 0.6084
Columns
11 through 12
0.6154 0.6179
图形结果见图2
![MATLAB中ODE的使用]( "MATLAB中ODE的使用")
图2
解微分方程的MATLAB命令
1、初值问题的求解
常微分方程求解MATLAB指令具有相同的格式，下面以最常用的ode45为例。
常用格式：[t,y]=ode45(odefun,tspan,y0)
参数说明：
odefun:用以表示f(t,y)的函数句柄或inline函数，t是标量，y是标量或向量。
tspan：如果是二维向量[t0,tf],表示自变量初值t0和终值tf;如果是高维向量[t0,t1,…,tn],则表示输出节点列向量。
y0:表示初始向量y0。
t：表示节点列向量(t0,t1,…,tn)T。
y: 表示数值解矩阵，每一列对应y的一个分量。
若无输出参数，则作出图形。
完整格式：[t,y]=ode45(odefun,tspan,y0,options,p1,p1,…)
options:
为计算参数（如精度要求）设置，默认可用空矩阵[]表示。
p1,p2,…:
为附加传递参数，这时的odefun表示f(t,y,p1,p2,…)。
注：ode45是最常用的求解微分方程的指令。它采用变步长四、五阶Runge-Kutta-Felhberg法，适合高精度问题，ode23与ode45类似，只是精度低一些。ode113是多步法，高低精度均可。这些指令对于刚性方程组不宜采用。ode23t,ode23s,ode23tb,ode15s都是求解刚性方程组的指令。
![MATLAB中ODE的使用]( "MATLAB中ODE的使用")
解：
（1）输入
>> dsolve('Dy=x+x*y','x')
得到含有一个任意常数的解析解：
ans =
C4*exp(x^2/2) –
1
（2）输入
>> r=dsolve('D2y=-a^2*y','y(0)=1','Dy(pi/a)=0')
r =
(1/exp(a*t*i))/2 +
exp(a*t*i)/2
>> % 化为最简式
>>
r=simple(r)
输出结果为:
r =
cos(a*t)
（3）输入
>> [x,y]=dsolve('Dx=y','Dy=-x')
x =
C10*cos(t) +
C9*sin(t)
y =
C9*cos(t) -
C10*sin(t)
例5 解微分方程
y’=y-2t/y,y(0)=1 (0
解：在指令窗口输入：
>> odefun=inline('y-2*t/y','t','y');
>>
[t,y]=ode45(odefun,[0,4],1);
>> [t,y]
ans =
0
1.0000
0.0502 1.0490
0.1005 1.0959
0.1507 1.1408
……
……
3.8507
2.9503
3.9005 2.9672
3.9502 2.9839
4.0000 3.0006
>> %
解函数图形表示，见图3
>>
plot(t,y,'o-')
>> ode45(odefun,[0,4],1);
>>
[t,y]=ode45(odefun,0:1:4,1);
>> [t,y]
ans =
0
1.0000
1.0000
1.7321
2.0000 2.2361
3.0000 2.6458
4.0000 3.0006
![MATLAB中ODE的使用]( "MATLAB中ODE的使用")
图3
事实上方程的准确解为y=sqrt(1+2t)。来比较一下几种方法的计算量和精度。下列结果中n为节点个数，反映计算量大小，e为每个节点均方误差。
>> odefun=inline('y-2*t/y','t','y');
>>
[t,y]=ode45(odefun,[0,4],1);
>>
n=length(t);
>> e=sqrt(sum((sqrt(1+2*t)-y).^2)/n);
>> [n,e]
ans =
45.0000 0.0002
>>
[t,y]=ode23(odefun,[0,4],1);
>>
n=length(t);
>> e=sqrt(sum((sqrt(1+2*t)-y).^2)/n);
>> [n,e]
ans =
13.0000 0.1905
>>
[t,y]=ode113(odefun,[0,4],1);
>>
n=length(t);
>> e=sqrt(sum((sqrt(1+2*t)-y).^2)/n);
>> [n,e]
ans =
18.0000 0.0097
>>
[t,y]=ode23t(odefun,[0,4],1);
>>
n=length(t);
>> e=sqrt(sum((sqrt(1+2*t)-y).^2)/n);
>> [n,e]
ans =
18.0000 0.0392
>>
[t,y]=ode23s(odefun,[0,4],1);
Warning: Failure
at t=3.925277e+000. Unable to
meet integration tolerances without
reducing the step size below
the smallest value
allowed
(1.394538e-014) at time t.
> In
ode23s at 478
>>
n=length(t);
>> e=sqrt(sum((sqrt(1+2*t)-y).^2)/n);
>>
[n,e]
ans =
81.0000 2.5437
>>
[t,y]=ode23tb(odefun,[0,4],1);
>>
n=length(t);
>> e=sqrt(sum((sqrt(1+2*t)-y).^2)/n);
>>
[n,e]
ans =
15.0000 0.2431
>>
[t,y]=ode15s(odefun,[0,4],1);
>>
n=length(t);
>> e=sqrt(sum((sqrt(1+2*t)-y).^2)/n);
>> [n,e]
ans =
22.0000 0.4551
可见，ode45精度高，但计算量比较大，ode23计算量小，但误差大。ode113适中。用刚性方程组解法解非刚性问题不太适合，特别是ode23s,计算量大且误差大。
![MATLAB中ODE的使用]( "MATLAB中ODE的使用")
解：将变量x,y合写成向量变量x,先写M函数eg6fun.m:
% eg6fun.m
function f=eg6fun(t,x)
f(1)=-x(1)^3-x(2);
f(2)=x(1)-x(2)^3;
% 保证f为列向量
f=f(:);
再在指令窗口中执行：
[t,x]=ode45(@eg6fun,[0,30],[1,0.5]);
subplot(1,2,1);
plot(t,x(:,1),t,x(:,2),':')
subplot(1,2,2)
plot(x(:,1),x(:,2));
图4 作出了解函数图和相平面图。
![MATLAB中ODE的使用]( "MATLAB中ODE的使用")
图4
![MATLAB中ODE的使用]( "MATLAB中ODE的使用")
先些M函数eg7fun.m:
% eg7fun.m
function f=eg7fun(t,y)
f=[y(2);y(3);-3*y(1)*y(3)+2*y(2)^2-y(4);y(5);-2.1*y(1)*y(5)];
再在命令窗口中执行：
>> y0=[0 0 0.68 1
-0.5];
>>
[t,y]=ode45(@eg7fun,[0,5],y0);
>>
plot(t,y(:,1),t,y(:,4),':r')
图5作出了f和T的图：
![MATLAB中ODE的使用]( "MATLAB中ODE的使用")
图5 ![MATLAB中ODE的使用]( "MATLAB中ODE的使用")
编写M函数文件eg8fun.m
% eg8fun.m
function f=eg8fun(t,x)
f=[x(2);7*(1-x(1)^2)*x(2)-x(1)];
再在命令窗口中输入：
>> y0=[1;0];
>> [t,x]=ode45('eg8fun',[0,40],y0);
>> y=x(:,1);
>>
plot(t,y)
![MATLAB中ODE的使用]( "MATLAB中ODE的使用")
图6
2、边值问题解法 ![MATLAB中ODE的使用]( "MATLAB中ODE的使用")
sinit=bvpinit(tinit,yinit)：由在粗略节点tinit的预估解yiniy生成粗略解网络sinit。
sol=bvp4c(odefun,bcfun,sinit)：odefun是微分方程组函数，bcfun为边值条件函数，sinit是由bvpinit得到的粗略解网络。求的边值问题解sol是一个结构，sol.x为求解节点，sol.y是y(t)的数值解。
sx=deval_r(sol,ti)：计算由bvp4c得到的解在ti的值。
例9 求解边值问题
z’’+|z|=0， z(0)=0， z(4)=-2
解：首先改写为标准形式。令y(1)=z,y(2)=z’，则方程为
y’(1)=y(2)，y’(2)=-|y(1)|
边界条件为
ya(1)=0，yb(1)+2=0
求解用下列M文件eg9fun.m。注意sol的域名。
clear;
close;
% 注意sinit的域名
sinit=bvpinit(0:4,[1;0])
odefun=inline('[y(2);-abs(y(1))]','t','y');
bcfun=inline('[ya(1),yb(1)+2]','ya','yb');
% 注意sol的域名
sol=bvp4c(odefun.bcfun,sinit)
t=linspace(0,4,101);
y=deval_r(sol,t);
plot(t,y(1,:),sol.x(1,:),'o',sinit.x,sinit.y(1,:),'s')
legend('解曲线','解点','粗略解')
