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

1、第三讲数值算法之一：单步法接下来的章节，我们将详细讲解常微分方程的数值求解方法，特别是编程实 现方程的初值问题，这也将成为解决实际问题的必要基础和课程设计的主要内 容.本章重点讲解一阶常微分方程的数值算法中最简单的一类方法：单步法.首 先介绍euler法、后退euler法与梯形法，并分析单步法的局部截断误差，然后 给出了改进的euler法.3.1简单的单步法及基本概念3.1.1 euler法、后退euler法与梯形法求初值问题(1. 2. 1)=a<x<b1 ）=乂）的一种最简单方法是将节点的导数>，u）用差商代替，于是 h（1.2. 1）的方程可近似写成xn+l） +n =

2、/o i ,从x0出发冲）=灿）”0,由（3. 1. 1）求得破什。，八）=乃再将 w（6）代入（3. 1. 1）右端，得至俨（心）的近似a =乃+wii）, 一般写成»+1，乃= o，i，".(3. 1.2)该数值解法称为解微分方程初值问题的euler法.euler法的几何意义如图3-1所示.初值问题（1. 2. 1）的解曲线y=y（x）过点 以八），从出发，以/（、）为斜率作一段直线，与直线x =交点于巧（心乃）， 显然有乃=a +舻(x。，八)，再从巧出发，以为斜率作直线推进到x = a上一点其余类推，这样得到解曲线的一条近似曲线，它就是折线.euler法也可利用2+

3、1的taylor展开式得到，由1 2y(xn + h) = y(xn) + hy1 m + y1'),e (“j(3. 1. 3)略去余项，以就得到近似计算公式(3. 1.2).另外，还可对(1.2.1)的方程两端由&到&+1积分得= p"(3. 1.4)若右端积分用左矩形公式，用+i"x+i),则得(3. 1.2).如果在 (3. 1. 4)的积分中用右矩形公式，则得人+1 =人 +v(+1，人+1)，乃=0，1，".(3. 1.5)该算法称为后退(隐式)euler法.若在(3. 1. 4)的积分中用梯形公式，则得h八+1 =人+/(+1

4、，人+1)，乃=0，v"(3. 1. 6)该算法称为梯形方法.上述三个公式(3. 1.2), (3. 1.5)及(3. 1.6)都是由八计算a+i，这种只用前一步即可算出h+1的公式，我们称之为单步法，其中(3. 1.2)可由力逐次求出>2，"" 的值，称为显式方法，而(3.1.5)及(3.1.6)右端含有八1，当亡对7非线性时它不能直接求出力化此时应把它看作一个方程，求解h+i，这类方法称为隐式 方法.此时可将(3. 1. 5)或(3. 1. 6)写成不动点形式的方程yn+1 = w(+l«+l)+ g这里对式(3. 1.5)有p = 1 s =

5、 ykj 对(3. 1.6)则a = g与y«+i无关，可构造迭代法g，厂 0丄(3. 1. 7)由丁-/(u7)对y满足lipschitz条件，即存在常数£0，使对e /?有 i f(x, y, )-f(x,y2)<ly-y2(1.2.2)故有ytf砂|/('+i，a5h) /(+1«+1)|当邸£1或迭代法(3. 1.7)收敛到a+i,因此只要步长h足够小，就 可保证迭代(3. 1.7)收敛.对后退euler法(3. 1.5)，当h j时迭代收敛，对梯形2法(3. 1.6)，当力时迭代序列收敛.例3. 1用euler法、隐式euler法

6、、梯形法解y' = - y + x +l，y(o) =1取h=0. 1，计算到x=0. 5，并与精确解比较.解直接用给出公式计算.由于0，7) = 7+1 + 1，11 = 0.1，;=0，八=1,euler法的计算公式为a+i = +kyv +1)=(1 一 h)yn +hxk +h = q.9yk + 0.1xs + 0.1n=0时，乃=0什。+ o uo +01= 1.00000.0.其余n=l, 2, 3, 4的计算结果见表3-1.对隐式euler法，计算公式为«+i =« +a(7«+1 + 兀+1+1)解出八+1 =y7icv« +如

7、+1 +) = 0« +0上 +o.n)当 n=o 时，乃=&(凡 +0.1xq +0.11) =1.009091.其余 n=l，2, 3, 4 的计算结果 见表3-1.表3-1例3. 1的三种方法及精确解的计算结果euler 法 yn隐式euler法yn梯形法yn精确解y()011110.11.0000001.0090911.0047621.0048370.21.0100001.0264461.0185941.0197310.31.0290001.0513151.0406331.0408180.41.0561001.0830141.0700971.0703200.51.09

8、04901.1209221.1062781.106531对梯形法，计算公式为1»+1+ 1)h人+i =八 + 5(一人 + w (一八+1 +解得火+1 =咖+ 啦 + a) + 2h2+ k= a(i.9+o.2xm+o.21)当 n:0 时，a = t，：1，9 + 0-21)= h(x)4762.其余 n=l, 2, 3, 4 的计算结果见表7-1.木题的精确解为x) =x +e'表3-1列出三种方法及精确解的计算结果.附:具体三种算法程序如下：首先定义一阶方程右端函数f如下 function y=f(x, y)y二-y+x+1;(1) euler 法function

9、 y = deeuler (f, h, a, b, yo, varvec)%这是euler法的函数命令 %阶常微分方程的一般表达式的右端函数：f这里可用inline %自变量下限a;上限b%函数初值yo %积分步长h%常微分方程的变量组varvecformat long;%数据显示方式，不影响计算和存储方式，%是指小数点后15位数字表示n = (b-a)/h; y = zeros (n+l, 1); x = a:h:b; y(l) = yo; for i=2:n+1y (i) = y(il) +h*funval(f, varvec, x(il), y (il);%简革euler公式迭代，%本题

10、单步法的格式为：y(i) = y(i-l)+h. *(-y(i-l)+x(i-l)+l); endformat short;其中，funval (f, varvec, x(il), y (i-1)的m函数力function fv = funval (f, varvec, varval)%varvec为deeuler中输入的变量 组var = findsym(f); %按字典顺序返冋符号函数f中的所有变量名(符号)if length (var) < 4 %if 语句求值，注意：2009a 版 mat lab 改写为 <3 if var(1)二=varvec (1)fv = subs(

11、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, yo, varvec) format long;n = (b-a)/h;y = zeros (n+l, 1);y(l) = y0；x = a:h:b;var = fin

12、dsym(f);for i:2:n+lfx = funval (f, var(l), x(i); gx = y(i-l)+h*fx - varvec(2); y (i) = newtonroot (gx, -10, 10, eps);endformat short;其中，n e w t o n r o o t函数的m文件为function root=newtonroot(f, a, b, eps)if(nargin=3) eps=l. 0e4;endfl=subs(sym(f)，findsym(sym(f)，a); f2=subs (sym(f), findsym(sym(f), b); if

13、(fl=0)root=a;endif(f2=0)root=b;endtol=l;fun=diff(sym(f);fa=subs(sym(f)，f indsym(sym(f)，a); fb=subs(sym(f), findsym(sym(f), b); dfa=subs(sym(fun)，findsym (sym (fun)，a); dfb=subs (sym(fun), findsym(sym(fun)，b); if(dfadfb)root=a-fa/dfa;elseroot=b-fb/dfb;endwhile(tol>eps) rl=root;fx=subs(sym(f)，finds

14、ym(sym(f), rl); dfx=subs(sym(fun)，findsym (sym (fun)，rl); root=rl-fx/dfx;tol=abs (root-rl);end(3)梯形法function y = deimpeulerl(f, h，a，b，yo) format 1ong;n = (b-a)/h;y = zeros (n+l，1);x = a:h:b-h;y(l) = yo；var = findsym(f);yd = diff (f，var(4) ;%关于f屮第二个变量y或者v的导数 for i=2:n+lfx = f - ydvar (4); dfy = subs(

15、yd, findsym (yd), x(il); ex = subs(fx，findsym(fx)，x(i-l); y(i) = (y(i-l)+h*cx)/(l-h*dfy);endformat short;3.1.2单步法的局部截断误差解初值问题（1.2. 1）的单步法可表示为其中分与/有关，称为增量函数，当於含有夂+1时，是隐式单步法，如(3. 1.5) 及(3. 1.6)均为隐式单步法，而当不含-1时，则为显式单步法，它表示为y肿1 =yx +a(xm,a),« = o,i,-(3. i 9)如euler法(3. 1. 2)，=7).为讨论方便，我们只对显式单步法 (3. 1

16、. 9)给出局部截断误差概念.定义3.1设y(x)是初值问题(7. 1. 1)的精确解，记tn+i =-(3.1.10)称7；+1为显式单步法(3. 1.9)在；的局部截断误差.这里，7；+1之所以称为局部 截断误差，可理解为用公式(3. 1.9)计算时，前面各步都没有误差，即只考虑由么计算到这一步的误差，此时由(3. 1. 10)有xxn+i)-yn+i =x«+i)-a +a(xs,a)=x»+i )-y(xk)-hxk, yxn )，a) = tn+1局部截断误差(3. 1. 10)实际上是将精确解yoo代入(3. 1.9)产生的公式误差，利用taylor展开式可得到

17、t«+i = °例如对euler法(3. 1. 2)有 於0,夂a) = /(x，y),故= x» +办)70«)吻(>«) =v.o«) + !什.0«) + =咐2)2 6它表明euler法(3. 1. 2)的局部截断误差为tn+1 =y() + 0什3)，称 2/ ()为局部截断误差主项.定义3. 2设y(4是初值问题(7. 1. 1)的精确解，若显式单步法(3. 1. 9)的局部截断误差1+1=+1 p是展开式的最大整数，称p为单步法(3. 1.9)的阶，含 的项称为局部截断误差主项.根据定义，euler法(3

18、. 1. 2)中的71故此方法为一阶方法.对隐式单步法(3. 1.8)也可类似求其局部截断误差和阶，如对后退euler法 (3. 1.5)有局部截断误差t«+i =池+1)-池)-wk+1，池+1)=+ 办)_+)=如()十/从)+ yd + -4乂()吻d + =-答z(x) + o(的故此方法的局部截断误差主项为=也是一阶方法.对梯形法(3. 1.6)同样有ht«+1 =池+1)-池)-$/(，x么)+ /(+1，池+1) h=-盖/.从)哪4)它的局部误差主项为= 2,方法是二阶的.3.1.3 改进 euler 法上述三种简单的单步法中，梯形法(3. 1.6)为二阶方

19、法，且局部截断误差最 小，但方法是隐式的，计算要用迭代法.为避免迭代，可先用euler法计算出a+1的近似h+i，将(3. 1.6)改为夂+i =y，+h_人+i =« +i/o«，y«) + /(+i，j+1)(3. 1. 11)称为改进euler法.改进的euler法是二阶显式方法，即a+i =« +4/(w«) + /(+i，人 +机，人)(3. 1.右端已不含可以证明t、+i =00, p=2,故方法仍为二阶的，与梯形法一样, 但用(3. 1.11)计算l+i不用迭代.例3. 2用改进euler法求例3. 1的初值问题并与euler法和

20、梯形法比较误 差的大小.解将改进euler法用于例3. 1的计算公式 hj+i =人 4- d + (一人一+ xk +1) +1)r, a(2-a), a(l-a) h,?a(2 -a)、=i一 。、b + )= 0.905人 +0.095xh +0.l当n=0时，乃=0905八+ 0 095 +ol=5000.其余结果见表3-2.表3-2改进euler法及三种方法的误差比较改进误差 |yn-y(xn)|euler方法|yn -y(xn)|梯形法 |yn-y(xn)|o.ll.005000l.6xl0.44.8xw37.5xl。-50.21.0190252.9xw48.7 x io-31.4

21、x10"40.31.0412184.0x101.2xw21.9xl0_40.4l.0708024.8xl0'41.4xl0-22.2xw40.51.1.070765.5xw41.6xl0-22.5x10从表3-2中看到改进euler法的误差数量级与梯形法大致相同，而比euler法小 得多，它优于euler法.附：改进的euler法程序function y = demodifeuler(f, h,a,b,yo,varvec) format long;n = (b-a) /h;y = zeros(n+l,1);y(1) = yo;x = a : h : b;var = findsym(f);for i=2:n+lvl = funval(f,varvec, x(i-l) y(i-l); t = y (i-1) + h*vl; v2 = funval(f,varvec,x(i) t); y(i) = y(i-1)+h*(vl+v2)/2;endformat short;讲解:求初值问题(1.2.1)的数值解就是在假定初值问题解存在唯一的前提下在 给定区间久利上的一组离散点<xi 么<6上求解析解7=犬幻 的一组近似八，乃，人，为此先要建立求数值解八的计算公式，通常称