第三讲(实践课)数值算法之单步法Euler法

1、第三讲数值算法之一：单步法接下来的章节，我们将详细讲解常微分方程的数值求解方法，特别是编程实现方程的初值问题，这也将成为解决实际问题的必要基础和课程设计的主要内 容.本章重点讲解一阶常微分方程的数值算法中最简单的一类方法：单步法 .首 先介绍Euler法、后退Euler法与梯形法，并分析单步法的局部截断误差，然后 给出了改进的Euler法.3.1简单的单步法及基本概念3.1.1 Euler法、后退Euler法与梯形法y' f(x, y), a x b(1.2.1)求初值问题y(a) y。的一种最简单方法是将节点Xn的导数y'(Xn)用差商y(xnh)一y(xn)代替，于是(1.

2、2.1)的方程可近似写成讯岳1)迨十炉亿丿)(G)" 7丄(3.1.1)从王。出发(切=肌无厂加,由(3.1.1)求得巩町厂*孙心小)=列再将k 舟】)代入(3.1.1)右端，得到pE)的近似卩厂几4财韜必)|, 一般写成儿+厂几十纱亿几)严=0丄(3.1.2)该数值解法称为解微分方程初值问题的Euler法.Euler法的几何意义如图3-1所示.初值问题(1.2.1)的解曲线y=y(x)过点%仏九),从丽出发，以盅坯凤)为斜率作一段直线，与直线"引交点于Pi(m”)，显然有”二儿+孙(知片),再从珂出发，以fZJ为斜率作直线推进到"呵上一点巴，其余类推，这样得到解

3、曲线的一条近似曲线，它就是折线R毘马-.Euler法也可利用卩尤时1)的Taylor展开式得到，由(3.1.3)略去余项，以叫*M兀)，就得到近似计算公式(3.1.2).另外，还可对(1.2.1)的方程两端由抵到Si积分得曲祕)-片=畑巩小沁*叫！(3.1.4)若右端积分用左矩形公式，用y宀E，儿/巩耳“(3.1.4)的积分中用右矩形公式，则得，则得(3.1.2).如果在+汕(心竝,儿沁)用= 0,1,-*(3.1.5)该算法称为后退(隐式)Euler法.若在(3.1.4)的积分中用梯形公式，则得凡+厂儿十才/(兀十川Xm几+L)”刃=0丄(3.1.6)该算法称为梯形方法.上述三个公式(3.1

4、.2) ,(3.1.5)及(3.1.6)都是由儿计算儿+1 ,这种只用前一 步即可算出儿口的公式，我们称之为单步法，其中(3.1.2)可由珂逐次求出尸1M宀 的值，称为显式方法，而(3.1.5)及(3.1.6)右端含有几L)当f对y非线性 时它不能直接求出 g 此时应把它看作一个方程，求解T 这类方法称为隐式 方法.此时可将(3.1.5)或(3.1.6)写成不动点形式的方程九1 "欽S畑)十呂这里对式(3.1.5)有0T呂=儿，对(3.1.6)贝叨=耳十彳了(耳必g与九门无关，可构造迭代法讯于=HffQz畑+呂,b Q乙(3.1.7)由于Y(2)对y满足Lipschitz 条件，即存

5、在常数L 0，使对y1, y2 R有|f(x,y1) f(x, y2)| L|y1 y? |(122)故有y防卩儿吐丈厶)山儿韵S04心-加I，因此只要步长h足够小，就当方朋< 1或X ，迭代法(3.1.7)收敛到X可保证迭代(3.1.7)收敛.对后退Euler法(3.1.5)，当 < -时迭代收敛，对梯形 2法(3.1.6)，当<2时迭代序列收敛.例3.1用Euler法、隐式Euler法、梯形法解W=-y + x*l,y(O) =1取h=0.1，计算到x=0.5，并与精确解比较.解 直接用给出公式计算.由于圮芯叭=- +量+1上=0丄誥0 = 0丿0 =Euler法的计算公

6、式为=(1 -册儿十饥十內=0 9兀十Og十0.1n=0时，卄79儿十05+011.0000).其余n=1,2,3,4的计算结果见表 3-1.对隐式Euler法，计算公式为儿41 "儿+風-畑+耳订*1)解出几+1 =耐山十抵+1十h)=订仇十0. g十011)当 n=0 时=订(必 + 0+ 0.11)= 1.009091.其余 n=1,2,3,4 的计算结果 见表3-1.表3-1例3.1的三种方法及精确解的计算结果夠Euler法丽隐式Euler袪Yn梯旳法为楕确解Jt和0J1110.11 DOOOOO1.0090911.0047d21 004337021.0100001.0364

7、41.0185P41.01P7310.31 03(>0001.0513151.0406331 D4D313041.05(51001.0830141.000071.0703200.51 0304PO1.1209221.10627S1 106531对梯形法，计算公式为= F. + 2 W-y，+心斗 1)+ (-儿+ 4 S +1)解得=肓(1 叽+ 0 2+0 21)当n=0时，”刃"TSr。汽其余 n=1,2,3,4的计算结果见表7-1.本题的精确解为-x-K表3-1列出三种方法及精确解的计算结果.附: 具体三种算法程序如下：首先定义一阶方程右端函数f如下fun ctio n

8、Y=f(x,y)丫二y+x+1;(1)Euler 法这是Euler法的函数命令f这里可用f=iniinefun ctio n y = DEEuler(f, h,a,b,yO,varvec)% 阶常微分方程的一般表达式的右端函数：%自变量下限a ；上限b%函数初值y0%积分步长h%常微分方程的变量组varvec format Iong;%数据显示方式，不影响计算和存储方式, %是指小数点后15位数字表示N = (b-a)/h;y = zeros(N+1,1);x = a:h:b;y(1) = y0;for i=2:N+1y(i) = y(i-1)+h*Fu nv al(f,varvec,x(i-

9、1), y(i-1);御单Euler公式迭代，%本题单步法的格式为：y(i) = y(i-1)+h.*(-y(i-1)+x(i-1)+1);endformat short;其中，Funval(f,varvec,x(i-1), y(i-1)的M函数为function fv = Fun val(f,varvec,varval)%varvec为 DEEuler中 输入的变量组var = findsym(f); %按字典顺序返回符号函数f中的所有变量名(符号)if length(var) < 4 %if语句求值，注意：2009a版matlab 改写为<3'if var(1) = v

10、arvec(1)fv = subs(f,varvec(1),varval(1);elsefv = subs(f,varvec(2),varval(2);endelsefv = subs(f,varvec,varval);end本题输入命令为 >> syms u,v>> Y = DEEuler(-v+u+1, 0.1,0,0.5,1,u v)2）隐式 Euler 法function y = DEimpEuler(f, h,a,b,y0,varvec) format long;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;

11、var = findsym(f);for i=2:N+1fx = Funval(f,var(1),x(i);gx = y(i-1)+h*fx - varvec(2);y(i) = NewtonRoot(gx,-10,10,eps);endformat short;其中，NewtonRoo®数的M文件为function root=NewtonRoot(f,a,b,eps) if(nargin=3)eps=1.0e-4; end f1=subs(sym(f),findsym(sym(f),a);f2=subs(sym(f),findsym(sym(f),b); if(f1=0)root=

12、a;endif(f2=0)root=b;endtol=1;fun=diff(sy m(f);fa=subs(sym(f),fi ndsym(sy m( f),a); fb=subs(sym(f),fi ndsym(sym(f),b); dfa=subs(sym (fun ),fi ndsym(sym (fun ),a); dfb=subs(sym (fun ),fi ndsym(sym (fun ),b); if(dfa>dfb)root=a-fa/dfa;elseroot=b-fb/dfb;endwhile(tol>e ps)r仁root;fx=subs(sy m(f),fin

13、dsym(sy m(f) ),r1); dfx=subs(sym (fun ),fi ndsym(sym (fun ),r1); root=r1-fx/dfx;tol=abs(root-r1);end(3)梯形法fun ctio n y = DEi mp Euler1(f, h,a,b,yO) format long;N = (b-a)/h;y = zeros(N+1,1);x = a:h:b-h;y(1) = yO;var = fin dsym(f);yd = diff(f,var(4);%关于f中第二个变量y或者v的导数for i=2:N+1fx = f - yd*var(4);dfy =

14、 subs(yd,fi ndsym(yd),x(i-1);cx = subs(fx,fi ndsym(fx),x(i-1);y(i) = (y(i-1)+h*cx)/(1-h*dfy);endformat short;3.1.2单步法的局部截断误差解初值问题(1.2.1)的单步法可表示为其中0与/有关，称为增量函数，当5含有"沁时，是隐式单步法，如(3.1.5)及(3.1.6)均为隐式单步法，而当不含儿*1时，则为显式单步法，它表示为儿1 =儿7敕m儿閻卫=0,1,(3.1.9)如Euler法(3.1.2)，呎町=(兀y).为讨论方便，我们只对显式单步法(3.1.9)给出局部截断误差

15、概念.定义3.1设y(x)是初值问题(7.1.1)的精确解，记(3.1.10)称Tn 1为显式单步法(3.1.9)在Xn 1的局部截断误差.这里，Tn 1之所以称为局部截断误差，可理解为用公式(3.1.9)计算时，前面各步都没有误差，即= ><£),只考虑由汇计算到这一步的误差，此时由(3.1.10)有就勒1+1)-珀+1 =丿0曲)-儿斗必曲兀，儿/) =y(眄¥1) /(召)-初呱务(心)出)=JT我门5护"丄).例如对Euler法(3.1.2)有局部截断误差(3.1.10)实际上是将精确解代入(3.1.9)产生的公式误 差，利用Taylor展开式

16、可得到$0耐-/(扎力，故入1 =- 7(G)-纱(r)皿 F-心)-妙幺)2心七+衲”区)+一0的它表明Euler法(3.1.2)的局部截断误差为*丫心)为局部截断误差主项.定义3.2设yx）是初值问题（7.1.1）的精确解，若显式单步法（3.1.9）的局部截断误差匚*1（护""），网是展开式的最大整数，称 国为单步法（3.1.9）的阶，含 血闷的项称为局部截断误差主项.根据定义，Euler法（3.1.2）中的|刘=1故此方法为一阶方法.对隐式单步法（3.1.8）也可类似求其局部截断误差和阶，如对后退Euler法（3.1.5）有局部截断误差/'g 尸-妙 0)+

17、T卄1 =7（心+1）-尹（召）-皆氐1 =心十Q - yg - V （耳十筒 =奶（讣分g詣故此方法的局部截断误差主项为兰TOJp = i，也是一阶方法.对梯形法2（3.1.6）同样有T* 仇+ 于(心+1 J(gi)=心砒)-畑)-細仏KqFl-八区）2（胪）它的局部误差主项为-冷气2，方法是二阶的.3.1.3 改进 Euler 法上述三种简单的单步法中，梯形法（3.1.6）为二阶方法，且局部截断误差最 小，但方法是隐式的，计算要用迭代法.为避免迭代，可先用Euler法计算出y殊討的近似兀将（3.1.6）改为A+1 =儿十妙（X片）(3.1.11)几十1 =儿+£了（和片）十十衣

18、眾）称为改进Euler法.改进的Euler法是二阶显式方法，即女儿十1 -儿近S心必）Tgi必+幼区必）(3.1.12)与梯形法一样,右端已不含4可以证明T" =0（沪），国=2,故方法仍为二阶的,但用（3.1.11）计算+1不用迭代.例3.2用改进Euler法求例3.1的初值问题并与Euler法和梯形法比较误 差的大小.解将改进Euler法用于例3.1的计算公式*地十1 A.十-【（-丿十心+ 1+ （-几-k（-人+耳+1 + %十1 + 1） “嘗十処宀”轨+山+竺矜=0.905/, 巩 +0 1.其余结果见表3-2.当 n=0 时 0 905)0 + 0 095心 +01 -

19、1 005000表3-2改进Euler法及三种方法的误差比较改进E必冲纵误差|仏-y（K J1抚肚方法1兀-曲J 1梯形袪1灯吨J0.11 00500016x10-*4&1尸75x10-0.21.0190252.98/7x101你10.31.0412184.0x10-*12x101.9x10-*0.41.07030248x10-14x1 尸2.2x1 旷0.51.1.070765.5x10-*1 6工沪2.5x10从表3-2中看到改进Euler法的误差数量级与梯形法大致相同，而比Euler法小 得多，它优于Euler法.附：改进的Euler法程序fun cti on y = DEMod

20、ifEuler(f, h,a,b,y0,varvec) format long;N = (b-a)/h;y = zeros(N+1,1);y=yO;x = a:h:b;var = fin dsym(f);for i=2:N+1v1 = Fu nval(f,varvec,x(i-1) y(i-1); t = y(i-1) + h*v1;v2 = Fu nval(f,varvec,x(i) t);y(i) = y(i-1)+h*(v1+v2)/2;endformat short;讲解:求初值问题（1.2.1）的数值解就是在假定初值问题解存在唯一的前提下在 给定区间|>上的一组离散点口二可*可U'心G+<力上求解析解y = 的一组近似几J1，儿1为此先要建立求数值解儿I的计算公式，通常称为差分 公式，简单的单步法就是由计算下一步坯，构造差分公式有三种