微分方程初值问题的数值解法

1、关于微分方程初值问题的数值解法第1页，共52页，2022年，5月20日，14点53分，星期三 包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中, 自变量的个数只有一个, 称为常微分方程。自变量的个数为两个或两个以上的微分方程叫偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数y及其各阶导数8.1 引 言8.1.1 微分方程知识回顾都是一次的,则称它是线性的,否则称为非线性的。 第2页，共52页，2022年，5月20日，14点53分，星期三 在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法： 一阶：可分离变量法、

2、齐次方程、一阶线性方程以及伯努利方程 高阶：可降阶方程、 常系数齐次线性方程的解法、常系数非齐次线性方程的解法 但能求解的常微分方程仍然是有限的，大多数的常微分方程是求不出解析解的。 这个一阶微分方程就不能用初等函数及其积分来表达它的解。 8.1.1 微分方程知识回顾 例如 第3页，共52页，2022年，5月20日，14点53分，星期三 从实际问题当中归纳出来的微分方程，通常主要依靠数值解法来解决。本章主要讨论一阶常微分方程初值问题 ( 8.1 ) 在区间 可以证明,如果函数在带形区域 R=axb,-y内连续，且关于y满足李普希兹(Lipschitz)条件，即存在常数L(它与x,y无关)使 对

3、R内任意两个 都成立,则方程( 8.1 )的解 在a, b上存在且唯一。 上的数值解法。 第4页，共52页，2022年，5月20日，14点53分，星期三 常微分方程初值问题(8.1)式的数值解法，首先要算出精确解y(x)在区间a,b上的一系列离散节点 处的函数值 的近似值: 相邻两个节点的间距 称为步长，本章总是假定h为定数，称为定步长，这时节点可表示为8.1.2 数值方法的基本思想 1、数值解法需要把连续性的问题加以离散化，从而求出离散节点处的数值解。第5页，共52页，2022年，5月20日，14点53分，星期三 描述这类算法，要求给出用已知信息 计算 的递推公式。 建立这类递推公式的基本方

4、法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法，对初值问题中的导数 进行不同的离散化处理。 2、数值解法的基本特点是采用“步进式”：即求解过程按照递推公式顺着节点排列的次序一步一步地向前推进。8.1.2 数值方法的基本思想第6页，共52页，2022年，5月20日，14点53分，星期三 递推公式通常有两类， 一类是计算yi+1时只用到xi+1, xi 和yi，即前一步的值，因此有了初值以后就可以逐步往下计算，此类方法称为单步法，其代表是龙格库塔法。 另一类是计算yi+1时，除用到xi+1,xi和yi以外，还要用到 ，即前面k步的值，此类方法称为多步法，其代表是亚当斯法。 8.1.2

5、数值方法的基本思想第7页，共52页，2022年，5月20日，14点53分，星期三一、Euler方法及其改进 将a , b n 等分, 记 微分法: 积分法: 积分项利用矩形公式计算 1. 显式Euler方法 ()第8页，共52页，2022年，5月20日，14点53分，星期三Taylor公式推导: 第9页，共52页，2022年，5月20日，14点53分，星期三Oyxy=y(x)(x1,y1)p1p0 x0 x1x2xixi+1xn-1xnpipi+1pn-1pnp1p2pipi+1pn-1pn切线p0p1的斜率为f (x0,y0)p2(x2,y2) 欧拉公式的几何意义：第10页，共52页，202

6、2年，5月20日，14点53分，星期三Euler法的求解过程是: 从初始点P0(即点(x0,y0)出发,作积分曲线y=y(x)在P0点上切线 (其斜率为 ),与直线x=x1相交于P1点(即点(x1,y1),得到y1作为y(x1)的近似值）这样就获得了P1点的坐标。 当 时,得 重复以上过程,就可获得一系列的点:p1,p2,pn，相应的可求出y1,y2,yn,取第11页，共52页，2022年，5月20日，14点53分，星期三从图形上看,就获得了一条近似于曲线y=y(x)的折线通常取 (常数),则Euler法的计算格式 i=0,1,n ( 7.2 ) 第12页，共52页，2022年，5月20日，1

7、4点53分，星期三2. 梯形法 称之为梯形公式.这是一个隐式公式,通常用迭代法求解.具体做法: 取 先用Euler法求出初值 ,即 ,将其代入梯形公式的右端,使之转化为显式公式,即 注: 当 f (x , y)关于y满足Lipschitz条件且步长h 满足 直至满足: 若采用梯形公式计算()中的积分项,则有类似地,可得() 第13页，共52页，2022年，5月20日，14点53分，星期三时,迭代格式 () 收敛 . 3. 改进的Euler方法 把Euler法作为预报(称为预估公式),把隐式的梯形公式作为校正(称为校正公式 ),则得改进的Euler方法:或也称为预估-校正法.第14页，共52页，

8、2022年，5月20日，14点53分，星期三有时为了方便,预估-校正格式也写成下面形式:第15页，共52页，2022年，5月20日，14点53分，星期三 改进的欧拉公式比欧拉公式精度高的原因是：改进欧拉公式用梯形面积代替曲边梯形面积，而欧拉公式用矩形面积代替曲边梯形面积。 数值积分的梯形公式比矩形公式的精度高。第16页，共52页，2022年，5月20日，14点53分，星期三二、单步法的局部截断误差及精度 Def 1: 先假设 ,再估计误差这种误差称为单步迭代法在 xk+1处的局部截断误差.Def 2: 若某种数值方法的局部截断误差为 ,则称该数值方法的精度为P 阶的.注: 通常情况下,P 越大

9、, h 越小,则截断误差越小,数值方法越精确. 第17页，共52页，2022年，5月20日，14点53分，星期三所以Euler方法为一阶方法. 而 设 10 .Euler方法是一阶方法.第18页，共52页，2022年，5月20日，14点53分，星期三20 . 梯形法是二阶方法.Taylor展开 第19页，共52页，2022年，5月20日，14点53分，星期三将 代入上式,得 而代入上式得:当h充分小时, 若 , 则可选取 h , 使得第20页，共52页，2022年，5月20日，14点53分，星期三故梯形法的精度为2 . 同样可以证明改进的Euler法也是二阶方法. 梯形法的局部截断误差为: 从

10、而第21页，共52页，2022年，5月20日，14点53分，星期三例1: 取步长 h = 2/10, 2/20, 2/30, 2/40, 分别用欧拉法、改进的欧拉法和梯形法求解.解: 记 f (x, y) = y x y2, xk= k h (k = 0, 1, 2, n )(1) . Euler法: yk+1 = yk + h( yk xk yk2 ) (k = 0, 1, ,n) y0 = 1当 h = 2/10时, n=10. 由Euler公式可得: k01234yk+11.21.38241.5061.535041.46503k56789yk+11.328771.170771.02113

11、0.891690.783788第22页，共52页，2022年，5月20日，14点53分，星期三(2) . 改进的Euler法: k01234yk+11.19121.343841.423481.419051.3473k56789yk+11.237261.114240.9941510.8847510.788666(3) . 梯形法(计算过程略) 第23页，共52页，2022年，5月20日，14点53分，星期三 n 10 20 30 40 h 0.2 0.1 0.0667 0.05误差 0.1059 0.0521 0.0342 0.0256Euler法误差:改进的Euler法误差: n 10 20

12、30 40 h 0.2 0.1 0.0667 0.05误差 0.0123 0.0026 0.0011 5.9612e-004第24页，共52页，2022年，5月20日，14点53分，星期三预-校方法, h=0.2时误差最大值: 0.0123欧拉方法, h=0.2时误差最大值: 0.1059解析解:第25页，共52页，2022年，5月20日，14点53分，星期三三、Runge -Kutta 方法1、Taylor 级数法 设初值问题 有解 y (x), 由Tayler公式得:令当 时, 有 . 此时为 p 阶Taylor方法. p=1时即为Euler公式.称之为Taylor级数法. 其中例2: 取

13、步长 h = 0.1, 用一阶、二阶和四阶Taylor方法求解下列初值问题第26页，共52页，2022年，5月20日，14点53分，星期三解: (1) 一阶Taylor法k01234yk+11.11.2211.370081.557791.80046(2) 二阶Taylor法k01234yk+11.111.246891.421751.652631.97088第27页，共52页，2022年，5月20日，14点53分，星期三(3) 四阶Taylor法k01234yk+11.11111.249961.428481.666441.99942第28页，共52页，2022年，5月20日，14点53分，星期三

14、记由得称为xk , xk+1上的平均斜率. 故2、Runge -Kutta方法只要对K*提供不同的算法, 就会得出不同的计算公式. 如取则得改进的Euler公式, 它是利用xk , xk+1两点的斜率值K1 , K2 的算术平均值作为K*, 精度比Euler法高.则得Euler公式; 取第29页，共52页，2022年，5月20日，14点53分，星期三Runge-Kutta法的基本思想: 设法在xk , xk+1内多预报几个点的斜率, 再将它们的加权平均值作为平均斜率K*一般显式Runge-Kutta公式为: 其中 为待定参数, 且 . 称为r 级Runge-Kutta方法计算公式. 第30页，

15、共52页，2022年，5月20日，14点53分，星期三即可得 p 个方程, 从而确定出待定参数. 代入表达式即可得到计算公式. 如果要求两个表达式的前p+1项完全重合, 即局部截断误差达到 , 则称式为 p 阶 r 级的Runge-Kutta方法. 常用的是 r =2,3,4 级的R-K方法, 且适当选取参数使得 p = r . 如要求: 注: 式中待定参数的确定: 先将式右端在(xk , yk ) 处展成h的幂级数(即将 yk+1 展成 h 的幂级数); 再将 y(xk+1) 作Taylor 级数展开;最后比较两式中hk ( k=0,1,2,)的系数, 以确定出所有待定参数.第31页，共52

16、页，2022年，5月20日，14点53分，星期三Runge-Kutta方法的推导(以r =2为例): 当r =2 时记第32页，共52页，2022年，5月20日，14点53分，星期三则又第33页，共52页，2022年，5月20日，14点53分，星期三这是一个四个参数三个方程的非线性方程组. 它有一个自由度. 称满足上述方程组的一族公式为二级二阶Runge-Kutta方法.为使局部截断误差为 ,比较上述两式右端同次幂系数,应取第34页，共52页，2022年，5月20日，14点53分，星期三(1) 常用的二阶Runge-Kutta方法: 预估-校正算法(2)中间点方法 第35页，共52页，2022

17、年，5月20日，14点53分，星期三注: 二级Runge-Kutta方法的精度最高是二阶的,不可能达到三阶. 要提高计算方法的阶, 就必须增加预报点. 常用的三阶Runge-Kutta方法(r =3): (1) Heun (休恩)方法 (3) 三阶Kutta方法 第36页，共52页，2022年，5月20日，14点53分，星期三(1) 三阶Heun方法 标准(经典)四阶Runge-Kutta方法 (2) 常用的四阶Runge-Kutta方法(r =4): 第37页，共52页，2022年，5月20日，14点53分，星期三(2) 称为Gill(吉尔)方法 注: 从理论上讲,可以构造任意高阶的计算方法

18、. 但事实上, 精度的阶数与预报点的个数之间并非等量关系. 预报点的个数 r123456789r 10精度的阶数123445667 r -2一般情况下, 四阶Runge-Kutta方法已可满足精度要求. 第38页，共52页，2022年，5月20日，14点53分，星期三例3: 用经典Runge-Kutta方法求解下列初值问题(取 h = 0.1) 解: 标准Runge-Kutta公式为: 计算结果见下表. 为比较在相同计算量条件下近似解的精度, 表中列出了Euler法(h =0.025)和改进的Euler法(h=0.05)在相应节点上的计算结果. 第39页，共52页，2022年，5月20日，14

19、点53分，星期三xiEuler法h=0.025改进Euler法h=0.05经典R-K法h=0.1准确解0.11.1114391.1153801.1155121.1155130.21.2552091.2639141.2642081.2642080.31.4346671.4490891.4495761.4495760.41.6535171.6747561.6754731.6754740.51.9158491.9451711.9461621.9461640.62.2261782.2650402.2663542.2663560.72.5894852.6395612.6412552.6412580.83

20、.0112713.0744793.0766193.0766230.93.4976063.5761443.5788043.5788091.04.0551924.1515734.1548394.154845注: 用表中每种方法计算 yi 都需要计算四次 f 的值, 即它们的计算量基本相等. 第40页，共52页，2022年，5月20日，14点53分，星期三四、单步法的进一步讨论收敛性、相容性与稳定性注: 由定义可知, 数值方法的收敛性并不涉及计算过程的舍入误差, 只与方法的截断误差有关. 若格式收敛, 则整体截断误差必趋于零. Def : (整体截断误差) 称 为某一数值方法在点 xk 处的整体截断

21、误差. 它不仅与 xk 有关, 也与xk -1, xk -2 , x1, x0 有关. 则称该单步法收敛. Def : 对满足解存在唯一性条件的初值问题(1), 如果一个显式单步法(3)产生的近似解对于任一固定的 , 均有 1. 收敛性第41页，共52页，2022年，5月20日，14点53分，星期三由于 , 且 关于 y 满足Lipschitz条件, 得 则存在常数 c 0 使得 且单步法中函数 关于 y 满足Lipschitz条件, 则 定理1: 若初值问题的一个单步法的局部截断误差为 记 证: 由局部截断误差的定义知第42页，共52页，2022年，5月20日，14点53分，星期三故 从而有

22、 故 若 y(x0) = y0, 则e0=0, 由不等式 得 第43页，共52页，2022年，5月20日，14点53分，星期三设单步法为 注: 定理表明, 数值方法的整体截断误差比局部截断误差低一阶. 收敛的方法至少是一阶方法. 在该定义条件下, Euler方法是一阶的, 预估-校正方法是二阶. 当f (x , y)关于 y 也满足Lipschitz条件, r 级Runge-Kutta方法中的 关于 y 也满足Lipschitz条件, 故定理中的条件得到满足, 解的收敛性得到保证.由于R n, h0(h0), 且 xn为任意点, 故该式相当于用近似方程 当x = xn+1固定时, , 所以有

23、2. 相容性第44页，共52页，2022年，5月20日，14点53分，星期三通过在 x = xn 处求解近似方程而获得原方程的近似解. 因此,必须要求当h0 时, 近似方程应逼近于原方程.来代替 因此, 要使 h0 时, 近似方程的极限状态为原微分方程, 需且只需下列极限成立:由于 由于假设 是连续函数, 故上式可表示为 Def : 如果当h0时, 近似方程逼近微分方程, 则称数值公式与原微分方程相容. 相容的充要条件: 第45页，共52页，2022年，5月20日，14点53分，星期三事实上: Remark: 可以证明若单步法的阶大于或等于1, 则单步法与微分方程相容; 反之, 如果单步法与微

24、分方程相容, 且 关于 y 满足Lipschitz条件, 则单步法至少为一阶方法. (h0) (1) 若单步法的阶大于或等于1, 由 知 即单步法与微分方程相容. 故有 (2) 如果单步法与微分方程相容, 且 关于y 满足Lipschitz条件, 则第46页，共52页，2022年，5月20日，14点53分，星期三关于 单步法的收敛性以及收敛性定理都是在计算过程中无任何舍入误差的前提条件下建立的, 但在实际计算时通常会有舍入误差及其积累, 数值求解微分方程的过程是一个递推公式, 必须考 即与微分方程相容的单步法至少为一阶方法. Remark: 在定理条件下, Euler方法、预估-校正方法以及Runge-Kutta方法都与原微分方程相容.中连续, 且关于变量y 满足Lipschitz条件, 则单步法收敛的充要条件为相容性条件成立. Th1. 设增量函数 在区域 3. 稳定性第47页，共52页，2022年，5月20日，14点53分，星期三 如果数值方法在计算过程中舍入误差的积累越来越大, 得不到有效控制, 则称其是不稳定的; 反之如果计算结果对初始数据的误差及计算过程中的误差不敏感, 即舍入误差不增长, 则称相应的算法是稳定的. 数值方法的稳定性有各种定义, 这里仅考虑