第六章36微分方程数值解法

1、 第6章 常微分方程数值解法 微分方程数值解一般可分为：微分方程数值解一般可分为：常微分方程数值解常微分方程数值解和和偏微分偏微分方程数值解方程数值解。自然界与工程技术中的许多现象，其数学表达式。自然界与工程技术中的许多现象，其数学表达式可归结为可归结为常微分方程（组）的定解问题常微分方程（组）的定解问题。一些。一些偏微分方程偏微分方程问题问题也可以也可以转化转化为为常微分方程常微分方程问题来（近似）求解。问题来（近似）求解。Newton最早采最早采用数学方法研究二体问题，其中需要求解的运动方程就是常微用数学方法研究二体问题，其中需要求解的运动方程就是常微分方程。许多著名的数学家，如分方程。许

2、多著名的数学家，如 Bernoulli（家族），（家族），Euler、Gauss、Lagrange和和Laplace等，都遵循历史传统，研究重要等，都遵循历史传统，研究重要的力学问题的数学模型，在这些问题中，许多是常微分方程的的力学问题的数学模型，在这些问题中，许多是常微分方程的求解。作为科学史上的一段佳话，求解。作为科学史上的一段佳话，海王星海王星的发现就是通过对常的发现就是通过对常微分方程的近似计算得到微分方程的近似计算得到的。本章主要介绍常微分方程数值解的。本章主要介绍常微分方程数值解的若干方法。的若干方法。1 1、常微分方程与解、常微分方程与解为为n n阶常微分方程阶常微分方程。0 )

3、, , ,()(nyyyyxF如果函数如果函数 在区间在区间a,ba,b内内n n阶可导，称方程阶可导，称方程)(xyy )(xyy 满足方程的函数满足方程的函数称为微分方程的称为微分方程的解解。则则如如为任意常数）为任意常数）xy2 CCxy(, 2一般称为方程的一般称为方程的通解通解。为方程的解为方程的解, ,12 xy如果如果则有则有10 )(y为方程满足定解条件的解。为方程满足定解条件的解。一一、初值问题的数值解法初值问题的数值解法 102)(yxy CCxy1212 xy方程的通解方程的通解满足定解条件的解满足定解条件的解微分关系（方程）微分关系（方程）解的图示解的图示本教材重点讨论

4、定解问题本教材重点讨论定解问题( (初值问题）初值问题）定解条件（初始条件）定解条件（初始条件） 00yxyyxfy)(),(),(yxf是否能够找到定解问题的解取决于是否能够找到定解问题的解取决于仅有极少数的方程可以通过仅有极少数的方程可以通过“常数变易法常数变易法”、“可分可分离变量法离变量法”等特殊方法求得初等函数形式的解，绝大等特殊方法求得初等函数形式的解，绝大部分方程至今无法理论求解。部分方程至今无法理论求解。如如xyxyeyxyxyy 212,),sin(sin等等等等2 2、数值解的思想、数值解的思想（1 1）将连续变量）将连续变量 离散为离散为,bax bxxxxank 10n

5、kxyykk,)(21 （2 2）用代数的方法求出解函数）用代数的方法求出解函数 在在 点的近似值点的近似值)(xyy kx)(kxy* *ky)(xyy 数学界关注数学界关注工程师关注工程师关注如果找不到解函数如果找不到解函数数学界还关注：数学界还关注：解的存在性解的存在性解的唯一性解的唯一性解的光滑性解的光滑性解的振动性解的振动性解的周期性解的周期性解的稳定性解的稳定性解的混沌性解的混沌性 求函数求函数 y(x) 在一系列节点在一系列节点 a = x0 x1 xn= b 处的近似值处的近似值 的方法称为微分方程的数值解法。的方法称为微分方程的数值解法。() (1,., )iiyy xin称

6、节点间距称节点间距 为步长，为步长，通常采用通常采用等距节点等距节点，即取，即取 hi = h (常数常数)。) 1,., 0(1 nixxhiii1,nyy称为微分方程的数值解。称为微分方程的数值解。所谓数值解法：所谓数值解法：称称 在区域在区域D上对上对 满足满足Lipschitz条件条件是指是指:1212120. .( ,)( ,), , , ( ),( )Ls tf x yf x yL yyxa byyy xy x ( , )f x yy( , ), ( )( )Dx y axb y xyy x记记3 3、相关定义、相关定义(2) 一般构造方法：一般构造方法：4、 迭代格式的构造迭代格

7、式的构造(1) 构造思想：构造思想：将连续的微分方程及初值条件离散为线性方程将连续的微分方程及初值条件离散为线性方程组加以求解。由于离散化的出发点不同，产生出各种不同的数组加以求解。由于离散化的出发点不同，产生出各种不同的数值方法。基本方法有：有限差分法（数值微分）、有限体积法值方法。基本方法有：有限差分法（数值微分）、有限体积法（数值积分）、有限元法（函数插值）等等。（数值积分）、有限元法（函数插值）等等。 (3) 如何保证迭代公式的稳定性与收敛性如何保证迭代公式的稳定性与收敛性?5、微分方程的数值解法需要解决的主要问题、微分方程的数值解法需要解决的主要问题(1) 如何将微分方程离散化，并建

8、立求其如何将微分方程离散化，并建立求其数值解的迭代公式？数值解的迭代公式？(2) 如何估计迭代公式的局部截断误差与整体误差？如何估计迭代公式的局部截断误差与整体误差？二、初值问题解的存在唯一性二、初值问题解的存在唯一性 考虑一阶常微分方程的初值问题考虑一阶常微分方程的初值问题 /* Initial-Value Problem */: 0)(,),(yaybaxyxfdxdy| ),(),(|2121yyLyxfyxf 则上述则上述IVP存在唯一解。存在唯一解。只要只要 在在 上连续上连续, 且关于且关于 y 满足满足 Lipschitz 条件，条件，( , )f x y1, a bR即存在与即

9、存在与 无关的常数无关的常数 L 使使, x y对任意定义在对任意定义在 上的上的 都成立，都成立，, a b 12,yxyx三三、初值问题的离散化方法初值问题的离散化方法 离散化方法的基本特点是依照某一递推公式，离散化方法的基本特点是依照某一递推公式，值值 ,取取 。按节点从左至右的顺序依次求出按节点从左至右的顺序依次求出 的近似的近似( )iy x(1,., )iyin0y 如果计算如果计算 ，只用到前一步的值，只用到前一步的值 ，则称这则称这类方法为类方法为单步方法单步方法。1iyiy如果计算如果计算 需用到前需用到前r步的值步的值 , ,则称这类方法为则称这类方法为r步方法步方法。1i

10、y11,ii ryy iy6.2 Euler6.2 Euler方法方法kp0p1p1npnpkx0 x1x1nxnx),(111212yxfxxyy 第一步：连续变量离散化第一步：连续变量离散化,nkxxxxx10第二步：用直线步进第二步：用直线步进),(000101yxfxxyy ),(),(nnnnnnnnnnyxhfyyyxfxxyy 111EulerEuler格式格式1 1、EulerEuler格式格式 00yxyyxfy)(),(l 18 18世纪最杰出的数学家之一，世纪最杰出的数学家之一，1313岁岁时入读巴塞尔大学，时入读巴塞尔大学，1515岁大学毕业，岁大学毕业，1616岁获得

11、硕士学位。岁获得硕士学位。l 17271727年年-1741-1741年（年（2020岁岁-34-34岁）在彼岁）在彼得堡科学院从事研究工作，在分析学、得堡科学院从事研究工作，在分析学、数论、力学方面均有出色成就，并应数论、力学方面均有出色成就，并应俄国政府要求，解决了不少地图学、俄国政府要求，解决了不少地图学、造船业等实际问题。造船业等实际问题。l 2424岁晋升物理学教授。岁晋升物理学教授。l 17351735年（年（2828岁）右眼失明。岁）右眼失明。l 1741 1741年年 - 1766- 1766（3434岁岁-59-59岁）任德国科学院物理数学所所岁）任德国科学院物理数学所所长，

12、任职长，任职2525年。在行星运动、刚体运动、热力学、弹道学、人年。在行星运动、刚体运动、热力学、弹道学、人口学、微分方程、曲面微分几何等研究领域均有开创性的工作。口学、微分方程、曲面微分几何等研究领域均有开创性的工作。l 17661766年应沙皇礼聘重回彼得堡，在年应沙皇礼聘重回彼得堡，在17711771年（年（6464岁）左眼失岁）左眼失明。明。l EulerEuler是数学史上最多产的数学家，平均以每年是数学史上最多产的数学家，平均以每年800800页的速页的速度写出创造性论文。他去世后，人们用度写出创造性论文。他去世后，人们用3535年整理出他的研究成年整理出他的研究成果果7474卷。

13、卷。 在假设在假设 yi = y(xi)，即第，即第 i 步计算是精确的前提步计算是精确的前提下，考虑的截断误差下，考虑的截断误差 Ri = y(xi+1) yi+1 称为称为局部截断局部截断误差误差 /* local truncation error */。定义定义2.2 若某算法的局部截断误差为若某算法的局部截断误差为O(hp+1)，则称该，则称该 算法有算法有p 阶精度。阶精度。定义定义2.12、欧拉法的局部截断误差、欧拉法的局部截断误差 欧拉法的局部截断误差：欧拉法的局部截断误差：11()iiiRy xy23()()2ihyxO hRi 的的主项主项/* leading term */

14、欧拉法具有欧拉法具有 1 阶精度阶精度。232 ( )( )( )() ( ,)hiiiiiiy xhy xy xO hyhf x y( )iiyy x( )( , ( )iiiy xf x y x2()O h例例1:1: 用欧拉公式求解初值问题用欧拉公式求解初值问题 2201.201yxyxy （）取步长取步长 。 0.1h 解解: : 应用应用EulerEuler公式于题给初值问题的具体形式为：公式于题给初值问题的具体形式为： 2120,1,.,1101iiiiyyhx yiy 其中其中 。0.1ixi计算结果列于下表：计算结果列于下表： iixiy iy xiiy xy 12345678

15、91011120.10.20.30.40.50.60.70.80.91.01.11.21.0000000.9800000.9415840.8883890.8252500.7571470.6883540.6220180.5601130.5036420.4529110.4077830.9900990.9615380.9174310.8630690.8000000.7352940.6711410.6097560.5524860.5000000.4524890.4098360.0099010.0184620.0241530.0263200.0252500.0218520.0172130.0122620

16、.0076260.0036420.0004220.002053可用来检验近似解的准确程度。可用来检验近似解的准确程度。 进行计算，数值解已达到了一定的精度。进行计算，数值解已达到了一定的精度。这个初值问题的准确解为这个初值问题的准确解为 ， 21 1y xx从上表最后一列，我们看到取步长从上表最后一列，我们看到取步长0.1h 3、 欧拉公式的改进：欧拉公式的改进： 隐式欧拉法隐式欧拉法 /* implicit Euler method */向后差商近似导数向后差商近似导数hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy 由于未知数由于未知数 yi+1 同时出现在

17、等式的两边，不能直接同时出现在等式的两边，不能直接得到，故称为得到，故称为隐式隐式 /* implicit */ 欧拉公式，而前者欧拉公式，而前者称为称为显式显式 /* explicit */ 欧拉公式。欧拉公式。111,0,1iiiiyhf xiyyn 一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：11()iiiRy xy23( )()2ihy xO h即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。 梯形公式梯形公式 / /* *trapezoid formula trapezoid formula

18、 * */ / 显、隐式两种算法的显、隐式两种算法的平均平均) 1,., 0(),(),(2111 niyxfyxfhyyiiiiii注：注：梯形公式的局部截断误差梯形公式的局部截断误差 ，311iiiRy xyO h即梯形公式即梯形公式具有具有2 阶精度阶精度，比欧拉方法有了进步。，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式公式隐式公式，计算时不得不用到，计算时不得不用到迭代法。迭代法。中点欧拉公式中点欧拉公式 /* midpoint formula */中心差商近似导数中心差商近似导数hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyx

19、y 1,., 1),(211 niyxfhyyiiii假设假设 , 则可以导出则可以导出即中点公式具有即中点公式具有 2 阶精度。阶精度。)(),(11iiiixyyxyy )()(311hOyxyRiii 方方 法法 显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式中点公式中点公式简单简单精度低精度低稳定性最好稳定性最好精度低精度低, 计算量大计算量大精度提高精度提高计算量大计算量大精度提高精度提高, 显式显式多一个初值多一个初值, 可能影响精度可能影响精度 改进欧拉法改进欧拉法 /* modified Eulers method */Step 1: 先用显式欧拉公式作预测，算出先用显式欧拉

20、公式作预测，算出Step 2: 再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1 ny1,nnnnyyhf xy111,2nnnnnnyyyyhf xf x11121211()2(,) (0,1,2,.)Euler, (,)(,(,)2,()nnnnnnnnnnnnnnyykkkhhyyf xyf xh yhf xf xynkhfyxh yk上式还常写成该式称为改进方法 亦可写成注注: :此法亦称为此法亦称为预测预测- -校正法校正法 / /* * predictor-corrector method predictor-corrector method * *

21、/ /可以证明该算法可以证明该算法具有具有 2 阶精度阶精度，同时可以看到它，同时可以看到它是个是个单步单步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。后面将看到，它的。后面将看到，它的稳定性高于稳定性高于显式欧拉法。显式欧拉法。改进的欧拉法改进的欧拉法111(,)(,)2nnnnnnhyyf xyf xy在实际计算时，可将欧拉法与梯形法则相结合，在实际计算时，可将欧拉法与梯形法则相结合，计算公式为计算公式为 ,.2, 1, 0),(),(2),()(11)1(1)0(1kyxfyxfhyyyxhfyyknnnnnknnnnn应用改进欧拉法应用改进欧拉法, ,

22、如果序列如果序列 收敛收敛, ,)1(1)0(1 nnyy它的极限便满足方程它的极限便满足方程111(,)(,)2nnnnnnhyyf xyf xy改进欧拉法的截断误差改进欧拉法的截断误差)(0)(311hyxynn 因此，改进欧拉法公式具有因此，改进欧拉法公式具有 2 2 阶精度阶精度例例2:2: 用改进用改进Euler公式求解例公式求解例1中的初值问题，中的初值问题， 取步长取步长 。0.1h 解：解：对此初值问题采用改进对此初值问题采用改进EulerEuler公式，公式， 其具体形式为其具体形式为 21( )211111111( ,)2(,)2()12piiiiiiicppiiiiiii

23、pciiiyyhf x yyhx yyyhf xyyhxyyyy 计算结果列于下表：计算结果列于下表：例例1:1: 用欧拉公式求解初值问题用欧拉公式求解初值问题 2201.201yxyxy （）0,1,.,11i 01y iixiy 1piy 1ciy iiy xy改进的改进的Euler法法 iiy xyEuler法法01234567891011120.00.10.20.30.40.50.60.70.80.91.01.11.21.0000000.9900000.9613660.9172460.8619540.8000340.7355270.6175870.6103990.5532890.500

24、9190.4534790.4108591.0000000.9703890.9243970.8667650.8025170.7360290.6706070.6084430.5507850.4981860.4507350.4082370.9800000.9523330.9100950.8571430.7975510.7350250.6725670.6123550.5557930.5036510.4562230.4134810.0000000.0000990.0001730.0001850.0001150.0000340.0002330.0004460.0006430.0008030.0009190

25、.0009900.0010230.0000000.0099010.0184620.0241530.0263200.0252500.0218520.0172130.0122620.0076260.0036420.0004220.0020532(0)101EulerEuler,12dyxydxyyxyx练习 设初值问题试分别用法和改进法求解 并与精确解进行比较。 xEuler法y改进的Euler法y精确解01.000000 1.0000001.0000000.11.000000 1.0959091.0954450.21.191818 1.1840971.1832160.31.277438 1.26

26、62011.2649110.41.358213 1.3433601.3416410.51.435133 1.4164021.4142140.61.508966 1.4859561.4832400.71.580338 1.5525141.5491930.81.649783 1.6164751.6124520.91.717779 1.6781661.6733201.01.784770 1.7378671.732051通过计算结果的比较可以看出，改进的通过计算结果的比较可以看出，改进的Euler方法方法的计算精度比的计算精度比Euler方法要高。方法要高。Euler方法、隐式方法、隐式Euler方法

27、、梯形方法与单步法计算公式的方法、梯形方法与单步法计算公式的显式单步法显式单步法对应关系对应关系),(1hyxhyykkkk 隐式单步法隐式单步法),(11hyyxhyykkkkk ),(),(kkkkyxfhyx ),(1kkkkkyxfhyy 显式显式 Euler方法方法隐式隐式Euler方法方法),(111 kkkkkyxfhyy),(),(2111 kkkkkkkyxfyxfhyy梯形方法梯形方法(隐式隐式),(),(111 kkkkkyxfhyyx),(21),(21),(111 kkkkkkkyxfyxfhyyx单步法小结单步法小结 收敛性收敛性 /* Convergency */

28、定义定义 若某算法对于任意固定的若某算法对于任意固定的 x = xi = x0 + i h，当，当 h0 ( 同时同时 i ) 时有时有 yi y( xi )，则称该算法是，则称该算法是收敛的。收敛的。 收敛性定理收敛性定理若某单步法满足以上条件，则该方法是收敛的若某单步法满足以上条件，则该方法是收敛的则该单步法的整体截断误差为：则该单步法的整体截断误差为：若单步法若单步法 具有具有 p 阶精度，阶精度，且增量函数且增量函数 关于关于 y 满足：满足：1(, )nnnnyyhxy h ( , , )x y h ()()pnny xyO h( , , )( , , )x y hx y hL yy

29、 Lipschitz 条件：条件： 初值初值 y0 是准确的是准确的 稳定性稳定性 /* Stability */例：例：考察初值问题考察初值问题 在区间在区间0, 0.5上的解上的解.分别用欧拉显、隐式格式和改进的欧拉格分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。式计算数值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精确解精确解改进欧拉法改进欧拉法 欧拉隐式欧拉隐式欧拉显式欧拉显式 节点节点 xixey30 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.2000 101 1.00002.5000 10 1 6.2500 1

30、0 21.5625 10 23.9063 10 39.7656 10 41.00002.50006.25001.5626 1013.9063 1019.7656 1011.00004.9787 10 22.4788 10 31.2341 10 46.1442 10 63.0590 10 7What is wrong ?! An Engineer complains: Math theorems are so unstable that a small perturbation on the conditions will cause a crash on the conclusions!定义定

31、义若某算法在计算过程中任一步产生的误差在若某算法在计算过程中任一步产生的误差在以后的计算中都以后的计算中都逐步衰减逐步衰减，则称该算法是，则称该算法是绝对稳定绝对稳定的的 /*absolutely stable */。一般分析时为简单起见，只考虑一般分析时为简单起见，只考虑试验方程试验方程 /* test equation */ 0yy 常数，可以常数，可以是复数是复数我们称我们称算法算法A 比算法比算法B 稳定稳定，就是指，就是指 A 的绝对稳定的绝对稳定区域比区域比 B 的的大大。当步长取为当步长取为 h 时，将某算法应用于上式，并假设只在时，将某算法应用于上式，并假设只在初值产生误差初值

32、产生误差 ，则若此误差以后逐步衰减，则若此误差以后逐步衰减,000yy 就称该算法相对于就称该算法相对于 绝对稳定，绝对稳定， 的全体构成的全体构成h h h绝对稳定区域绝对稳定区域。yy 例：例：考察显式欧拉法考察显式欧拉法110(1)iiiiyyhyhy000yy 011)1 (yhyii 01111)1 ( iiiihyy0-1-2ReImg由此可见，要保证初始误差由此可见，要保证初始误差 0 以后逐步衰减，以后逐步衰减，必须满足：必须满足：hh 1|1| h例：例：考察隐式欧拉法考察隐式欧拉法11 iiiyhyy iiyhy 11101111 iih可见绝对稳定区域为：可见绝对稳定区域

33、为：1|1 | h210ReImg注：注：一般来说，隐式欧拉法的绝对稳定性比同阶一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。的显式法的好。定理定理 某单步法用于解模型方程某单步法用于解模型方程y =y，若得到的解，若得到的解yn+1=E(h)yn，满足，满足|E(h)|1，则称该单步法是，则称该单步法是绝对绝对稳定稳定的的. 在平面上在平面上, 使使|E(h)|1的变量围成的区域，的变量围成的区域，称为称为绝对稳定区域绝对稳定区域，它与实轴的交称为，它与实轴的交称为绝对稳定区绝对稳定区间间.稳定性定理稳定性定理6.3 龙格库塔方法 对许多实际问题来说，欧拉公式与改进欧拉对许多实际问题来

34、说，欧拉公式与改进欧拉公式精度还不能满足要求，为此从另一个角度来分公式精度还不能满足要求，为此从另一个角度来分析这两个公式的特点，从而探索一条构造高精度方析这两个公式的特点，从而探索一条构造高精度方法的途径法的途径. 受改进的Euler方法启发，更一般算式可设为1121211()2(,)(,)(0,1,2,.)nnnnnnyykkkhf xykhf xh ykn123111()-()()nnnnnTy xyO hyy x适当选择参数 ， ， ，使局部截断误差,这里仍假定。改进欧拉法改进欧拉法11 122121(,)(0,1,2,.)(,)nnnnnnyykkkhf xynkhf xh yk11

35、22232()()()(,)(,)(,)()nnnxnnnnynnyy xhy xfxyf xyfxyhO h)(),(),(),()()(),(),(),(:Taylor323122hOyxfyxfyxfhxyhhOyxfhkyxfhyxhfknnynnnnxnnnynnxnn展开式由二元函数11 122121 (,) (,)(0,1,2,.)nnnnnnnyykkkhf xyhykhf xh ykn由于四个参数，三个方程，因此有一个自由参数，即解答不唯一。1122232()()()(,)(,)(,)()nnnxnnnnynnyy xhy xfxyf xyfxyhO hTaylor:与展式相

36、比较得1212122122312()()()()()hnnnny xy xhy xy xO h 这是改进的Euler方法。1211(1),1,22取可得此时算式为1121211()2(,) (,)nnnnnnyykkkhf xykhf xh yk-R K这是二阶方法.121(2)0,1,2取可得此时算式为12121(,)11(,) 22nnnnnnyykkhf xykhf xh yk R-K这也是二阶方法。12132(3),443取可得又有算式11212113)4(,)22(,)33nnnnnnyykkkhf xykhf xh yk（三阶龙格-库塔方法三阶龙格-库塔方法是用三个值 k1, k2, k3 的线性组合112312123221331332112(,)(,)(,)nnnnnnnnyykkkkhf xykhf xhcccababykkhf xh ykb k 要使三阶龙格-库塔方法具有三阶精度，必须使其局部截断误差为 O(h4)将 k1, k2, k3 代入 yn+1 的