计算方法word教案第六章常微分方程的数值解法.doc

第六章 常微分方程的数值解法6.0 引言6.1 算法构造的主要途径 6.2 Runge-Kutta Method算法6.3 线性多步法6.4 线性多步法的一般形式6.5 算法的稳定性、收敛性6.6 小结376.0 引 言1主要考虑如下的一阶常微分方程初值问题的求解:微分方程的解就是求一个函数，使得该函数满足微分方程并且符合初值条件。2. 例如微分方程:初始条件.可得一阶常微分方程的初始问题。显然函数y(x)=x2-4x满足以上条件,因而是该初始问题的微分方程的解。3. 但是,只有一些特殊类型的微分方程问题能够得到用解析表达式表示的函数解，而大量的微分方程问题很难得到其解 析解，有的甚至无法用解析表达式来表示。因此，只能依赖于数值方法去获得微分方程的数值解。4 微分方程的数值解：设微分方程问题的解y(x)的存在区间是a,b，初始点x0=a，将a,b进行划分得一系列节点x0 , x1 ,.,xn，其中a= x0 x1 xn =b。y(x)的解析表达式不容易得到或根本无法得到，我们用数值方法求得y(x)在每个节点xk的近似值y(xk)，即，这样y0 , y1 ,.,yn称为微分方程的数值解。如图所示:ab x0 x1 x2 . xn-1 xn6.1 算法构造的主要途径1 欧拉公式1.1 构造的思想：微分方程初值问题: 利用差商代替一阶导数，即，则。于是，可求出y(x1)的近似值y1，同样地，可利用x1处的微分方程可得：一般地，利用在xn处的微分方程可得： 此式称为欧拉公式。1.2 几何意义：对于微分方程y=2(x+1),其通解是y=(x+1)2+c,是一个曲线族，当给定初值条件y(0)=2，其特解为y=(x+1)2+1。如图所示：由y(0)=2，过该曲线上一点(0,2)作曲线的切线，其斜率，切线为：,因此可计算出y1,如此，可根据：，故欧拉法又称欧拉折线法。1.3 算例：例：在区间0,1上以h=0.1为步长，用欧拉法求初值问题的数值解?解：该方程为贝努利方程，其精确解为欧拉公式的具体形式为取步长 , 计算结果见下表： 欧拉法准确解近似解与准确解比较，欧拉法的结果大致只有两位有效数字，而预估校正法的结果则 有3位有效数字。1.4 Euler法的特点和误差迭代格式：特点：(1)单步方法；(2)显式格式；(3)局部截断误差为 。局部截断误差：当时，由按照欧拉方法计算来的的误差称为局部截断误差。即，是局部截断误差。如：欧拉法得：因此，局部截断误差是。2改进Euler法2.1方法构造在微分方程初值问题,对其从到进行定积分得:将右端的定积分用梯形公式来进行近似计算。用和来分别代替和得计算格式:这就是改进欧拉方法。2.2 显式格式和隐式格式在欧拉式中每一步计算已知,直接用格式可以计算出,此类格式称为显式格式。而在改进欧拉方法中在每一步计算中是未知,待求的,未知量在中这是一个方程,如是非线形或超越函数,此方程是无法直接解出来(要依靠迭代法才能解出)。这类格式称为隐式格式。2.3 算例例: 用改进欧拉方法求解。解：解得:注意：由于,是线形函数可以从隐式格式中解出。问题的精确解是欧拉方法误差例：在区间0,1上以h=0.1为步长，分别用欧拉法与预估 校正法求初值问题的数值解?解：该方程为贝努利方程，其精确解为欧拉公式的具体形式为预估校正公式的具体形式为取步长 , 计算结果见下表： 欧拉法预估-校正法准确解近似解与准确解比较，欧拉法的结果大致只有两位有效数字，而预估校正法的结果则 有3位有效数字。3预测校正方法由于改进Euler法是隐式格式，无法从格式中直接求出必须要解方程。下面用预测校正方法来求隐式格式中的。 预测值: 校正值: 此式相当于对隐式格式求时采用迭代的方法,用欧拉格式得到的作为初始值迭代公式迭代一次而已,此公式代入后得:如改写成平均的形式为:6.2 龙格-库塔法Runge-Kutta Method1 龙格-库塔法的思想1.1 考虑微分方程的初值问题:根据微分中值定理有：，其中01。于是即我们称为y(x)在区间xk, xk+1上的平均斜率，记作K，其中，是存在但是未知的。因此，如何对平均斜率K进行近似计算，相应地就得到一种近似公式，或称为微分方程的一种计算格式。1.2 例如：用f(xk,yk)作为平均斜率K的近似值就得到欧拉格式：。用作为平均斜率K的近似值就得到比欧拉格式高一阶精度的格式，即，改进欧拉格式的预测-校正方法:1.3 启发(Motivative)?能否在二维平面中xxk, xk+1, yyk, yk+1上多找一些f(x,y)点，在这些点上作函数值的平均数并以此作为平均斜率K的近似值。由于自由度的增加，使得的p能够提高，从而达到提高精度的目标，这就是龙格-库塔法的基本思想。1.4 R-K的一般形式 在二维平面区域上取N个点，得到公式:，其中Kj是二元函数f(x,y)在这N个点上的值。其中是待定常数。2 二阶龙格-库塔法2.1计算格式: 应用Taylor公式求参数：。将f(x,y)在(xk , yk)上展开又因 以上两式代入中，又由得比较两边得四个未知量,三个方程,有无穷多组解，但局部截断误差均为,是二阶精度的。例如:(1) 取,则,就是改进欧拉公式的预测-校正方法(或二阶龙格-库塔法)。(2) 取，则，则,得二阶中点法：3 三阶龙格库塔方法3.1 计算格式: 参数有,共有八个参数。将f和均在点Tayor展开,有 利用,按展开式代入(高于的不计)得:以上展开式代入中比较得:(1)(2)(3)(4)(5)(6)(7)(8)上述八个方程有六个独立方程。八个未知量,六个方程有无穷多组解,但其截断误差均为,具有三阶精度。例如:，得:此方法称为三阶Kutta方法。例如:,得:此方法称为三阶Heun法。4 四阶龙格库塔方法4.1 显式四阶龙格-库塔方法的一般形式是:利用与三阶龙格-库塔类似的处理方法，经过推导，可得古典四阶龙格-库塔算法:其特点是：截断误差为，阶数为四阶精度，为显式四阶龙格库塔方法。计算均要计算四个f的值。一般一阶常微分方程初值问题的解均用四阶龙格-库塔方法进行计算，其精度满足实际问题的精度要求。数值例子：例：初值问题用四阶古典RungeKutta方法，。0.21.18322931.18321600.41.34168031.34164080.61.48328381.48323970.81.61251721.61245151.01.73214631.7320508例： 6.3 线性多步法1 单步方法和多步方法单步方法：欧拉方法、改进欧拉方法、龙格-库塔方法均是单步方法，即在每一步计算时，只要前面一个值已知的条件下就可以计算出。单步方法特点：(1)可以自成系统进行直接计算，因为初始条件只有一个已知，由可以计算，不必借助于其它方法，这种单步方法是自开始的。(2)如果格式简单如欧拉方法，则只有一阶精度，如果提高精度，则计算很复杂，如RungeKutta方法。(3)公式的构造推导也很复杂。多步方法：利用前面已知计算出来的，由已经计算好的个值来计算，这样可以提高算法的精度，该方法称为多步方法，利用k个值计算，称为k步方法。多步方法的特点：(1)因初始条件只有一个，运用多步方法要借助高阶的单步方法来开始。例如，已知用单步的四阶Runge-Kutta方法计算，再计算，再由计算，用单步方法有后运用四阶的四步方法，由计算；由计算；由计算；一直下去， 可以用多步方法，并且始终达到四阶精度。(2)多步方法比较简单，只要在这四个点的函数值的线性组合，而且每步中后三个函数值下一步还可使用。2 显式Adams方法：2.1 构造的思想考虑微分方程初值问题将微分方程在上积分，2.2 显式Adams方法若已知来计算，简记，用的拉格朗日插值多项式代替f， 截断部分，用等距步长，上面积分很简单，得到的方法就是显式四阶Adams方法。以上可以看到该方法的局部截断误差是因而是四阶精度的。例如：解：取，首先用四阶RungeKutta方法来起步，计算出，,下面不必用RungeKutta方法，而开始用四阶Adams方法。(1)、求(2)、求只要补算(3)、求只要补算现列表看用Adams方法求出的误差，精解为0.81.61142311.61245151.0210-31.01.72984031.73205082.2110-31.21.84066161.84390893.2410-32.3 隐式Adams方法用作为插值结点，由于也是插值结点，必带来从而导致是隐式格式。用插值多项式来代替积分中的得：截掉得近似公式：得：，从而得四阶隐式Adams方法。因，而是未知的，故这是隐式格式。隐式格式的解法用预测-校正法：用显式格式作为预测值，再用隐式格式来校正。预测值：校正值：6.4 线性多步法的一般形式微分方程初值问题: 1 k步线性法的一般结构是常数线性多步法 , k步法不同时为零 k步法 ,隐式结构;否则为显式结构。2 k步p阶线性法的推导由: 定义算子：由Taylor展开右端有代入上式：其中，因此，若有次可微，令,则因而算子对应的线性多步法为p阶k步方法。3 算例例：考虑2步法记于是：故一般的二步法为：同时可以算出：当时，上述公式为三阶的。当时，上述公式为四阶的。 6.5 收敛性与稳定性1 定义收敛性： 若某算法对于任意固定的 x = xn= x0 + n h，当 h0 ( 同时 n ) 时,有 yn y( xn)，则称该算法是收敛的。 1.2 算例例：对于初值问题 考察欧拉显式格式的收敛性。解：该问题的精确解为 欧拉公式为即：对任意固定的 x = xn = n h ，有由 2 稳定性 例：考察初值问题 在区间0, 0.5上的解。分别用欧拉显式、隐式格式和改进的欧拉格式计算数值解。节点 xi欧拉显式欧拉隐式改进欧拉法精确解0.00.10.20.30.40.51.0000-2.00004.0000-8.0000 1.6000101-3.20001011.00002.500010-16.250010-21.562510-23.906310-39.765610-41.00002.50006.25001.56261013.90631019.76561011.00004.978710-22.478810-31.234110-46.144210-63.059010-7 误差包括截断误差(算法理论误差)和 舍入误差两个部分。后者由计算机字长等决定，属于稳定性问题。若某算法在计算过程中任一步产生的误差在以后的计算中都逐步衰减，则称该算法是绝对稳定的 讨论数值方法的稳定性，通常只考虑试验方程 当步长取为 h 时，将某算法应用于上式，并假设只在初值产生误差 ，则若此误差以后逐步衰减，就称该算法相对于 绝对稳定，的全体构成绝对稳定区域。有时，我们称算法A 比算法B 稳定，就是指 A 的绝对稳定区域比 B 的绝对稳定区域大。2.1 算例：考察显式欧拉法算法的稳定性解：由 得：于是：因此，要保证初始误差e0 以后逐步衰减