第十章 常微分方程(组)求解

第三篇 第十章 常微分方程（组）求解 170． 一、 求微分方程的解 （一） 相关函数（命令）及简介 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之一。

（2）odefun是显式常微分方程：00

(,),().dyftydtyty







（3）在积分区间0[,]ftspantt上，从0t到ft，用初始条件0y求解。 （4）要获得问题在其他指定时间点012,,,...ttt上的解，则令012[,,,...,]ftspantttt（要求是单调的） （5）因为没有一种算法可以有效地解决所有的ODE问题，为此，Matlab提供了多种求解器solver，对于不同的ODE问题，采用不同的solver。

求解器solver ODE类型 特点 说明 ode45 非刚性 单步算法：4，5阶的Runge-Kutta 大部分场合的首选算法 方程；累计截断误差达3()x

ode23 非刚性 单步算法：2，3阶的Runge-Kutta 使用于精度较低的情形

方程；累计截断误差达3()x

Matlab常微分方程（组）求解 第三篇 第十章 常微分方程（组）求解 171． ode113 非刚性 多步法：Adams算法；高低精度均 计算时间比ode45短 可到3610~10。 ode23t 适度刚性 采用梯形算法 适度刚性情形 ode15s 刚性 多步法：Gears反向数值微分 若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).

例1：求解微分方程22xdyxyxedx，并加以验证。

求解本问题的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：求解微分方程0xxyye在初始条件(1)2ye下的特解，并画出解函数的图形。 求解本问题的Matlab程序： syms x y y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') ezplot(y) 第三篇 第十章 常微分方程（组）求解 172． 微分方程的特解为y=1/x*exp(x)+1/x*exp(1)(Matlab格式)，即xeeyx，解函数图形（略）。

例3：求解微分方程组5,30tdxxyedtdyxydt在初始条件01tx，

00ty

下的特解，并画出解函数的图形。 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等求解非刚性的标准形式的一阶常微分方程（组）的初值问题的数值解（近似解）

例4：求解微分方程初值问题2222,(0)1dyyxxdxy







的数值解，求解范围

为[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折线法求解微分方程初值问题2

2,(0)1dyxydxyy





的数

值解（步长h取0.4），求解范围为[0，2]。 解：本问题的差分方程为

001210,1,0.4,,(,)((,)2/),0,1,2,...,1kkkkkkxyhxxhyyhfxyfxyyxykn





其中：

Matlab程序： clear f=sym('y+2*x/y^2'); a=0; 第三篇 第十章 常微分方程（组）求解 173． 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，求微分方程2(1)2sin0xyxyx的通解。 2，求微分方程25sinxyyyex的通解。

3，求微分方程组0,0dxxydtdyxydt







在初值条件001,0ttxy下的特解，并画

出函数()yfx的图形。 4，分别用23,45odeode求上述第3题中的微分方程初值问题的数值解（近似解），求解区间为[0,2]t。利用画图来比较两种求解器之间的差异。

5，用Euler折线法求解微分方程初值问题2312,(0)1x

yyyy





的数值

解（步长h取0.1），求解范围为[0，2]。 6，用四阶RungeKutta法求解微分方程初值问题cos,(0)1xyyexy的数值解（步长h取0.1），求解范围为[0，3]。 四阶RungeKutta法的迭代公式为（Euler法实为一阶RungeKutta第三篇 第十章 常微分方程（组）求解

174． 法）：001112341213243(),,(22),6(,),0,1,2,...,1(,),22(,),22(,),kkkkkkkkkkkkyyxxxhhyyLLLLLfxyknhhLfxyLhhLfxyLLfxhyhL 试用该方法求解第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求常微分方程：()(,,,,...,)0nFxyyyy （10.1）