常微分方程的Euler解法32162.pdf

《常微分方程的Euler解法32162.pdf》由会员分享，可在线阅读，更多相关《常微分方程的Euler解法32162.pdf（14页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、 毕 业 论 文 题 目：常微分方程的 Euler 解法 及其程序设计 学 院：数学与信息科学学院 专 业：数学与应用数学 毕业年限：2011 年 6 月 学生：学 号：指导教师：.常微分方程的 Euler 解法及其程序设计 摘要 本文总结了常微分方程的 Euler 解法，对各种格式给出了误差估计，设计了这些格式的计算程序 关键词 常微分方程；Euler 解法；误差分析；程序设计 Euler Method of Ordinary Differential Equation and Its Programming Abstract Euler method of ordinary differe

2、ntial equation is summarized，the error of each format is analyzed and its programming is designed in this paper.Keywords Ordinary differential equation;Euler method;Error analysis;Programming .科学技术中常常需要求解常微分方程的定解问题，这类问题最简单的形式，即为微分方程 (,)dyf x ydx （1）的初值问题 00(,),().dyf x ydxy xy （2）定理 （存在与唯一性定理）如果方程（1

3、）的右端函数(,)f x y在闭矩形域 000000:,R xaxxa ybyyb 上满足如下条件：（1）在R上连续；（2）在R上关于变量y满足利普希茨（Lipschitz）条件，即存在常数L，使 对于R上任何一对点(,)x y和(,)x y有不等式：(,)(,)f x yf x yL yy,则初值问题（2）在区间00000 xhxxh上存在唯一解 00(),()yy xy xy,其中0(,)min(,),max(,)x yRbhaMf x yM 根据存在与唯一性定理，只要(,)f x y关于y满足 Lipschitz 条件 (,)(,)f x yf x yL yy，即可保证其解()yy x存

4、在并唯一 然而解析方法只能用来求解少数较简单和典型的常微分方程，例如线性常系数微分方程等，对于变系数常微分方程的解析求解就比较困难，而一般的非线性常微分方程就更不用说，因此，在大多数情况下，实际问题中归结出来的微分方程主要靠数值解法求解 .所谓数值解法，就是寻找()y x在一系列离散节点121nnxxxx上的近似值121,nny yyy相邻两个节点的间距1nnhxx称为步长，假定h为定数，节点为0,0,1,2,nxxnh n.1 Euler 解法（1）Euler 格式 Euler 格式的计算公式为 1(,)nnnnyyhf xy (0,1,2)n.（3）下面用 4 种方法推导公式（3）a 泰勒

5、展开法 在nx处展开1()()nny xy xh有 211()()()()2!nnnny xy xhy xh y 1()nnnxx,（4）略去余项，得 1()()()()(,()nnnnnny xy xhy xy xhf xy x.用近似值ny代替()ny x，把上式右端所得结果记为1ny，得 1(,)nnnnyyhf xy (0,1,2)n.b 数值差商 由导数定义知 1()()()(,()nnnnny xy xy xf xy xh，所以 1()()(,()nnnny xy xhf xy x.用ny代替()ny x，把上式右端所得结果记为1ny，即得公式（3）.c 数值积分 在区间1,nnx

6、 x上对微分方程（1）进行积分，得 11()()(,()nnxnnxy xy xf x y x dx,（5）.利用左矩形公式得 1()()(,()nnnny xy xhf xy x.用ny代替()ny x，把上式右端所得结果记为1ny，即得公式（3）.d 几何方法 在Oxy平面中，微分方程 00.(,),()yf x yy xy 的解()yy x称为它的积分曲线积分曲线上的一点(,)x y的切线斜率等于函数(,)f x y的值如果按函数(,)f x y在Oxy平面上建立一个方向场，那么，积分曲线上的每一点的切线方向均与方向场在该点的方向一致 基于上述对微分方程解的几何解释，从初始点000(,)

7、pxy出发，依方向场在该 点的方向上推进到1xx上一点1p，然后再从1p依方向场的方向推进到2xx上 的一点2p循环前进，可作出一条折线012p p p（如图 1）图 1 一般地，设已作出该折线的极点np，过(,)nnnpxy依方向场的方向再推进到111(,)nnnpxy，显然，两个极点np，1np的坐标有下列关系：11(,)nnnnnnyyf xyxx,1(,)nnnnyyhf xy.若初值0.y已知，则依上式可逐步算出 1000(,)yyhf xy，2111(,)yyhf x y，1(,)nnnnyyhf xy.即得 Euler 格式的计算公式（3）.（2）梯形格式 如果用梯形公式计算（5

8、）式右端的积分，得 111()()(,()(,()2nnnnnnhy xy xf xy xf xy x (0,1,2)n.对于上式右端，用ny代替()ny x，把上式右端所得结果记为1ny，即得梯形公式 111(,)(,)2nnnnnnhyyf xyf xy (0,1,2)n.（6）梯形公式（6）是关于1ny的隐式方程，因此称梯形法为隐式方法，而公式（3）是1ny的显式形式，因此称欧拉法为显式方法（3）改进的 Euler 格式 先用 Euler 格式求得一个初步的近似值1ny，称之为预测值，预测值1ny的精度可能很差，再用梯形公式将它校正一次，得1ny，这个结果称之为校正值 而这样建立的预测校

9、正系统通常称为改进的 Euler 格式 预测 1(,)nnnnyyhf xy 校正 111(,)(,)2nnnnnnhyyf xyf xy 这一计算格式亦可表示为：11(,)(,(,)2nnnnnnnnhyyf xyf xyhf xy,或表示为下列平均化形式：11(,)(,)1()2pnnncnnpnpcyyhf xyyyhf xyyyy 2 误差与精度（1）截断误差 在实际计算中，即使初始值是准确的，但所得的数值解往往与初始值问题的.准确解有一定的误差.用 111()nnny xy,表示在节点1nx处准确解与数值解之差，称之为整体截断误差，它不仅与1nxx这步计算有关，而且与121,nnxx

10、x这些步计算有关.为方便估计，假定()iy in为准确的，即在()()iiyy xin的前提下估计误差 111()nnnEy xy,这种误差称为局部截断误差.如果局部截断误差11(),pnEO h这里h是等距节点的步长，则称所用的数值解的方法的阶数是p，或称该方法有p阶精确度.显然，步长h越小，阶数p越高，则局部截断误差越小，计算结果的精度也越高.（2）Euler 格式是一阶方法 首先，可以通过几何直观来考察 Euler 格式的精度 假设()nnyy x，即顶点np落在积分曲线()yy x上，那么，按 Euler 格式作出的折线1nnp p便是()yy x过点np的切线（如图 2）从图形上看，

11、这样定出的顶点1np显著地偏离了原来的积分曲线，可见 Euler格式是相当粗糙的 图 2 首先考虑 Euler 格式的局部截断误差1nE，由（4）式得 211()()()()2!nnnny xy xhy xh y 1()nnnxx.由 Euler 格式，有.1(,)()nnnnnnyyhf xyyhy x,上面两式相减，则Euler格式的局部截断误差显然为 22111()()()22nnnnhhEy xyyyx.若令 max()a x bMyx,则有 221()2nMEhO h,这表明 Euler 格式是一阶方法.下面考虑 Euler 格式的整体截断误差n.由（3）式和（4）式，有 1111(

12、)()(,()nnnnny xy xhf xy xE,111(,)nnnnyyhf xy.两式相减得 111111(,()(,)nnnnnnnh f xy xf xyE.由 Lipschitz 条件及22iMEh，有 22111(1)22nnnnMMhLhhLh.反复使用这个不等式，得 21(1)2nnMhLh 12222200(1)(1)(1)(1)222nniniMMMhLhhLhhLhhL 20(1)1(1)2nnMhLhLhhL 0(1)(1)12nnMhhLhLL.又1(),(1)(1)nnhLb a LhLnhbahLhLe，.所以()()0(1)2b a Lb a LnMheeL

13、.上式表明，Euler 格式的整体截断误差由初始误差0和局部截断误差界22Mh决定.又00()y xy，所以00，则()nO h.这表明当0h时，0n，也就是说当h充分小时近似值ny能和)(nxy充分接近，即数值解是收敛的.综上所述，Euler 格式的整体截断误差比局部截断误差低一阶,这表明 Euler格式的精度是比较差的，为了提高精度，减小误差，可采用梯形格式.（3）梯形格式是二阶方法 首先考察梯形格式的局部截断误差.由泰勒展式，有 231()()()()()2nnnnhy xy xhy xyxO h，（7）21()()()()nnny xy xhyxO h.（8）由（8）式，有 1()()

14、()()nnny xy xyxO hh.将上式代入（7）式，得 311()()()()()2nnnnhy xy xy xy xO h,（9）又由梯形公式（6），有 111(,)(,)2nnnnnnhyyf xyf xy.（10）由微分中值定理，有 1111111(,)(,()(,)()nnnnynnnf xyf xy xfxyy x 1111()(,)()nynnny xfxyy x,这里是介于1ny与1()ny x之间的值，将上式代入（10）式，有 11111()()(,)()22nnnnynnnhhyyy xy xfxyy x.将上式与（9）式相减，得梯形格式的局部截断误差为 311111

15、()(,)()()2nnynnnhyy xfxyy xO h,即有 31111()()1(,)2nnynyy xO hhfx.因此梯形方法是二阶方法.梯形方法是隐式的，可用迭代法求解，其迭代公式如下：(0)1(1)()111(,)(,)(,)2nnnnkknnnnnnyyhf xyhyyf xyf xy (0,1,2,)k.下面分析迭代过程的收敛性：因为(1)()111111(,)(,)2kknnnnnnhyyf xyf xy，于是有(1)()11112kknnnnhLyyyy，式中L为(,)f x y对y满足 Lipschitz 常数，如果h选取充分小，使得12hL，则当k 时有(1)11k

16、nnyy，这说明迭代过程是收敛的 可见，梯形格式虽提高了精度，但在应用迭代公式进行实际计算时，每迭代一次，都要重新计算函数 f 的值，而迭代又要反复进行若干次，计算量很大，为减小计算量，便可采用改进的 Euler 格式 3 程序设计（1）Euler 格式的算法及其程序设计 应用 Euler 方法计算初值问题 (,),()yf t yy a atb 的解()y t在区间,a b上的 N 个等距点的近似值的算法 输入 端点,a b;区间等分数 N;初值 输出 ()y t在t的 N 个点处的近似值y.1step ()/;hbaNtay 2step 对1,2,iN做34stepstep 3step (

17、,);yyhf t ytaih 4step 输出(,)t y 5step 停机 例：求解初值问题 2/(0)1yyx yy 01x.源程序如下：#include#include void main()int i;double a,b,c,h,t,y,z;a=0;b=1;c=1;h=(b-a)/10;t=a;y=c;for(i=1;i=10;i+)y=y+h(y-2*t/y);t=a+i*h;z=sqrt(1+2*t);printf(%f%f%fn,t,y,z);.运行结果如下：（2）改进的 Euler 格式的算法及其程序设计 应用改进的 Euler 方法计算初值问题 (,),()yf t yy

18、 a atb 的解()y t在区间,a b上的 N+1 个等距点的近似值的算法 输入 端点,a b;区间等分数 N;初值 输出 ()y t在t的 N+1 个点处的近似值y 1step ()/;hbaNtay 2step 输出(,)t y 3step 对1,2,iN做45stepstep 4step (,);(,);1().2pcppcyyhf t ytaihyyhf t yyyy 5step 输出(,)t y 6step 停机 下面用改进的 Euler 格式计算上面的例题：求解初值问题 2/(0)1yyx yy 01x.源程序如下：#include#include void main().in

19、t i;double a,b,c,h,t,y,y1,y2,z;a=0;b=1;c=1;h=(b-a)/10;t=a;y=c;for(i=1;i=10;i+)y1=y+h(y-2*t/y);t=a+i*h;y2=y+h(y1-2*t/y1);y=(y1+y2)/2;z=sqrt(1+2*t);printf(%f%f%fn,t,y,z);运行结果如下：从上面的运行结果可以看出，该常微分方程欧拉解法的精确解ny和其数值解()ny x很接近由此可知，常微分方程欧拉解法对大量计算数据来说精度是很高的，可以满足绝大多数工程上的要求由于欧拉解法的稳定性和收敛性，可以.通过减少h的值来进一步提高其精度 参考文

20、献 1 庆扬，王能超，易大义数值计算方法M：华中科技大学，2006 2 林成森数值计算方法M：科学，2005 3 信真，车刚明，欧阳洁，封建湖计算方法M：西北工业大学，2000 4 朱长青数值计算方法及其应用M：科学，2006 5 曾喆昭，文卉数值计算M：清华大学，2005 6 丽娟常微分方程的 Euler 解法及其计算机实现J师学院学报（自然科学版），2005，24（6）：11-14 7 马燕，根耀，杜利锋，竹林常微分方程的数值解法及其在计算机中的实现J大学学报（自然科学版），2000，19（2）：23-25 8 均亮，邓易冬，徐宏坤基于 C+的常微分方程欧拉解法的误差分析J 信息技术，2007，7：92-94