微分方程数值解法

《微分方程数值解法》

【摘要】自然界与工程技术中的很多现象，可以归结为微分方程定解问题。其中，常微分方程求解是微分方程的重要基础内容。但是，对于许多的微分方程，往往很难得到甚至不存在精确的解析表达式，这时候，数值解提供了一个很好的解决思路。，针对于此，本文对常微分方程数值解法进行了简单研究，主要讨论了一些常用的数值解法，如欧拉法、改进的欧拉法、Runge —Kutta 方法、Adams 预估校正法以及勒让德谱方法等，通过具体的算例，结合MA TLAB 求解画图，初步给出了一般常微分方程数值解法的求解过程。同时，通过对各种方法的误差分析，让大家对各种方法的特点和适用范围有一个直观的感受。 【关键词】 常微分方程 数值解法 MA TLAB 误差分析

引言

在我国高校，《微分方程数值解法》作为对数学基础知识要求较高且应用非常广泛的一门课程，不仅

在数学专业，其他的理工科专业的本科及研究生教育中开设这门课程．近四十年来，《微分方程数值解法》不论在理论上还是在方法上都获得了很大的发展．同时，由于微分方程是描述物理、化学和生物现象的数学模型基础，且它的一些最新应用已经扩展到经济、金融预测、图像处理及其他领域 在实际应用中，通过相应的微分方程模型解决具体问题，采用数值方法求得方程的近似解，使具体问题迎刃而解。

2 欧拉法和改进的欧拉法

2.1 欧拉法

2.1.1 欧拉法介绍

首先，我们考虑如下的一阶常微分方程初值问题 ???==0

0)()

,('y x y y x f y

（2--1）

事实上，对于更复杂的常微分方程组或者高阶常微分方程，只需要将x 看做向量，（2--1）就成了一个一阶常微分方程组，而高阶常微分方程也可以通过降阶化成一个一阶常微分方程组。

欧拉方法是解常微分方程初值问题最简单最古老的一种数值方法，其基本思路就是把（2--1）中的导数项'y 用差商逼近，从而将一个微分方程转化为一个代数方程，以便求解。

设在[]b a ,中取等距节点h ，因为在节点n x 点上，由（2--1）可得：

))(()(',n n n x y x f x y =，

（2--2）

又由差商的定义可得：

h

x y x y x y n n n )

()(('1-≈

+）

（2--3） 所以有 ))(,()()(1n n n n x y x hf x y x y +≈+ （2--4）

用)(k x y 的近似值k y )1,(+=n n k 代入（2--4），则有计算1+n y 的欧拉公式

))(,(1n n n n x y x hf y y +=+ （2--5）

2.1.2欧拉法误差分析

从欧拉公式中可以看出，右端的n y 都是近似的，所以用它计算出来的

1+n y 会有累计误

差，累计误差比较复杂，为简化分析，我们考虑局部截断误差，即认为n y 是精确的前提下来估计11)(++-n n y x y ，记为1+n ε，泰勒展开有

)()()(''2

)(')()(132

1++<<+++=n n n n n x x h O y h x hy x y x y ξξ

（2--6）

联立（2--5），（2--6）即得1+n ε=)(''2

2

ξy h +)(3h O =)(2h O ，根据数值算法精度的定义，如果一个数值方法的局部截断误差1+n ε=)(1

+p h

O 则称这个算法具有P 阶精度，所以，欧拉

方法具有一阶精度或者称欧拉方法为一阶方法。 2.2 改进的欧拉方法

2.2.1 改进的欧拉法介绍

用数值积分离散化问题（1），两边做积分有： dx x y x f x y x y n n

x x n n ?

+=-+1

))(,()()(1

（2--7）

对右端积分使用梯形积分公式可得：

[]))(,())(,(2

))(,(111

+++≈

?

+n n n n x x x y x f x y x f h

dx x y x f n n

（2--8）

同欧拉方法，用)(k x y 的近似值k y )1,(+=n n k 代入（2--7），联立（2--8）得到改进的欧拉方法计算公式： ),(),((2

111+++++=n n n n n n y x f y x f h

y y （2--9）

一般来说，如果求解1+n y 时，算法右端不包含1+n y ，称为显性计算公式，如果包含，则求解时还需要解方程，这种称为隐式计算公式。显然公式（2--9）是一个隐式计算公式，事实上，改进的欧拉方法是用欧拉方法先求一个预测值1+n y ，再用这个预测值来计算1+n y ，即：

??

?

??++=+=++++)

,(),((2),(1111n n n n n n n n n n y x f y x f h

y y y x hf y y （2--10）

2.2.2改进的欧拉法误差分析

同欧拉法误差分析类似，用泰勒展开容易知道改进的欧拉方法具有二阶精度，证明略。 2.3 算例

2.3.1（一阶常微分方程）

求解初值问题

??

?

??∈-==]

1,0[21)0('

x y x y y y

解析：在MATLAB 中求解这个方程

y=dsolve('Dy=y-2*x/y','y(0)=1','x') 得y =(2*x+1)^(1/2)

它的解析解为x y 21+=，下面我们分别用欧拉方法和改进的欧拉方法来求其数值解。 欧拉方法：创建M 文件euler1.m,内容如下： function [x,y]=euler1(fun,x0,xfinal,y0,n) if nargin<5,n=50; end

h=(xfinal-x0)/n; x(1)=x0;y(1)=y0;

for i=1:n

x(i+1)=x(i)+h;

y(i+1)=y(i)+h*feval(fun,x(i),y(i));

End

再定义函数方程组中的函数f1，创建f1.m文件，内容如下：

function f=f1(x,y)

f=y-2*x/y

在MATLAB中输入：[x,y]=euler1('f1',0,1,1,20)输出f,[x,y]的值，将数值解跟精确解画图表示，输入：plot(x,y,'r*-',x,sqrt(1+2*x),'g+--')；xlabel('x');ylabel('y');title('y‘=y-2x/y');legend('数值解','精确解')得到图形，保存为euler1.fig，图形如下：

改进的欧拉方法：创建M文件eulerprove1.m,内容如下：

function [x,y]=eulerprove1(fun,x0,xfinal,y0,n)

if nargin<5,n=50;

end

h=(xfinal-x0)/n;

y(1)=y0,x(1)=x0;

for i=1:n

x(i+1)=x(i)+h;

y1=y(i)+h*feval(fun,x(i),y(i))/2;

y2=h*feval(fun,x(i+1),y1)/2;

y(i+1)=y1+y2

End

在MATLAB的command窗口输入: [x,y1]=eulerprove1('f1',0,1,1,20)返回f,x1,y1的值，作图，输入：plot(x,y1,'r*-',x,sqrt(1+2*x),'g+--')；xlabel('x');ylabel('y');title('y‘=y-2x/y');legend('数值解','精确解')，将图片保存为eulerprove1.fig,图形如下：

为了便于比较两种方法的误差，将两者的误差作到同一个图上，继续输入：

plot(x,abs(y-sqrt(1+2*x)),'y*-',x,abs(y1-sqrt(1+2*x)),'g+--')；

xlabel('x');ylabel('y');title('误差曲线');legend('欧拉方法','改进的欧拉方法')将图片保存为error1.fig，图形如下：

从该图形来看，改进的欧拉方法与欧拉方法误差接近，欧拉方法误差稍微大些，将x 的取值扩宽，n 取值增大时，可以发现改进的欧拉方法相比欧拉方法有更高的精度。

2.3.2（高阶微分方程）

对于二阶常微分方程）

（]1,0[1)0(2)0('''∈?

??==-=x x x x

x ，求数值解 解析：先算出其解析解，在MATLAB 中输入：y=dsolve('D2x=-x','x(0)=1','Dx(0)=2')得到解为：

Y=2*sin(t)+cos(t)，前面已经分别给出过欧拉方法和改进的欧拉方法的算例跟误差比较，这里我们就用精度更高的改进欧拉法进行数值求解。

改进的欧拉方法：先换元，令y x ='，则原方程可以转化为])1,0[(2)0(1)0(''∈??

?

??===-=x y x y

x x

y ，

现在，二阶常微分方程转化为了一个一阶常微分方程组，同2.3.2的方法，建立M 文件eulerprove3.m,内容如下：

function [t,x,y]=eulerprove2(t0,tfinal,x0,y0,n) f1=inline('y'); f2=inline('-x') if nargin<5,n=50; end

h=(tfinal-t0)/n; t(1)=t0,x(1)=x0;y(1)=y0; for i=1:n t(i+1)=t(i)+h;

x1=x(i)+h*feval(f1,y(i)); y1=y(i)+h*feval(f2,x(i)); x2=x(i)+h*feval(f1,y1); y2=y(i)+h*feval(f2,x1); x(i+1)=(x1+x2)/2; y(i+1)=(y1+y2)/2;

end

在command窗口输入[t,x,y]=eulerprove3(0,1,1,2,10)，得到t,x,y的值，其中x就是我们要求的数值解，作图,输入：plot(t,x,'r*-',t,2*sin(t)+cos(t),'b+-')；

xlabel('x');ylabel('y‘）；legend('数值解','精确解')，将图形保存为eulerprove3.fig,图形如下：

上面，我们已经通过例子看出，改进的欧拉法相比于欧拉法，在每一个节点处的误差值更下，下面，我们来讨论节点的多少（步长大小）对误差的影响，创建erro3r.m文件，内容如下：function [N,Y]=error3(n0,nfinal)

N(1)=n0,m=fix((nfinal-n0)/4);

for i=1:m

N(i+1)=N(i)+4;

[t,x1,y1]=eulerprove3(0,1,1,2,N(i));

Y(i)=log10(max(x1-2*sin(t)-cos(t)));

[t,x2,y1]=eulerprove3(0,1,1,2,N(i));

Y(i+1)=log10(max(x2-2*sin(t)-cos(t)));

end

输入[N,Y]=error3(4,100)，返回节点个数值和Y值，Y代表在N个节点时，数值解与精确解差的绝对值的最大值的对数（10为

底）。:plot(N,Y,'d-.'),xlabel('N'),ylabel('log_1_0max|error|) title('误差曲线'),将图片保存为error3.fig,图形如下:

从这个图的误差曲线，随着插值节点个数增多，误差也是越来越小的，所以，为了得到尽可能准确的解，使用改进的欧拉法时插值节点数应当越大越好。但是，随着节点个数的增多，计算量也会随之增大，所以，在实际应用时，应当结合具体要求灵活选取。

3龙格-库塔法

3.1 龙格-库塔法与泰勒展开

3.1.1泰勒（Taylor ）展开 首先，考虑考虑微分方程

),('y x f y = （3--1） 引入算子符号][

y

f x ??

+??， y x ff f f y

f x +=??

+??][ （3--2） 于是有

f

y

f x y f y f x f y f x y f x y f y

f x ff f y f f y m m y x y x 1)(2][,][)]]([['''][

'''-??

+??=??

+??=??+????+??=??+??=+=+=

（3--3）

设不带括号的f 及其各阶微商都在))(,(n n x y x 处取值，于是有泰勒展开式：

)

(]435333[!

4]

2[!

3][!2)()(][!4][!3][!2)()(53

22324

22325342321h O ff f f f f f f ff f ff f f f f f f f f ff f h ff f f f f ff f h ff f h hf x y h O f y

f x h f

y

f x h f y f x h hf x y x y y y x yy y yy x xy

y xy x xx y yyy xyy xxy xxx y y x yy xy xx y x n n n ++++++++++++++++++++=+??

+??+??

+??+??+??++=+ （3--4）

3.1.2 龙格-库塔（R-K ）方法介绍

龙格-库塔方法的基本思路是想办法计算),(y x f 在某些点上的函数值，然后对这些函数值做数值线性组合，构造出一个近似的计算公式；再把近似的计算公式和解的泰勒展开式相比较，使得前面的若干项相吻合，从而达到较高的精度，一般的显示R-K 方法的形式如下：

????

?

????

≤≤++==+=∑∑-==+)2(),,(),,(,

1

1111r i k y h x hf k y x hf k k y y i j j ij n n i n n r i i i n n βαω

（3--5）

当式（3--5）的局部截断误差达到方法段阶）为时称公式（K R r p h

O p ---+53)(1

。

龙格库塔方法是应用最广的求解常微分方程初值问题的单步法，下面我们给出几个推导公式。

3.2 龙格-库塔法公式与ode 函数

3.2.1 二阶二段R-K 法公式推导

由式子（3--5），二阶二段的R-K 方法可以写成

???

??++==++=+),

,(),,(1212

122111k y h x hf k y x hf k k k y y n n n n n n βαωω

（3--6）

不带括号的f 及其各阶微商都在))(,(n n x y x 处取值，由泰勒展开式（3--4）及截断误差的定义可知：

)

()2

1

()21()1()

()()()(2

)()))

(,()(,

())(,()()(32212222213212212

2122111h O ff h f h hf h O hff hf f h hf x y ff f h hf x y x y x hf x y h x hf x y x hf x y x y y x y x n y x n n n n n n n n n n +--+-+--=+++---+++=++---=++βωαωωωβαωωβαωωε （3--7）

由于上式是二阶二段的R-K 法，即必须有)(31h O n =+ε，所以

???

??===+，，，

2/12/1121

22221βωαωωω

（3--8）

式（3--8）有四个未知元，三个方程，故有无穷组解。

若令12=α，由式子（3--9）可解得，

，1,2/12121===βωω于是 ???

?

???

++==++=+)

,(),(,21211

21211k y h x hf k y x hf k k k y y n n n n n n

这个公式称为预估-校正公式。常用的二阶公式还有中点公式（取2/12=α），这里不再列举。

3.2.2 其他常见公式和ode 函数

从二阶二段R-K 公式的推导可以看出，二段方法每一步需要计算两次函数值，而它也只能达到2阶精度，如果我们提高函数的计算次数，就可以得到精度更高的计算公式，高阶的

R-K 法的推导与二阶方法的推导完全相同，只是随着阶数的增高计算量会逐渐增大。下面给出几个常用的计算公式。

（1）库塔公式（3阶3段）

?????????

+-+=++==+++=+)

2,()2/,2/()

,()4(6

12131213211

k k y h x hf k k y h x hf k y x hf k k k k y y n n n n n n n n

（2）标准四阶R-K 公式

?????

?

??

???++=++=++==++++=+),()

2/,2/()2/,2/(),()22(6

1342

3

12143211

k y h x hf k k y h x hf k k y h x hf k y x hf k k k k k y y n n n n n n n n n n

上述三阶公式具有三阶精度，四阶公式具有四阶精度，要验证这一点，我们观察如下泰勒展开：

...

)33(!

31

)2(!21)

,()()!1(1),()(!

1),(32232210++++++++++=??+??++??+??=+++=∑yyy xyy xxy xxx

yy xy xx y x n i n

i f k f hk kf h f h f k hkf f h kf hf f y x f y k x h n y x f y k x h i h y h x f (3--9)

其中y x ,落在平面上连接),(),(k y h x y x ++和的线段上。

只需要按照式（3.9）将k 泰勒展开，代入截断误差的计算公式中即可证明。(详细证明见参考文献[2]P253)

由于R-K 法是基于泰勒展开的数值方法，所以，如果解具有较好的光滑性，则能够得到较好的计算精度，反之，则可能四阶R-K 法的数值精度还不如低阶的数值方法，这个需要根据具体情况自行选择数值方法。一般而言，R-K 法在求解范围较大而精度要求较高时是比较好的方法，四阶的R-K 法也可以用作阿达姆斯预估校正法的前几步计算，以保持四阶精度（下一节阿达姆斯预估校正法中也提到了这一点）。

在MATLAB 中，专门提供了求解微分方程的ode 函数，如ode45，ode23，ode113，ode15s ，ode23s 等等。ode 求解器分为变步长和固定步长两种求解模式，其中，ode45是采用R-K 法的变步长求解器，和它一样的还有ode23。目前，ode45是使用最多的求解函数，主要用于求解非刚性常微分方程。（如果求解时长时间没有响应，则方程是刚性的，可以换用ode23求解），ode 函数的调用方式大都类似，以ode45为例，常用的语法格式有：[T,Y] = ode45(odefun,tspan,y0) ； [T,Y] = ode45(odefun,tspan,y0,options) ；sol =

ode45(odefun,[t0tf],y0...) ，其中odefun 是函数句柄，可以是函数文件名或内联函数名，tspan 是求解区间 [t0 tf] 或者一系列散点[t0,t1,...,tf]，y0 是初始值向量，T 返回列向量的时间点，Y 返回对应T 的求解列向量 ，options 是求解参数设置,可以用odeset 在计算前设定误差,输出参数,事件等 。 3.3算例

用自带ODE45求解高阶常微分方程组

将常微分方程组初值问题?????????=====++=++--0

.1)0(',25.0)0(75.0)0(',5.0)0(0))(''0)(''23

2

22

3

22y y x x y x y y y x x x 转化为一阶常微分方程组并用

ode45求解).1(),1(y x

解析：首先我们尝试用MATLAB 求出该方程组的解析解，在command 窗口输入：

[x,y]=dsolve('D2x+x*(x^2+y^2)^(-3/2)=0','D2y+y*(x^2+y^2)^(-3/2)=0','x(0)=0.5','y(0)=0.25','Dx(0)=0.75','Dy(0)=1')，运行大概五分钟，提示：Warning: Explicit solution could not be found.x =[ empty sym ]，y =[]，产生该问题的原因可能在于该非线性方程无法用solve 求解，或该方程的解析解不存在。

下面，我们用数值方法求解该方程组，调用的是函数库中的ode45函数。首先，我们可以做变量替换，令b y a x ==','，则该方程组可以化为一个一阶常微分方程组初值问题：

?????????=====+-==+-=--25

.0)0('5.0)0('0.1)0()

)('75.0)0()('23

2

22

3

22y b

y x a x b y x y b a y x x a 现在，可以直接调用ode45函数进行数值计算。先建立一个名为fun45.m 的M 文件来定

义方程组，内容如下：

function y=fun45(t,x)

y(1)=-x(3)*(x(3)^2+x(4)^2)^(-3/2);

y(2)=-x(4)*(x(3)^2+x(4)^2)^(-3/2);

y(3)=x(1);

y(4)=x(2);

y=y(:)

在command窗口输入： x0=[0.75,1,0.5,0.25];[t,x]=ode45('fun45',[0,1],x0),输出x 的值，x是一个四列的矩阵，四列分别代表a,b,x,y在区间[0,1]上的值，可以读出：

x(1)=0.5499, y(1)=0.7407，或者作图：plot(t,x)；legend('a','b','x','y')，xlabel('t')将图片保存为fun45.fig.图形如下：

通过Data Cursor一样可以读出x(1)=0.5499, y(1)=0.7407。

结语

掌握并灵活运用上述数值解法对于求常微分方程模型中的数值解有着重要的意义。本文通过对常微分方程的数值解法进行比较分析其中的优缺点，并使用实例证明。数值求解常微分方程有着重要意义，我们应该在不断探索其解法的前提下，创新出更合适软件来解决实际问题。有待于进一步的研究。

参考文献

[1] 李荣华,刘播. 微分方程数值解法[M]. 第四版. 北京：高等教育出版社.,2009.1