常微分方程初值问题数值解法

桂林理工大学理学院2007级信息与计算科学专业·生产实习 张玉娟：常微分方程初值问题数值解法 7/6/2010 1 常微分方程初值问题数值求解 学生：张玉娟 学号：3070942232 指导老师：梁鹏

摘 要：数值求解是指在没有办法知道未知函数的解析表达式的情况下，近似计算未知函数在其定义域中的某些离散点上的函数值。本文利用欧拉方法、梯形方法和龙格_库塔方法三种单步法基于Matlab求解常微分方程初值问题。 关键词：常微分初值问题；数值求解；欧拉方法；梯形方法；龙格_库塔方法；Matlab实现

1 引言 很多科学技术和工程问题常用微分方程的形式建立数学模型，因此微分方程的求解是很有意义的。建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些典型的方程，而对于绝大多数的微分方程问题，很难或者根本不可能得到它的解析解，实际问题终归结出来的微分方程主要靠数值解法。因此，研究微分方程求解的数值方法是非常有意义的。本文介绍了欧拉法、梯形法和四阶龙格_库塔方法三种单步法，通过Matlab的平台运行。

2 建模 常微分方程初值问题的数学模型是：求y()yx，使之满足

0'()(,)()yxfxyxbyay，(a),, (1) 其中0y，a，b是已知的常数，(,)fxy是已知函数，且满足条件： (,)(,')'fxyfxyLyy， 式中的L是不依赖于y，'y的常数，称为利普希茨（Lipschitz）常数。 由常微分方程理论知识可知，上述问题存在唯一解()yx。现在的目标就是计算区间,ab上等节

点ixaih处该未知函数的函数值()iyx，其中()/hbaN，0,1,2,...,iN。为此，可以用求解常微分方程问题的单步法，即欧拉方法、梯形方法和龙格_库塔方法求解。

3 算法求解和程序设计 3.1 欧拉方法

将微分方程离散化，用向前商1()()nnyxyxh代替微分()nyx，代入（1）中的微分方程，可得

1()()(,())1,2,3,...nnnnyxyxfxyxnh，（）

化简可得 桂林理工大学理学院2007级信息与计算科学专业·生产实习 张玉娟：常微分方程初值问题数值解法 7/6/2010 2 1()()(,())1,2,3,...nnnnyxyxfxyxhn，（） 如果用ny近似()nyx代入上式便可得到1()nyx的近似值1ny，计算式为：

1(,)1,2,3,...nnnnyyfxyhn，（） (2) 这样问题(1)的近似解课通过求解下面的差分初值问题：

10(,)1,2,3,...()nnnnyyfxyhnyay，（） (3)

得到，按(3)式由初值0y可逐次求出12y,,...y。 欧拉方法就是用差分方程初值问题(3)的解来近似微分方程初值问题(1)的解。即由公式(3)算出()nyx的近似值。这组公式求问题(1)的数值解就是著名的欧拉(Euler)公式。欧拉法是数值求解微分方程最简单的方法，它特别适合快速编程，尽管精度不是很高。 Euler方法的局部截断误差为

11()()(,())nnnnnTyxyxhfxyx

()()'()nnnyxhyxhyx 232''()/2()()nhyxhh

即局部截断误差是2h阶的。而数值算法的精度定义为：若一种算法的局部截断误差为1()ph，则称该算法具有p阶精度。显然p越大，方法的精度越高。上式说明，欧拉方法是一阶方法，因此它的精度不高。 Euler方法的整体截断误差为

()nnneyxy

((1)1)/((1)1)nnehLhLT ()(1)/()()LbaehLTh

即整体截断误差是h阶的。 欧拉法流程图： 桂林理工大学理学院2007级信息与计算科学专业·生产实习

张玉娟：常微分方程初值问题数值解法 7/6/2010 3 程序代码见附录。 3.2 梯形方法

在一阶微分方程111()()'()(,())nnnnxxnnxxyxyxyxdxfxyxdx右端积分用梯形求积公式近

似，并用ny代替()nyx，1ny代替1()nyx 111((,)(,))2nnnnnnhyyfxyfxy 梯形法的迭代公式为： (0)1(1)()111(,)((,)(,))2nnnkknnnnnnyfxyhhyyfxyfxy





，(k=0,1,2,...)

梯形方法也是隐式方法，要迭代求解，迭代收敛的条件是h12L。 梯形法的误差估计的局部截断误差是3()h，整体截断误差是2()h。即梯形方法具有2阶精度，是二阶方法，它优于欧拉方法。

x0+h=>x1 y0+h*f(x0,y0)=>y1

n=1

输出x1,y1 n=1+n x1=> x0 y1=> y0

结束 n=N ?

读入x0,y0,b,h 开始 计算N=fix((b-x0)/h) 桂林理工大学理学院2007级信息与计算科学专业·生产实习

张玉娟：常微分方程初值问题数值解法 7/6/2010 4 梯形法的流程图：

程序代码见附录。 3.3 龙格_库塔方法 由Lagrange微分中值定理，

11()()'()()()(,())nnnnnyxyxyxxyxhfy

记*(,())khfy，则得到 *1()()nnyxyxk

这样，给出*k的一种算法，就得到求解微分方程初值问题的一种计算公式。

n=1 x0=>x1 y0=>y (1) y0=>y1

n=2 x1+h=>x2 y1+h*f(x1,y1)=>y (2) y1+h/2*(f(x2,y2)+f(x1,y1))=>y2

输出x1,y1,x2,y2 n=n+1 x2=>x1 y2=>y1 abs(y2-y (2))<1e-3

n=N ? 结束

开始 计算N=fix((b-x0)/h) 读入x0,y0,b,h 桂林理工大学理学院2007级信息与计算科学专业·生产实习

张玉娟：常微分方程初值问题数值解法 7/6/2010 5 四阶龙格_库塔法是用1k，2k，3k和4k的加权平均值来近似*k。最经典的四阶龙格_库塔公式为： 121324311234(,)(,)22(,)22(,)y(22)2nnnn

nnnnnnKfxyhhKfxyKhhKfxyKKfxhyhKhyKKKK

















四阶龙格_库塔法的误差估计局部截断误差为5()Oh。 四阶龙格_库塔法流程图：

程序代码见附录。 值得指出的是，龙格_库塔方法的推导基于泰勒展开方法，因而它要求的解具有较好的光滑性质。

计算N=fix((b-x0)/h) n=1 x0+h=>x1 f(x0,y0)=>K1 f(x0+h/2,y0+h/2*K1)=>K2 f(x0+h/2,y0+h/2*K2)=>K3 f(x0+h,y0+h*K3)=>K4 y0+h/6*(K1+2*K2+2*K3+K4)=>y1

输出x1,y1 n=n+1 x1=>x0 y1=>y0 n=N ?

结束

开始 读入x0,y0,b,h 桂林理工大学理学院2007级信息与计算科学专业·生产实习

张玉娟：常微分方程初值问题数值解法 7/6/2010 6 反之，如果解得光滑性差，那么，使用四阶龙格_库塔法方法求得的数值解，其精度可能反而不如梯形方法。实际计算中，应当针对问题的具体特点选择合适的算法。 对欧拉法、梯形法和龙格_库塔法进行比较：在相同步长的情况下，欧拉法每步只计算一个函数值，梯形法每步计算两个函数值，四阶龙格_库塔法每步需计算四个函数值，就是说，四阶龙格_库塔法的计算量差不多是欧拉法的四倍，是梯形法的两倍；为了比较它们的计算精度，可以将欧拉法的步长取为h，将梯形法的步长取为2h,将四阶龙格_库塔法的步长取为4h。这样，如果用三种方法求解同一初值问题，则它们的计算量相当。在计算量相当的条件下，比较它们的计算结果，就能够看出它们的精度差异。 4 实例 用数值解法求解如下初值问题：

y'()1,01(0)1xxyxy



分别用欧拉法、梯形法和龙格_库塔法求解。 解： 为了比较三种方法的计算精度，将欧拉法的步长取为0.05，梯形法的步长取为0.1四阶龙格_库塔法的步长取为0.2，则三种方法在变量每增加0.1时，都需要计算4个函数值。从下面的图像界面显示可以看出，在0.2x处，欧拉法、梯形法和四阶龙格_库塔法的y值分别是1.01451，

1.01902，1.01873。 运行结果显示：

5 小结 对于一些不能求出解析解的方程通过用欧拉法、梯形法、四阶龙格_库塔法解决，只要输入相关给出的初始值和步长就可以相应的近似解。