-
1 文档
-
2 视频
-
3 PPT
1.常微分方程求解方法
在Matlab中可以计算显式、线性隐式和完全隐式常微分方程初值问题的数值解。任何高阶常微分方程都可以变换成一阶微分方程，这是利用龙格-库塔法求解微分方程的前提。有高阶微分方程：
2.常微分方程求解函数
Matlab提供了多个求解常微分方程初始值问题的函数，见表4-3。
表4-3 常微分方程求解函数列表
| 函数 | 采用算法 | 特点 |
| ode45 | 四阶/五阶龙格-库塔 | 单步解法，速度较快，应作为你试解的首选。 |
| ode23 | 二阶/三阶龙格-库塔 | 单步解法，对误差限较大的稍带刚性问题的求解比前者有效。 |
| ode113 | 可变阶的Adams PECE算法 | 变阶次多步解法，在相同精度下，比ode45、ode23更快些，更适合求解高阶或需要大量计算的问题，不适合不连续的系统。 |
| ode15s | 可变阶的NDFS算法 | 刚性系统的变阶次多步解法，采用数值差分方法，速度较慢。 |
| ode23s | 基于改进的Rosenbrock公式 | 刚性系统固定阶次的单步解法，低阶方法。 |
| ode23t | 龙格-库塔公式采用梯形规则 | 单步解法，采用梯形法，用于不带无数值阻尼的适度刚性问题。 |
| ode23tb | 龙格-库塔公式的第一级采用TR梯形规则，第二级采用阶数为2的BDF反向微公式 | 单步解法，比ode15s更适合用于误差容许范围较宽的情况。 |
可将表中函数分成两大类：一类适用非刚性系统，一般的常微分方程都是非刚性系统，选用ode45、ode23、ode113来求解，使用一般的ode45、ode23、ode113来求解刚性系统，由于稳定性的要求可能会使计算步长变得很小，导致计算时间过长。另一类适用刚性系统，求解刚性系统函数包括ode15s、ode23s、ode23t、ode23tb。刚性问题指同时包含了快变环节和慢变环节的系统，例如
。
常微分方程求解函数格式基本相同，下面以ode45为例说明其使用。
[T,Y] = ode45(odefun,tspan,y0)
[T,Y] = ode45 (odefun,tspan,y0,options)
[T,Y,TE,YE,IE] = ode45 (odefun,tspan,y0,options)
sol = ode45 (odefun,[t0 tf],y0...)
说明：odefun代表为常微分方程组，其格式是
，
tspan是由两个元素组成的向量[t0,tfinal]表示求取t0~tfinal范围内的常微分方程组解，y0是初始值。
（2）然后，将方程组表示成函数形式。
function dydt=odefun1(t,y)
dydt=[0;0];
dydt(1)=y(2);
dydt(2)=-y(1)+1-t^2/(2*pi);
或
function dydt=odefun1(t,y)
u=1-(t.^2)/(pi*2);
dydt=[0,1;-1,0]*y+[0,1]'*u;
（3）应用ode45求解。
t0=0,tf=3*pi;
[t,y]=ode45('odefun1',[t0,tf], [0;0])
注意：本题求得的对应时间序列t的计算结果y有两列，第一列是数值积分解
，第二列是值对应的一阶导数值。
（2）然后，将方程组表示成函数odefun2。
function dydt=odefun2(t,y)
dydt=[0;0];
dydt(1)=y(2);
dydt(2)=-(y(1) ^2-1)*y(2)-y(1) ;
（3）应用ode45求解。
[t,y]=ode45('odefun2',[0,20],[0.25;0]);
plot(t,y(:,1));
figure(2);plot(t,y(:,2));
图4-9 Van del Pol方程解y 图4-10 Van del Pol方程解y的一阶导数
注意，调用微分函数，通常将微分方程化为一阶微分方程组的标准形式：
2.带边界条件常微分方程求解函数
Matlab提供了bvp4c、bvpinit，用于解常微分方程的边值问题。格式如下：
sol = bvp4c(odefun,bcfun,solinit)
sol = bvp4c(odefun,bcfun,solinit,options)
solinit = bvpinit(x, yinit, params)
说明：odefun为常微分方程函数，bcfun是边界条件函数，solinit为求解的初始猜测值，可以由函数bvpinit计算取得，solinit是一个结构体，包括x，y，parameters 三种属性，并且满足边界条件a= solinit.x(1)和b = solinit.x(end)。options结构体来设定解法器的参数。
（2）然后，将方程组表示成函数mat4ode。
function dydx=mat4ode(x,y,lambda)
q=5;
dydx=[y(2);-(lambda-2*q*cos(2*x))*y(1)];
（3）创建边界条件函数
function res=mat4bc(ya,yb,lambda)
res=[ya(2);yb(2);ya(1)-1];
表示边界条件时，必须写成等于零的形式。
（4）生成初始猜测值
function yinit=mat4init(x)
yinit=[cos(4*x);-4*sin(4*x)];
设未知参数猜测解等于15，调用该函数生成初始猜测值。
lambda=15;
solinit=bvpinit(linspace(0,pi,10),@mat4init,lambda);
（5）应用bvp4c求解并输出显示结果。
>> sol=bvp4c(@mat4ode,@mat4bc,solinit)
sol =
solver: 'bvp4c'
x: [1x37 double]
y: [2x37 double]
yp: [2x37 double]
parameters:17.0973
stats: [1x1 struct]
绘图显示结果。4阶Mathieu方程
的特征值等于17.0973。
plot(sol.x,sol.y(1,:));
图4-11 Mathieu方程在边界条件下的解
