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

第六章常微分方程的数值解法本章内容§6.1引言§6.2欧拉方法§6.3龙格—库塔方法§6.4边值问题的数值方法一.问题提出有一个或多个导数及其函数的方程式称为微分方程，在工程中常遇到求解微分方程的问题。§6.1引言)()(],[),(dd00xyyyxybaxyxfxy解函数初值问题如，一阶常微分方程的§6.1引言二.两类定解问题◆常微分方程的定解问题有两种基本类型类:初值问题和边值问题◆定解指已知因变量和/或其导数在某些点上是已知的(约束条件)。●1.边值问题约束条件为已知，在自变量的任一非初值上，已知函数值和/或其导数值，如例如,受连续分布横向荷载的变截面简支梁弯曲问题ybyaybaxyyxfy求解)(,)(],[),,(0)()0()()(dd22lwwxEIxMxw—简支梁跨度弯曲刚度—挠度，—其中，lxEIxw,)()(q(x)xwOl●2.初值问题00)(),(yxyyxfy例如，单自由度系统的非线性受迫振动0000)(,)(),,(yxyyxyyyxfyxm()Pt023200dd()ddd(0)dxxxxmckxxPtttxxxxt，唯一。初值问题的解必存在且由常微分方程理论知：，使得即存在常数条件，满足李普希兹且关于适当光滑连续，只要函数yyLyxfyxfLLipschitzyyxf),(),()(),(实际问题中还存在初边值混合问题，如梁在横向激励下的弯曲振动。高阶常微分方程可以化成一阶的常微分方程组00)(),(yxyyxfy很多微分方程的解不能用初等函数来表示，有时即使能够用解析式表示其解，但计算量太大而不实用(表达式过于复杂)。需要用数值方法来求解，一般只要求得到若干个点上的近似值或者解的简单的近似表达式(精度要求满足即可)。计算）数值解（适合于计算机近似解析解，求近似解求精确解一般较为困难§6.1引言),,2,1,0(),,2,1,0()(,:1)(),(011000niihxxxxhniyxybxxxayxyyxfyiiiiin，有通常取等距节点，即。的近似值得到在数值解式§6.1引言yx()yfx0y0x0hxxxxx)()()()()()(1211211210nnnnnnnnxyyyyyyxyxyxyxyxxxxxxy的近似解上的值在一些离散点即解函数§6.1引言初值问题的常见解法单步法：利用前一个单步的信息(一个点)，在y=f(x)上找下一点yi，有欧拉法，龙格－库格法。预测校正法：多步法，利用一个以上的前点信息求f(x)上的下一个yi，常用迭代法，如改进欧拉法，阿当姆斯法。§6.1引言得到。，这可以通过递推公式，求出yy,,y,知的y进，由已点：按节点顺序依次推初值问题的数值解法特1ii10§6.2欧拉方法及其改进Euler’sMethod内容一.欧拉格式二.Euler预估—校正法三.误差估计、收敛性和稳定性6.2.1欧拉公式：/*Euler’sMethod*/向前差商近似导数hxyxyxy)()()(010),()()()(000001yxfhyxyhxyxy1y记为)1,...,0(),(1niyxfhyyiiii§6.2欧拉方法及其改进1111111()(,)(,)-(,),(,)-()nnnnnnnnnnnnnnnnnnnyyxPPxyPxyPPyyfxyyyhfxyxxEuler几何意义：折线逼近解曲线。设已做出折线的顶点，过依方向场的方向再推进，显然两个顶点，的坐标有关系即这就是欧拉公式。xP0P1P2P3P4PnyO()yyx，计算结果如下表：分别取步长公式为解：此时的。确解是的数值解，此问题的精，方法求初值问题例：利用05.0,1.0,2.02,1,0,0)211()1/()(0)0(202110221222hiyyxhyyEulerxxxyyxyxyEuleriiii§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越大越好。——阶方法则称该方法为。误差为若某种方法的局部截断：定义的误差，记为步严格算出逐，由—去掉假设—）整体截断误差（这一步的误差，记为严格算出，按某方法由—假设—）局部截断误差（：定义计数值解法的截断误差估pphORyxyyyxyyyxyRyyxyyPiiiiiiiiiiiiii)(2)()(2)()(11111111011116.2.2§6.2欧拉方法及其改进计算方法——09计11、61§7.2欧拉方法，一阶方法相减得与式式显格式为：，则假设积分法）展开法（也可使用数值解：使用显格式的局部截断误差例：求)()(2)()()()())(,()(),()()()(2))(,()(2)()()()(2121111221hOyhyxyRbabxyxhfxyyxhfyyEulerxyyayhxyxhfxyhyhxyxyxyTaylorEulerniiiiiiiiiiiiiiiiiiii§6.2欧拉方法及其改进计算方法——09计11、61§7.2欧拉方法（具有普遍意义）低一阶。比局部截断误差，则＝若初始误差为（证略）为求解范围。，（则足估计式方法的整体截断误差满有界，则且条件有关于设：定理系—整体”截断误差的关—“局部●iibxaababLiRbaxyMieLhMyxyeEulerxyyyLyxfyxfLipschitzyyxf0],[)(max),2,112)()(),(),(),(102)(200)(212211§6.2欧拉方法及其改进6.2.3隐式欧拉法/*implicitEulermethod*/向后差商来近似导数hxyxyxykkk)()()(11)1,...,0(),(111niyxfhyyiiii由于未知数yi+1同时出现在等式的两边，不能直接得到，故称为隐式/*implicit*/欧拉公式，而前者称为显式/*explicit*/欧拉公式。一般先用显式计算一个初值，再用隐式法（迭代）求解。)(,)(111kkkxyxfxy)211(2,1,0,0)211(0)0(202112121)0(1)0(1102121122iiiiiiiiiiyxhyyEuleryyiyyxhyyEuleryxyxy公式得到可由显式，迭代初值迭代求解公式为此时的隐式，例：求初值问题§6.2欧拉方法及其改进1221111221111221111()()()()()()()2()()()(,())()()2()(,)()()2()()iiiiiiiiiiiiiiiiiiiiiiiiyxxTayloryyxyxyxxxxxyyxyxfxyxhxxayyyfxyhxxbabEulerR将在作展开即又由式和得隐式法的截断误差为22111()()2iiiiyyxyh§6.2欧拉方法及其改进显格式确定。，初值用、用迭代法求情况）、化为显格式（仅小数是隐格式。实用上（（高一阶期望截断误差的阶数会若将两者做平均，可以正负号法的截断误差相差一个法与隐式显式EulerynyxfyxfhyyEulerEuleriiiiiii111121)1,0,1,2,(i)],),[2显、隐式两种算法的平均。需要迭代求解，能否不迭代？§6.2欧拉方法及其改进6.2.4梯形格式§6.2欧拉方法及其改进校正法—预估Euler6.2.5/*predictor-correctormethod*/Step1:先用显式欧拉公式作预测，算出),(1iiiiyxfhyyStep2:再将代入隐式梯形公式的右边作校正，得到1iy)],(),([2111iiiiiiyxfyxfhyy)1,...,0(),(,),(211niyxfhyxfyxfhyyiiiiiiii1,2,1,0,),(),()(20121211niyhKyhxfKyxfKKKhyyEuleriiiiii方法也可以写成：改进的§7.2欧拉方法及其改进§7.2欧拉方法9,2,1,0,11.0)1.0(21.0/2)(05.0)2(9,2,1,0,1/2.01.1)1()21()(1.01)0(10201121211012/1iyKyxKyKyxyKKKyyEuleriyyxyyEulerxxyhyxyxyyiiiiiiiiiiii方法利用改进方法利用解：）。（精确解为的数值解，取步长，例：求初值问题§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龙格—库塔方法内容一.2阶龙格—库塔格式三.高阶龙格—库塔格式单步法：即利用前一个节点的函数值yi，计算后一个节点的函数值yi+1。目的：建立高精度的单步递推格式。单步递推法的基本思想是从(xi,yi)点出发，以某一斜率沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形所能达到的最高精度为2阶。§6.3龙格—库塔方法二.2阶龙格—库塔格式斜率一定取K1K2的平均值吗？步长一定是一个h吗？)(),(),(2121)],(),([2),(00121211111111xyyhKyhxfKyxfKKKhyyyxfyxfhyyyxhfyyiiiiiikkkkkkkk