微分方程数值解

1、第6章 常微分方程初值问题数值解法本章探讨常微分方程特解的常用数值方法的构造和原理，主要介绍求常微分方程初值问题的常用方法和有关知识。重点论述Euler方法、Runge-Kutta方法和线性多步法的原理、构造、局部截断误差和稳定性等内容。6.1 实际案例工程技术里某些振动问题可以表示为单摆的运动，其运动规律的微分方程为：怎样求出其特解？该微分方程不能用通常的解析方法来求解！怎样解不能用解析方法求解的微分方程特解问题是本章要解决的问题。6.2 问题的描述和基本概念1、常微分方程初值问题l 一般形式 （6.1）式中已知，称为初值条件。l 初值问题的数值方法和数值解求函数在若干离散点上的近似值的方法

2、称为初值问题的数值方法，而称为初值问题的数值解。2. 建立数值解法的思想与方法微分方程初值问题的数值解法是用离散化方法将初值问题化为差分方程后再求解的方式。设区间a,b上的一组节点为距离称为步长。求数值解一般是从开使逐次顺序求出。初值问题的解法有单步法和多步法两种：l 单步法：计算时只用到一个值；l 多步法：计算时要用多个值。数值解法还有显格式和隐格式之分。l 微分方程离散化方法主要有数值微分法，数值积分法和Taylor展开法。l 初值问题化为差分方程的方法1. 用离散方法去掉方程中的导数得到近似离散化方程；2. 在近似离散化方程中用代替；3. 在近似离散化方程中将近似号“”用等号“=”代替。

3、1) 数值微分法由初值问题（6.1）有，用数值微分的2点前差公式代替，得近似离散化方程记，做，“”，得差分方程写容易计算的形式 （6.2）（Euler公式）由初值条件及式（6.2）可求出（6.1）的数值解。公式（6.2）是显式单步法。2）数值积分法在上对两边取定积分，得对右端积分采用梯形公式（数值积分公式）得近似离散化方程：于是得到求初值问题（6.1）的梯形方法该公式是隐式单步法。3）Taylor展开法因为初值问题中函数是已知函数，由，可以计算，于是有函数在处的Taylor展式取上式右端前若干项，得近似离散化方程。例如取前两项有于是又得到Euler公式：6.3数值解法的误差、阶与绝对稳定性单步

4、法数学描述为显式： （6.4）隐式：（6.5）其中称为增量函数，与有关。l 显式单步法的一些概念定义6.1 设是问题（6.1）的解，是经过式（6.4）求出的的计算解，则称为单步法（6.4）在节点的整体截断误差，而称 （6.6）为在点的局部截断误差。表示解在的值，是准确值，没有误差；表示由数值解公式得出的近似值，是数值解，有截断误差；表示用计算机计算给出的计算解，有舍入误差。l 局部截断误差的理解假设在计算时没有误差（）下，由式（6.4）计算出的（）与的误差。整体截断误差还要加上与的误差。考察初值问题解法的优劣，引入阶的概念。定义6.2 如果初值问题（6.1）对某种数值解法的局部截断误差为则称该

5、方法具有p阶精度或该方法是p阶方法。方法的阶越高，方法越好。l 主局部截断误差或局部截断误差的主项如果某方法是p阶方法，按可展为则称为主局部截断误差。在同阶方法中，主局部截断误差越小，方法越好。例6.1常微分方程初值问题的单步法为试求其局部截断误差主项并回答它是几阶精度的？解 该单步公式的局部截断误差是故局部截断误差主项是，方法是一阶的。l 求阶p的另一方法因为去掉下标，有若将在x点展开有则知该方法的阶是p。例如，对Euler方法，有那么将在x点展开，有故有因此，Euler方法是一阶方法。定义6.3 设用某种数值方法求初值问题（6.1）在任意节点的数值解时，满足，则称该数值方法是绝对稳定的。这

6、里是计算机计算时得出的计算解的舍入误差，。通常用试验方程 （为复数）来讨论求解初值问题的数值方法绝对稳定性；对具体初值问题，可取。稳定性常与步长h有关。定义6.4 某方法在试验方程中绝对稳定的复平面范围称为该方法的绝对稳定域，它与复平面实轴的交称为该方法的绝对稳定区间。绝对稳定域包含复平面左半平面的方法称为是A-稳定的。绝对稳定域越大，对应的方法绝对稳定性越好。6.4 Euler方法的有关问题用Euler公式，求解初值问题（6.1）数值解的方法称为Euler方法。1 )Euler方法的几何意义Euler方法常称为折线法。2) Euler方法的误差设为的计算解，满足其中。则有Euler方法的局部

7、截断误差；Euler方法的总体截断误差由，,有因为对任意m都有，可得由k的任意性，可知Euler方法的总体截断误差为说明当h足够小时，由Euler方法计算所得数值解可以很好地逼近准确解，从而Euler方法是收敛的。3) Euler方法稳定性将Euler公式用于试验方程，得到 （6.8）设计算时有舍入误差，则有得要想，只须，因此Euler方法在时是绝对稳定的，其绝对稳定域为复平面上以-1为中心的单位圆盘。由图形可知Euler方法的绝对稳定区间为。-1-2图6.3图6.4若指定是负实数，则有当步长h满足时，可保证Euler方法的计算绝对稳定。例6.2 ，取用Euler方法求其数值解并与其准确解作比

8、较。 解 本题准确解为，Euler法计算公式为直接计算结果与对应准确值的误差有如下表xk0.0250.050.0750.1yk-1.52.25-3.3755.0625y(xk)- yk1.58-2.243.37-5.06可见计算结果波动大，不稳定，计算结果不可用。原因：，h的稳定性范围，而本题的h=0.0250.02，故计算结果不稳定。4) 改进的Euler方法易证梯形方法（6.3）的绝对稳定域为复平面上的部分，它是复平面的整个左半平面（图6.4）。因此是A-稳定的。其绝对稳定域比Euler方法大，是二阶方法，故梯形方法比Euler方法好，但美中不足的是梯形方法是隐式方法。为避免迭代，实用中采

9、取如下两次计算格式由来计算（改进的Euler方法）：预测 （6.9）校正 （6.10）预测校正公式。易证它也是二阶方法。6.5 Runge-Kutta方法理论上，由Taylor展开法可以构造出解初值问题（6.1）的高阶数值方法，但这涉及到要计算的高阶导数，因此很不方便。本节的Runge-Kutta方法从函数本身着手，避开计算的高阶导数，来构造高阶数值方法。基本思想将微分方程初值问题转化为积分方程问题，再对积分方程中的定积分使用待定的m点插值型求积公式构造高阶的函数展开模式以获得高阶数值方法。1. 构造原理将初值问题的方程在上积分对定积分用m个点的插值型求积公式取，并对中的做展开，得解初值问题的

10、Runge-Kutta一般公式形式（6.11）利用Taylor展开式对这些参数适当选择就可以构造高阶方法了。Runge-Kutta方法（6.11）的增量函数是要构造p阶数值方法，可选择参数使局部截断误差即可。2 构造过程以来说明Runge-Kutta方法的构造方法和过程，对一般的Runge-Kutta方法可类似处理。的Runge-Kutta公式为 （6.12）式中 ，。其增量函数为引进符号， ，由，可得，在处做Taylor展开，有对在做二元Taylor展开，有由，有选 （6.13）有局部截断误差，这样从式（6.13）中解出，代入到式（6.12）中，即得到一组二阶Runge-Kutta公式。由于

11、式（6.13）中有3个方程4个参数，若令一个参数为自由变量，如取，则式（6.13）的解为，取不同的可得出不同的二阶Runge-Kutta公式。如取时，得到改进的Euler公式3. Runge-Kutta方法的阶与级的关系计算函数值的次数称为级。l Runge-Kutta方法的级与阶对应关系：计算的次数 1 2 3 4 5 6 7对应方法的最高阶 1 2 3 4 4 5 6l 经典Runge-Kutta公式 （6.14）例 设初值问题为分别用Euler方法（），改进Euler方法（）和经典Runge-Kutta方法（）计算。解 Euler方法计算格式（）为改进的Euler方法计算格式（）为经典R

12、unge-Kutta方法计算格式（）为它们的初值，计算结果及准确解列于下表Euler方法改进Euler方法经典R-K法 000000.10.096 3120.095 1230.095 162500.095 162 580.20.183 3480.181 1930.181 269 100.181 269 250.30.262 0010.259 0850.259 181 580.259 181 780.40.333 0790.329 5630.329 679 710.329 679 95例6.5给定初值问题1）分析求解公式的局部截断误差，指出它是几阶公式；2）证明用上面求解公式计算初值问题的数值解

13、成立极限 本题中的节点是等距节点，h为步长，n为由节点分割区间a,b的份数。解 由题意有1）局部截断误差将在点做Taylor展开到项，将在点做二元Taylor展开到项，则有得所用公式是2阶的。2）显然所给初值问题的准确解为。由给出的数值解计算公式有故式中6.6线性多步法线性多步法的一般计算格式为 （6.15）式中均为常数,，为等距节点，步长为h。若不同时为零，计算一个需要用到的前个值,式（6.15 ）称为线性步法。当n>1时就称为线性多步方法。构造解初值问题（6.1）的线性多步法有基于数值积分方法和Taylor展开方法两种手段。1. 基于数值积分的构造方法1) 构造原理与分析将初值问题在

14、积分，m为正整数，得然后选择的不同插值函数代替并积分之，即可得到计算公式。l Adams方法给出构造过程。已初值问题的解在个点上的数值解为，记，考虑数表构造的次Lagrange插值多项式，在式（6.16）中取=1，有因为式中 。记，有于是可得n步Adams显示（外插）公式 （6.17）其局部截断误差为式中的定积分计算利用了变换。常用的Adams公式1、4阶Adams显式公式局部截断误差 2、4阶Adams隐式公式局部截断误差 例6.7 构造3阶Adams显示公式。解 取三个相邻的等距节点对应的数据做二次Lagrange插值多项式注意到，作变换，有故有于是得到3阶Adams外推公式其局部截断误差

15、为2. 基于Taylor展开的构造方法l 先给出线性多步法的计算公式模式l 然后对局部截断误差表达式在处作Taylor展开以确定公式的系数。例 6.8 设求初值问题的线性3步公式具有如下形式h为步长，试求系数a,b,c使该公式的阶数尽可能高，并写出其局部截断误差。解 局部截断误差为因为有3个待定系数，选择截断误差的前3项的系数为零，得关于待定系数的线性方程组该方程组有唯一解 ，用此值代入第四项的系数，有故有当时，所给公式的阶数达到最高，其值为3，对应的局部截断误差为6.7步长的自动选取问题数值解法步长h不变，会产生求出的各节点上的数值的解具有不同的精度的情况！要使所有数值解都满足给定的精度，必

16、须选用较小的步长h，但这又会产生计算量的增加和舍入误差的累积！因此最好的方式是步长可以变化。通常计算机上多采用根据误差估计来自动选步长值的方式处理这个问题。l 自动选步长方法的原理和过程设某种方法的阶为p>0, 设为步长。记表示由计算一步算出的在的数值解，有如将步长折半为，用表示由计算两个半步长算出的在的数值解，有假设，则有整理后，有记，则在给定精度后，若有，就有说明用步长折半前后在的两次数值解之差的大小可判别在处数值解的精度！步长自动选取的算法1. 先选定某个步长，然后计算；2. 若，（说明精度足够，但为防止精度太高导致求数值解的计算量太大，做步长扩大处理）将加倍为，计算，若还有，再将加倍为，类似做下去直到某个，则选定步长；3. 若，（说明精度不够，做步长缩小处理，以提高精度）将折半为，计算，若还有，再将折半为，类似做下去知道某个，则选定步长。6.8 一阶微