第六章常微分方程的数值解法演示

（优选）第六章常微分方程的数值解法当前1页，总共60页。本章内容§6.1引言§6.2欧拉方法§6.3龙格—库塔方法§6.4边值问题的数值方法当前2页，总共60页。一.问题提出有一个或多个导数及其函数的方程式称为微分方程，在工程中常遇到求解微分方程的问题。§6.1引言当前3页，总共60页。§6.1引言二.两类定解问题◆常微分方程的定解问题有两种基本类型类:初值问题和边值问题◆定解指已知因变量和/或其导数在某些点上是已知的(约束条件)。●1.边值问题约束条件为已知，在自变量的任一非初值上，已知函数值和/或其导数值，如例如,受连续分布横向荷载的变截面简支梁弯曲问题当前4页，总共60页。q(x)xwOl●2.初值问题例如，单自由度系统的非线性受迫振动当前5页，总共60页。实际问题中还存在初边值混合问题，如梁在横向激励下的弯曲振动。高阶常微分方程可以化成一阶的常微分方程组当前6页，总共60页。很多微分方程的解不能用初等函数来表示，有时即使能够用解析式表示其解，但计算量太大而不实用(表达式过于复杂)。需要用数值方法来求解，一般只要求得到若干个点上的近似值或者解的简单的近似表达式(精度要求满足即可)。§6.1引言当前7页，总共60页。§6.1引言当前8页，总共60页。§6.1引言当前9页，总共60页。初值问题的常见解法单步法：利用前一个单步的信息(一个点)，在y=f(x)上找下一点yi，有欧拉法，龙格－库格法。预测校正法：多步法，利用一个以上的前点信息求f(x)上的下一个yi，常用迭代法，如改进欧拉法，阿当姆斯法。§6.1引言当前10页，总共60页。§6.2欧拉方法及其改进

Euler’sMethod

内容一.欧拉格式二.Euler预估—校正法三.误差估计、收敛性和稳定性当前11页，总共60页。6.2.1欧拉公式：/*Euler’sMethod*/向前差商近似导数记为§6.2欧拉方法及其改进当前12页，总共60页。xP0P1P2P3P4PnyO当前13页，总共60页。§6.2欧拉方法及其改进当前14页，总共60页。§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当前15页，总共60页。§6.2欧拉方法及其改进当前16页，总共60页。§7.2欧拉方法§6.2欧拉方法及其改进当前17页，总共60页。§7.2欧拉方法当前18页，总共60页。§6.2欧拉方法及其改进6.2.3隐式欧拉法/*implicitEulermethod*/向后差商来近似导数)1,...,0(),(111-=+=+++niyxfhyyiiii由于未知数yi+1

同时出现在等式的两边，不能直接得到，故称为隐式

/*implicit*/

欧拉公式，而前者称为显式

/*explicit*/欧拉公式。一般先用显式计算一个初值，再用隐式法（迭代）求解。当前19页，总共60页。§6.2欧拉方法及其改进当前20页，总共60页。§6.2欧拉方法及其改进当前21页，总共60页。显、隐式两种算法的平均。需要迭代求解，能否不迭代？§6.2欧拉方法及其改进6.2.4梯形格式当前22页，总共60页。

§6.2欧拉方法及其改进/*predictor-correctormethod*/Step1:先用显式欧拉公式作预测，算出),(1iiiiyxfhyy+=+Step2:再将代入隐式梯形公式的右边作校正，得到1+iy)],(),([2111+++++=iiiiiiyxfyxfhyy当前23页，总共60页。§7.2欧拉方法及其改进当前24页，总共60页。§7.2欧拉方法当前25页，总共60页。§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计算结果如下:当前26页，总共60页。§6.3龙格—库塔方法内容一.2阶龙格—库塔格式三.高阶龙格—库塔格式当前27页，总共60页。单步法：即利用前一个节点的函数值yi，计算后一个节点的函数值yi+1。目的：建立高精度的单步递推格式。单步递推法的基本思想是从(xi,yi)点出发，以某一斜率沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形所能达到的最高精度为2阶。§6.3龙格—库塔方法二.2阶龙格—库塔格式当前28页，总共60页。斜率一定取K1K2的平均值吗？步长一定是一个h吗？§6.3龙格—库塔方法当前29页，总共60页。§7.3龙格—库塔方法首先希望能确定系数1、2、p，使得到的算法格式有2阶精度，即在的前提假设下，使得

Step1:

将K2在(xi,yi)点作Taylor展开将改进欧拉法推广为：),(),(][12122111phKyphxfKyxfKKKhyyiiiiii++==++=+ll当前30页，总共60页。Step3:

将y(xi+1)在xi点的泰勒展开并与yi+1作比较要求，则必须有：这里有3个未知数，2个方程。存在无穷多个解。所有满足上式的格式统称为2阶龙格-库塔格式。注意到，就是改进的欧拉法。Step2:

将K2代入第1式，得到当前31页，总共60页。其中i

(i=1,…,m)，i

(i=2,…,m)

和ij

(i=2,…,m;j=1,…,i1

)

均为待定系数，确定这些系数的步骤与前面相似。)...,(......),(),(),(]...[1122112321313312122122111--++++++=+++=++==++++=mmmmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyybbbabbaballl问题:为获得更高的精度，应该如何进一步推广？§6.3龙格—库塔方法三.高阶龙格—库塔格式当前32页，总共60页。3阶龙格-库塔法§6.3龙格—库塔方法当前33页，总共60页。§6.3龙格—库塔方法最常用为四级4阶经典龙格-库塔法当前34页，总共60页。注：

龙格-库塔法的主要运算在于计算Ki

的值，即计算f

的值。Butcher于1965年给出了计算量与可达到的最高精度阶数的关系：753可达到的最高精度642每步须算Ki的个数由于龙格-库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算法而将步长h取小。§6.3龙格—库塔方法当前35页，总共60页。§7.3龙格—库塔方法例：使用高阶R-K方法计算初值问题解:（1）使用三阶R-K方法(公式)§6.3龙格—库塔方法当前36页，总共60页。§7.3龙格—库塔方法其余结果如下:

ixiK1K2K3

yi1.00000.10001.00001.10251.25551.11112.00000.20001.23451.37551.59451.24993.00000.30001.56241.76372.09221.42844.00000.40002.04042.34232.86581.66645.00000.50002.77683.25874.16341.9993§6.3龙格—库塔方法当前37页，总共60页。§7.3龙格—库塔方法（2）如果使用四阶R-K方法(公式)§6.3龙格—库塔方法当前38页，总共60页。§7.3龙格—库塔方法其余结果如下:

ixiK1K2K3K4yi1.00000.10001.00001.10251.11331.23511.11112.00000.20001.23461.37561.39211.56331.25003.00000.30001.56251.76391.79082.04231.42864.00000.40002.04082.34282.38922.78051.66675.00000.50002.77773.26003.34764.00572.0000§6.3龙格—库塔方法当前39页，总共60页。步长过大，达不到精度要求；步长过小，虽然局部截断误差小，加大了计算工作量，舍入误差的累积增大。解决途径之一——引入变步长技术，常用的有Richardson外推法。从结点xi出发，先以h为步长，通过一步计算出y(xi+1)的近似值§6.4步长的自动选择当前40页，总共60页。当前41页，总共60页。当前42页，总共60页。当前43页，总共60页。§7.2欧拉方法/*Convergency*/§6.5收敛性与稳定性当前44页，总共60页。例：就初值问题考察欧拉显式格式的收敛性。解：该问题的精确解为欧拉公式为对任意固定的x

=xi=ih，有§6.5收敛性与稳定性当前45页，总共60页。§7.2欧拉方法例：考察初值问题在区间[0,0.5]上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。0.00.10.20.30.40.5精确解改进欧拉法

欧拉隐式欧拉显式

节点xi1.00002.00004.00008.00001.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.1§6.5收敛性与稳定性当前46页，总共60页。§7.2欧拉方法§6.5收敛性与稳定性当前47页，总共60页。一般分析时为简单起见，只考虑试验方程常数，可以是复数§6.5收敛性与稳定性当前48页，总共60页。§6.5收敛性与稳定性当前49页，总共60页。Euler法的绝对稳定区域当前50页，总共60页。§7.3龙格—库塔方法例：隐式龙格-库塔法而显式1~4阶方法的绝对稳定区域为其中2阶方法的绝对稳定区域为0ReImk=1k=2k=3k=4-1-2-3---123ReIm无条件稳定当前51页，总共60页。§6.6一阶常