第六章—常微分方程的数值解法.ppt

第六章第六章常微分方程的数值解法常微分方程的数值解法.本章内容本章内容n6.1 引言引言n6.2 欧拉方法欧拉方法n6.3 龙格龙格库塔方法库塔方法n6.4 边值问题的数值方法边值问题的数值方法.n一一.问题提出问题提出 有一个或多个导数及其函数的方程式称为微分方有一个或多个导数及其函数的方程式称为微分方程，在工程中常遇到求解微分方程的问题。程，在工程中常遇到求解微分方程的问题。6.1引言引言.6.1引言引言n二二.两类定解问题两类定解问题常微分方程的定解问题有两种基本类型类常微分方程的定解问题有两种基本类型类:初值问题初值问题和和边值问题边值问题定解指已知因变量和定解指已知因变量和/或其导数在某些点上是已知的或其导数在某些点上是已知的(约束条件约束条件)。1.边值问题边值问题 约束条件为已知，在自变量的任一非初值上，已知约束条件为已知，在自变量的任一非初值上，已知函数值和函数值和/或其导数值，如或其导数值，如 例如例如,受连续分布横向荷载的变截面简支梁弯曲问题受连续分布横向荷载的变截面简支梁弯曲问题.q(x)xwOl2.初值问题初值问题例如，单自由度系统的非线性受迫振动例如，单自由度系统的非线性受迫振动.实际问题中还存在初边值混合问题，如梁在横向激励下的弯曲振动。高阶常微分方程可以化成一阶的常微分方程组.很多微分方程的解不能用初等函数来表示，有时很多微分方程的解不能用初等函数来表示，有时即使能够用解析式表示其解，但计算量太大而不实即使能够用解析式表示其解，但计算量太大而不实用用(表达式过于复杂表达式过于复杂)。需要用数值方法来求解，一般只要求得到若干个需要用数值方法来求解，一般只要求得到若干个点上的近似值或者解的简单的近似表达式点上的近似值或者解的简单的近似表达式(精度要求精度要求满足即可满足即可)。6.1引言引言.6.1引言引言.6.1引言引言.初值初值问题问题的的常见常见解法解法单步法：单步法：利用前一个单步的信息利用前一个单步的信息(一个点一个点)，在，在y=f(x)上找下一点上找下一点yi，有欧拉法，龙格库格法。有欧拉法，龙格库格法。预测校正法：预测校正法：多步法，利用一个以上的前点信息求多步法，利用一个以上的前点信息求f(x)上上的下一个的下一个yi，常用迭代法，如改进欧拉法，阿当姆斯法。常用迭代法，如改进欧拉法，阿当姆斯法。6.1引言引言.6.2欧拉方法及其改进欧拉方法及其改进EulersMethodn内容内容一一.欧拉格式欧拉格式二二.Euler预估预估校正法校正法三三.误差估计、收敛性和稳定性误差估计、收敛性和稳定性.6.2.1 欧拉公式：欧拉公式：/*EulersMethod*/向前差商近似导数向前差商近似导数记为记为6.2欧拉方法及其改进欧拉方法及其改进.xP0P1P2P3P4PnyO.6.2欧拉方法及其改进欧拉方法及其改进.7.2欧拉方法欧拉方法hxiyi真值真值y(xi)误差误差y(xi)-yih=0.20.000.400.801.201.602.000.000000.376310.542280.527090.466320.406820.000000.344830.487800.491800.449440.400000.00000-0.03148-0.05448-0.03529-0.01689-0.00682h=0.10.000.400.801.201.602.000.000000.360850.513710.509610.458720.404190.000000.344830.487800.491800.449440.400000.00000-0.01603-0.02590-0.01781-0.00928-0.00419h=0.050.000.400.801.201.602.000.000000.352870.500490.500730.454250.402270.000000.344830.487800.491800.449440.400000.00000-0.00804-0.01268-0.00892-0.00481-0.00227.6.2欧拉方法及其改进欧拉方法及其改进.计算方法09计11、617.2欧拉方法欧拉方法6.2欧拉方法及其改进欧拉方法及其改进.计算方法09计11、617.2欧拉方法欧拉方法.6.2欧拉方法及其改进欧拉方法及其改进6.2.3隐式欧拉法隐式欧拉法/*implicitEulermethod*/向后差商来向后差商来近似导数近似导数)1,.,0(),(111=+=+niyxfhyyiiii由于未知数由于未知数 yi+1 同时出现在等式的两边，不能直接得同时出现在等式的两边，不能直接得到，故称为到，故称为隐式隐式/*implicit*/欧拉公式，而前者称为欧拉公式，而前者称为显显式式/*explicit*/欧拉公式。欧拉公式。一般先用显式计算一个初值，再用隐式法（一般先用显式计算一个初值，再用隐式法（迭代）迭代）求解。求解。.6.2欧拉方法及其改进欧拉方法及其改进.6.2欧拉方法及其改进欧拉方法及其改进.显、隐式两种算法的显、隐式两种算法的平均。需要迭代求解，能否不迭代？平均。需要迭代求解，能否不迭代？6.2欧拉方法及其改进欧拉方法及其改进6.2.4 梯形格式梯形格式.6.2欧拉方法及其改进欧拉方法及其改进/*predictor-correctormethod*/Step 1:先用先用显式显式欧拉公式作欧拉公式作预测预测，算出，算出),(1iiiiyxfhyy+=+Step 2:再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1+iy),(),(2111+=iiiiiiyxfyxfhyy.7.2欧拉方法及其改进欧拉方法及其改进.7.2欧拉方法欧拉方法.7.2欧拉方法欧拉方法ixiEuler方法方法yi改进改进Euler法法yi精确解精确解y(xi)01234567891000.10.20.30.40.50.60.70.80.9111.11.1918181.2774381.3582131.4351331.5089661.5803381.6497831.7177791.78477011.0959091.1840961.2662011.3433601.4164021.4859561.5525151.6164761.6781681.73786911.0954451.1832161.2649911.3416411.4142141.4832401.5491931.6124521.6733201.732051计算结果如下计算结果如下:.6.3龙格龙格库塔方法库塔方法n内容内容一一.2阶龙格阶龙格 库塔格式库塔格式三三.高阶龙格高阶龙格 库塔格式库塔格式.单步法：单步法：即利用前一个节点的函数值即利用前一个节点的函数值yi，计算后计算后一个节点的函数值一个节点的函数值yi+1。目的：目的：建立高精度的单步递推格式。建立高精度的单步递推格式。单步递推法的单步递推法的基本思想基本思想是从是从(xi,yi)点出发，以点出发，以某一斜率某一斜率沿直线达到沿直线达到(xi+1,yi+1)点。欧拉法及其点。欧拉法及其各种变形所能达到的最高精度为各种变形所能达到的最高精度为2 阶。阶。6.3龙格龙格库塔方法库塔方法n二二.2阶龙格阶龙格 库塔格式库塔格式.斜率斜率一定取一定取K1 K2的的平均值平均值吗？吗？步长一定是一个步长一定是一个h 吗？吗？6.3龙格龙格库塔方法库塔方法.7.3龙格龙格库塔方法库塔方法首先希望能确定系数首先希望能确定系数 1、2、p，使得到的算法格式有，使得到的算法格式有2阶阶精度，即在精度，即在 的前提假设下，使得的前提假设下，使得 Step 1:将将 K2 在在(xi,yi)点作点作 Taylor 展开展开将改进欧拉法推广为：将改进欧拉法推广为：),(),(12122111phKyphxfKyxfKKKhyyiiiiii+=+=+.Step 3:将将 y(xi+1)在在 xi 点的点的泰勒泰勒展开并与展开并与yi+1作比较作比较要求要求 ，则必须有：，则必须有：这里有这里有3个未知数，个未知数，2个方程。个方程。存在存在无穷多个解无穷多个解。所有满足上式的格式统称为。所有满足上式的格式统称为2阶龙阶龙格格-库塔格式库塔格式。注意到，注意到，就是就是改进的欧拉法改进的欧拉法。Step 2:将将 K2 代入第代入第1式，得到式，得到.其中其中 i (i=1,m)，i (i=2,m)和和 ij(i=2,m;j=1,i 1)均为待定系数，确定这些系数的步骤均为待定系数，确定这些系数的步骤与前面相似。与前面相似。).,(.),(),(),(.1122112321313312122122111 +=+=+=+=mm mmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy 问题问题:为获得更高的精度，应该如何进一步推广？为获得更高的精度，应该如何进一步推广？6.3龙格龙格库塔方法库塔方法n三三.高阶龙格高阶龙格 库塔格式库塔格式.3阶阶龙格龙格-库塔法库塔法6.3龙格龙格库塔方法库塔方法.6.3龙格龙格库塔方法库塔方法最常用为四级最常用为四级4阶阶经典龙格经典龙格-库塔法库塔法.注：注：龙格龙格-库塔法的主要运算在于计算库塔法的主要运算在于计算 Ki 的值，即计算的值，即计算 f 的的值。值。Butcher 于于1965年给出了计算量与可达到的最高精年给出了计算量与可达到的最高精度阶数的关系：度阶数的关系：753可达到的最高精度可达到的最高精度642每步须算每步须算Ki 的个数的个数 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精度主要受库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好解函数的光滑性影响。对于光滑性不太好的解，最好采用采用低阶算法低阶算法而将步长而将步长h 取小取小。6.3龙格龙格库塔方法库塔方法.计算方法09计11、617.3龙格龙格库塔方法库塔方法例：使用高阶例：使用高阶R-K方法计算初值问题方法计算初值问题解解:（1）使用三阶）使用三阶R-K方法方法(公式公式)6.3龙格龙格库塔方法库塔方法.计算方法09计11、617.3龙格龙格库塔方法库塔方法其余结果如下其余结果如下:i xi K1 K2 K3 yi 1.0000 0.1000 1.0000 1.1025 1.2555 1.1111 2.0000 0.2000 1.2345 1.3755 1.5945 1.2499 3.0000 0.3000 1.5624 1.7637 2.0922 1.4284 4.0000 0.4000 2.0404 2.3423 2.8658 1.6664 5.0000 0.5000 2.7768 3.2587 4.1634 1.9993 6.3龙格龙格库塔方法库塔方法.计算方法09计11、617.3龙格龙格库塔方法库塔方法（2）如果使用四阶）如果使用四阶R-K方法方法(公式公式)6.3龙格龙格库塔方法库塔方法.计算方法09计11、617.3龙格龙格库塔方法库塔方法其余结果如下其余结果如下:i xi K1 K2 K3 K4 yi 1.0000 0.1000 1.0000 1.1025 1.1133 1.2351 1.1111 2.0000 0.2000 1.2346 1.3756 1.3921 1.5633 1.2500 3.0000 0.3000 1.5625 1.7639 1.7908 2.0423 1.4286 4.0000 0.4000 2.0408 2.3428 2.3892 2.7805 1.6667 5.0000 0.5000 2.7777 3.2600 3.3476 4.0057 2.00006.3龙格龙格库塔方法库塔方法.步长过大，达不到精度要求；步长过小，虽然局部截断误差小，加大了计算工作量，舍入误差的累积增大。解决途径之一引入变步长技术，常用的有Richardson外推法。从结点 xi出发，先以h为步长，通过一步计算出y(xi+1)的近似值6.4步长的自动选择步长的自动选择.计算方法09计11、617.2欧拉方法欧拉方法/*Convergency*/6.5收敛性与稳定性收敛性与稳定性.例：例：就初值问题就初值问题 考察欧拉显式格式的收敛性。考察欧拉显式格式的收敛性。解：解：该问题的精确解为该问题的精确解为 欧拉公式为欧拉公式为对任意固定的对任意固定的 x=xi=i h，有，有 6.5收敛性与稳定性收敛性与稳定性.计算方法09计11、617.2欧拉方法欧拉方法例：例：考察初值问题考察初值问题 在区间在区间0,0.5上的解。上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。0.00.10.20.30.40.5精确解改进欧拉法 欧拉隐式欧拉显式 节点 xi 1.00002.0000 4.00008.0000 1.6000101 3.2000101 1.00002.5000101 6.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.90631019.76561011.00004.97871022.47881031.23411046.14421063.0590107原因原因?!取步长取步长 h=0.16.5收敛性与稳定性收敛性与稳定性.计算方法09计11、617.2欧拉方法欧拉方法6.5收敛性与稳定性收敛性与稳定性.一般分析时为简单起见，只考虑一般分析时为简单起见，只考虑试验方程试验方程常数，可常数，可以是复数以是复数6.5收敛性与稳定性收敛性与稳定性.6.5收敛性与稳定性收敛性与稳定性.Euler法的绝对稳定区域.7.3龙格龙格库塔方法库塔方法例：例：隐式龙格隐式龙格-库塔法库塔法而而显式显式 1 4 阶方法的绝对稳定阶方法的绝对稳定区域为区域为其中其中2阶方法阶方法 的绝对稳定区域为的绝对稳定区域为0ReImk=1k=2k=3k=4-1-2-3-123ReIm无条件稳定无条件稳定.6.6 一阶常微分方程组与高阶方程一阶常微分方程组与高阶方程 可以把单个方程 中的 f 和 y 看作向量来处理，这样就可把前面介绍的各种差分算法推广到求一阶方程组初值问题中来。6.6.1 一阶常微分方程组对于一阶常微分方程组的初值问题对于一阶常微分方程组，可类似得到各种解法，而高阶常微分方对于一阶常微分方程组，可类似得到各种解法，而高阶常微分方程可转化为一阶常微分方程组来求解。程可转化为一阶常微分方程组来求解。.设 为节点上的近似解，则有改进的Euler格式为 预报：校正：校正：又，相应的四阶龙格库塔格式（经典格式）为.式中式中.例 用改进的Euler法求解初值问题 取步长h=0.1，保留六位小数。解:改进的Euler法公式为.由初值由初值 ,计算得计算得.6.6.1 一阶常微分方程组 高阶微分方程(或方程组)的初值问题,原则上都可以归结为一阶方程组来求解。例如,有二阶微分方程的初值问题 在引入新的变量 后,即化为一阶方程组初值问题:此为一个一阶方程组的初值问题，对此可用前面中介绍的方法来求解。例如应用四阶龙格-库塔公式得.例例 求解下列二阶微分方程的初值问题求解下列二阶微分方程的初值问题 取步长取步长h=0.1 h=0.1 解解:先作变换：令先作变换：令 ，代入上式，得一阶方程组，代入上式，得一阶方程组 用用四四阶阶龙龙格格-库库塔塔方方法法求求解解,按按式式(9.37)(9.37)及及(9.38)(9.38)进行计算：进行计算：取步长取步长 ,时时 .然后计算然后计算 时的时的y y2 2和和z z2 2；依此类推，直到；依此类推，直到i=9i=9时的时的y y1010和和z z1010，即可得到，即可得到其数值解。其数值解。.6.7 边值问题的数值解法边值问题的数值解法 打靶法与有限差分法.