微分方程数值解法

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

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

3、程与解为为n n阶常微分方程阶常微分方程。0 ), , ,()(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、数值解的思想

5、、数值解的思想（1 1）将连续变量）将连续变量 离散为离散为,bax bxxxxank 10nkxyykk,)(21 （2 2）用代数的方法求出解函数）用代数的方法求出解函数 在在 点的近似值点的近似值)(xyy kx)(kxy* *ky)(xyy 数学界关注数学界关注工程师关注工程师关注如果找不到解函数如果找不到解函数数学界还关注：数学界还关注：解的存在性解的存在性解的唯一性解的唯一性解的光滑性解的光滑性解的振动性解的振动性解的周期性解的周期性解的稳定性解的稳定性解的混沌性解的混沌性微分方程数值解法 求函数求函数 y(x) 在一系列节点在一系列节点 a = x0 x1 xn= b 处的近似值

6、处的近似值 的方法称为微分方程的数值解法。的方法称为微分方程的数值解法。() (1,., )iiyy xin称节点间距称节点间距 为步长，为步长，通常采用通常采用等距节点等距节点，即取，即取 hi = h (常数常数)。) 1,., 0(1 nixxhiii1,nyy称为微分方程的数值解。称为微分方程的数值解。所谓数值解法：所谓数值解法：微分方程数值解法称称 在区域在区域D上对上对 满足满足Lipschitz条件条件是指是指:1212120. .( ,)( ,), , , ( ),( )Ls tf x yf x yL yyxa byyy xy x ( , )f x yy( , ), ( )(

7、)Dx y axb y xyy x记记3 3、相关定义、相关定义微分方程数值解法(2) 一般构造方法：一般构造方法：4、 迭代格式的构造迭代格式的构造(1) 构造思想：构造思想：将连续的微分方程及初值条件离散为线性方程将连续的微分方程及初值条件离散为线性方程组加以求解。由于离散化的出发点不同，产生出各种不同的数组加以求解。由于离散化的出发点不同，产生出各种不同的数值方法。基本方法有：有限差分法（数值微分）、有限体积法值方法。基本方法有：有限差分法（数值微分）、有限体积法（数值积分）、有限元法（函数插值）等等。（数值积分）、有限元法（函数插值）等等。 微分方程数值解法(3) 如何保证迭代公式的稳

8、定性与收敛性如何保证迭代公式的稳定性与收敛性?5、微分方程的数值解法需要解决的主要问题、微分方程的数值解法需要解决的主要问题(1) 如何将微分方程离散化，并建立求其如何将微分方程离散化，并建立求其数值解的迭代公式？数值解的迭代公式？(2) 如何估计迭代公式的局部截断误差与整体误差？如何估计迭代公式的局部截断误差与整体误差？微分方程数值解法二、初值问题解的存在唯一性二、初值问题解的存在唯一性 考虑一阶常微分方程的初值问题考虑一阶常微分方程的初值问题 /* Initial-Value Problem */: 0)(,),(yaybaxyxfdxdy| ),(),(|2121yyLyxfyxf 则上

9、述则上述IVP存在唯一解。存在唯一解。只要只要 在在 上连续上连续, 且关于且关于 y 满足满足 Lipschitz 条件，条件，( , )f x y1, a bR即存在与即存在与 无关的常数无关的常数 L 使使, x y对任意定义在对任意定义在 上的上的 都成立，都成立，, a b 12,yxyx微分方程数值解法三、初值问题的离散化方法三、初值问题的离散化方法 离散化方法的基本特点是依照某一递推公式，离散化方法的基本特点是依照某一递推公式，值值 ,取取 。按节点从左至右的顺序依次求出按节点从左至右的顺序依次求出 的近似的近似( )iy x(1,., )iyin0y 如果计算如果计算 ，只用到

10、前一步的值，只用到前一步的值 ，则称这则称这类方法为类方法为单步方法单步方法。1iyiy如果计算如果计算 需用到前需用到前r步的值步的值 , ,则称这类方法为则称这类方法为r步方法步方法。1iy11,ii ryy iy微分方程数值解法6.2 Euler6.2 Euler方法方法kp0p1p1npnpkx0 x1x1nxnx),(111212yxfxxyy 第一步：连续变量离散化第一步：连续变量离散化,nkxxxxx10第二步：用直线步进第二步：用直线步进),(000101yxfxxyy ),(),(nnnnnnnnnnyxhfyyyxfxxyy 111EulerEuler格式格式1 1、Eul

11、erEuler格式格式 00yxyyxfy)(),(微分方程数值解法l 18 18世纪最杰出的数学家之一，世纪最杰出的数学家之一，1313岁岁时入读巴塞尔大学，时入读巴塞尔大学，1515岁大学毕业，岁大学毕业，1616岁获得硕士学位。岁获得硕士学位。l 17271727年年-1741-1741年（年（2020岁岁-34-34岁）在彼岁）在彼得堡科学院从事研究工作，在分析学、得堡科学院从事研究工作，在分析学、数论、力学方面均有出色成就，并应数论、力学方面均有出色成就，并应俄国政府要求，解决了不少地图学、俄国政府要求，解决了不少地图学、造船业等实际问题。造船业等实际问题。l 2424岁晋升物理学教

12、授。岁晋升物理学教授。l 17351735年（年（2828岁）右眼失明。岁）右眼失明。微分方程数值解法l 1741 1741年年 - 1766- 1766（3434岁岁-59-59岁）任德国科学院物理数学所所岁）任德国科学院物理数学所所长，任职长，任职2525年。在行星运动、刚体运动、热力学、弹道学、人年。在行星运动、刚体运动、热力学、弹道学、人口学、微分方程、曲面微分几何等研究领域均有开创性的工作。口学、微分方程、曲面微分几何等研究领域均有开创性的工作。l 17661766年应沙皇礼聘重回彼得堡，在年应沙皇礼聘重回彼得堡，在17711771年（年（6464岁）左眼失岁）左眼失明。明。l Eu

13、lerEuler是数学史上最多产的数学家，平均以每年是数学史上最多产的数学家，平均以每年800800页的速页的速度写出创造性论文。他去世后，人们用度写出创造性论文。他去世后，人们用3535年整理出他的研究成年整理出他的研究成果果7474卷。卷。 微分方程数值解法在假设在假设 yi = y(xi)，即第，即第 i 步计算是精确的前提步计算是精确的前提下，考虑的截断误差下，考虑的截断误差 Ri = y(xi+1) yi+1 称为称为局部截断局部截断误差误差 /* local truncation error */。定义定义2.2 若某算法的局部截断误差为若某算法的局部截断误差为O(hp+1)，则称

14、该，则称该 算法有算法有p 阶精度。阶精度。定义定义2.12、欧拉法的局部截断误差、欧拉法的局部截断误差微分方程数值解法 欧拉法的局部截断误差：欧拉法的局部截断误差：11()iiiRy xy23()()2ihyxO hRi 的的主项主项/* leading term */欧拉法具有欧拉法具有 1 阶精阶精度度。232 ( )( )( )() ( ,)hiiiiiiy xhy xy xO hyhf x y( )iiyy x( )( , ( )iiiy xf x y x2()O h微分方程数值解法例例1:1: 用欧拉公式求解初值问题用欧拉公式求解初值问题 2201.201yxyxy （）取步长取步

15、长 。 0.1h 解解: : 应用应用EulerEuler公式于题给初值问题的具体形式为：公式于题给初值问题的具体形式为： 2120,1,.,1101iiiiyyhx yiy 其中其中 。0.1ixi计算结果列于下表：计算结果列于下表： 微分方程数值解法iixiy iy xiiy xy 1234567891011120.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

16、.9615380.9174310.8630690.8000000.7352940.6711410.6097560.5524860.5000000.4524890.4098360.0099010.0184620.0241530.0263200.0252500.0218520.0172130.0122620.0076260.0036420.0004220.002053微分方程数值解法可用来检验近似解的准确程度。可用来检验近似解的准确程度。 进行计算，数值解已达到了一定的精度。进行计算，数值解已达到了一定的精度。这个初值问题的准确解为这个初值问题的准确解为 ， 21 1y xx从上表最后一列，我们看

17、到取步长从上表最后一列，我们看到取步长0.1h 微分方程数值解法3、 欧拉公式的改进：欧拉公式的改进： 隐式欧拉法隐式欧拉法 /* implicit Euler method */向后差商近似导数向后差商近似导数hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy 微分方程数值解法由于未知数由于未知数 yi+1 同时出现在等式的两边，不能直接同时出现在等式的两边，不能直接得到，故称为得到，故称为隐式隐式 /* implicit */ 欧拉公式，而前者欧拉公式，而前者称为称为显式显式 /* explicit */ 欧拉公式。欧拉公式。111,0,1iiiiyhfxiy

18、yn微分方程数值解法一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：11()iiiRy xy23( )()2ihy xO h即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶精度。微分方程数值解法 梯形公式梯形公式 / /* *trapezoid formula trapezoid formula * */ / 显、隐式两种算法的显、隐式两种算法的平均平均) 1,., 0(),(),(2111 niyxfyxfhyyiiiiii注：注：梯形公式的局部截断误差梯形公式的局部截断误差 ，311iiiRy xyO h即

19、梯形公式即梯形公式具有具有2 阶精度阶精度，比欧拉方法有了进步。，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式公式隐式公式，计算时不得不用到，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。迭代法，其迭代收敛性与欧拉公式相似。微分方程数值解法中点欧拉公式中点欧拉公式 /* midpoint formula */中心差商近似导数中心差商近似导数hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,., 1),(211 niyxfhyyiiii假设假设 , 则可以导出则可以导出即中点公式具有即中点公式具有 2 阶精度。阶精度。)(),(

20、11iiiixyyxyy )()(311hOyxyRiii 微分方程数值解法方方 法法 显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式中点公式中点公式简单简单精度低精度低稳定性最好稳定性最好精度低精度低, 计算量大计算量大精度提高精度提高计算量大计算量大精度提高精度提高, 显式显式多一个初值多一个初值, 可能影响精度可能影响精度微分方程数值解法 改进欧拉法改进欧拉法 /* modified Eulers method */Step 1: 先用显式欧拉公式作预测，算出先用显式欧拉公式作预测，算出Step 2: 再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1 n

21、y1,nnnnyyhf xy111,2nnnnnnyyyyhfxfx微分方程数值解法注注: :此法亦称为此法亦称为预测预测- -校正法校正法 / /* * predictor-corrector method predictor-corrector method * */ /可以证明该算法可以证明该算法具有具有 2 阶精度阶精度，同时可以看到它，同时可以看到它是个是个单步单步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。后面将看到，它的。后面将看到，它的稳定性高于稳定性高于显式欧拉法。显式欧拉法。改进的欧拉法改进的欧拉法111(,)(,)2nnnnnnhyyf

22、xyf xy微分方程数值解法在实际计算时，可将欧拉法与梯形法则相结合，在实际计算时，可将欧拉法与梯形法则相结合，计算公式为计算公式为 ,.2, 1, 0),(),(2),()(11)1(1)0(1kyxfyxfhyyyxhfyyknnnnnknnnnn应用改进欧拉法应用改进欧拉法, ,如果序列如果序列 收敛收敛, ,)1(1)0(1 nnyy它的极限便满足方程它的极限便满足方程111(,)(,)2nnnnnnhyyf xyf xy微分方程数值解法改进欧拉法的截断误差改进欧拉法的截断误差)(0)(311hyxynn 因此，改进欧拉法公式具有因此，改进欧拉法公式具有 2 2 阶精度阶精度微分方程数

23、值解法例例2:2: 用改进用改进Euler公式求解例公式求解例1中的初值问题，中的初值问题， 取步长取步长 。0.1h 解：解：对此初值问题采用改进对此初值问题采用改进EulerEuler公式，公式， 其具体形式为其具体形式为 21( )211111111( ,)2(,)2()12piiiiiiicppiiiiiiipciiiyyhf x yyhx yyyhf xyyhxyyyy 计算结果列于下表：计算结果列于下表：例例1:1: 用欧拉公式求解初值问题用欧拉公式求解初值问题 2201.201yxyxy （）0,1,.,11i 01y 微分方程数值解法iixiy 1piy 1ciyiiy xy改

24、进的改进的Euler法法iiy xyEuler法法01234567891 01 11 20.00.10.20.30.40.50.60.70.80.91.01.11.21.0000000.9900000.9613660.9172460.8619540.8000340.7355270.6175870.6103990.5532890.5009190.4534790.4108591.0000000.9703890.9243970.8667650.8025170.7360290.6706070.6084430.5507850.4981860.4507350.4082370.9800000.9523330

25、.9100950.8571430.7975510.7350250.6725670.6123550.5557930.5036510.4562230.4134810.0000000.0000990.0001730.0001850.0001150.0000340.0002330.0004460.0006430.0008030.0009190.0009900.0010230.0000000.0099010.0184620.0241530.0263200.0252500.0218520.0172130.0122620.0076260.0036420.0004220.002053微分方程数值解法通过计算结

26、果的比较可以看出，改进的通过计算结果的比较可以看出，改进的Euler方法方法的计算精度比的计算精度比Euler方法要高。方法要高。微分方程数值解法欧拉法误差概述)()(2)O()(y2Euler式对于)O(h)(2Euler对于显式则精度相对的越高。越大，误差阶，一个算法，局部截断一般431321321hOxyhThxhTxyhTnnnnnn ：对于梯形公式：法隐：法说来微分方程数值解法6.3 龙格库塔方法 对许多实际问题来说，欧拉公式与改进欧拉对许多实际问题来说，欧拉公式与改进欧拉公式精度还不能满足要求，为此从另一个角度来分公式精度还不能满足要求，为此从另一个角度来分析这两个公式的特点，从而

27、探索一条构造高精度方析这两个公式的特点，从而探索一条构造高精度方法的途径法的途径. 微分方程数值解法受改进的Euler方法启发，更一般算式可设为1121211()2(,)(,)(0,1,2,.)nnnnnnyykkkhf xykhf xh ykn123111()-()()nnnnnTy xyO hyy x适当选择参数 ， ， ，使局部截断误差,这里仍假定。改进欧拉法改进欧拉法11 122121(,)(0,1,2,.)(,)nnnnnnyykkkhf xynkhf xh yk微分方程数值解法 这是改进的Euler方法。1211(1),1,22取可得此时算式为1121211()2(,) (,)nn

28、nnnnyykkkhf 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 的表达式中，在 (xn, yn) 处用二元泰勒公式展开，与 y(xn+1) 在 xn 处的泰勒展开式比较微分方程数值解法 类似二阶龙格-库塔方法的推导过程，8 个待定系数 c1, c2, c3, a2, a3, b21, b31, b32 应满足： 123221331322233222233232