第三篇 第十章 常微分方程（组）求解
Matlab常微分方程（组）求解 一、 求微分方程的解
（一） 相关函数（命令）及简介
1， dsolve('equ1','equ2',…):Matlab求微分方程的解析解。
equ1,equ2,…为方程（或条件）。写方程（或条件）时用Dy表示y关于自变量的一阶导数，用D2y表示y关于自变量的二阶导数，依次类推。
2， simplify(s):对表达式s使用maple的化简规则进行化简。 例如： syms x
simplify(sin(x)^2+cos(x)^2) ans=1
3,[r,how]=simple(s):由于Matlab提供了多种化简规则，simple命令就是对表达式s用各种规则进行化简，然后用r返回最简形式，how返回形成这种形式所用的规则。 例如： syms x
[r,how]=simple(cos(x)^2-sin(x)^2) r=cos(2*x) how=combine
4,[T,Y]=solver(odefun,tspan,y0),求微分方程的数值解。 （1）其中的solver为命令
ode45,ode23,ode113,ode15s,ode23s,ode23t,ode23tb之一。
?dy?f(t,y),（2）odefun是显式常微分方程：? ?dt??y(t0)?y0.（3）在积分区间tspan?[t0,tf]上，从t0到tf，用初始条件y0求解。
（4）要获得问题在其他指定时间点t0,t1,t2,...上的解，则令tspan?[t0,t1,t2,...,tf]（要求是单调的）
（5）因为没有一种算法可以有效地解决所有的ODE问题，为此，Matlab提供了多种求解器solver，对于不同的ODE问题，采用不同的solver。
求解器solver ODE类型 特点 说明
ode45 非刚性 单步算法：4，5阶的Runge-Kutta 大部分场合的首选算法 方程；累计截断误差达(?x)3
ode23 非刚性 单步算法：2，3阶的Runge-Kutta 使用于精度较低的情形
3 方程；累计截断误差达(?x)
170．
第三篇 第十章 常微分方程（组）求解
ode113 非刚性 多步法：Adams算法；高低精度均 计算时间比ode45短 可到10~10。
ode23t 适度刚性 采用梯形算法 适度刚性情形
ode15s 刚性 多步法：Gear?s反向数值微分 若ode45失败时，可尝
精度中等 试使用
ode23s 刚性 单步法：2阶Rosebrock算法； 当精度较低时， 计算
低精度 时间比ode15s短
ode23tb 刚性 梯形算法；低精度 当精度较低时， 计算
时间比ode15s短
（6）要特别说明的是：ode23,ode45是极其常用的用来求解非刚性的标准形式的一阶常
微分方程（组）的初值问题的解的Matlab的常用程序，其中：
ode23采用龙格-库塔2阶算法，用3阶公式作误差估计来调节步长，具有低等的精度。 ode45则采用龙格-库塔4阶算法，用5阶公式作误差估计来调节步长，具有中等的精度。 5，ezplot(x,y,[tmin,tmax]):符号函数的作图命令，x,y为关于参数t的符号函数，[tmin,tmax]
为t的取值范围。
6，inline():建立一个内联函数。格式：inline('expr','var1','var2',…)，注意括号里的表达式要
加引号。
如：Q=dblquad(inline('y*sin(x)'),pi,2*pi,0,pi).
?3?6例1：求解微分方程dy?2xy?xe?x，并加以验证。
2dx求解本问题的Matlab程序： syms x y
y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') simplify(diff(y,x)+2*x*y-x*exp(-x^2)) 说明：
（1） 行1是用命令定义x，y为符号变量。这里可以
不写，但为确保正确性，建议写上；
（2） 行2是用命令求出的微分方程的解：
1/2*exp(-x^2)*x^2+exp(-x^2)*C1
（3） 行3使用所求得的解。这里是将解代入原微分
方程，结果应该为0，但这里给出：-x^3*exp(-x^2)-2*x*exp(-x^2)*C1+2*x*(1/2*exp(-x^2)*x^2+exp(-x^2)*C1)
（4） 行4用simplify()函数对上式进行化简，结果为
0，表明y=y(x)的确是微分方程的解。
例2：求解微分方程xy??y?ex?0在初始条件y(1)?2e下的特解，并画出解函数的图形。
求解本问题的Matlab程序： syms x y
y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') ezplot(y)
171．
第三篇 第十章 常微分方程（组）求解
微分方程的特解为y=1/x*exp(x)+1/x*exp(1)(Matlab格式)，
e?ex即y?，解函数图形（略）。
x?dxt?5x?y?e,??dt例3：求解微分方程组?在初始条件xt?0?1，yt?0?0dy??x?3y?0??dt下的特解，并画出解函数的图形。
Matlab程序： syms x y
[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t')
simple(x); simple(y);
ezplot(x,y,[0,1.3]);axis auto 2，用ode23,ode45等求解非刚性的标准形式的一阶常微分方程（组）的初值问题的数值解（近似解） 例
?dy2???2y?2x?2x,4：求解微分方程初值问题?dx的数值解，求解范围
??y(0)?1为[0，0.5]。
Matlab程序：
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-') 例5：用Euler
2x?dy??y?2,折线法求解微分方程初值问题?dxy的数
?y(0)?1?值解（步长h取0.4），求解范围为[0，2]。 解：本问题的差分方程为
?x0?0,y0?1,h?0.4,? ?xk?1?xk?h,?2?yk?1?yk?hf(xk,yk)(其中：f(x,y)?y?2x/y),k?0,1,2,...,n?1Matlab程序： clear
f=sym('y+2*x/y^2'); a=0;
172．
第三篇 第十章 常微分方程（组）求解
b=2; h=0.4;
n=(b-a)/h+1; x=0; y=1;
szj=[x,y]; for i=1:n-1
y=y+h*subs(f,{'x','y'},{x,y}); x=x+h;
szj=[szj;x,y]; end szj
plot(szj(:,1)szj(:,2)) 练习：
1，求微分方程(x2?1)y??2xy?sinx?0的通解。 2，求微分方程y???2y??5y?exsinx的通解。
?dx?x?y?0,?dt3，求微分方程组?在初值条件xt?0?1,yt?0?0下的特解，并画??dy?x?y?0??dt出函数y?f(x)的图形。
4，分别用ode23,ode45求上述第3题中的微分方程初值问题的数值解（近似解），求解区间为t?[0,2]。利用画图来比较两种求解器之间的差异。
5，用Euler
?12x2?y??y?3,折线法求解微分方程初值问题?y的数值
?y(0)?1?解（步长h取0.1），求解范围为[0，2]。
?y??y?excosx,6，用四阶Runge?Kutta法求解微分方程初值问题?的数
?y(0)?1值解（步长h取0.1），求解范围为[0，3]。
四阶Runge?Kutta法的迭代公式为（Euler
法实为一阶Runge?Kutta 173．
第三篇 第十章 常微分方程（组）求解
??y0?y(x0),??xk?1?xk?h,?h?yk?1?yk?(L1?2L2?2L3?L4),6??k?0,1,2,...,n?1 法）：?L1?f(xk,yk),?hh?L2?f(xk?,yk?L1),22??hh?L3?f(xk?,yk?L2),22???L4?f(xk?h,yk?hL3),试用该方法求解第5题中的初值问题。
Matlab程序： clear
f=sym('y-exp(x)*cos(x)'); a=0; b=3; h=0.1;
n=(b-a)/h+1; x=0; y=1;
szj=[x,y]; for i=1:n-1
l1=subs(f,{'x','y'},{x,y});
l2=subs(f,{'x','y'},{x+h/2,y+l1*h/2}); l3=subs(f,{'x','y'},{x+h/2,y+l2*h/2}); l4=subs(f,{'x','y'},{x+h,y+l3*h}); y=y+h*(l1+2*12+2*13+14)/6; x=x+h;
szj=[szj;x,y]; end szj
plot(szj(:,1)szj(:,2))
7,用ode45方法求上述第6题的常微分方程初值问题的数值解（近似解），从而利用画图来比较两者的差异。
二、常微分方程（组）的MATLAB符号求解
10.1.1 用MATLAB求常微分方程（组）的通解
用Matbal函数dsolve求常微分方程：F(x,y,y?,y??,...,y(n))?0 （10.1）
174．