甘肃中医药大学《数值分析》课件-第9章 常微分方程数值解法.ppt

数值分析1第第9 9章章 常微分方程初值问题数值解法常微分方程初值问题数值解法/*Numerical Method for Ordinary Differential Equations*/9.1 9.1 引言引言9.2 9.2 简单的数值方法简单的数值方法 9.3 9.3 龙格龙格- -库塔方法库塔方法9.4 9.4 单步法的收敛性与稳定性单步法的收敛性与稳定性 9.5 9.5 线性多步法线性多步法 9.6 9.6 线性多步法的收敛性与稳定性线性多步法的收敛性与稳定性9.7 9.7 一阶方程组与刚性方程一阶方程组与刚性方程 数值分析29.1 引引 言言待求解的问题待求解的问题（1.1）（1.2）（1.3） 即存在实数即存在实数 ，使得，使得称为称为 的的利普希茨常数利普希茨常数（简称简称Lips.Lips.常数）常数）. .解的存在唯一性定理解的存在唯一性定理1 设设 在区域在区域 上连续，关于上连续，关于 满足利普希茨条件，满足利普希茨条件，/* Initial-Value Problem */:则对任意则对任意 ，常微分方程，常微分方程(1.1),(1.2)(1.1),(1.2)式式当当 时存在唯一的连续可微解时存在唯一的连续可微解 . .一阶常微分方程初值问题一阶常微分方程初值问题 数值分析3解对扰动的敏感性结论解对扰动的敏感性结论 定理定理2 设设 在区域在区域 （如定理（如定理1 1所定义）上连续，且所定义）上连续，且关于关于 满足利普希茨条件，设初值问题满足利普希茨条件，设初值问题的解为的解为 ，则，则 这个定理表明解对初值的敏感性，它与右端函数这个定理表明解对初值的敏感性，它与右端函数 有有关，当关，当 的的Lips.Lips.常数常数 比较小时，解对初值和右端函数相比较小时，解对初值和右端函数相对不敏感，可视为对不敏感，可视为好条件好条件. .若若 较大则可认为较大则可认为坏条件坏条件，即为，即为病态问题病态问题. .数值分析4 如果右端函数可导，由中值定理有如果右端函数可导，由中值定理有 若假定若假定 在域在域 内有界，设内有界，设 ，则，则 它表明它表明 满足满足利普希茨条件，且利普希茨条件，且 的大小反映了右端函的大小反映了右端函数数 关于关于 变化的快慢，刻画了初值问题变化的快慢，刻画了初值问题(1.1)(1.1)和和(1.2)(1.2)式是否式是否是好条件是好条件. .数值分析5解析解法解析解法:如何求解如何求解计算解函数计算解函数 y(x) 在一系列节点在一系列节点 a = x0 x1 xn= b 处的近似值处的近似值节点间距节点间距 为步长，通常采用为步长，通常采用等距节等距节点点，即取，即取 hi = h (常数常数)。

求解所有的常微分方程求解所有的常微分方程步进式步进式：根据已知的或已求出的节点上的函数值计算：根据已知的或已求出的节点上的函数值计算当前节点上的函数值，一步一步向前推进因此只需当前节点上的函数值，一步一步向前推进因此只需建立由已知的或已求出的节点上的函数值求当前节点建立由已知的或已求出的节点上的函数值求当前节点函数值的递推公式即可函数值的递推公式即可常数常数变易变易法、法、Lapalace变换等变换等分离变量分离变量法、变量法、变量代换代换、数值解法数值解法:数值分析6数值分析7数值分析8 一类是计算一类是计算 时只用到前一点的值时只用到前一点的值 ，称为称为单步法单步法. . 另一类是用到另一类是用到 前面前面 点的值点的值 ，称为称为 步法步法. . 其次其次: :要研究公式的局部截断误差和阶，数值解要研究公式的局部截断误差和阶，数值解 与与精确解精确解 的误差估计及收敛性，还有递推公式的计算的误差估计及收敛性，还有递推公式的计算稳定性等问题稳定性等问题. . 首先首先: :对方程对方程 离散化，建立求数值解的递推离散化，建立求数值解的递推公式公式. .数值分析9 欧拉（欧拉（Euler）方法是解初值问题的最简单的数值方法。

方法是解初值问题的最简单的数值方法9.2 简单的数值方法简单的数值方法 9.2.1 欧拉法与后退欧拉法欧拉法与后退欧拉法 初值问题初值问题 的解的解y=y(x)代表通过点代表通过点 的一条称之为微分方程的的一条称之为微分方程的积分积分曲线曲线 积分曲线上每一点积分曲线上每一点 的切线的斜率的切线的斜率 等于函数等于函数 在这点的值在这点的值 导数用均差近似，即导数用均差近似，即 Euler方法的导出方法的导出1然后用然后用 代替代替 ，即得，即得（2.1）数值分析10从初始点从初始点P0(即点即点(x0,y0)出发出发,作积作积分曲线分曲线y=y(x)在在P0点上切线点上切线 (其其斜率为斜率为 ),若初值若初值 已知，则依公式可逐步算出已知，则依公式可逐步算出(即点即点(x1,y1),得到得到y1作为作为y(x1)的近似值的近似值,如上图所示如上图所示过点过点(x0,y0),以以f(x0,y0)为斜率的切线方程为为斜率的切线方程为获得了获得了P1点的坐标点的坐标 当当 时时,得得 Euler法的求解过程法的求解过程与与x=x1直线相交于直线相交于P1点点数值分析11同样同样, 过点过点P1(x1,y1),作作当当 时时,得得 由此获得了由此获得了P2的坐标。

重复以上过程的坐标重复以上过程,就可获得一系列的点就可获得一系列的点:P1,P1,Pn对已求得点对已求得点 ，以，以 = 为为斜率作直线斜率作直线 当当 时时,得得取取获得了一条近似于曲线获得了一条近似于曲线y=y(x)的折线的折线这样这样,从从x0逐个算出逐个算出 对应的数值解对应的数值解 交直线交直线x=x2于于P2点点,斜率斜率=数值分析12 例例1 用欧拉公式用欧拉公式求解初值问题求解初值问题 （2.2）欧拉公式为欧拉公式为 取步长取步长 ，计算结果见表计算结果见表9-1.9-1. 解解 欧拉方法的精度很差欧拉方法的精度很差数值分析13误差的几何直观误差的几何直观 假设假设 ，即顶点即顶点 落在积分曲线落在积分曲线 上，上，那么，按欧拉方法做出的折线那么，按欧拉方法做出的折线 便是便是 过点过点 的切线（图的切线（图9-29-2）. .图图9-29-2将将 在在 处展开，则有处展开，则有 EulerEuler法的误差分析法的误差分析在在 的前提下，的前提下，误差误差 （2.3） 显然，这种近似有一定误差，显然，这种近似有一定误差，而且步长越大，误差越大，而且步长越大，误差越大，如何估计这种误差如何估计这种误差y(xn+1) yn+1 ？在在假假设设 yn = y(xn)，即即第第 n 步步计计算算是是精精确确的的前前提提下下，考考虑虑公公式式或或方方法法本本身身带带来来的的误误差差: Rn+1 = y(xn+1) yn+1 , 称称为为局局部部截截断断误误差差 /* local truncation error */。

数值分析14截断误差截断误差: 实际上，实际上，y(xn) yn， yn 也有误差，它对也有误差，它对yn+1的误的误差也有影响，见下图但这里不考虑此误差的影响，仅考虑差也有影响，见下图但这里不考虑此误差的影响，仅考虑方法或公式本身带来的误差，因此称为方法误差或截断误差方法或公式本身带来的误差，因此称为方法误差或截断误差局部截断误差的分析局部截断误差的分析：由于假设：由于假设yn = y(xn) ，即，即yn准确，因此准确，因此分析局部截断误差时将分析局部截断误差时将y(xn+1) 和和 yn+1都用点都用点xn上的信息来表上的信息来表示，工具：示，工具：Taylor展开 欧拉法的局部截断误差：欧拉法的局部截断误差：Tn+1 的的主项主项/* leading term */若若某某算算法法的的局局部部截截断断误误差差为为O(hp+1)，则则称称该该算算法法有有p 阶精度 欧拉法具有欧拉法具有 1 阶精度数值分析15（2.4） 对方程对方程 从从 到到 积分，得积分，得 左矩形公式左矩形公式 近似近似 再以再以 代替代替代替代替 也得到欧拉法也得到欧拉法 后退的后退的Euler法法 如果在（如果在（2.42.4）中右端积分用）中右端积分用右矩形公式右矩形公式 （2.5）近似，近似， Euler方法的导出方法的导出2显式显式隐式隐式数值分析16逐步显示化逐步显示化隐式公式不能直接求解，一般需要用隐式公式不能直接求解，一般需要用Euler显式公式得到显式公式得到初值，然后用初值，然后用Euler隐式公式迭代求解。

因此隐式公式较隐式公式迭代求解因此隐式公式较显式公式计算复杂，但显式公式计算复杂，但稳定性好稳定性好（后面分析）后面分析）2.6）迭代的收敛性迭代的收敛性由于由于 对对 满足利普希茨条件满足利普希茨条件由此可知，只要由此可知，只要 迭代法迭代法(2.6)(2.6)就收敛到解就收敛到解 . .迭代公式两端取极限即得迭代公式两端取极限即得数值分析17后退欧拉法的局部误差后退欧拉法的局部误差这里这里 ，是是1阶方法，局部截断误差主项为阶方法，局部截断误差主项为 . . 数值分析18 分别用分别用 代替代替 得到计算公式得到计算公式（2.7） 梯形方法梯形方法隐式单步法，用迭代法求解隐式单步法，用迭代法求解. . （2.8）迭代公式迭代公式9.2.2 梯形方法梯形方法数值分析19 迭代公式的收敛性迭代公式的收敛性如果选取如果选取 充分小，使得充分小，使得 则当则当 时有时有 ，迭代过程迭代过程(2.8)(2.8)是收敛的是收敛的. . 式中式中 为为 关于关于 的利普希茨常数的利普希茨常数. . 数值分析20梯形法的局部误差估计梯形法的局部误差估计梯形方法是梯形方法是2阶的，其局部误差主项为阶的，其局部误差主项为数值分析21迭代一两次就转入下一步的计算迭代一两次就转入下一步的计算. .预测值预测值精度可能很差精度可能很差 改进的欧拉公式改进的欧拉公式校正一次校正一次迭代一次得迭代一次得 ，这个结果称，这个结果称校正值校正值. .迭代计算量大迭代计算量大Step 1: 先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),(1nnnnyxfhyy+ += =+ +Step 2: 再将再将 代入代入隐式隐式梯形公式梯形公式1+ +ny9.2.3 改进欧拉公式改进欧拉公式梯形法的缺点梯形法的缺点改进改进数值分析22 预测预测- -校正系统：校正系统： 预测预测校正校正（2.9） 平均化形式平均化形式 几何解释几何解释xnxn+1ABPn+1=(A+B)/2欧拉法欧拉法后退欧拉法后退欧拉法梯形法梯形法数值分析23 例例2 用改进的欧拉方法求解初值问题（用改进的欧拉方法求解初值问题（2.22.2）. .（2.2） 解解 改进的欧拉公式为改进的欧拉公式为 仍取仍取 ，计算计算结果见表结果见表9-2. 9-2. 改进欧拉法明显改进欧拉法明显改善了精度改善了精度. . 数值分析24单步法的一般形式单步法的一般形式（2.10）其中多元函数其中多元函数 与与 有关，有关，当当 含有含有 时，方法是时，方法是隐式的隐式的，若不含，若不含 则为则为显式显式方法，方法，（2.11） 称为称为增量函数增量函数，?显式单步法显式单步法局部截断误差局部截断误差（1.1）（1.2）欧拉法欧拉法9.2.4 单步法的局部截断误差与阶单步法的局部截断误差与阶数值分析25 定义定义1 设设 是初值问题是初值问题(1.1)(1.1)，(1.2)(1.2)的准确解，的准确解，称称（2.12）为显式单步法（为显式单步法（2.112.11）的）的局部截断误差局部截断误差. . 之所以称为局部的，是假设在之所以称为局部的，是假设在 前各步没有误差前各步没有误差. . 当当 时，计算一步，则有时，计算一步，则有 定义定义2 设设 是初值问题是初值问题(1.1)(1.1)，(1.2)(1.2)的准确解，的准确解，若存在最大整数若存在最大整数 。