第7章-常微分方程数值解法课件

1第七章第七章常微分方程的数值解法常微分方程的数值解法1引引言言2欧拉方法欧拉方法3龙格龙格-库塔方法库塔方法2 1引引言言 在工程和科学技术的实际问题中，常需要解常微分方程。但常微分方程组中往往只有少数较简单和典型的常微分方程（例如线性常系数常微分方程等）可求出其解析解。对于变系数常微分方程的解析求解就比较困难，而一般的非线性常微分方程就更不用说了。在大多数情况下，常微分方程只能用近似法求解。这种近似解法可分为两大类：近似解法可分为两大类：一类是近似解析法近似解析法，如级数解法、逐次逼近法等；另一类则是数值解法数值解法，它给出方程在一些离散点上的近似解。在具体求解微分方程时，需要具备某种定解条件定解条件，微分方程和定解条件合在一起组成定解问题定解问题。定解条3件有两种：一种是给出积分曲线在初始点的状态，称为初始条件初始条件，相应的定解问题称为初值问题初值问题；另一种是给出积分曲线首尾两端的状态，称为边界条件边界条件，相应的定解问题则称为边值问题边值问题。例如，例如，弹簧-质量系统的振动问题（图图7-1），作一定的 简化后，可用一个二阶常微分方程来描述。式中，x是质量m离平衡位置（0点）的距离；t 是时间；c是弹簧常数。4当弹簧在振动开始时刻t=t0 时的初始位置x(t0)=x0和初速度 确定时，弹簧的振动规律x(t)也就唯一确定。这就是一个常微分方程的初值问题，可写成：5 m x x oc图图7-1 本章先从一阶常微分方程的初值问题：（1.1）出发进行讨论。6 由常微分方程的理论知，只要上式中的函数f(x,y)在区域 内连续，且关于y满足李普希李普希兹兹(Lipschitz)条件条件，即存在常数L(它与x,y无关)使对 内任意两个y1和 y2 都成立，则方程的解必定存在且唯一。下面的分析均假定满足上述条件。初值问题（1.1）的数值解法，常采用差分方法，即把一个连续的初值问题离散化为一个差分方程来求解。即将（1.1）离散化后，求找其解y=y(x)在一系列离散节点7上的近似值y0,y1,yn。两相邻节点间的距离称为步长步长。当 （常值）时称为等步长等步长，有 或 a=x0 x1 xi xn=bhi=xi+1 -xi (i=0,1,2,n-1)8因为初值问题中的初始条件 为已知，即可利用已知的 来求出下一节点处 的近似值 再从 来求出 如此继续，直到求出 为止。这种用按节点的排列顺序一步一步地向前推进的方式求解的差分算法称为“步进式步进式”或“递进式递进式”算法。它是初值问题数值解法的各种差分格式的共同特点。因此，只要能写出由前几步已知信息 来计算的递推公式（即差分格式差分格式），即可完全表达该种算法。92欧拉方法欧拉方法2.1欧拉格式欧拉格式 对于初值问题（1.1）式，先将其离散化，即把a,b 作 n等分，得各离散节点 xi=a +ih (i=0,1,2,n-1)式中 h=(b-a)/n 设y=y(x)为方程（1.1）的解，则 y=y(xi+1)在(xi,yi)点处的泰勒展开式为：10 当 有界且 h 充分小时，可忽略高阶无穷小量 将上式写成(2.1)或11 若将 和 的近似值分别记为 和 ，则得这就是欧拉（欧拉（Euler）公式）公式，又称欧拉格式欧拉格式。利用它可由已知的初值 出发，逐步算出 。这类形式的方程也称为差分方程差分方程。当假定 为准确，即在 的前提下来估计误差 ，这种截断误差称为局部截断误差局部截断误差。由（2.1）和（2.2）可知，欧拉格式在节点 处的局部截断误差显然为：(2.2)12(2.3)如果局部截断误差为 ，则称这种数值算法的精度为P阶。故欧拉格式的精度为一阶。从几何意义上来看欧拉格式，可如图图7-2中所示。由方程（1.1）知，其积分曲线 y=y(x)上任意一点(x,y)的切线斜率dy/dx 都等于函数 f(x,y)的值。从初值点P0（即点(x0,y0)）出发，作积分曲线y=y(x)在P0 点上的切线 （其斜率为f(x0,y0)），与直线 x=x1 相交于点 P1（即点(x1,y1)），得到y1作为y(x1)的近似值，则有13 相比较可知，这时是用切线段 近似代替了曲线段 ；P1 点近似代替了 点；y1 近似代替了y(x1)；hf(x0,y0)近似代替了递推继续从P1点出发，作一斜率为f(x1,y1)的直线 ，与直线x=x2 相交于P2 点（即点(x2,y2)）14与yox图图7-215得到 y2 作为y(x2)的近值；如此继续，直到Pn 点。这样，得出一条折线P0P1P2PiPn 近似代替积分曲线P0P1P2PiPn。当步数越多时，由于误差的积累，折线P0P1P2PiPn可能会越来越偏离真解 P0P1P2PiPn 曲线。差分是微分的近似计算，所以欧拉格式也可用差商近似代替导数的离散方法得到。在节点xi处有（2.4）用向前差商16近似代替上式中的导数项y(xi),即（2.5）代入（2.4）,可得：y(xi+1)和 y(xi)用其近似值 yi+1 和 yi 代入，则得17此即（2.2）（欧拉格式）（欧拉格式）。显然，欧拉格式具有递推性，在计算yi+1时只要用到前一步所得结果yi一个信息就够了，因此是一种单单步步格格式式或称一一步步格格式式。其中（2.5）是一个数值微分公式。故用其他数值微分公式也可导出略异于（2.2）的其他形式的算式来。例如，用向后差商表示的两点数值微分公式代入18（2.6）可得y(xi+1)和 y(xi)用其近似值 yi+1 和 yi 代入，则得19其局部截断误差为：(2.7)(2.7)(2.6)(2.6)称为向后欧拉格式向后欧拉格式，又称为隐式欧拉格式隐式欧拉格式。因为在此式的右边也包含未知的yi+1，所以是yi+1的一个函数方程，故称它为隐式格式隐式格式。（2.22.2）的右边则没有未知的yi+1，因此是一种显式格式显式格式。隐式格式的计算当然比显式格式要困难得多，一般情况下，只能用迭代法求解，计算量比较大。20不论是显式欧拉格式（2.22.2），还是隐式欧拉格式（2.62.6），它们都是单步格式单步格式或称为一步格式一步格式。因为它们在计算yi+1时只用到前一步所得结果yi一个信息；而格式（2.82.8）则除了yi外，还需用到更前一步所得信息yi-1，即需调用前两步的信息，因此（2.82.8）称为两步欧拉格式两步欧拉格式，或称为中点欧拉格式中点欧拉格式。比较（2.32.3）,（2.72.7）和（2.92.9）可知，两步欧拉格式比显式或隐式欧拉格式具有更高的精度，因为它的局部截断误差是O(h3)。由（2.32.3）和(2.7)(2.7)可见，显式欧拉格式与隐式欧拉格式的局部截断误差的符号正好相反，因此可以设想取（2.22.2）和(2.6)(2.6)的平均，即两式相加除以2，得21 （2.10）使误差相互抵消，提高误差阶数，从而提高算法的精度。事实上，格式(2.10)(2.10)的局部截断误差为O(h3)，即其精度为二阶。(2.10)(2.10)称为梯形格式梯形格式，是一种隐式格式。从几何角度上来看，梯形格式是取xi,xi+1的两端点的平均斜率。从图图7-37-3可见，用梯形格式得到的P点比用显式欧拉格式得到的A点和用隐式欧拉格式得到的B点都要合理。22 APBxyO图图7-37-3232.2改进的欧拉格式改进的欧拉格式 显式欧拉格式计算工作量比较小，但精度低。梯形格式虽提高了精度，但为隐式格式，需用迭代法求解，计算工作量大。综合这两种格式可得到改进的欧拉格式。先用欧拉格式（2.2）求出一个初步的近似值 ，称为预测值，它的精度不高，接着用梯形格式（2.10）对它校正一次，即迭代一次，求得yi+1，称为校正值。这种预测-校正方法称为改进的欧拉格式改进的欧拉格式：预测（2.11）校正24可以证明，（2.11）格式的精度为二阶。它是一种一步显式格式。（2.11）也可写成 （2.12）或将它写成下列平均化的形式：（2.13）25例例1 试分别用欧拉格式和改进的欧拉格式求解下 列初值问题,并比较两法所得结果的精度:解解取步长h=0.1.这样,欧拉格式和改进的欧拉格 式的具体算式分别为及26 两式格式的计算结果分别列于表表7-17-1中。由此表可见，与精确解 相比，改进的欧拉格式的精度较欧拉格式有明显的提高。27表表7-1283龙格龙格库塔方法库塔方法 要进一步提高求解的精度，可用一种高精度的单步法龙格龙格库塔（库塔（RungeKutta）方法）方法，简称RK法法。它采用了间接使用泰勒级数法的技术。3.1龙格龙格库塔公式的导出库塔公式的导出 对于一阶常微分方程（1.1）的解 y=y(x)，利用微分中值定理得即29 也即(3.1)(3.2)式中 K 可看作是y=y(x)在区间xi,xi+1 上的平均平均斜率斜率。这样，欧拉格式(2.2)相当于取(xi,yi)点上斜率 作为平均斜率K 的近似值，当然30与（3.1）比较知，它相当于把(xi,yi)和(xi+1,yi+1)两个点上的斜率K1和K2 的算术平均值作为（3.1）中的平均斜率K的近似值。其中K2 是通过已知信息yi来近似地预测的。是十分粗糙的，因此精度必然很低。而改进的欧拉格式（2.12）可改写成31 由此可以设想，取区间xi,xi+1 内的某一个节点上的斜率K2 与(xi,yi)点上的斜率K1作线性组合（即加权平均），作为平均斜率K的近似值，即为了要得到(xi+p,yi+p)点上的斜率 K2，需先预测yi+p=yi+ph K132 由此构成的计算格式为(3.3)根据预测值 yi+p 再来算出K2：上式含有三个待定参数:1,2 和 p。适当选定其值可使算法的局部截断误差为O(h3)，即有二阶精度。33 将它与y(x)在xi+1 处的二阶泰勒展开式假定yi=y(xi)，分别将K1 和K2 作泰勒展开：代入（3.3），得34进行比较系数后可知，只要（3.4）成立，格式（3.3）的局部截断误差就等于O(h3)，从而能具有二阶精度。35（3.4）中有三个待定系数:1,2 和 p，但却只有两个方程式，因此还有一个自由度一个自由度。凡满足条件（3.4）的一族格式（3.3）统称为二阶龙格二阶龙格-库塔格式库塔格式。当 p=1,1=2=1/2时，二阶RK格式（3.3）即为改进的欧拉公式（2.12）。如取p=1/2，则 1=0,2=1，（3.3）就成为（3.5）36（3.5）称为变形的欧拉格式变形的欧拉格式。由于（3.5）中的 是欧拉格式预测出来的区间xi,xi+1的中点xi+1/2的近似解,K2 就近似地等于此中点的斜率 ，因此（3.5）就相当于用中点xi+1/2的斜率作为（3.1）中的平均斜率K的近似值，故格式（3.5）也称为中中点点格格式式。粗看起来，yi+1=yi+hK2中只含有一个斜率值K2，但实际上K2 是通过K1 才能算出来的，因此，式37中还隐含着K1。这样，每完成一步仍需计算函数 f 值两次，其计算工作量仍与改进的欧拉格式一样。3.2 高阶龙格高阶龙格-库塔公式库塔公式 要提高精度，可在区间xi,xi+1内取二个节点：上的斜率K2，K3与点(xi,yi)上的斜率K1 加权平均，作为平均斜率K的近似值，即 38其中K1和K2仍如（3.3）。利用区间xi,xi+q 内的两个斜率K1和K2，加权平均作为其平均斜率K，来预测 y(xi+q)：从而得到由此构成的计算格式为39(3.6)类似于二阶龙格-库塔格式的导出过程，运用泰勒展开的方法，可找出格式（3.6）的局部截断误差为 ，从而具有三阶精度所必须满足的条件为：40(3.7)其中共有七个待定系数：，但只有五个方程式，因此还有两个自由度两个自由度。凡满足条件（3.7）的一族格式（3.6）统称为三阶龙格三阶龙格-库库塔格式塔格式。41当待定系数取为(3.8)时的三阶龙格-库塔格式称为库塔格式库塔格式，其具体形式为：42 继续推广这种处理过程，可得四阶龙格-库塔格式。最常用的一种经典龙格-库塔格式的具体形式如下：(3.9)43例例2 试分别用欧拉方法 ，改进的欧拉 方法 及经典RK方法 ，求解下列初值问题，并比较三种方法所得结 果的精度：解解 三种方法的具体算式分别如下：欧拉格式：44改进的欧拉公式：经典RK格式：45精确解表表7-246三种格式的 计算结果分别列于表表72中。上例中对三种方法采用了不同的步长h值，是为了使它们所耗的计算工作量大致相同，以便于比较。从表表72可见，经典RK方法的精度较改进的欧拉方法又有了很大的提高。关于这一结论也可以从理论上大致分析出来：（1）欧拉方法的局部截断误差为计算四步后截断误差为47 （2）经典RK方法的局部截断误差为 为大致相同数量级下的常数，故有 注意：注意：RK方法的导出利用了泰勒展开，因此要求所求的解有较好的光滑性，如果解的光滑性差，则采用经典 RK 法所得的数值解，其精度有可能反而不及改进的欧拉法。因此在实际计算中，应根据问题的具体情况来选择合适的算法。