微分方程数值解法课件.ppt

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

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

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

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

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

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

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

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

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

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

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

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

13、理出他的研究成年整理出他的研究成果果7474卷。卷。15ppt课件在假设在假设 yi=y(xi)，即第，即第 i 步计算是精确的前提步计算是精确的前提下，考虑的截断误差下，考虑的截断误差 Ri=y(xi+1)yi+1 称为称为局部截断局部截断误差误差/*local truncation error*/。定义定义2.2 若某算法的局部截断误差为若某算法的局部截断误差为O(hp+1)，则称该，则称该 算法有算法有p 阶精度。阶精度。定义定义2.12、欧拉法的局部截断误差、欧拉法的局部截断误差16ppt课件 欧拉法的局部截断误差：欧拉法的局部截断误差：11()iiiRy xy23()()2ihyxO

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

15、xiy 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.9615380.9174310.8630690.8000000.7352940.6711410.6097560.5524860.5000000.4524890.4098360.0099010.0184620.0241530.0263200.0252500.

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

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

18、 /*trapezoid formula trapezoid formula*/显、隐式两种算法的显、隐式两种算法的平均平均)1,.,0(),(),(2111 niyxfyxfhyyiiiiii注：注：梯形公式的局部截断误差梯形公式的局部截断误差 ，311iiiRy xyO h即梯形公式即梯形公式具有具有2 阶精度阶精度，比欧拉方法有了进步。，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式公式隐式公式，计算时不得不用到，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。迭代法，其迭代收敛性与欧拉公式相似。24ppt课件中点欧拉公式中点欧拉公式/*midpoint formula*/

19、中心差商近似导数中心差商近似导数hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,.,1),(211 niyxfhyyiiii假设假设 ,则可以导出则可以导出即中点公式具有即中点公式具有 2 阶精度。阶精度。)(),(11iiiixyyxyy )()(311hOyxyRiii 25ppt课件方方 法法 显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式中点公式中点公式简单简单精度低精度低稳定性最好稳定性最好精度低精度低,计算量大计算量大精度提高精度提高计算量大计算量大精度提高精度提高,显式显式多一个初值多一个初值,可能影响精度可能影响精度26pp

20、t课件 改进欧拉法改进欧拉法/*modified Eulers method*/Step 1:先用显式欧拉公式作预测，算出先用显式欧拉公式作预测，算出Step 2:再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1 ny1,nnnnyyhf xy111,2nnnnnnyyyyhf xf x27ppt课件1121211Euler,1()2(,)(,)(0,1,2,.(,)(,(,).)2nnnnnnnnnnnnnnyykkkhf xykhf xh yhyyf xyf xh yhf xykn上式还常写成该式称为改进方法 亦可写成28ppt课件注注:此法亦称为此法亦称为

21、预测预测-校正法校正法 /*predictor-corrector method predictor-corrector method*/可以证明该算法可以证明该算法具有具有 2 阶精度阶精度，同时可以看到它，同时可以看到它是个是个单步单步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。后面将看到，它的。后面将看到，它的稳定性高于稳定性高于显式欧拉法。显式欧拉法。改进的欧拉法改进的欧拉法111(,)(,)2nnnnnnhyyf xyf xy29ppt课件在实际计算时，可将欧拉法与梯形法则相结合，在实际计算时，可将欧拉法与梯形法则相结合，计算公式为计算公式为 ,.2

22、,1,0),(),(2),()(11)1(1)0(1kyxfyxfhyyyxhfyyknnnnnknnnnn应用改进欧拉法应用改进欧拉法,如果序列如果序列 收敛收敛,)1(1)0(1 nnyy它的极限便满足方程它的极限便满足方程111(,)(,)2nnnnnnhyyf xyf xy30ppt课件改进欧拉法的截断误差改进欧拉法的截断误差)(0)(311hyxynn 因此，改进欧拉法公式具有因此，改进欧拉法公式具有 2 2 阶精度阶精度31ppt课件例例2:2:用改进用改进Euler公式求解例公式求解例1中的初值问题，中的初值问题，取步长取步长 。0.1h 解：解：对此初值问题采用改进对此初值问题

23、采用改进EulerEuler公式，公式，其具体形式为其具体形式为 21()211111111(,)2(,)2()12piiiiiiicppiiiiiiipciiiyyhf x yyhx yyyhf xyyhxyyyy 计算结果列于下表：计算结果列于下表：例例1:1:用欧拉公式求解初值问题用欧拉公式求解初值问题 2201.201yxyxy （）0,1,.,11i 01y 32ppt课件iixiy 1piy 1ciyiiy xy改进的改进的Euler法法iiy xyEuler法法01234567891011120.00.10.20.30.40.50.60.70.80.91.01.11.21.000

24、0000.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.9100950.8571430.7975510.7350250.6725670.6123550.5557930.5036510.4562230.4134810.0000000

25、.0000990.0001730.0001850.0001150.0000340.0002330.0004460.0006430.0008030.0009190.0009900.0010230.0000000.0099010.0184620.0241530.0263200.0252500.0218520.0172130.0122620.0076260.0036420.0004220.00205333ppt课件通过计算结果的比较可以看出，改进的通过计算结果的比较可以看出，改进的Euler方法方法的计算精度比的计算精度比Euler方法要高。方法要高。34ppt课件欧拉法误差概述)()(2)O()(

26、y2Euler式对于)O(h)(2Euler对于显式则精度相对的越高。越大，误差阶，一个算法，局部截断一般431321321hOxyhThxhTxyhTnnnnnn ：对于梯形公式：法隐：法说来35ppt课件6.3 龙格库塔方法 对许多实际问题来说，欧拉公式与改进欧拉对许多实际问题来说，欧拉公式与改进欧拉公式精度还不能满足要求，为此从另一个角度来分公式精度还不能满足要求，为此从另一个角度来分析这两个公式的特点，从而探索一条构造高精度方析这两个公式的特点，从而探索一条构造高精度方法的途径法的途径.36ppt课件受改进的Euler方法启发，更一般算式可设为1121211()2(,)(,)(0,1,

27、2,.)nnnnnnyykkkhf xykhf xh ykn123111()-()()nnnnnTy xyO hyy x适当选择参数，使局部截断误差,这里仍假定。改进欧拉法改进欧拉法11 122121(,)(0,1,2,.)(,)nnnnnnyykkkhf xynkhf xh yk37ppt课件1122232()()()(,)(,)(,)()nnnxnnnnynnyy xhy xfxyf xyfxyhO h)(),(),(),()()(),(),(),(:Taylor323122hOyxfyxfyxfhxyhhOyxfhkyxfhyxhfknnynnnnxnnnynnxnn展开式由二元函数11

28、 122121 (,)(,)(0,1,2,.)nnnnnnnyykkkhf xyhykhf xh ykn38ppt课件由于四个参数，三个方程，因此有一个自由参数，即解答不唯一。1122232()()()(,)(,)(,)()nnnxnnnnynnyy xhy xfxyf xyfxyhO hTaylor:与展式相比较得1212122122312()()()()()hnnnny xy xhy xy xO h39ppt课件 这是改进的Euler方法。1211(1),1,22取可得此时算式为1121211()2(,)(,)nnnnnnyykkkhf xykhf xh yk40ppt课件-R K这是二阶

29、方法.121(2)0,1,2取可得此时算式为12121(,)11(,)22nnnnnnyykkhf xykhf xh yk41ppt课件 R-K这也是二阶方法。12132(3),443取可得又有算式11212113)4(,)22(,)33nnnnnnyykkkhf xykhf xh yk（42ppt课件三阶龙格-库塔方法三阶龙格-库塔方法是用三个值 k1,k2,k3 的线性组合112312123221331332112(,)(,)(,)nnnnnnnnyykkkkhf xykhf xhcccababykkhf xh ykb k 要使三阶龙格-库塔方法具有三阶精度，必须使其局部截断误差为 O(h

30、4)将 k1,k2,k3 代入 yn+1 的表达式中，在(xn,yn)处用二元泰勒公式展开，与 y(xn+1)在 xn 处的泰勒展开式比较43ppt课件 类似二阶龙格-库塔方法的推导过程，8 个待定系数 c1,c2,c3,a2,a3,b21,b31,b32 应满足：123221331322233222233232211 21 31 6cccababbc ac ac ac ac b a 8 个未知参数，个未知参数，6 个个方程，有无穷多组解方程，有无穷多组解14166611211231213122(,)(,)(,)2nnnnnnnnyykkkkhf xykhf xh ykkhf xh ykk 三

31、阶龙格库塔公式三阶龙格库塔公式44ppt课件四阶四阶Runge-Kutta方法方法1123412132335RK1(22)6(,)11(,)2211(,)22(,),()nnnnnnnnnnyykkkkkhf xykhf xh ykkhf xh ykkhf xh ykO h类似以上推导过程，可求得经典四阶算式该公式结构整齐 局部截断误差为。45ppt课件附注：附注：二阶二阶Runge-Kutta方法的局部截断误差方法的局部截断误差 只能达到只能达到3();O h 五阶五阶Runge-Kutta方法的局部截断误差方法的局部截断误差 只能达到只能达到 四阶四阶Runge-Kutta方法的局部截断误

32、差方法的局部截断误差 只能达到只能达到 三阶三阶Runge-Kutta方法的局部截断误差方法的局部截断误差 只能达到只能达到 4();O h5();O h6().O h46ppt课件附注：附注：龙格龙格-库塔法的主要运算在于计算库塔法的主要运算在于计算 的值，即计的值，即计算算 的值。的值。Butcher 于于1965年给出了计算量与可达年给出了计算量与可达到的最高精度阶数的关系：到的最高精度阶数的关系：iKf753可达到的最高精度可达到的最高精度642每步须算每步须算Ki 的个数的个数3()O h4()O h5()O h7()O h8()O h6()O h)(2nhO8n 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精库塔法的导出基于泰勒展开，故精太好的解，最好采用太好的解，最好采用低阶算法低阶算法而将步长而将步长h 取小取小。度主要受解函数的光滑性影响。对于光滑性不度主要受解函数的光滑性影响。对于光滑性不47ppt课件感谢您的听讲！感谢您的听讲！48ppt课件