第十章 常微分方程（组）求解

第三篇 第十章 常微分方程（组）求解

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．

第三篇 第十章 常微分方程（组）求解

的通解的主要调用格式如下： 调用格式一： S=dsolve ('eqn'，'var')

输入的量：eqn是常微分方程（10.1）改用符号方程表示的常微分方程F(x,y,Dy,D2y,...,Dny)?0。var表示自变量，默认的自变量为t。 输出的量：S是常微分方程（10.1）的通解。

?,y1??,...,y1(n1))?0,?F1(x,y1,y1?(n2)?,y2??,...,y2)?0,?F2(x,y2,y2求常微分方程组? （10.2）

?..............................?F(x,y,y?,y??,...,y(nm))?0,mmmm?m的通解的主要调用格式如下：

调用格式二： S=dsolve ('eqn1','eqn2', ... ,'eqnm'，'var')

输入的量：eqn1,eqn2, ... ,eqnm分别是常微分方程组（10.2）中用符号方程表示的m个常微分方程。默认的自变量为t，也可将自变量t变为其他的符号变量var。 输出的量：S是常微分方程组（10.2）的通解。 例1：求下列常微分方程的通解：

dyd2ydy(1)??at;(2)2?a?b(sinx?cosx);dtdxdxdyd2xadxbF(t)2(3)?y(1?y);(4)2??x?;

dxdtmdtmmd2y(5)2?siny.dx解：输入程序：

??y1=dsolve('Dy=-a*t'),

y2=dsolve('D2y+a*Dy=b*(sin(x)+cos(x))','x'), y3=dsolve('Dy=y^2*(1-y)'),

x4=dsolve('D2x+a*Dx/m+b*x/m=F(t)/m'), y5=dsolve('D2y=sin(y)'),pretty(y5) 例2：求下列常微分方程组的通解：

?df(t)?dy?af(t)?bg(t),?4z?sinx,???dt?dx(1)?(2)?

dg(t)dydz???af(t)?bg(t);??3?cosx.??dx?dt?dx解：输入程序：

??S1=dsolve('Df=a*f+b*g','Dg=-a*f+b*g'), f=S1.f,g=S1.g,de1='Dy+4*z=sin(x)'; de2='Dy+3*Dz=cos(x)';

S2=dsolve(de1,de2,'x'),Y=S2.y,z=S2.z

10.1.2 用MATLAB求常微分方程（组）的特解

如果给定常微分方程（10.1）的初始条件

y(x0)?a0,y?(x0)?a1,...,y(n)(x0)?an, （10.3）

则求方程（10.1）的特解的主要调用格式如下：

调用格式三： S=dsolve ('eqn'，'condition1'，…，'conditionn'，'var')

输入的量：eqn是常微分方程（10.1）改用符号方程表示的常微分方程F(x,y,Dy,D2y,...,Dny)?0；condition1,condition2,...,conditionn是初始条件（10.3）；var表示自变量，默认的自变量为t。 输出的量：S是常微分方程（10.1）的特解。

同理，如果给定常微分方程组（10.2）的初始条件，则求方程（10.2）的特解的主要调用格式如下：

175．

第三篇 第十章 常微分方程（组）求解

调用格式四：S=dsolve('eqn1','eqn2', ... ,'eqnm'，'condition1'，'condition2'，…，'var') 输入的量：eqn1,eqn2, ... ,eqnm分别是常微分方程组（10.2）中用符号方程表示的m个常微分方程；condition1,condition2,...是常微分方程组（10.2）的初始条件；默认的自变量为t，也可将自变量t变为其他的符号变量var。 输出的量：S是常微分方程组（10.2）的特解。 例3：

?dy?(1)???y2?1,y(0)?0;?dx?d3wdwd2w(2)3??w,w(0)?1,?0,2?0;

dtdtt?0dtt?0d2f(x)df(x)2df(x)(3)??a,f(0)?1,?0.2?dxdxdxx?a2解：输入程序：

??y=dsolve('(Dy)^2+y^2=1','y(0)=0')

w=dsolve('D3w=-w','w(0)=1,Dw(0)=0,D2w(0)=0') f=dsolve('D2f=-a^2*Df','f(0)=1,Df(pi/a)=0') 例4：

?df(t)?af(t)?bg(t),??dt(1)?f(0)?1,g(0)?2;dg(t)???af(t)?bg(t);??dt?dydydz?2?y??z?0,?dy?dx2dxdx(2)?y(0)?0,2dx?dy?y?dz?2dz?z?ex,?dx2dx?dx解：输入程序：

2

?2,z(5)?1,x?3dzdx?0.x?0??S1=dsolve('Df=a*f+b*g','Dg=-a*f+b*g','f(0)=1','g(0)=2'),

f=S1.f;fs=simple(f),g=S1.g;

gs=simple(g),de1='D2y+2*Dy+y+Dz+z=0'; de2='Dy+y+D2z+2*Dz+z=exp(x)';

S2=dsolve(de1,de2,'y(0)=0,Dy(3)=2,z(5)=1,Dz(0)=0') Y=S2.y;ys=simple(y),z=S2.z;zs=simple(z),

10.1.3 线性常微分方程组的解法

求解齐次线性常微分方程组Y?AX的MATLAB主程序

输入量：齐次线性常微分方程组Y=AX的矩阵A.

输出量：X是Y=AX的解，E是矩阵A的特征值，V是矩阵A的特征向量。

名为qcxxcwz.m的求解齐次线性常微分方程组Y?AX的MATLAB主程序如下：

function [X,E,V]=qcxxcwz(A) syms x E=eig(A);

[V,n]=eig(A); X=exp(E*x)'*V;

例1：用两种方法求下面常微分方程组的通解：

176．

第三篇 第十章 常微分方程（组）求解

?df(t)?2f(t)?3g(t),??dt ?dg(t)???2f(t)?3g(t).??dt解：输入程序： ??syms C1 C2

A=[2,3;-2,3];[X,E,V]=qcxxcwz(A),S=X*[C1;C2]

S1=dsolve('Df=2*f+3*g','Dg=-2*f+3*g','x'),f=S1.f,g=S1.g 习题一：

1， 求下列常微分方程的通解：

(1)d6yd27ydydt6??at?cos5t;(2)dx2?5dx?7x(sin3x?cos3x);(3)d2xdxF(t)d2ydy dt2??dt??x?m;(4)dx2?2dx?y?sin3x.2，求下列常微分方程组的通解：

?df(t)?5f(t)?dy(1)???7g(t),3?8z?sin5x,?dt?dg(t)(2)???dxdydz

??dt??5f(t)?7g(t);???dx?3dx?cos5x.3， 求下列常微分方程在给定初始条件下的特解：

2(1)??dy??dx???y2?cosx,y90)?0;d3wdwd2(2)wdt3??w,w(0)?1,dt?1,2?0;

t?0dtt?0(3)d2f(x)df(dx2?4x)dx,f(0)?1,df(x)dx0.x???24， 求下列常微分方程组满足所给初始条件下的特解：

?df(t)?2f(t)?3g(t),(1)???dtdg(t)f(0)?1,g(0)?2;???dt??2f(t)?3g(t);?d2y?2dy?y?dz?z?(2)???dx2dxdx0,dy?2,z(5)?1,dz?dy2y(1)?0,?y?dz?2dz?z?exdxx?3dx?3.x?2??dxdx2dx,5，用两种方法求下面常微分方程组的通解：

??df(t)??dt?5f(t)?7g(t),?dg(t) ??dt??5f(t)?7g(t).10.3 欧拉（Euler）方法的MATLAB程序

10.3.1 向前欧拉公式及其误差估计

177．

第三篇 第十章 常微分方程（组）求解

（一）欧拉方法

如果f(x,y)中的x取小区间?xn,xn?1?的左端点xn，则向前欧拉公式

?yn?1?yn?hf(xn,yn),n?0,1,2,..., ??xn?x0?nh.欧拉方法就是用差分方程初值问题

?yn?1?yn?hf(xn,yn)(n?0,1,2,...), (10.8) ??y0?y(x0)?dy??f(x,y),的解作为常微分方程初值问题?dx (10.5)

??y0?y(x0).的数值解，即将初值y0代入yn?1?yn?hf(xn,yn)(n?0,1,2,...)，逐次求出y1,y2,...

例10.3.1 用欧拉方法求初值问题

?dy??x?y,0?x?1, ?dx??yx?0?10.0075，并计算误差，画出精确解和数值解的图形. 的数值解，分别取h?0.0750，解 （1）先求精确解，输入程序：

??y=dsolve('(Dy)+y-x=0','y(0)=1','x') 对此问题的前欧拉公式具体形式为：

?yn+1?yn?h(xn?yn),n?0,1,2,..., ??xn?nh,(2)编写并保存名为Eulerli1.m的MATLAB计算和画图的主程序如下

function P=Eulerli1(x0,y0,b,h) n=(b-x0)/h; X=zeros(n,1); Y=zeros(n,1); k=1; X(k)=x0; Y(k)=y0; for k=1:n

X(k+1)=X(k)+h;

Y(k+1)=Y(k)+h*(X(k)-Y(k)); k=k+1; end

y=X-1+2*exp(-X); plot(X,Y,'mp',X,y,'b-') grid

xlabel('自变量 X'), ylabel('因变量 Y') title('用向前欧拉公式求dy/dx=x-y，y(0)=1在[0,1]上的数值解和精确

解y=x-1+2 exp(-x)')

legend('h=0.075时，dy/dx=x-y，y(0)=1在[0,1]上的数值解','精确

解y=x-1+2 exp(-x)')

jwY=y-Y;xwY=jwY./y; k1=1:n;k=[0,k1];

P=[k',X,Y,y,jwY,xwY]; 在MATLAB工作窗口输入下面的程序

>>x0=0;y0=1;b=1;h=0.0750; P=Eulerli1(x0,y0,b,h) 在MATLAB工作窗口输入下面的程序

>>h1=0.0075; P1=Eulerli1(x0,y0,b,h1)

legend('h1=0.0075时，dy/dx=x-y，y(0)=1在[0,1]上的数值解','精

确解y=x-1+2 exp(-x)')

178．

第三篇 第十章 常微分方程（组）求解

10.3.2 向前欧拉方法的三种MATLAB程序

(一) 向前欧拉公式的MATLAB主程序1

用向前欧拉公式（10.8）求解常微分方程初值问题（10.5）的数值解及其截断误差公式的MATLAB主程序1

输入量：funfcn是函数f(x,y)，x0和y0是初值问题y(x0)?y0,b是自变量x的最

大值，n是自变量x取值的个数，tol是精度。

输出量：数组X的元素是由自变量x的取值组成，数组Y的元素是利用向前欧拉公式(10.8)求出常微分方程初值问题(10.5)在X的元素处的数值解，REn是在每个子

h是步长，区间[xx,xk?1]上的局部截断误差公式，其中dy2是y(x)的2阶导数y??(x)。

function [h,k,X,Y,P,REn]=Qeuler1(funfcn,x0,y0,b,n,tol)

x=x0; h=(b-x)/n; X=zeros(n,1); y=y0; Y=zeros(n,1); k=1; X(k)=x; Y(k)=y'; for k=2:n+1

fxy=feval(funfcn,x,y); delta=norm(h*fxy,'inf');

wucha=tol*max(norm(y,'inf'),1.0); if delta>=wucha

x=x+h; y=y+h*fxy; X(k)=x;Y(k)=y'; end

plot(X,Y,'rp') grid

xlabel('自变量 X'), ylabel('因变量 Y')

title('用向前欧拉（Euler）公式计算dy/dx=f(x,y)，y(x0)=y0

在[x0,b]上的数值解')

end

P=[X,Y]; syms dy2,

REn=0.5*dy2*h^2;

例10.3.2 用向前欧拉公式（10.8）求解初值问题

dy??3y?8x?7,y(0)?1， dx分别取n?10,100，并将计算结果与精确解作比较，写出在每个子区间[xk,xk?1]上的局部截断误差公式，画出数值解与精确解在区间[0,1]上的图形.

解 (1)建立并保存名为funfcn.m的M文件函数。 function f=function(x,y) f=8*x-3*y-7;

(2)建立并保存名为Qeuler1.m的M文件函数。 (3)输入程序

>> S1=dsolve('Dy=8*x-3*y-7','y(0)=1','x') (4)输入程序

>> subplot(2,1,1)

x0=0;y0=1;b=1-1.e-4;

n=100;tol=1.e-4;

[h1,k1,x1,Y1,P1,Ren1]=QEuler1(@funfcn,x0,y0,b,n,tol) hold on

S1= 8/3*x1-29/9+38/9*exp(-3*x1), plot(x1,S1,'b-') title('用向前欧拉公式计算dy/dx=8x-3y-7，y(0)=1在[0,1]上的数值解

')

legend('n=100时，dy/dx=8x-3y-7，y(0)=1在[0,1]上的数值解','

dy/dx=8x-3y-7，y(0)=1在[0,1]上的精确解')

179．

第三篇 第十章 常微分方程（组）求解

hold off jdwuc1=S1-Y1; jwY1=S1-Y1;

xwY1=jwY1./S1;k1=1:n;k=[0,k1]; P1=[k',x1,Y1,S1,jwY1,xwY1] subplot(2,1,2) n1=10;

[h2,k2,x2,Y2,P2,Ren2]=QEuler1(@funfcn,x0,y0,b,n1,tol)

hold on

S1 = 8/3*x2-29/9+38/9*exp(-3*x2), plot(x2,S1,'b-')

legend('n=10时，dy/dx=8x-3y-7，y(0)=1在[0,1]上的数值解','

dy/dx=8x-3y-7，y(0)=1在[0,1]上的精确解')

hold off

jwY2=S1-Y2;xwY2=jwY2./S1;k1=1:n1;k=[0,k1]; P2=[k',x2,Y2,S1,jwY2,xwY2]

运行后屏幕显示分别取n?10,100时，所给的初值问题在[0,1]上的自变量x的值构成的数组Xi (i=1,2)，利用向前欧拉公式（10.8）求出的与Xi (i=1,2)对应的数值解Yi (i=1, 2)，Yi的绝对误差jwYi (i=1,2)和相对误差xwYi (i=1,2)（略），每个子区间[xk,xk?1]上的局部截断误差公式Ren2 =1/200*dy2， 近似解Y与精确解的图形.

例10.3.3 分别用向前欧拉公式(10.8)、一、四阶泰勒逼近及用Matlab函数double求解

?dy??1?ycosx,0?x?2,初值问题?dx分别取n?40,80时，并将计算结果与精确解做比

??y(0)?0.较，写出在每个子区间?xn,xn?1?上的局部截断误差公式，画出数值解与精确解在区间[0，

2]上的图形。

解：(1)建立并保存名为funfcn.m的M文件函数。 function f=function(x,y) f=1-y*cos(x); (2)输入程序

>> S1=dsolve('Dy=1-y*cos(x)','y(0)=0','x')

取一阶泰勒逼近

y?e?sinx(1?x?cosx).

同理，取四阶泰勒逼近，输入程序 >> syms X

y=exp(-sin(X))*int(1+sin(u)+((sin(u))^2)/2+((sin(u))^3)/6+((sin(u))^4)/24,u,0,X)

根据运行后，将屏幕显示结果整理得四阶泰勒逼近为

?108111?1017??y?e?sinx??x?cosx??sinx?sin2x?sin3x??.

1896?964???964（3）取n?40时，输入程序

>> syms u

h=2/40,X=0:h:2;

S=exp(-sin(X))*int(exp(sin(u)),u,0,X);S1=double(S); subplot(4,1,1)

plot(X,S1,'bo'),grid

legend('用double计算dy/dx=1-ycos(x),y(0)=0在[0,2]上的精确解的数值解')

subplot(4,1,2)

y=exp(-sin(X)).*(1+X-cos(X));plot(X,y,'g-'),grid

legend('y=exp(-sin(X)).*(1+X-cos(X))在[0,2]上的图形') subplot(4,1,3)

y1=exp(-sin(X)).*81/64*X-10/9*cos(X)-17/64*cos(X).*sin(X)-1/18*sin(X).^2.*cos(X)-1/96*sin(X).^3.*cos(X)+10/9);plot(X,y1,'m-'),grid

legend('四阶泰勒逼近的解在[0,2]上的图形') subplot(4,1,4)

180．

第三篇 第十章 常微分方程（组）求解

x0=0;y0=0;b=2-1.e-4;n=40;tol=1.e-4; [h,k,X,Y,p,Ren]=QEuler1(@funfcn,x0,y0,b,n,tol) legend('n=40时，用向前欧拉公式计算dy/dx=1-ycos(x),y(0)=0在[0,2]上的数值解')

syms u

X=0:0.05:2;S=exp(-sin(X))*int(exp(sin(u)),u,0,X); S1=double(S);jwS1y1=y1'-S1',

xwS1y1=jwS1y1./S1',jwYy1=y1'-Y,xwYy1=jwYy1./Y, k1=1:n;k=[0,k1];p=[k',X',Y,y1',jwYy1,xwYy1], pS=[k',X',S1',jwS1y1],xwS1y1

运行后屏幕显示取n=40时，分别用向前欧拉公式(10.8)和一、四阶泰勒多项式逼近法及double求解此初值问题在[0,2]上的自变量X对应的数值解Y，y，y1，S1及其图形，Y与y1，S1与y1的绝对误差向量jwYy1，xwS1y1和相对误差向量xwYy1，xwS1y1（略）.向前欧拉公式在每个子区间

?xk,xk?1?上的局部截断误差公式为

Ren?0.0012y??(?k),?k??xk,xk?1?,k?1,2,...,40。

取n?40时，输入程序 >> syms u

X=0:0.05:2;

S=exp(-sin(X))*int(exp(sin(u)),u,0,X);S1=double(S); y=exp(-sin(X)).*(1+X-cos(X));plot(X,y,'g-'),grid x0=0;y0=0;b=2-1.e-4;n=40;tol=1.e-4;

[h,k,X,Y,p,Ren]=QEuler1(@funfcn,x0,y0,b,n,tol)

y4=exp(-sin(X)).*81/64*X-10/9*cos(X)-17/64*cos(X).*sin(X)-1/18*sin(X).^2.*cos(X)-1/96*sin(X).^3.*cos(X)+10/9);plot(X,S1,'mo',X,y,'g-.',X,Y,'rp',X,y4,'b-'),grid

legend('用double计算dy/dx=1-ycos(x),y(0)=0在[0,2]上的精确解的数值解','y=exp(-sin(X)).*(1+X-cos(X))在[0,2]上的图形','n=40时，用向前欧拉公式计算dy/dx=1-ycos(x),y(0)=0在[0,2]上的数值解

','y4=exp(-sin(X))(10/9+81/64X-cos(X)(10/9+17/64sin(X)+1/18sin(X)^2+1/96sin(X)^3))在[0,2]上的图形')

运行后屏幕显示取n=40时，分别用向前欧拉公式(10.8)和一、四阶泰勒多项式逼近法及double求解此初值问题在[0,2]上的自变量X对应的数值解Y，y，y1，S1及其图形。 取n=40时，向前欧拉公式和四阶泰勒多项式逼近法所求初值问题的数值解的曲线几乎重合，二者的绝对误差和相对误差最小。一阶泰勒多项式逼近法所求初值问题的数值解的曲线与前两条曲线较近。而用函数double所求初值问题的数值解的曲线与前三条曲线很远，且绝对误差和相对误差很大，尤其是相对误差，当x=0时，Y，y1，和y与真值0相等，而S1=0.6912e-093不等于真值0，由此推测，用向前欧拉公式和一、四阶泰勒多项式逼近法所求初值问题的值都接近与真值，是有实际意义的，而用函数double计算的值严重失真，没有实际意义。 下面再取n=80验证。

（4）取n=80时，输入程序 >> syms u

X=0:0.0250:2;

S=exp(-sin(X))*int(exp(sin(u)),u,0,X);S1=double(S); subplot(4,1,1)

plot(X,S1,'bo'),grid

legend('用double计算dy/dx=1-ycos(x),y(0)=0在[0,2]上的精确解的数值解')

subplot(4,1,2)

y=exp(-sin(X)).*(1+X-cos(X));plot(X,y,'g-'),grid

legend('y=exp(-sin(X)).*(1+X-cos(X))在[0,2]上的图形') subplot(4,1,3)

y1=exp(-sin(X)).*81/64*X-10/9*cos(X)-17/64*cos(X).*sin(X)-1/18*sin(X).^2.*cos(X)-1/96*sin(X).^3.*cos(X)+10/9);plot(X,y1,'m-'),grid

legend('四阶泰勒逼近的解在[0,2]上的图形') subplot(4,1,4)

x0=0;y0=0;b=2-1.e-4;n=80;tol=1.e-4;

[h,k,X,Y,p,Ren]=QEuler1(@funfcn,x0,y0,b,n,tol)

181．

第三篇 第十章 常微分方程（组）求解

legend('n=80时，用向前欧拉公式计算dy/dx=1-ycos(x),y(0)=0在[0,2]上的数值解')

syms u

X=0:0.0250:2;S=exp(-sin(X))*int(exp(sin(u)),u,0,X); S1=double(S);jwS1y=S1'-y',

xwS1y=jwS1y./S1',jwYy=Y-y',xwYy=jwYy./y', k1=1:n;k=[0,k1];p=[k',X',Y,y',jwYy,xwYy], pS=[k',X',S1',jwS1y],xwS1y

运行后屏幕显示取n=80时，分别用向前欧拉公式(10.8)和一、四阶泰勒多项式逼近法及double求解此初值问题在[0,2]上的自变量X对应的数值解Y，y，y1，S1及其图形，Y与y，S1与y的绝对误差向量jwYy，xwS1y和相对误差向量xwYy，xwS1y（略）.向前欧拉公式在每个子区间

?xk,xk?1?上的局部截断误差公式为

Ren?0.0003y??(?k),?k??xk,xk?1?,k?1,2,...,80。

?xk,xk?1?上的局部截断误差公式为

由运行的结果可见，当x=0时，Y，y1，和y与真值0相等，但是S1却随着n的改变而发生变化，但是不等于真值0，由此推测，用向前欧拉公式和一、四阶泰勒多项式逼近法所求初值问题的值都接近与真值，是有实际意义的，而用函数double计算的值严重失真，没有实际意义。下面再取n=100时，向前欧拉公式在每个子区间

Ren?1.9998e?004y??(?k),?k??xk,xk?1?,k?1,2,...,100。

这说明，节点数n越大，向前欧拉公式在每个子区间?xk,xk?1?上的局部截断误差越小，所以函数double计算的值严重失真。

(二) 向前欧拉公式的MATLAB主程序2

用向前欧拉公式（10.8）求解常微分方程初值问题（10.5）的数值解及其截断误差公式的MATLAB主程序2

输入量：funfcn是函数f(x,y)，x0和y0是初值问题y(x0)?y0,b是自变量x的最

大值，h是步长。

输出量：k是自变量x取值的个数，向量X的元素是自变量x的取值，数组Y的元素是利用向前欧拉公式(10.8)求出常微分方程初值问题(10.5)在X的元素处的数值解，REn是在每个子区间[xx,xk?1]上的局部截断误差公式，h是步长，其中dy2是

y(x)的2阶导数y??(x)。

function [k,X,Y,P,REn]=Qeuler2(funfcn,x0,y0,b,h) x=x0; n=fix((b-x)/h); X=zeros(n+1,1); y=y0; Y=zeros(n+1,1); k=1; X(k)=x; Y(k)=y'; for k=2:n+1

X(k)=x+(k-1)*h,

fxy=feval(funfcn,x,y), Y(k)=y+h*fxy

y=Y(k), k=k+1,plot(X,Y,'rp'),

grid,xlabel('自变量 X'), ylabel('因变量 Y') title('用向前欧拉（Euler）公式计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]

上的数值解')

end

k1=1:n;k=[0,k1]; P=[k',X,Y]; syms dy2,

REn=0.5*dy2*h^2;

例10.3.4 用向前欧拉公式（10.8）求解初值问题

dyx?y?,y(0)?1，取dx3y 182．

第三篇 第十章 常微分方程（组）求解

n?10,100，并将计算结果与精确解作比较，写出在每个子区间[xk,xk?1]上的局部截断误差公式，画出数值解与精确解在区间[0,2]上的图形.

解 (1)建立并保存名为funfcn.m的M文件函数。 function f=function(x,y) f=y-x/(3*y);

(2)建立并保存名为Qeuler2.m的M文件函数。 (3)输入程序

>> S1=dsolve('Dy=y-x/(3*y)','y(0)=1','x') (4)输入程序

>> subplot(2,1,1)

x0=0;y0=1;b=2;n=10;h=2/10;

[k,X,Y,P,REn]=Qeuler2(@funfcn,x0,y0,b,h) hold on

S1=1/6*(6+12*X+30*exp(2*X)).^(1/2); plot(X,S1,'b-')

title('用向前欧拉公式计算dy/dx=y-x/(3y)，y(0)=1在[0,2]上的数值

解')

legend('n=10时，dy/dx=y-x/(3y)，y(0)=1在[0,2]上的数值解','

dy/dx =y-x/(3y)，y(0)=1在[0,2]上的精确解')

hold off

jdwucY=S1-Y; jwY=S1-Y;

xwY=jwY./Y;k1=1:n;k=[0,k1]; P1=[k',X,Y,S1,jwY,xwY] subplot(2,1,2) n1=100;h1=2/100;

[k,X1,Y1,P1,Ren1]=Qeuler2(@funfcn,x0,y0,b,h1) hold on

S2=1/6*(6+12*X1+30*exp(2*X1)).^(1/2); plot(X1,S2,'b-') legend('n=100时，dy/dx=y-x/(3y)，y(0)=1在[0,2]上的数值解','

dy/dx =y-x/(3y)，y(0)=1在[0,2]上的精确解')

hold off

jwY1=S2-Y1;xwY1=jwY1./Y1;k1=1:n1;k=[0,k1]; P2=[k',X1,Y1,S2,jwY1,xwY1]

运行后屏幕显示取n?10,100时，此初值问题在[0,2]上的自变量x的向量X，X1，利用向前欧拉公式（10.8）求出的与X，X1对应的数值解Y，Y1及其绝对误差向量jwY，jwY1和相对误差向量xwY，xwY1（略），每个子区间[xk,xk?1]上的局部截断误差公式Ren2 =1/200*dy2， 数值解Y与精确解的图形.

(三) 自适应向前欧拉公式的MATLAB主程序

用自适应向前欧拉公式求解常微分方程初值问题（10.5）的数值解的MATLAB主程序2

输入量：funfcn是函数f(x,y)，x0和y0是初值问题y(x0)?y0,b是自变量x的最大值，tol是精度。

输出量：向量X的元素是由自变量x的取值组成，数组Y的元素是利用向前欧拉公式(10.8)求出常微分方程初值问题(10.5)在X的元素处的数值解，向量H是步长h序列，n是自变量x取值的序号。并画出数值解向量Y的图形。

function [H,X,Y,k,h,P]=QEuler(funfcn,x0,b,y0,tol) %初始化. pow=1/3;

if nargin < 5 | isempty(tol), tol = 1.e-6; end;

if nargin < 6 | isempty(trace),

183．

第三篇 第十章 常微分方程（组）求解

trace = 0; end;

x=x0; h=0.0078125*(b-x); y=y0(:);p=128; H=zeros(p,1); X=zeros(p,1); Y=zeros(p,length(y)); k=1; X(k)=x;

Y(k,:)=y';

% 绘图.

clc,x,h,y end

% 主循环.

while (xx) if x+h>b

h=b-x; end

%计算斜率.

fxy=feval(funfcn,x,y); fxy=fxy(:); %计算误差，设定可接受误差.

delta=norm(h*fxy,'inf');

wucha=tol*max(norm(y,'inf'),1.0);

% 当误差可接受时重写解. if delta<=wucha

x=x+h; y=y+h*fxy; k=k+1; if k>length(X)

X=[X;zeros(p,1)]; Y=[Y;zeros(p,length(y))]; H=[H;zeros(p,1)];

end

H(k)=h;X(k)=x;Y(k,:)=y'; plot(X,Y,'rp'), grid xlabel('自变量 X'), ylabel('因变量 Y') title('用向前欧拉（Euler）公式计算dy/dx=f(x,y)，y(x0)=y0

在[x0,b]上的数值解')

end

%更新步长.

if delta~=0.0

h=min(h*8,0.9*h*(wucha/delta)^pow); end end

if (x

disp('Singularity likely.'), x end

H=H(1:k); X=X(1:k); Y=Y(1:k,:);

n=1:k; P=[n',H,X,Y]

例 10.3.5 用向前欧拉公式（10.8）求解在区间[0,2]上的初值问题

dy??3y?8x?7,y(0)?1，取精度10?1并与精确解作比较，并在同一个坐标系中作出dx图形。

解：(1)建立并保存名为funfcn.m的M文件函数。 function f=function(x,y) f=8*x-3*y-7;

(2)建立并保存名为Qeuler2.m的M文件函数。 (3)输入程序

>> S1=dsolve('Dy=8*x-3*y-7','y(0)=1','x') (4)输入程序

>> subplot(2,1,1)

>> x0=0;y0=1;b=2;tol=1.e-1;

[H,X,Y,k,h,P]=Qeuler2(@funfcn,x0,b,y0,tol)

184．

本文来源：https://www.bwwdw.com/article/at5t.html