《《常微分方程课程设计》指导书 3》由会员分享，可在线阅读，更多相关《《常微分方程课程设计》指导书 3（23页珍藏版）》请在金锄头文库上搜索。
1、第第 3 3 章章 数值算法之一：单步法数值算法之一：单步法接下来的章节，我们将详细讲解常微分方程的数值求解方法，特别是编程 实现方程的初值问题，这也将成为解决实际问题的必要基础和课程设计的主要 内容。 本章重点讲解一阶常微分方程的数值算法中最简单的一类方法：单步法。 首先介绍 Euler 法、后退 Euler 法与梯形法，并分析单步法的局部截断误差， 然后给出了改进的 Euler 法。3.1 简单的单步法及基本概念3.1.13.1.1 EulerEuler 法、后退法、后退 EulerEuler 法与梯形法法与梯形法求初值问题(1.2.1)的一种最简单方法是将节点的导数用差商nx()ny x
2、代替，于是(1.2.1)的方程可近似写成()()nny xhy x h(3.1.1)从出发，由(3.1.1)求得再将 代入(3.1.1)右端，得到的近似，一般写成(3.1.2)该数值解法称为解微分方程初值问题的解微分方程初值问题的 EulerEuler 法法. .Euler 法的几何意义如图 3-1 所示。初值问题(1.2.1)的解曲线 y=y(x)过点，从出发，以为斜率作一段直线，与直线交点于，显然有，再从出发，以为斜率作直线推进到上一点，其余类推，这样得到解曲线的一条近似曲线，它就是折线.图图 3-13-1 EulerEuler 法的几何意义显示法的几何意义显示Euler 法也可利用的 T
3、aylor 展开式得到，由(3.1.3)略去余项，以，就得到近似计算公式(3.1.2).另外，还可对(1.2.1)的方程两端由到积分得(3.1.4)若右端积分用左矩形公式，用，则得(3.1.2).如果在(3.1.4)的积分中用右矩形公式，则得(3.1.5)该算法称为后退后退( (隐式隐式)Euler)Euler 法法。若在(3.1.4)的积分中用梯形公式，则得(3.1.6)该算法称为梯形方法梯形方法。上述三个公式(3.1.2)，(3.1.5)及(3.1.6)都是由计算，这种只用前这种只用前 一步即可算出一步即可算出的公式，我们称之为的公式，我们称之为单步法单步法，其中(3.1.2)可由逐次求出
4、的值，称为显式方法，而(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.13.1 用 Euler 法、隐式 Euler
5、法、梯形法解取 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-13-1 例例 3.13.1 的三种方法及精确解的计算结果的三种方法及精确解的计算结果对梯形法，计算公式为解得当 n=0 时，.其余 n=1,2,3,4 的计算结果见表7-1.本题的精确解为，表 3-1 列出三种方法及精确解的计算结果。附：具体三种算法程序如下：附：具体三种算法程序如下： 首先定义
6、一阶方程右端函数 f 如下function Y=f(x,y) Y=-y+x+1; （1）Euler 法 function y = DEEuler(f, h,a,b,y0,varvec)%这是Euler法的函数命令%一阶常微分方程的一般表达式的右端函数：f 这里可用 f=inline% 自变量下限 a；上限 b% 函数初值 y0% 积分步长 h% 常微分方程的变量组 varvecformat long;%数据显示方式，不影响计算和存储方式，%是指小数点后 15 位数字表示 N = (b-a)/h;y = zeros(N+1,1);x = a:h:b;y(1) = y0;for i=2:N+1y(
7、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); endformat short;本题输入命令为 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;var = findsym(
8、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; （3）梯形法 function y = DEimpEuler1(f, h,a,b,y0)format long;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);for i=2:N+1fx = f - yd*var(4);dfy = su
9、bs(yd,findsym(yd),x(i-1);cx = subs(fx,findsym(fx),x(i-1);y(i) = (y(i-1)+h*cx)/(1-h*dfy);endformat short;3.1.23.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.12
10、.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.22.2 设是初值问题(7.1.1)的精确解，若显式单步法(3.1.9)的局部截断误差， 是展开式的最
11、大整数，称 为单步法(3.1.9)的阶，含的项称为局部截断误差主项.根据定义，EulerEuler 法法(3.1.2)(3.1.2)中的中的 =1=1 故此方法为一阶方法故此方法为一阶方法. .对隐式单步法(3.1.8)也可类似求其局部截断误差和阶，如对后退 Euler 法 (3.1.5)有局部截断误差故此方法的局部截断误差主项为，也是一阶方法.对梯形法(3.1.6)同样有它的局部误差主项为，方法是二阶的.3.1.3 改进改进 Euler 法法上述三种简单的单步法中，梯形法(3.1.6)为二阶方法，且局部截断误差最 小，但方法是隐式的，计算要用迭代法。为避免迭代，可先用 Euler 法计算出
12、的近似，将(3.1.6)改为(3.1.11)称为改进 Euler 法。改进的改进的 EulerEuler 法是二阶显式方法。法是二阶显式方法。即(3.1.12)右端已不含.可以证明， =2,故方法仍为二阶的，与梯形法一样，但用(3.1.11)计算不用迭代.例例 3.23.2 用改进 Euler 法求例 3.1 的初值问题并与 Euler 法和梯形法比较误 差的大小.解解 将改进 Euler 法用于例 3.1 的计算公式当 n=0 时，.其余结果见表 3-2.表表 3-23-2 改进改进 EulerEuler 法及三种方法的误差比较法及三种方法的误差比较从表 3-2 中看到改进 Euler 法的
13、误差数量级与梯形法大致相同，而比 Euler 法 小得多，它优于 Euler 法。附：改进的 Euler 法程序function y = DEModifEuler(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=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;en
14、dformat short;讲解：讲解： 求初值问题（1.1.1）的数值解就是在假定初值问题解存在唯一的前提下在 给定区间上的一组离散点上求解析解 的一组近似为此先要建立求数值解的计算公式，通常 称为差分公式，简单的单步法就是由计算下一步，构造差分公式有三种方 法，一是用均差（即差商）近似，二是用等价的积分方程（3.1.4）用数值积分 方法，三是用函数的 Taylor 展开，其中 Taylor 展开最有普遍性，可以得到任何数值解的计算公式及其局部截断误差。计算公式是微分方程的一种 近似，局部截断误差的概念就是刻划这种逼迫的好坏。当为微分方程的解，即，而用，定义局部截断误差，它表示用精确解代入计
15、算公式（3.1.9）产生的公式误差为越大表明公式逼近微分方程的精度越高，因此 就定义为公式的阶，通常的公式才能用于计 算初值问题（1.1.1）的数值解。利用 Taylor 展开时，只要将 的表达式在 处展开成 Taylor 公式就可得到不同公式的局部截断误差。如 3.1.2 所给出的 Euler 法。后退 Euler 法和梯形法，它们只需用一元函数的 Taylor 展开，与后面第 4 章的多步法完全一致，而通常单步法（3.1.9）的一般情况则需要用二元函数的 Taylor 展开，才能得到公式的具体形式和局部截断误差。例如对改进 Euler 法， 其局部截断误差由（3.1.12）可得 要求出它的
16、结果就要用到二元函数的 Taylor 展开，将在 3.3 节再作 介绍。3.2 Runge-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 级的显
