第七章(微分方程-3)

第七章常微分方程的数值解法1引言1、一阶常微分方程初值问题（微分方程加初值条件）(1)特征：为已知的二元函数，只有一个变元，一阶导数，要求满足微分方程且过点的曲线，一个未知函数（微分方程的解函数）。 2、解的存在性 若连续且对满足Lipschitz条件，即存在常数，使对，有，则初值问题(1)的解存在唯一，本章均假设对满足Lipschitz条件。3、数值解法 问题：求精确解（解析表达式）极其困难，实际应用中只要求数值解。数值解：在取一系列等距离散节点，步长，求的近似值。方法：建立的递推公式，从而按节点排列顺序（步进式）。2Euler方法1、 显示、隐式Euler法与梯形法 对（1）的方程两端由到定积分得：从而有：（*）（积分中值定理）（1）若取，得（左矩形公式）由（*）得：， ,则有： 称式为解初值问题的Euler法.（显格式）几何意义：初值问题(1)的解曲线过点，从出发，以为斜率作一段直线，与直线交点于，显然有，再从出发，以为斜率作直线推进到上一点，其余类推，这样得到解曲线的一条近似曲线，它就是折线（某点导数值的几何意义为该点切线的斜率）。（2）若取，得（右矩形公式） 类似可得：称式为隐式Euler法（向后Euler法，隐格式）（3）取（梯形公式）则得：称式为梯形法.（隐格式）例1 用Euler法、隐式Euler法、梯形法及改进Euler法求解 ，取，计算到，并与精确解比较.解：由于,Euler法： 时，.隐式Euler法：解出 当时，梯形法： 解出 当时，改进Euler法: 即 当时， 精确解：，(梯形法效果最好，改进Euler较好：具有相同的误差数量级，其它不好!)Euler法隐式Euler法梯形法改进Euler法精确解0.11.000 0001.009 0911.004 7621.005 0001.004 8370.21.010 0001.026 4461.018 5941.019 0251.019 7310.31.029 0001.051 3151.040 6331.041 2181.040 8180.41.056 1001.083 0141.070 0971.070 8021.070 3200.51.090 4901.120 9221.106 2781.107 0761.106 5312、隐式法的计算（Euler法及梯形法）方法1：显示化（对线性时，以为未知量的一元线性方程，见上例） 方法2：迭代法（对非线性时，可看作一个关于方程，利用迭代法求解）Euler法： （*1）（以为未知量的一元非线性方程）,n （*2）梯形法：（*3），（*4）计算步骤：（1）初值 ；迭代 （当）（2）；迭代收敛性：只要步长足够小，就可保证迭代收敛。（a）当（*2）收敛 （b）当（*4）收敛.证明：（b）迭代序列收敛当 有；由（*3）、（*4）式知： （满足李氏条件），则当，（）时，有，从而收敛得证。3、 改进Euler法（梯形法隐式迭代式仅限迭代一次）,为避免迭代，可先用Euler法计算出的近似 显式预报，隐式校正 或 （6） 称以上两式为改进Euler法（校正就是迭代一次） 右端已不含，显式方法。例2：用梯形法的迭代格式求的数值解，计算到解：梯形公式迭代式，解为：3、 Euler法的局部截断误差定义1 在的前提下，称为在的局部截断误差。定义2若一种数值方法的局部截断误差：，则称这种方法是精度为阶的。（含的项称为局部截断误差主项.）注：按某种方法由算出这一步的误差。（1）Euler法（*）（*）由(*)（*）式有：根据定义2，Euler法中的故此方法为一阶方法.（2）隐式Euler法（分析时隐式公式右边的）（*）由（*）及（*）式：故局部截断误差主项为，也是一阶方法.（3）梯形法 类似有，方法是二阶的.（4）改进Euler法 则 =,方法是二阶的.3 Runge-Kutta方法一、 显式 Runge-Kutta法的一般形式 一般形式：二、 2级显式R-K方法（r=2）设想构造RK公式：（*）试确定，使上式为二阶格式，即：假设，对(7.3.6)式在处按Taylor公式展开，由于 ；将上述结果代入(*)得：要使公式为二阶,必须：（三个方程四个求知数，无穷多解） （1）取，有（改进欧拉法） （2）取，有（中点法）取其他数时，也可得到其他公式，但系数较复杂，一般不再给出.三、 四阶R-K方法及步长的自动选择 类似可得到，经典的四阶R-K方法是：、 它的局部截断误差，故p=4,这是最常用的四阶R-K方法，数学库中都有用此方法求解初值问题的软件.这种方法的优点是精度较高，缺点是每步要算4个右端函数值，计算量较大.例7.3 用经典四阶R-K方法解例7.1的初值问题 ，取，计算到，并与改进Euler法、梯形法在处比较其误差大小.解 用四阶R-K方法公式(7.3.12)，此处，于是当n=0时 于是，按公式(7.3.12)可算出此方法误差：改进Euler法误差：梯形法误差：可见四阶R-K方法的精度比二阶方法高得多.用四阶R-K方法求解初值问题(1)精度较高，但要从理论上给出误差的估计式则比较困难.那么应如何判断计算结果的精度以及如何选择合适的步长h?通常是通过不同步长在计算机上的计算结果近似估计.设在处的值,当时，的近似为，于是由四阶R-K方法有若以为步长，计算两步到，则有 于是得： 即或(7.3.13)它给出了误差的近似估计.如果 (为给定精度)，则认为以为步长的计算结果满足精度要求，若，则还可放大步长.因此(7.3.13)提供了自动选择步长的方法.4 线性多步法一、 线性多步法的一般公式前面给出了求解初值问题(1)的单步法，其特点是计算时只用到 的值（要计算多点处的函数值，计算量较大），此时 的值均已算出.如果在计算时除用的值外，还用到的值，这就是多步法.若记,为步长， ，则线性多步法可表示为（一般形式）：其中为常数，若，称(7.5.1)为线性步法.计算时用到前面已算出的个值.当时，为显式方法，当则称为隐式多步法.隐式方法与梯形方法一样，计算时要用迭代法求多步法(7.5.1)的局部截断误差定义也与单步法类似.线性多步法在处的局部截断误差定义为：假定时若，称其为是p阶的。二、 Adams显式与隐式方法 （即）的步法称为 Adams 方法，当时为 Adams 显式方法，当时，称为Adams隐式方法.对初值问题(1)的方程两端从到积分得选取不同的节点，对上式中被积函数进行不同的多项式插值，从而得到不同的的线性多步法。1、 Adams四步显式公式（外插公式）选 ，四节点作三次多项式插值：从而取与为插值节点做的插值多项式：用代替，令，并有用令代替，从而得到四步Admas（外插）公式则有：（注：记） 步的Adams方法计算时必须先用其他方法求出前面个初值才能按给定公式算出后面各点的值，它每步只需计算一个新的f值，计算量少，但改变步长时前面的也要跟着重算，不如单步法简便.局部截断误差：当时，由上知：从而： 保号令，并有，从而因此是四阶方法2、 Adams三步隐式公式（内插公式）选 ，四节点作三次多项式插值：类似可得：用代替则有：局部截断误差：，三步四阶隐式方法三、Adams预测-校正方法Adams显式方法计算简单，但精度比隐式方法差，而隐式方法由于每步要做迭代，计算不方便.为了避免迭代，通常可将同阶的显式Adams方法与隐式Adams方法结合，组成预测-校正方法.以四阶方法为例，可用显式方法(7.5.8)计算初始近似，这个步骤称为预测(Predictor)，以P表示，接着计算f值(Evaluation),，这个步骤用E表示，然后用隐式公式(7.5.10)计算，称为校正(Corrector)，以C表示，最后再计算，为下一步计算做准备.整个算法如下：上式称为四阶Adams预测-校正方法(PECE).例7.6 用四阶显式Adams方法及四阶隐式Adams方法解初值问题,步长h=0.1用到的初始值由精确解计算得到.解 本题直接由公式(7.5.8)及(7.5.10)计算得到.对于显式方法，将直接代入式(7.5.8)得到其中.对于隐式方法，由式(7.5.10)可得到 直接求出，而不用迭代，得到 计算结果如表所示.5 单步法的收敛性与绝对稳定性一、 单步法的收敛性定义4.1 设y(x)是初值问题(1)的精确解，是单步法(7.3.2)在处产生的近似解，若则称方法(7.3.2)产生的数值解收敛于.，实际上，定义中是一固定点，当h0时n，n不是固定的.因显然方法收敛，则在固定点处的整体误差，当p1时.下面定理给出方法(7.3.2)收敛的条件.定理4.1设初值问题(1)的单步法(7.3.2)是p阶方法(p1)，且函数对y满足Lipschitz条件，即存在常数L0,使对,均有则方法(7.3.2)收敛，且.定理证明略.可见3.二、 绝对稳定性用单步法(7.3.2)求数值解，由于原始数据及计算过程舍入误差影响，实际得到的不是而是，其中是误差，再计算下一步得到以Euler法为例，若令,则(7.4.1)如果，则从计算到误差不增长，它是稳定的.但如果条件不满足就不稳定.例7.4 y=-100y,y(0)=1,精确解为，用Euler法求解得 若取h=0.025，则,当，而，显然计算是不稳定的.如果用后退Euler法(7.2.5)解此例，仍取h=0.025,则,即显然当，计算是稳定的.由此看到稳定性与方法有关，也与有关，在此例中.在研究方法的稳定性时，通常不必对一般的f(x,y)进行讨论，而只针对模型方程(7.4.2)这里可能为复数.规定是因为时微分方程(7.4.2)本身是不稳定的，而讨论数值方法(7.3.2)的稳定性，必须在微分方程本身稳定的前提下进行.另一方面，对初值问题(1)，若将f(x,y)在处线性展开，可得于是方程(1)可近似表示为 它表明用模型方程(7.4.2)是合理的，至于模型方程(7.4.2)中所以用复数是因为初值问题(1)如果是方程组，即，则是(mm)阶矩阵，其特征值可能是复数.当然对单个方程，就是实数，此时只要规定0即可.用单步法(7.3.2)解模型方程(7.4.2)可得到(7.4.3)其中依赖所选方法，如用Euler法则(7.4.4)此时由(7.4.1)看到误差方程也为，与(7.4.4)是一样的.因此对一般单步法(7.3.2)误差方程也与(7.4.3)一致.下面再考虑二阶R-K方法有 对四阶R-K方法，可得 定义4.2将单步法(7.3.2)用于解模型方程(7.4.2)，若得到(7.4.3)中的 则称方法是绝对稳定的.在复平面上复变量满足 的区域，称为方法(7.3.2)的绝对稳定域，它与实轴的交点称为绝对稳定区间.例如对Euler法， 在复平面上是以(-1，0)为圆心，以1为半径的单位圆域内部，当为实数时，则得绝对稳定区间为，因0，故有.在例7.4中 时方法稳定，而例中取h=0.025故不稳定.对后退Euler法(7.2.5)， 因0，故，其绝对稳定域是以(1，0)为圆心的单位圆外部，绝对稳定区间为，即对任何h0方法都是绝对稳定的.二阶R-K方法的绝对稳定区间为.三阶R-K方法的绝对稳定区间为.四阶R-K方法的绝对稳定区间为.例7.5 用经典四阶R-K方法计算初值问题 步长取h=0.1及0.2,给出计算误差并分析其稳定性.解 本题直接按R-K方法(7.3.12)的公式计算.因精确解为，其计算误差如表所示.从计算结果看到，h=0.2时误差很大，这是由于在=-20,h=0.2时h=-4，而四阶R-K方法的绝对稳定区间为-2.785,0,故h=0.2时计算不稳定，误差很大.而h=0.1时=-2，其值在绝对稳定区间-2.785,0内，计算稳定，故结果是可靠的.讲解：由于微分方程初值问题数值解公式求出的解是一个逐次递推的过程，因此原始数据误差及计算过程舍入误差对解的影响就是数值方法绝对稳定性研究的问题，如果由计算误差不增长，方法就是绝对稳定的。为使问题得到简化通常就是将方法用于解模型方程（7.4.2），对于单步法得到的差分方程为，由于模型方程的，代入Euler法，得，对二阶R-K方法，例如，用改进Euler法 于是对三阶R-K方法有对四阶R-K方法有 只要方法，就是绝对稳定的，这时的值当n增大式是减少的，故计算稳定。这时舍入误差影响可忽略不计，而当，则增大，方法不稳定，计算结果是不可靠的。因此用显式单步法必须使，也就是步长选择要满足这一要求。对于隐式的梯形公式将模型方程，即代入得于是注意，于是有：，对成立。这就表明对任意步长h，梯形法都是绝对稳定的。16