《计算机数值方法教学课件》第二章-常微分方程数值解法

1、第二章 常微分方程数值解法Chapter 2 Numerical Solution of Ordinary Differential Equation(s)第二章 常微分方程数值解法Chapter 2 Numer22.1 引 言 Fm假设牵引力F为恒定值为了确定待定常数，可以给定初始条件：或者给定边界条件：假设牵引力不恒定呢？求速度22.1 引 言 Fm假设牵引力F为恒定值为了确定待定常3虽然求解微分方程有许多解析方法，但解析方法只能够求解一些特殊类型的方程。还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值。利用计算机解微分方程主要使用数值方法。 解析方法与数值方法3虽然求解微

2、分方程有许多解析方法，但解析方法只能够求解一些特4主要研究对象：初值问题求数值解：求y(x)在离散数据点xk处的近似值yk 。 y=y(x)xyx0=ax1x2x3xk-1xn-1xn=b xk4主要研究对象：初值问题求数值解：求y(x)在离散数据点xk5则称 f (x,y) 对y 满足李普希兹条件，此时初值问题在a,b上存在唯一的连续可微的解。定理1：设 f (x,y) 是定义在区域 G=(x,y)|axb, yR上的连续函数，若存在正的常数 L 使：(Lipschitz)条件有解条件：5则称 f (x,y) 对y 满足李普希兹条件，此时初值问题6则Lipschitz条件成立：在 f (x,

3、y) 对y可微的情况下，若偏导数有界：有解条件的判断：6则Lipschitz条件成立：在 f (x,y) 对y可微7定理2：如果f(x,y)在G=(x,y)|axb, yR上满足Lipschitz条件，则初值问题是适定的。适定性：指初值问题中，初始值y0及微分方程的右端函数f (x,y) 有微小变化时，只能引起解的微小变化。适定性条件：7定理2：如果f(x,y)在G=(x,y)|axb, 8在解的存在区间 a, b上取n + 1个节点 这里把 称为由xi到xi+1的步长 一般取成等间距的： 求解方法：步进法（分为单步法和多步法） 数值方法的基本思想8在解的存在区间 a, b上取n + 1个节点

4、 这里把 9本章规定：在 处初值问题的理论解用 表示，数值解法的近似解用 表示。记 ，它和 是不同的，后者等于 。9本章规定：在 处初值问题的理论解用 102.2 几种简单的数值方法 （一） 欧拉（Euler）法102.2 几种简单的数值方法 （一） 欧拉（Euler）11、泰勒公式解释 、求导的两点公式解释、积分公式解释欧拉公式的的分析解释11、泰勒公式解释 欧拉公式的的分析解释12泰勒公式解释其中：可以得到：12泰勒公式解释其中：可以得到：13求导的两点公式解释可以得到：13求导的两点公式解释可以得到：14对微分方程(1.1)两端从进行积分积分公式解释14对微分方程(1.1)两端从进行积分

5、积分公式解释15右端积分用左矩形数值求积公式：即：15右端积分用左矩形数值求积公式：即：16欧拉公式的的几何描述yxx0 x1 x2 x3 x4y=y(x)16欧拉公式的的几何描述yxx0 x1 x17例题1：（取步长h=0.1）用Euler方法求满足条件 的y(t) 数值解。 解：17例题1：（取步长h=0.1）用Euler方法求满足条件 18ntnyny(tn)y(tn)-yn01234101.01.11.21.31.43.00.00.271830.684761.276983.0935515.398240.00.345920.866641.607223.6203618.6830.00.07

6、4090.181880.330240.526813.28486数值解列表为 18ntnyny(tn)y(tn)-yn01.00.00.019欧拉方法的误差估计 通过数值方法进行计算时，考虑每一步产生的误差，从x0开始一步步累积到xn，称 为该数值方法在xn点处的整体截断误差，该误差与xn 及之前的各步计算误差都有关系。19欧拉方法的误差估计 通过数值方法进行计算时，考虑每一步产20欧拉方法的误差估计 为了简化分析，着重分析xn点单步计算产生的误差，即把xn点之前的计算当作无误差： 称该误差为数值方法在xn+1点处的局部截断误差。局部截断误差的第一个非零项为局部截断误差主项。20欧拉方法的误差估

7、计 为了简化分析，着重分析xn点单步计算21欧拉方法的误差估计 如果求解公式的局部截断误差为R (h)=O(hp+1)，则称该求解公式具有p阶精度，称该方法为p阶方法。定义：欧拉方法：具有1阶精度。21欧拉方法的误差估计 如果求解公式的局部截断误差为R (h22（二）向后欧拉法 （1）方法 其公式为：22（二）向后欧拉法 （1）方法 其公式为：23（2）局部截断误差23（2）局部截断误差24例题2：用向后Euler法解初值问题 向后Euler法的公式为 解：x0=0, y0=1, 取h=0.124例题2：用向后Euler法解初值问题 向后Euler法的25方法比较及推广： Euler方法 显式

8、公式向后Euler方法 隐式公式 解一个非线性方程 难求解 显式和隐式相结合 隐式的显化 25方法比较及推广： Euler方法 显式公式26计算公式为： 由显式得到，称为预估值；yn+1由隐式得到，称为校正值。这种求解方法统称为预估校正方法。其求解过程为：26计算公式为： 由显式得到，称为预估值；这种求27例3 用预估校正方法求解微分方程（取h=0.1）：解： 27例3 用预估校正方法求解微分方程（取h=0.1）：解28（三）梯形公式 28（三）梯形公式 29梯形公式局部截断误差29梯形公式局部截断误差30预估-校正方法：称为改进的Euler求解公式或改进Euler法。30预估-校正方法：称为

9、改进的Euler求解公式或改进Eul31为了表示方便，可以改写为：31为了表示方便，可以改写为：32（四）欧拉方法的收敛性分析由初值问题的单步法产生的近似解 ，如果对于任一固定的 均有 ，则称该方法是收敛的。定义：局部截断误差：若初值问题的一个单步法的局部截断误差为： 定理：且增量函数关于y满足Lipschitz条件，则整体截断误差： 整体截断误差比局部截断误差低1阶32（四）欧拉方法的收敛性分析由初值问题的单步法产生的近似解33证明：存在常数c，使得33证明：存在常数c，使得34当 固定时，所以34当 固定时，所以35（五）欧拉方法的稳定性分析问题：常微分方程初值问题数值解的每步计算都是在前

10、一步计算的结果上进行的，所以必须考虑前面的误差对以后计算结果的影响，误差的积累会不会盖过真解呢？选用代表性试验方程： y=y (Re()0，在计算yn时引入了误差n 。若这个误差在计算后面的yn+k(k=1,2,.)中所引的误差n+k按绝对值均不增加，就说这个数值方法对于这个步长h和复数是绝对稳定的。 若在区域R内数值方法是绝对稳定的，则R为该数值方法的绝对稳定区域(区间)。 35（五）欧拉方法的稳定性分析问题：常微分方程初值问题数值解36把欧拉方法用于试验方程： y=y误差方程 ：要求误差不增加： O-2-1Re(h)Im(h)36把欧拉方法用于试验方程： y=y误差方程 ：要求误差37把向

11、后欧拉方法用于试验方程： y=yO21Re(h)Im(h)要求误差不增加： 可见隐式的向后Euler方法比显式的Euler方法的绝对稳定域要大得多。同阶精度的数值方法，往往隐式方法比显式方法的绝对稳定域大。 37把向后欧拉方法用于试验方程： y=yO21Re(h38例4：以 y=y为例判断梯形公式的稳定性：解出：这种稳定性称为无条件稳定！38例4：以 y=y为例判断梯形公式的稳定性：解出：这种39例5：用Euler法、向后Euler法、改进的欧拉法（梯形公式）解初值问题 取步长h=0.2，小数点后至少保留4位。 解：Euler法向后Euler法39例5：用Euler法、向后Euler法、改进的

12、欧拉法（梯40改进的Euler法梯形公式法40改进的Euler法梯形公式法41h=0.23.00003.30083.46593.55653.60623.63351.01.21.41.61.83.0真解改进的欧拉法梯形公式向后EulerEulerxi 3.0000 3.0000 3.0000 3.0000 3.4000 3.2500 3.3077 3.2800 3.5600 3.4063 3.4734 3.4424 3.6240 3.5039 3.5626 3.5366 3.6496 3.5649 3.6106 3.5912 3.6598 3.6031 3.6365 3.6229如何得到高精度的

13、求解公式？41h=0.23.00001.0真解改进的欧拉法梯形公式向后422.3 Runge-Kutta法 如何构造更高精度的求解公式？一种思路是：对微分方程右端积分采用高次多项式近似，如用二次多项式近似可得到422.3 Runge-Kutta法 如何构造更高精度43该公式称为Simpson公式。该公式的局部截断误差为：43该公式称为Simpson公式。44Taylor展开：另一种得到高阶方法的想法是直接利用泰勒级数展开。如果能计算得到y的高阶微商,则可写出r阶的计算公式44Taylor展开：另一种得到高阶方法的想法是直接利用泰勒45例题：取h=0.1，求解初值问题 解： 一阶公式：二阶公式：

14、四阶公式：45例题：取h=0.1，求解初值问题 解： 一阶公式：46真解： p0.10.20.30.40.511.100001.221001.370081.557791.8004621.110001.246891.421741.662621.9208741.111101.249661.428481.666451.99942y(xn)1.111111.250021.428571.666672.00000 xn问题： 1、求微商麻烦； 2、计算量大。46真解： p0.10.20.30.40.511.100047R-K方法是通过对不同点上的函数值做线性组合，构造近似公式，把近似公式和泰勒展开相比较，

15、使前面的若干项相吻合，从而使近似公式达到一定的阶数。 问题：如何组合函数值？其中，47R-K方法是通过对不同点上的函数值做线性组合，构造近似公48选择参数ci, ai, bij 的原则是，要求的Runge-Kutta法的右端项在(xn,yn)处泰勒展开后按h的幂次重新整理得到的结果与微分方程的解y(xn+1)在xn处的Taylor展开式 有尽可能多的项重合。 48选择参数ci, ai, bij 的原则是，要求的Rung49以计算两个函数值为例说明选择c1、c2、a2、b21 阶最高为了计算49以计算两个函数值为例说明选择c1、c2、a2、b21 阶50于是若要求局部截断误差达到,则要求有 选取

16、50于是若要求局部截断误差达到,则要求有 选取51二阶R-K方法的公式 它的局部阶段误差为O(h3)。这是计算两次函数值的情况下所能达到的最高阶。 51二阶R-K方法的公式 它的局部阶段误差为O(h3)。这是52中间点法： 二阶休恩(Heun)法： 52中间点法： 二阶休恩(Heun)法： 53经典的R-K方法是一个四阶的方法,公式为 1、一步法，可以自开始；特点：2、精度较高；3、便于改变步长；4、计算量较小。53经典的R-K方法是一个四阶的方法,公式为 1、一步法，可54例2：用二阶R-K方法和四阶R-K方法(h=0.1)求解解：二阶R-K方法四阶R-K方法54例2：用二阶R-K方法和四阶

17、R-K方法(h=0.1)求解55xn二阶R-K误差四阶R-K误差y=1-e-5x20.10.20.30.40.50.60.70.80.91.0 0.0500 0.0012 0.0488 0.0000 0.0488 0.1830 0.0017 0.1813 0.0000 0.1813 0.3627 0.0004 0.3624 0.0000 0.3624 0.5475 -0.0031 0.5507 -0.0000 0.5507 0.7059 -0.0076 0.7134 -0.0001 0.7135 0.8235 -0.0112 0.8346 -0.0001 0.8347 0.9012 -0.01

18、25 0.9135 -0.0002 0.9137 0.9476 -0.0116 0.9590 -0.0002 0.9592 0.9733 -0.0093 0.9823 -0.0002 0.9826 0.9866 -0.0066 0.9931 -0.0002 0.9933说明：1、高阶方法比低阶方法精度高，但需要y(x)光滑性好；2、步长h并非越小越好；步长选择从大到小。55xn二阶R-K误差四阶R-K误差y=1-e-5x20.156数学家Butcher于1965年证明了运算量与可以达到的最高精度阶数关系如下：每步计算kR的个数2345678910可以达到最高精度阶数23445666R-2可见，

19、通常只选用R 4的公式是合适的。56数学家Butcher于1965年证明了运算量与可以达到的57自动选步长方法 设在xn处之前误差不计：误差：57自动选步长方法 设在xn处之前误差不计：误差：58隐式R-K方法类似于显式R-K公式，稍加改变，就得到隐式R-K方法。显式公式中对系数求和的上限是i-1。而在隐式公式中对系数求和的上限是L，需要用迭代法求出Ki。推导隐式公式的思路和方法与显式R-K法类似。通常，同级的隐式公式可以获得比显式公式更高的阶。58隐式R-K方法类似于显式R-K公式，稍加改变，就得到隐式59通常，同级的隐式公式获得比显式公式更高的阶。常用的隐式R-K法有:1级2阶中点公式 :

20、2级2阶梯形公式： 2级4阶R-K公式： 59通常，同级的隐式公式获得比显式公式更高的阶。常用的隐式R60隐式龙格-库塔法其中2阶方法 的绝对稳定区域为0ReImg而显式 1 4 阶方法的绝对稳定区域为k=1k=2k=3k=4-1-2-3-123ReImg无条件稳定龙格-库塔法稳定区域60隐式龙格-库塔法其中2阶方法 612.4 线性多步法 单步法的优缺点：简单，可以自开始；提高精度时需要增加中间函数值；没有充分利用前几步得到的信息。多步法： yn-p , yn-p+1 , , yn-1 , yn yn +1612.4 线性多步法 单步法的优缺点：简单，可以自开始62yn-p, yn-p+1,

21、 yn-1, yn yn +1的方法考虑如下形式的求解公式此计算公式为线性的，所以称为线性多步法。当 时公式含有 ，这时公式是隐式的，而当 时公式是显式的62yn-p, yn-p+1, yn-1, yn 63常微分初值问题与积分公式等价基本思想：对f(x,y(x)进行多项式插值，利用插值公式计算右边积分，可以得到常微分初值问题求解公式。63常微分初值问题与积分公式等价基本思想：对f(x,y(x)64三次多项式插值求积分插值余项为令 ，假设已知 ，取插值节点 ，对函数 作三次Lagrange多项式插值64三次多项式插值求积分插值余项为令 65四阶Adams内插公式内插公式局部截断误差： 65四阶

22、Adams内插公式内插公式局部截断误差： 66令 ，假设已知 ，取插值节点 ，对函数 作三次Lagrange多项式插值外插公式局部截断误差： 66令 67内插公式外插公式通常计算公式为（预估校正系统）： 特点：精度高，但须与同精度的单步法配合使用。67内插公式外插公式通常计算公式为（预估校正系统）： 特点68例3：用Adams预测校正系统求初值问题其中，E=200, L=3, =100, =50解：由于Adams预测校正系统为四阶公式，选用四阶R-K公式计算开始值。a=t0=0; b=0.05I0=0;h=0.002f(t,I)=(200-50I3-100I)/368例3：用Adams预测校正

23、系统求初值问题其中，E=2069再利用Adams预测校正系统n45625t0.0080.0100.0120.050I0.4646910.5592060.644933.1.16486169再利用Adams预测校正系统n45625t0.008702.5 常微分方程组、高阶微分方程及边值问题的数值解 考虑如下含二个未知函数的方程组：一、常微分方程组的数值解 则：可依照常微分方程数值方法来求解。令：702.5 常微分方程组、高阶微分方程及边值问题的数值解71则其改进的Euler格式具有预估形式校正公式为如对于方程组71则其改进的Euler格式具有预估形式校正公式为如对于方程72两方程微分方程组的四阶R-K公式72两方程微分方程组的四阶R-K公式73二、高阶微分方程的数值解 对于高阶微分方程，可以把它化成微分方程