第七章 微分方程问题的解法
? 微分方程的解析解方法
? 微分方程问题的数值解法
– 微分方程问题算法概述
– 四阶定步长 Runge-Kutta算法及 MATLAB 实现
– 一阶微分方程组的数值解
– 微分方程转换
? 特殊微分方程的数值解
7.1 微分方程的解析解方法
? 格式,
y=dsolve(f1,f2,…,f m)
? 格式：指明自变量
y=dsolve(f1,f2,…,f m,’x’)
fi即可以描述微分方程，又可描述初始条件
或边界条件。如,
描述微分方程时
描述条件时
( 4 ) ( ) 7 4 7y t D y? ? ?
( 2 ) 3 2 ( 2 ) 3y D y? ? ?&&
例,
>> syms t; u=exp(-5*t)*cos(2*t+1)+5;
>> uu=5*diff(u,t,2)+4*diff(u,t)+2*u
uu =
87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10
>> syms t y;
>> y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',..,
'87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10'])
>> y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',..,
'87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1),.,
+10'],'y(0)=3','Dy(0)=2','D2y(0)=0','D3y(0)=0')
分别处理系数，如,
>> [n,d]=rat(double(vpa(-445/26*cos(1)-51/13*sin(1)-69/2)))]
ans =
-8704 185 % rat()最接近有理数的分数
判断误差,
>> vpa(-445/26*cos(sym(1))-51/13*sin(1)-69/2+8704/185)
ans =
.114731975864790922564144636e-4
>> y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',..,
'87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1) +,.,
10'],'y(0)=1/2','Dy(pi)=1','D2y(2*pi)=0','Dy(2*pi)=1/5');
如果用推导的方法求 Ci的值，每个系数的解析解至少要写出 10
数行,故可采用有理式近似 的方式表示,
>> vpa(y,10) %有理近似值
ans =
1.196361839*exp(-5.*t)+.4166666667-
.4785447354*sin(t)*cos(t)*exp(-5.*t)-.4519262218e-
1*cos(2.*t)*exp(-5.*t)-2.392723677*cos(t)^2*exp(-
5.*t)+.2259631109*sin(2.*t)*exp(-5.*t)-473690.0893*exp(-
3.*t)+31319.63786*exp(-2.*t)-219.1293619*exp(-
1.*t)+442590.9059*exp(-4.*t)
? 例,
>> syms t x
>> x=dsolve('Dx=x*(1-x^2)')
x =
[ 1/(1+exp(-2*t)*C1)^(1/2)]
[ -1/(1+exp(-2*t)*C1)^(1/2)]
>> syms t x; x=dsolve('Dx=x*(1-x^2)+1')
Warning,Explicit solution could not be found; implicit solution
returned,
> In D:\MATLAB6p5\toolbox\symbolic\dsolve.m at line 292
x =
t-Int(1/(a-a^3+1),a=``..x)+C1=0
故只有部分非线性微分方程有解析解。
7.2 微分方程问题的数值解法
7.2.1 微分方程问题算法概述
微分方程求解的误差与步长问题,
7.2.2 四阶定步长 Runge-Kutta算法
及 MATLAB 实现
function [tout,yout]=rk_4(odefile,tspan,y0) ％ y0初值列向量
t0=tspan(1); th=tspan(2);
if length(tspan)<=3,h=tspan(3); % tspan=[t0,th,h]
else,h=tspan(2)-tspan(1); th=tspan(end); end ％等间距数组
tout=[t0:h:th]'; yout=[];
for t=tout'
k1=h*eval([odefile ‘(t,y0)’]); % odefile是一个字符串变
量，为表示微分方程的文件名。
k2=h*eval([odefile '(t+h/2,y0+0.5*k1)']);
k3=h*eval([odefile '(t+h/2,y0+0.5*k2)']);
k4=h*eval([odefile '(t+h,y0+k3)']);
y0=y0+(k1+2*k2+2*k3+k4)/6;
yout=[yout; y0'];
end ％由效果看，该算法不是一个较好的方法。
7.2.3 一阶微分方程组的数值解
7.2.3.1 四阶五级 Runge-Kutta-Felhberg算法
通过误差向量调节步长，此为自动变步长方法。
四阶五级 RKF算法有参量系数表。
7.2.3.2 基于 MATLAB 的微分方程
求解函数
格式 1,直接求解
[t,x]=ode45(Fun,[t0,tf],x0)
格式 2,带有控制参数
[t,x]=ode45(Fun,[t0,tf],x0,options)
格式 3,带有附加参数
[t,x]=ode45(Fun,[t0,tf],x0,options,p1,p2,… )
[t0,tf]求解区间,x0初值问题的初始状态变量 。
描述需要求解的微分方程组,
不需附加变量的格式
function xd=funname(t,x)
可以使用附加变量
function xd=funname(t,x,flag,p1,p2,…)
％ t是时间变量或自变量（必须给）,x为状态向量,
xd为状态向量的导数。 flag用来控制求解过程,指定初
值，即使初值不用指定，也必须有该变量占位。
修改变量,options唯一结构体变量，用 odeset( )修改。
options=odeset(‘RelTol’,1e-7);
options= odeset; options,RelTol= 1e-7;
? 例,
自变函数
function xdot = lorenzeq(t,x)
xdot=[-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3);…
-x(1)*x(2)+28*x(2)-x(3)];
>> t_final=100; x0=[0;0;1e-10]; % t_final为设定的仿真终止时间
>> [t,x]=ode45('lorenzeq',[0,t_final],x0); plot(t,x),
>> figure; % 打开新图形窗口
>> plot3(x(:,1),x(:,2),x(:,3));
>> axis([10 42 -20 20 -20 25]); % 根据实际数值手动设
置坐标系
? 可采用 comet3( )函数绘制动画式的轨迹。
>> comet3(x(:,1),x(:,2),x(:,3))
? 描述微分方程是常微分方程初值问题数值求解
的关键。
>> f1=inline(['[-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3);',..,
'-x(1)*x(2)+28*x(2)-x(3)]'],'t','x');
>> t_final=100; x0=[0;0;1e-10];
>> [t,x]=ode45(f1,[0,t_final],x0);
>> plot(t,x),figure;
>> plot3(x(:,1),x(:,2),x(:,3)); axis([10 42 -20 20 -20 25]);
得出完全一致的结果。
7.2.3.3 MATLAB 下带有附加参数的微
分方程求解
? 例,
? 编写函数
function xdot=lorenz1(t,x,flag,beta,rho,sigma)
％ flag变量是不能省略的
xdot=[-beta*x(1)+x(2)*x(3);
-rho*x(2)+rho*x(3);
-x(1)*x(2)+sigma*x(2)-x(3)];
求微分方程,
>> t_final=100; x0=[0;0;1e-10];
>> b2=2; r2=5; s2=20;
>> [t2,x2]=ode45('lorenz1',[0,t_final],x0,[],b2,r2,s2);
>> plot(t2,x2),％ options位置为 []，表示不需修改控制选项
>> figure; plot3(x2(:,1),x2(:,2),x2(:,3)); axis([0 72 -20 22
-35 40]);
f2=inline(['[-beta*x(1)+x(2)*x(3); -rho*x(2)+rho*x(3);',..,
'-x(1)*x(2)+sigma*x(2)-x(3)]'],…
't','x','flag','beta','rho','sigma');
％ flag变量是不能省略的
7.2.4 微分方程转换
7.2.4.1 单个高阶常微分方程处理方法
? 例,
? 函数描述为,
function y=vdp_eq(t,x,flag,mu)
y=[x(2); -mu*(x(1)^2-1)*x(2)-x(1)];
>> x0=[-0.2; -0.7]; t_final=20;
>> mu=1; [t1,y1]=ode45('vdp_eq',[0,t_final],x0,[],mu);
>> mu=2; [t2,y2]=ode45('vdp_eq',[0,t_final],x0,[],mu);
>> plot(t1,y1,t2,y2,':')
>> figure; plot(y1(:,1),y1(:,2),y2(:,1),y2(:,2),':')
>> x0=[2;0]; t_final=3000;
>> mu=1000; [t,y]=ode45('vdp_eq',[0,t_final],x0,[],mu);
由于变步长所采用的步长过小，所需时间较长，导致输出的 y
矩阵过大，超出计算机存储空间容量。所以不适合采用
ode45()来求解，可用刚性方程求解算法 ode15s( )。
7.2.4.2 高阶常微分方程组的变换方法
? 例,
? 描述函数,
function dx=apolloeq(t,x)
mu=1/82.45; mu1=1-mu;
r1=sqrt((x(1)+mu)^2+x(3)^2);
r2=sqrt((x(1)-mu1)^2+x(3)^2);
dx=[x(2);
2*x(4)+x(1)-mu1*(x(1)+mu)/r1^3-mu*(x(1)-mu1)/r2^3;
x(4);
-2*x(2)+x(3)-mu1*x(3)/r1^3-mu*x(3)/r2^3];
? 求解,
>> x0=[1.2; 0; 0; -1.04935751];
>> tic,[t,y]=ode45('apolloeq',[0,20],x0); toc
elapsed_time =
0.8310
>> length(t),
>> plot(y(:,1),y(:,3))
ans =
689
得出的轨道不正确,
默认精度 RelTol设置
得太大，从而导致的
误差传递，可减小该
值。
? 改变精度,
>> options=odeset; options.RelTol=1e-6;
>> tic,[t1,y1]=ode45('apolloeq',[0,20],x0,options); toc
elapsed_time =
0.8110
>> length(t1),
>> plot(y1(:,1),y1(:,3)),
ans =
1873
>> min(diff(t1))
ans =
1.8927e-004
>> plot(t1(1:end-1),…
diff(t1))
? 例,
>> x0=[1.2; 0; 0; -1.04935751];
>> tic,[t1,y1]=rk_4('apolloeq',[0,20,0.01],x0); toc
elapsed_time =
4.2570
>> plot(y1(:,1),y1(:,3))
% 绘制出轨迹曲线
显而易见，这样求解
是错误的，应该采用
更小的步长。
>> tic,[t2,y2]=rk_4('apolloeq',[0,20,0.001],x0); toc
elapsed_time =
124.4990 ％计算时间过长
>> plot(y2(:,1),y2(:,3))
% 绘制出轨迹曲线
严格说来某些点仍不
满足 10－ 6的误差限,
所以求解常微分方程
组时建议采用变步长
算法，而不是定步长
算法。
? 例,
用 MATLAB符号工具箱求解,
令 ％
>> syms x1 x2 x3 x4
>> [dx,dy]=solve('dx+2*x4*x1=2*dy','dx*x4+ …
3*x2*dy+x1*x4-x3=5','dx,dy')
dx =
-2*(3*x4*x1*x2+x4*x1-x3-5)/(2*x4+3*x2)
dy =
(2*x4^2*x1-x4*x1+x3+5)/(2*x4+3*x2)
对于更复杂的问题来说，手工变换的难度将很
大，所以如有可能，可采用计算机去求解有关方程，
获得解析解。如不能得到解析解，也需要在描写一
阶常微分方程组时列写出式子，得出问题的数值解。
1 2 3 4,,,,,x x x x x y x y d x x d y y? ? ? ? ? ?& & && &&
7.3特殊微分方程的数值解
7.3.1 刚性微分方程的求解
? 刚性微分方程
一类特殊的常微分方程，其中一些解变
化缓慢，另一些变化快，且相差悬殊，这类
方程常常称为刚性方程。
MATLAB采用求解函数 ode15s()，该函数的调用
格式和 ode45()完全一致。
[t,x]=ode15s(Fun,[t0,tf],x0,options,p1,p2,… )
? 例,
％计算
>> h_opt=odeset; h_opt.RelTol=1e-6;
>> x0=[2;0]; t_final=3000;
>> tic,mu=1000;
[t,y]=ode15s('vdp_eq',[0,t_final],x0,h_opt,mu);
toc
elapsed_time =
2.5240
％作图
>> plot(t,y(:,1)); figure; plot(t,y(:,2))
y(:,1)曲线变化较平滑,y(:,2)变化在某些点上
较快。
? 例,
定义函数
function dy=c7exstf2(t,y)
dy=[0.04*(1-y(1))-(1-y(2))*y(1)+0.0001*(1-y(2))^2;
-10^4*y(1)+3000*(1-y(2))^2];
? 方法一
>> tic,[t2,y2]=ode45('c7exstf2',[0,100],[0;1]); toc
elapsed_time =
229.4700
>> length(t2),plot(t2,y2)
ans =
356941
? 步长分析,
>> format long,[min(diff(t2)),max(diff(t2))]
ans =
0.00022220693884 0.00214971787184
>> plot(t2(1:end-1),diff(t2))
? 方法二，用 ode15s()代替 ode45()
>> opt=odeset; opt.RelTol=1e-6;
>> tic,[t1,y1]=ode15s('c7exstf2',[0,100],[0;1],opt); toc
elapsed_time =
0.49100000000000
>> length(t1),
>> plot(t1,y1)
ans =
169
7.3.2 隐式微分方程求解
? 例,
? 编写函数,
function dx=c6ximp(t,x)
A=[sin(x(1)) cos(x(2)); -cos(x(2)) sin(x(1))];
B=[1-x(1); -x(2)]; dx=inv(A)*B;
求解,
>> opt=odeset; opt.RelTol=1e-6;
>> [t,x]=ode45('c7ximp',[0,10],[0; 0],opt); plot(t,x)
? 微分方程的解析解方法
? 微分方程问题的数值解法
– 微分方程问题算法概述
– 四阶定步长 Runge-Kutta算法及 MATLAB 实现
– 一阶微分方程组的数值解
– 微分方程转换
? 特殊微分方程的数值解
7.1 微分方程的解析解方法
? 格式,
y=dsolve(f1,f2,…,f m)
? 格式：指明自变量
y=dsolve(f1,f2,…,f m,’x’)
fi即可以描述微分方程，又可描述初始条件
或边界条件。如,
描述微分方程时
描述条件时
( 4 ) ( ) 7 4 7y t D y? ? ?
( 2 ) 3 2 ( 2 ) 3y D y? ? ?&&
例,
>> syms t; u=exp(-5*t)*cos(2*t+1)+5;
>> uu=5*diff(u,t,2)+4*diff(u,t)+2*u
uu =
87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10
>> syms t y;
>> y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',..,
'87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10'])
>> y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',..,
'87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1),.,
+10'],'y(0)=3','Dy(0)=2','D2y(0)=0','D3y(0)=0')
分别处理系数，如,
>> [n,d]=rat(double(vpa(-445/26*cos(1)-51/13*sin(1)-69/2)))]
ans =
-8704 185 % rat()最接近有理数的分数
判断误差,
>> vpa(-445/26*cos(sym(1))-51/13*sin(1)-69/2+8704/185)
ans =
.114731975864790922564144636e-4
>> y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',..,
'87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1) +,.,
10'],'y(0)=1/2','Dy(pi)=1','D2y(2*pi)=0','Dy(2*pi)=1/5');
如果用推导的方法求 Ci的值，每个系数的解析解至少要写出 10
数行,故可采用有理式近似 的方式表示,
>> vpa(y,10) %有理近似值
ans =
1.196361839*exp(-5.*t)+.4166666667-
.4785447354*sin(t)*cos(t)*exp(-5.*t)-.4519262218e-
1*cos(2.*t)*exp(-5.*t)-2.392723677*cos(t)^2*exp(-
5.*t)+.2259631109*sin(2.*t)*exp(-5.*t)-473690.0893*exp(-
3.*t)+31319.63786*exp(-2.*t)-219.1293619*exp(-
1.*t)+442590.9059*exp(-4.*t)
? 例,
>> syms t x
>> x=dsolve('Dx=x*(1-x^2)')
x =
[ 1/(1+exp(-2*t)*C1)^(1/2)]
[ -1/(1+exp(-2*t)*C1)^(1/2)]
>> syms t x; x=dsolve('Dx=x*(1-x^2)+1')
Warning,Explicit solution could not be found; implicit solution
returned,
> In D:\MATLAB6p5\toolbox\symbolic\dsolve.m at line 292
x =
t-Int(1/(a-a^3+1),a=``..x)+C1=0
故只有部分非线性微分方程有解析解。
7.2 微分方程问题的数值解法
7.2.1 微分方程问题算法概述
微分方程求解的误差与步长问题,
7.2.2 四阶定步长 Runge-Kutta算法
及 MATLAB 实现
function [tout,yout]=rk_4(odefile,tspan,y0) ％ y0初值列向量
t0=tspan(1); th=tspan(2);
if length(tspan)<=3,h=tspan(3); % tspan=[t0,th,h]
else,h=tspan(2)-tspan(1); th=tspan(end); end ％等间距数组
tout=[t0:h:th]'; yout=[];
for t=tout'
k1=h*eval([odefile ‘(t,y0)’]); % odefile是一个字符串变
量，为表示微分方程的文件名。
k2=h*eval([odefile '(t+h/2,y0+0.5*k1)']);
k3=h*eval([odefile '(t+h/2,y0+0.5*k2)']);
k4=h*eval([odefile '(t+h,y0+k3)']);
y0=y0+(k1+2*k2+2*k3+k4)/6;
yout=[yout; y0'];
end ％由效果看，该算法不是一个较好的方法。
7.2.3 一阶微分方程组的数值解
7.2.3.1 四阶五级 Runge-Kutta-Felhberg算法
通过误差向量调节步长，此为自动变步长方法。
四阶五级 RKF算法有参量系数表。
7.2.3.2 基于 MATLAB 的微分方程
求解函数
格式 1,直接求解
[t,x]=ode45(Fun,[t0,tf],x0)
格式 2,带有控制参数
[t,x]=ode45(Fun,[t0,tf],x0,options)
格式 3,带有附加参数
[t,x]=ode45(Fun,[t0,tf],x0,options,p1,p2,… )
[t0,tf]求解区间,x0初值问题的初始状态变量 。
描述需要求解的微分方程组,
不需附加变量的格式
function xd=funname(t,x)
可以使用附加变量
function xd=funname(t,x,flag,p1,p2,…)
％ t是时间变量或自变量（必须给）,x为状态向量,
xd为状态向量的导数。 flag用来控制求解过程,指定初
值，即使初值不用指定，也必须有该变量占位。
修改变量,options唯一结构体变量，用 odeset( )修改。
options=odeset(‘RelTol’,1e-7);
options= odeset; options,RelTol= 1e-7;
? 例,
自变函数
function xdot = lorenzeq(t,x)
xdot=[-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3);…
-x(1)*x(2)+28*x(2)-x(3)];
>> t_final=100; x0=[0;0;1e-10]; % t_final为设定的仿真终止时间
>> [t,x]=ode45('lorenzeq',[0,t_final],x0); plot(t,x),
>> figure; % 打开新图形窗口
>> plot3(x(:,1),x(:,2),x(:,3));
>> axis([10 42 -20 20 -20 25]); % 根据实际数值手动设
置坐标系
? 可采用 comet3( )函数绘制动画式的轨迹。
>> comet3(x(:,1),x(:,2),x(:,3))
? 描述微分方程是常微分方程初值问题数值求解
的关键。
>> f1=inline(['[-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3);',..,
'-x(1)*x(2)+28*x(2)-x(3)]'],'t','x');
>> t_final=100; x0=[0;0;1e-10];
>> [t,x]=ode45(f1,[0,t_final],x0);
>> plot(t,x),figure;
>> plot3(x(:,1),x(:,2),x(:,3)); axis([10 42 -20 20 -20 25]);
得出完全一致的结果。
7.2.3.3 MATLAB 下带有附加参数的微
分方程求解
? 例,
? 编写函数
function xdot=lorenz1(t,x,flag,beta,rho,sigma)
％ flag变量是不能省略的
xdot=[-beta*x(1)+x(2)*x(3);
-rho*x(2)+rho*x(3);
-x(1)*x(2)+sigma*x(2)-x(3)];
求微分方程,
>> t_final=100; x0=[0;0;1e-10];
>> b2=2; r2=5; s2=20;
>> [t2,x2]=ode45('lorenz1',[0,t_final],x0,[],b2,r2,s2);
>> plot(t2,x2),％ options位置为 []，表示不需修改控制选项
>> figure; plot3(x2(:,1),x2(:,2),x2(:,3)); axis([0 72 -20 22
-35 40]);
f2=inline(['[-beta*x(1)+x(2)*x(3); -rho*x(2)+rho*x(3);',..,
'-x(1)*x(2)+sigma*x(2)-x(3)]'],…
't','x','flag','beta','rho','sigma');
％ flag变量是不能省略的
7.2.4 微分方程转换
7.2.4.1 单个高阶常微分方程处理方法
? 例,
? 函数描述为,
function y=vdp_eq(t,x,flag,mu)
y=[x(2); -mu*(x(1)^2-1)*x(2)-x(1)];
>> x0=[-0.2; -0.7]; t_final=20;
>> mu=1; [t1,y1]=ode45('vdp_eq',[0,t_final],x0,[],mu);
>> mu=2; [t2,y2]=ode45('vdp_eq',[0,t_final],x0,[],mu);
>> plot(t1,y1,t2,y2,':')
>> figure; plot(y1(:,1),y1(:,2),y2(:,1),y2(:,2),':')
>> x0=[2;0]; t_final=3000;
>> mu=1000; [t,y]=ode45('vdp_eq',[0,t_final],x0,[],mu);
由于变步长所采用的步长过小，所需时间较长，导致输出的 y
矩阵过大，超出计算机存储空间容量。所以不适合采用
ode45()来求解，可用刚性方程求解算法 ode15s( )。
7.2.4.2 高阶常微分方程组的变换方法
? 例,
? 描述函数,
function dx=apolloeq(t,x)
mu=1/82.45; mu1=1-mu;
r1=sqrt((x(1)+mu)^2+x(3)^2);
r2=sqrt((x(1)-mu1)^2+x(3)^2);
dx=[x(2);
2*x(4)+x(1)-mu1*(x(1)+mu)/r1^3-mu*(x(1)-mu1)/r2^3;
x(4);
-2*x(2)+x(3)-mu1*x(3)/r1^3-mu*x(3)/r2^3];
? 求解,
>> x0=[1.2; 0; 0; -1.04935751];
>> tic,[t,y]=ode45('apolloeq',[0,20],x0); toc
elapsed_time =
0.8310
>> length(t),
>> plot(y(:,1),y(:,3))
ans =
689
得出的轨道不正确,
默认精度 RelTol设置
得太大，从而导致的
误差传递，可减小该
值。
? 改变精度,
>> options=odeset; options.RelTol=1e-6;
>> tic,[t1,y1]=ode45('apolloeq',[0,20],x0,options); toc
elapsed_time =
0.8110
>> length(t1),
>> plot(y1(:,1),y1(:,3)),
ans =
1873
>> min(diff(t1))
ans =
1.8927e-004
>> plot(t1(1:end-1),…
diff(t1))
? 例,
>> x0=[1.2; 0; 0; -1.04935751];
>> tic,[t1,y1]=rk_4('apolloeq',[0,20,0.01],x0); toc
elapsed_time =
4.2570
>> plot(y1(:,1),y1(:,3))
% 绘制出轨迹曲线
显而易见，这样求解
是错误的，应该采用
更小的步长。
>> tic,[t2,y2]=rk_4('apolloeq',[0,20,0.001],x0); toc
elapsed_time =
124.4990 ％计算时间过长
>> plot(y2(:,1),y2(:,3))
% 绘制出轨迹曲线
严格说来某些点仍不
满足 10－ 6的误差限,
所以求解常微分方程
组时建议采用变步长
算法，而不是定步长
算法。
? 例,
用 MATLAB符号工具箱求解,
令 ％
>> syms x1 x2 x3 x4
>> [dx,dy]=solve('dx+2*x4*x1=2*dy','dx*x4+ …
3*x2*dy+x1*x4-x3=5','dx,dy')
dx =
-2*(3*x4*x1*x2+x4*x1-x3-5)/(2*x4+3*x2)
dy =
(2*x4^2*x1-x4*x1+x3+5)/(2*x4+3*x2)
对于更复杂的问题来说，手工变换的难度将很
大，所以如有可能，可采用计算机去求解有关方程，
获得解析解。如不能得到解析解，也需要在描写一
阶常微分方程组时列写出式子，得出问题的数值解。
1 2 3 4,,,,,x x x x x y x y d x x d y y? ? ? ? ? ?& & && &&
7.3特殊微分方程的数值解
7.3.1 刚性微分方程的求解
? 刚性微分方程
一类特殊的常微分方程，其中一些解变
化缓慢，另一些变化快，且相差悬殊，这类
方程常常称为刚性方程。
MATLAB采用求解函数 ode15s()，该函数的调用
格式和 ode45()完全一致。
[t,x]=ode15s(Fun,[t0,tf],x0,options,p1,p2,… )
? 例,
％计算
>> h_opt=odeset; h_opt.RelTol=1e-6;
>> x0=[2;0]; t_final=3000;
>> tic,mu=1000;
[t,y]=ode15s('vdp_eq',[0,t_final],x0,h_opt,mu);
toc
elapsed_time =
2.5240
％作图
>> plot(t,y(:,1)); figure; plot(t,y(:,2))
y(:,1)曲线变化较平滑,y(:,2)变化在某些点上
较快。
? 例,
定义函数
function dy=c7exstf2(t,y)
dy=[0.04*(1-y(1))-(1-y(2))*y(1)+0.0001*(1-y(2))^2;
-10^4*y(1)+3000*(1-y(2))^2];
? 方法一
>> tic,[t2,y2]=ode45('c7exstf2',[0,100],[0;1]); toc
elapsed_time =
229.4700
>> length(t2),plot(t2,y2)
ans =
356941
? 步长分析,
>> format long,[min(diff(t2)),max(diff(t2))]
ans =
0.00022220693884 0.00214971787184
>> plot(t2(1:end-1),diff(t2))
? 方法二，用 ode15s()代替 ode45()
>> opt=odeset; opt.RelTol=1e-6;
>> tic,[t1,y1]=ode15s('c7exstf2',[0,100],[0;1],opt); toc
elapsed_time =
0.49100000000000
>> length(t1),
>> plot(t1,y1)
ans =
169
7.3.2 隐式微分方程求解
? 例,
? 编写函数,
function dx=c6ximp(t,x)
A=[sin(x(1)) cos(x(2)); -cos(x(2)) sin(x(1))];
B=[1-x(1); -x(2)]; dx=inv(A)*B;
求解,
>> opt=odeset; opt.RelTol=1e-6;
>> [t,x]=ode45('c7ximp',[0,10],[0; 0],opt); plot(t,x)
