《常微分方程课程设计》指导书 3

第3章数值算法之一：单步法接下来的章节，我们将详细讲解常微分方程的数值求解方法，特别是编程实现方程的初值问题，这也将成为解决实际问题的必要基础和课程设计的主要内容。本章重点讲解一阶常微分方程的数值算法中最简单的一类方法：单步法。首先介绍Euler法、后退Euler法与梯形法，并分析单步法的局部截断误差，然后给出了改进的Euler法。3.1简单的单步法及基本概念3.1.1Euler法、后退Euler法与梯形法求初值问题(1.2.1)的一种最简单方法是将节点nx的导数'()nyx用差商()()nnyxhyxh代替，于是(1.2.1)的方程可近似写成(3.1.1)从出发，由(3.1.1)求得再将代入(3.1.1)右端，得到的近似，一般写成(3.1.2)该数值解法称为解微分方程初值问题的Euler法.Euler法的几何意义如图3-1所示。初值问题(1.2.1)的解曲线y=y(x)过点，从出发，以为斜率作一段直线，与直线交点于，显然有，再从出发，以为斜率作直线推进到上一点，其余类推，这样得到解曲线的一条近似曲线，它就是折线.图3-1Euler法的几何意义显示Euler法也可利用的Taylor展开式得到，由(3.1.3)略去余项，以，就得到近似计算公式(3.1.2).另外，还可对(1.2.1)的方程两端由到积分得(3.1.4)若右端积分用左矩形公式，用，，则得(3.1.2).如果在(3.1.4)的积分中用右矩形公式，则得(3.1.5)该算法称为后退(隐式)Euler法。若在(3.1.4)的积分中用梯形公式，则得(3.1.6)该算法称为梯形方法。上述三个公式(3.1.2)，(3.1.5)及(3.1.6)都是由计算，这种只用前一步即可算出的公式，我们称之为单步法，其中(3.1.2)可由逐次求出的值，称为显式方法，而(3.1.5)及(3.1.6)右端含有当f对y非线性时它不能直接求出，此时应把它看作一个方程，求解，这类方法称为隐式方法。此时可将(3.1.5)或(3.1.6)写成不动点形式的方程这里对式(3.1.5)有，对(3.1.6)则，g与无关，可构造迭代法(3.1.7)由于对y满足Lipschitz条件(1.1.2)，故有当或，迭代法(3.1.7)收敛到，因此只要步长h足够小，就可保证迭代(3.1.7)收敛。对后退Euler法(3.1.5)，当时迭代收敛，对梯形法(3.1.6)，当时迭代序列收敛。例3.1用Euler法、隐式Euler法、梯形法解取h=0.1，计算到x=0.5，并与精确解比较。解直接用给出公式计算。由于，Euler法的计算公式为n=0时，.其余n=1,2,3,4的计算结果见表3-1.对隐式Euler法，计算公式为解出当n=0时，.其余n=1,2,3,4的计算结果见表3-1.表3-1例3.1的三种方法及精确解的计算结果对梯形法，计算公式为解得当n=0时，.其余n=1,2,3,4的计算结果见表7-1.本题的精确解为，表3-1列出三种方法及精确解的计算结果。附：具体三种算法程序如下：首先定义一阶方程右端函数f如下functionY=f(x,y)Y=-y+x+1;（1）Euler法functiony=DEEuler(f,h,a,b,y0,varvec)%这是Euler法的函数命令%一阶常微分方程的一般表达式的右端函数：f这里可用f=inline%自变量下限a；上限b%函数初值y0%积分步长h%常微分方程的变量组varvecformatlong;%数据显示方式，不影响计算和存储方式，%是指小数点后15位数字表示N=(b-a)/h;y=zeros(N+1,1);x=a:h:b;y(1)=y0;fori=2:N+1y(i)=y(i-1)+h*Funval(f,varvec,[x(i-1),y(i-1)]);%简单Euler公式迭代，%本题单步法的格式为：y(i)=y(i-1)+h.*(-y(i-1)+x(i-1)+1);endformatshort;本题输入命令为Y=DEEuler(-v+u+1,0.1,0,0.5,1,[uv])（2）隐式Euler法functiony=DEimpEuler(f,h,a,b,y0,varvec)formatlong;N=(b-a)/h;y=zeros(N+1,1);y(1)=y0;x=a:h:b;var=findsym(f);fori=2:N+1fx=Funval(f,var(1),x(i));gx=y(i-1)+h*fx-varvec(2);y(i)=NewtonRoot(gx,-10,10,eps);endformatshort;（3）梯形法functiony=DEimpEuler1(f,h,a,b,y0)formatlong;N=(b-a)/h;y=zeros(N+1,1);x=a:h:b-h;y(1)=y0;var=findsym(f);yd=diff(f,var(4));fori=2:N+1fx=f-yd*var(4);dfy=subs(yd,findsym(yd),x(i-1));cx=subs(fx,findsym(fx),x(i-1));y(i)=(y(i-1)+h*cx)/(1-h*dfy);endformatshort;3.1.2单步法的局部截断误差解初值问题(1.1.1)的单步法可表示为(3.1.8)其中与有关，称为增量函数，当含有时，是隐式单步法，如(3.1.5)及(3.1.6)均为隐式单步法，而当不含时，则为显式单步法，它表示为(3.1.9)如Euler法(3.1.2)，.为讨论方便，我们只对显式单步法(3.1.9)给出局部截断误差概念.定义2.1设y(x)是初值问题(7.1.1)的精确解，记(3.1.10)称为显式单步法(3.1.9)在的局部截断误差.之所以称为局部截断误差，可理解为用公式(3.1.9)计算时，前面各步都没有误差，即，只考虑由计算到这一步的误差，此时由(3.1.10)有局部截断误差(3.1.10)实际上是将精确解代入(3.1.9)产生的公式误差，利用Taylor展开式可得到.例如对Euler法(3.1.2)有，故它表明Euler法(3.1.2)的局部截断误差为，称为局部截断误差主项.定义2.2设是初值问题(7.1.1)的精确解，若显式单步法(3.1.9)的局部截断误差，是展开式的最大整数，称为单步法(3.1.9)的阶，含的项称为局部截断误差主项.根据定义，Euler法(3.1.2)中的=1故此方法为一阶方法.对隐式单步法(3.1.8)也可类似求其局部截断误差和阶，如对后退Euler法(3.1.5)有局部截断误差故此方法的局部截断误差主项为，也是一阶方法.对梯形法(3.1.6)同样有它的局部误差主项为，方法是二阶的.3.1.3改进Euler法上述三种简单的单步法中，梯形法(3.1.6)为二阶方法，且局部截断误差最小，但方法是隐式的，计算要用迭代法。为避免迭代，可先用Euler法计算出的近似，将(3.1.6)改为(3.1.11)称为改进Euler法。改进的Euler法是二阶显式方法。即(3.1.12)右端已不含.可以证明，=2,故方法仍为二阶的，与梯形法一样，但用(3.1.11)计算不用迭代.例3.2用改进Euler法求例3.1的初值问题并与Euler法和梯形法比较误差的大小.解将改进Euler法用于例3.1的计算公式当n=0时，.其余结果见表3-2.表3-2改进Euler法及三种方法的误差比较从表3-2中看到改进Euler法的误差数量级与梯形法大致相同，而比Euler法小得多，它优于Euler法。附：改进的Euler法程序functiony=DEModifEuler(f,h,a,b,y0,varvec)formatlong;N=(b-a)/h;y=zeros(N+1,1);y(1)=y0;x=a:h:b;var=findsym(f);fori=2:N+1v1=Funval(f,varvec,[x(i-1)y(i-1)]);t=y(i-1)+h*v1;v2=Funval(f,varvec,[x(i)t]);y(i)=y(i-1)+h*(v1+v2)/2;endformatshort;讲解：求初值问题（1.1.1）的数值解就是在假定初值问题解存在唯一的前提下在给定区间上的一组离散点上求解析解的一组近似为此先要建立求数值解的计算公式，通常称为差分公式，简单的单步法就是由计算下一步，构造差分公式有三种方法，一是用均差（即差商）近似，二是用等价的积分方程（3.1.4）用数值积分方法，三是用函数的Taylor展开，其中Taylor展开最有普遍性，可以得到任何数值解的计算公式及其局部截断误差。计算公式是微分方程的一种近似，局部截断误差的概念就是刻划这种逼迫的好坏。当为微分方程的解，即，而用，定义局部截断误差，它表示用精确解代入计算公式（3.1.9）产生的公式误差为越大表明公式逼近微分方程的精度越高，因此就定义为公式的阶，通常的公式才能用于计算初值问题（1.1.1）的数值解。利用Taylor展开时，只要将的表达式在处展开成Taylor公式就可得到不同公式的局部截断误差。如3.1.2所给出的Euler法。后退Euler法和梯形法，它们只需用一元函数的Taylor展开，与后面第4章的多步法完全一致，而通常单步法（3.1.9）的一般情况则需要用二元函数的Taylor展开，才能得到公式的具体形式和局部截断误差。例如对改进Euler法，其局部截断误差由（3.1.12）可得要求出它的结果就要用到二元函数的Taylor展开，将在3.3节再作介绍。3.2Runge-Kutta方法3.2.1显式Runge-Kutta法的一般形式上节已给出与初值问题(1.1.1)等价的积分形式(3.2.1)只要对右端积分用不同的数值求积公式近似就可得到不同的求解初值问题(1.1.1)的数值方法，若用显式单步法(3.2.2)当，即数值求积用左矩形公式，它就是Euler法(3.1.2)，方法只有一阶，若取(3.2.3)就是改进Euler法，这时数值求积公式是梯形公式的一种近似，计算时要用二个右端函数f的值，但方法是二阶的。若要得到更高阶的公式，则求积分时必须用更多的f值，根据数值积分公式，可将(3.2.1)右端积分表示为注意，右端f中还不能直接得到，需要像改进Euler法(3.1.11)一样，用前面已算得的f值表示为(3.2.3)，一般情况可将(3.2.2)的表示为(3.2.4)其中这里均为待定常数，公式(3.2.2)，(3.2.4)称为r级的显式Runge-Kutta法，简称R-K方法。它每步计算r个f值(即)，而ki由前面(i-1)个已算出的表示，故公式是显式的。例如当r=2时，公式可表示为(3.2.5)其中.改进Euler法(3.1.11)就是一个二级显式R-K方法。参数取不同的值，可得到不同公式。3.2.2二、三级显式R-K方法对r=2的显式R-K方法(3.2.5)，要求选择参数，使公式的阶p尽量高，由局部截断误差定义(3.2.6)令，对(3.2.6)式在处按Taylor公式展开，由于将上述结果代入(3.2.6)得要使公式(3.2.5)具有的阶p=2,即，必须(3.2.7)即由此三式求的解不唯一.因r=2，故，于是有解(3.2.8)它表明使(3.2.5)具有二阶的方法很多，只要都可得到二阶R-K方法.若取，则，则得改进Euler法(3.1.11)，若取，则得，此时(3.2.5)为(3.2.9)其中称为中点公式.后退Euler法(3.1.11)及中点公式(3.2.9)是两个常用的二级R-K方法，注意二级R-K方法只能达到二阶，而不可能达到三阶.因为r=2只有4个参数，要达到p=3则在(3.2.6)的展开式中要增加3项，即增加三个方程.加上(3.2.7)的三个方程求4个待定参数是无解的.当然r=2,p=2的R-K方法(3.2.5)当取其他数时，也可得到其他公式，但系数较复杂，一般不再给出.对r=3的情形，要计算三个k值，即其中将按二元函数在处按Taylor公式展开，然后代入局部截断误差表