微分方程数值解法

会计学1微分方程数值解法微分方程数值解法1 1、常微分方程与解、常微分方程与解为为n n阶常微分方程阶常微分方程。0 ), , ,()(nyyyyxF如果函数如果函数 在区间在区间a,ba,b内内n n阶可导，称方程阶可导，称方程)(xyy )(xyy 满足方程的函数满足方程的函数称为微分方程的称为微分方程的解解。则则如如为任意常数）为任意常数）xy2 CCxy(, 2一般称为方程的一般称为方程的通解通解。为方程的解为方程的解。12 xy如果如果则有则有10 )( y为方程满足定解条件的解。为方程满足定解条件的解。一一、初值问题的数值解法初值问题的数值解法第2页/共48页 102)(yxy CCxy1212 xy方程的通解方程的通解满足定解条件的解满足定解条件的解微分关系（方程）微分关系（方程）解的图示解的图示第3页/共48页本教材重点讨论定解问题本教材重点讨论定解问题( (初值问题）初值问题）定解条件（初始条件定解条件（初始条件） 00yxyyxfy)(),(),(yxf是否能够找到定解问题的解取决于是否能够找到定解问题的解取决于仅有极少数的方程可以通过仅有极少数的方程可以通过“常数变易法常数变易法”、“可分可分离变量法离变量法”等特殊方法求得初等函数形式的解，绝大等特殊方法求得初等函数形式的解，绝大部分方程至今无法理论求解。部分方程至今无法理论求解。如如xyxyeyxyxyy 212,),sin(sin等等等等第4页/共48页2 2、数值解的思想、数值解的思想（1 1）将连续变量）将连续变量 离散为离散为,bax bxxxxank 10nkxyykk,)(21 （2 2）用代数的方法求出解函数）用代数的方法求出解函数 在在 点的近似值点的近似值)(xyy kx)(kxy* *ky)(xyy 数学界关注数学界关注工程师关注工程师关注如果找不到解函数如果找不到解函数数学界还关注：数学界还关注：解的存在性解的存在性解的唯一性解的唯一性解的光滑性解的光滑性解的振动性解的振动性解的周期性解的周期性解的稳定性解的稳定性解的混沌性解的混沌性第5页/共48页 求函数求函数 y(x) 在一系列节点在一系列节点 a = x0 x1 xn= b 处的近似值处的近似值 的方法称为微分方程的数值解法。的方法称为微分方程的数值解法。() (1,., )iiyy xin称节点间距称节点间距 为步长，为步长，通常采用通常采用等距节点等距节点，即取，即取 hi = h (常数常数)。) 1,., 0(1 nixxhiii1,nyy称为微分方程的数值解。称为微分方程的数值解。所谓数值解法：所谓数值解法：第6页/共48页称称 在区域在区域D上对上对 满足满足Lipschitz条件条件是指是指:1212120. .( ,)( ,), , , ( ),( )Ls tf x yf x yL yyxa byyy xy x ( , )f x yy( , ), ( )( )Dx y axb y xyy x记记3 3、相关定义、相关定义第7页/共48页(2) 一般构造方法：一般构造方法：4、 迭代格式的构造迭代格式的构造(1) 构造思想：构造思想：将连续的微分方程及初值条件离散为线性方程将连续的微分方程及初值条件离散为线性方程组加以求解。由于离散化的出发点不同，产生出各种不同的数组加以求解。由于离散化的出发点不同，产生出各种不同的数值方法。基本方法有：有限差分法（数值微分）、有限体积法值方法。基本方法有：有限差分法（数值微分）、有限体积法（数值积分）、有限元法（函数插值）等等。（数值积分）、有限元法（函数插值）等等。 第8页/共48页(3) 如何保证迭代公式的稳定性与收敛如何保证迭代公式的稳定性与收敛性性?5、微分方程的数值解法需要解决的主要问题、微分方程的数值解法需要解决的主要问题(1) 如何将微分方程离散化，并建立求其如何将微分方程离散化，并建立求其数值解的迭代公式数值解的迭代公式？(2) 如何估计迭代公式的局部截断误差与整体误差？如何估计迭代公式的局部截断误差与整体误差？第9页/共48页二、初值问题解的存在唯一性二、初值问题解的存在唯一性 考虑一阶常微分方程的初值问题考虑一阶常微分方程的初值问题 /* Initial-Value Problem */: 0)(,),(yaybaxyxfdxdy| ),(),(|2121yyLyxfyxf 则上述则上述IVP存在唯一解。存在唯一解。只要只要 在在 上连续上连续, 且关于且关于 y 满足满足 Lipschitz 条件，条件，( , )f x y1, a bR即存在与即存在与 无关的常数无关的常数 L 使使, x y对任意定义在对任意定义在 上的上的 都成立，都成立，, a b 12,yxyx第10页/共48页三三、初值问题的离散化方法初值问题的离散化方法 离散化方法的基本特点是依照某一递推公式，离散化方法的基本特点是依照某一递推公式，值值 ,取取 。按节点从左至右的顺序依次求出按节点从左至右的顺序依次求出 的近似的近似( )iy x(1,., )iyin0y 如果计算如果计算 ，只用到前一步的值，只用到前一步的值 ，则称则称这类方法为这类方法为单步方法单步方法。1iyiy如果计算如果计算 需用到前需用到前r步的值步的值 , ,则称这类方法为则称这类方法为r步方法步方法。1iy11,ii ryy iy第11页/共48页6.2 Euler6.2 Euler方法方法kp0p1p1npnpkx0 x1x1nxnx),(111212yxfxxyy 第一步：连续变量离散化第一步：连续变量离散化,nkxxxxx10第二步：用直线步进第二步：用直线步进),(000101yxfxxyy ),(),(nnnnnnnnnnyxhfyyyxfxxyy 111EulerEuler格式格式1 1、EulerEuler格式格式 00yxyyxfy)(),(第12页/共48页l 18 18世纪最杰出的数学家之一，世纪最杰出的数学家之一，1313岁岁时入读巴塞尔大学，时入读巴塞尔大学，1515岁大学毕业，岁大学毕业，1616岁获得硕士学位。岁获得硕士学位。l 17271727年年-1741-1741年（年（2020岁岁-34-34岁）在彼岁）在彼得堡科学院从事研究工作，在分析学得堡科学院从事研究工作，在分析学、数论、力学方面均有出色成就，并、数论、力学方面均有出色成就，并应俄国政府要求，解决了不少地图学应俄国政府要求，解决了不少地图学、造船业等实际问题。、造船业等实际问题。l 2424岁晋升物理学教授。岁晋升物理学教授。l 17351735年（年（2828岁）右眼失明。岁）右眼失明。第13页/共48页l 1741 1741年年 - 1766- 1766（3434岁岁-59-59岁）任德国科学院物理数学所所岁）任德国科学院物理数学所所长，任职长，任职2525年。在行星运动、刚体运动、热力学、弹道学、年。在行星运动、刚体运动、热力学、弹道学、人口学、微分方程、曲面微分几何等研究领域均有开创性的人口学、微分方程、曲面微分几何等研究领域均有开创性的工作。工作。l 17661766年应沙皇礼聘重回彼得堡，在年应沙皇礼聘重回彼得堡，在17711771年（年（6464岁）左眼失岁）左眼失明。明。l EulerEuler是数学史上最多产的数学家，平均以每年是数学史上最多产的数学家，平均以每年800800页的速页的速度写出创造性论文。他去世后，人们用度写出创造性论文。他去世后，人们用3535年整理出他的研究年整理出他的研究成果成果7474卷。卷。 第14页/共48页在假设在假设 yi = y(xi)，即第，即第 i 步计算是精确的前提步计算是精确的前提下，考虑的截断误差下，考虑的截断误差 Ri = y(xi+1) yi+1 称为称为局部截断局部截断误差误差 /* local truncation error */。定义定义2.2 若某算法的局部截断误差为若某算法的局部截断误差为O(hp+1)，则称该，则称该 算法有算法有p 阶精度。阶精度。定义定义2.12、欧拉法的局部截断误差、欧拉法的局部截断误差第15页/共48页 欧拉法的局部截断误差：欧拉法的局部截断误差：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第16页/共48页例例1:1: 用欧拉公式求解初值问题用欧拉公式求解初值问题 2201.201yxyxy （）取步长取步长 。 0.1h 解解: : 应用应用EulerEuler公式于题给初值问题的具体形式为：公式于题给初值问题的具体形式为： 2120,1,.,1101iiiiyyhx yiy 其中其中 。0.1ixi计算结果列于下表：计算结果列于下表： 第17页/共48页iixiy iy x iiy 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.0218520.0172130.0122620.0076260.0036420.0004220.002053第18页/共48页可用来检验近似解的准确程度。可用来检验近似解的准确程度。 进行计算，数值解已达到了一定的精度。进行计算，数值解已达到了一定的精度。这个初值问题的准确解为这个初值问题的准确解为 ， 21 1y xx从上表最后一列，我们看到取步长从上表最后一列，我们看到取步长0.1h 第19页/共48页3、 欧拉公式的改进：欧拉公式的改进： 隐式欧拉法隐式欧拉法 /* implicit Euler method */向后差商近似导数向后差商近似导数hxyxyxy)()()(011 x0 x1)(,()(1101xyxfhyxy 第20页/共48页由于未知数由于未知数 yi+1 同时出现在等式的两边，不能直接同时出现在等式的两边，不能直接得到，故称为得到，故称为隐式隐式 /* implicit */ 欧拉公式，而前者欧拉公式，而前者称为称为显式显式 /* explicit */ 欧拉公式。欧拉公式。111,0,1iiiiyhf xiyyn 第21页/共48页一般先用显式计算一个初值，再一般先用显式计算一个初值，再迭代迭代求解。求解。隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：11()iiiRy xy23( )()2ihy xO h即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度阶精度。第22页/共48页 梯形公式梯形公式 / /* *trapezoid formula trapezoid formula * */ / 显、隐式两种算法的显、隐式两种算法的平均平均) 1,., 0(),(),(2111 niyxfyxfhyyiiiiii注：注：梯形公式的局部截断误差梯形公式的局部截断误差 ，311iiiRy xyO h即梯形公式即梯形公式具有具有2 阶精度阶精度，比欧拉方法有了进步。，比欧拉方法有了进步。但注意到该公式是但注意到该公式是隐式公式隐式公式，计算时不得不用到，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。迭代法，其迭代收敛性与欧拉公式相似。第23页/共48页中点欧拉公式中点欧拉公式 /* midpoint formula */中心差商近似导数中心差商近似导数hxyxyxy2)()()(021 x0 x2x1)(,(2)()(1102xyxfhxyxy 1,., 1),(211 niyxfhyyiiii假设假设 , 则可以导出则可以导出即中点公式具有即中点公式具有 2 阶精度。阶精度。)(),(11iiiixyyxyy )()(311hOyxyRiii 第24页/共48页方方 法法 显式欧拉显式欧拉隐式欧拉隐式欧拉梯形公式梯形公式中点公式中点公式简单简单精度低精度低稳定性最好稳定性最好精度低精度低, 计算量大计算量大精度提高精度提高计算量大计算量大精度提高精度提高, 显式显式多一个初值多一个初值, 可能影响精度可能影响精度第25页/共48页 改进欧拉法改进欧拉法 /* modified Eulers method */Step 1: 先用显式欧拉公式作预测，算出先用显式欧拉公式作预测，算出Step 2: 再将再将 代入代入隐式隐式梯形公式的右边作梯形公式的右边作校正校正，得到，得到1 ny1,nnnnyyhf xy111,2nnnnnnyyyyhf xf x第26页/共48页1121211Euler,1()2(,)(,)(0,1,2,.(,)(,(,).)2nnnnnnnnnnnnnnyykkkhf xykhf xh yhyyf xyf xh yhf xykn上式还常写成该式称为改进方法 亦可写成第27页/共48页注注: :此法亦称为此法亦称为预测预测- -校正法校正法 / /* * predictor-corrector method predictor-corrector method * */ /可以证明该算法可以证明该算法具有具有 2 阶精度阶精度，同时可以看到它，同时可以看到它是个是个单步单步递推格式，比隐式公式的迭代求解过程递推格式，比隐式公式的迭代求解过程简单简单。后面将看到，它的。后面将看到，它的稳定性高于稳定性高于显式欧拉法显式欧拉法。改进的欧拉法改进的欧拉法111(,)(,)2nnnnnnhyyf xyf xy第28页/共48页在实际计算时，可将欧拉法与梯形法则相结合，在实际计算时，可将欧拉法与梯形法则相结合，计算公式为计算公式为 ,.2, 1, 0),(),(2),()(11)1(1)0(1kyxfyxfhyyyxhfyyknnnnnknnnnn应用改进欧拉法应用改进欧拉法, ,如果序列如果序列 收敛收敛, ,)1(1)0(1 nnyy它的极限便满足方程它的极限便满足方程111(,)(,)2nnnnnnhyyf xyf xy第29页/共48页改进欧拉法的截断误差改进欧拉法的截断误差)(0)(311hyxynn 因此，改进欧拉法公式具有因此，改进欧拉法公式具有 2 2 阶精度阶精度第30页/共48页例例2:2: 用改进用改进Euler公式求解例公式求解例1中的初值问题，中的初值问题， 取步长取步长 。0.1h 解：解：对此初值问题采用改进对此初值问题采用改进EulerEuler公式，公式， 其具体形式为其具体形式为 21( )211111111( ,)2(,)2()12piiiiiiicppiiiiiiipciiiyyhf x yyhx yyyhf xyyhxyyyy 计算结果列于下表：计算结果列于下表：例例1:1: 用欧拉公式求解初值问题用欧拉公式求解初值问题 2201.201yxyxy （）0,1,.,11i 01y 第31页/共48页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.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.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第32页/共48页通过计算结果的比较可以看出，改进的通过计算结果的比较可以看出，改进的Euler方法方法的计算精度比的计算精度比Euler方法要高。方法要高。第33页/共48页)()(2)O()(y2Euler式对于)O(h)(2Euler对于显式则精度相对的越高。越大，误差阶，一个算法，局部截断一般431321321hOxyhThxhTxyhTnnnnnn ：对于梯形公式：法隐：法说来第34页/共48页 对许多实际问题来说，欧拉公式与改进欧拉对许多实际问题来说，欧拉公式与改进欧拉公式精度还不能满足要求，为此从另一个角度来公式精度还不能满足要求，为此从另一个角度来分析这两个公式的特点，从而探索一条构造高精分析这两个公式的特点，从而探索一条构造高精度方法的途径度方法的途径. 第35页/共48页受改进的Euler方法启发，更一般算式可设为1121211()2(,)(,)(0,1,2,.)nnnnnnyykkkhf xykhf xh ykn123111()-()()nnnnnTy xyO hyy x适当选择参数 ， ， ，使局部截断误差,这里仍假定。改进欧拉法改进欧拉法11 122121(,)(0,1,2,.)(,)nnnnnnyykkkhf xynkhf xh yk第36页/共48页1122232()()()(,)(,)(,)()nnnxnnnnynnyy xhy xfxyf xyfxyhO h)(),(),(),()()(),(),(),(:Taylor323122hOyxfyxfyxfhxyhhOyxfhkyxfhyxhfknnynnnnxnnnynnxnn展开式由二元函数11 122121 (,) (,)(0,1,2,.)nnnnnnnyykkkhf xyhykhf xh ykn第37页/共48页由于四个参数，三个方程，因此有一个自由参数，即解答不唯一。1122232()()()(,)(,)(,)()nnnxnnnnynnyy xhy xfxyf xyfxyhO hTaylor:与展式相比较得1212122122312()()()()()hnnnny xy xhy xy xO h第38页/共48页这是改进的Euler方法。1211(1),1,22取可得此时算式为1121211()2(,) (,)nnnnnnyykkkhf xykhf xh yk第39页/共48页-R K这是二阶方法.121(2)0,1,2取可得此时算式为12121(,)11(,) 22nnnnnnyykkhf xykhf xh yk第40页/共48页R-K这也是二阶方法。12132(3),443取可得又有算式11212113)4(,)22(,)33nnnnnnyykkkhf xykhf xh yk（第41页/共48页112312123221331332112(,)(,)(,)nnnnnnnnyykkkkhf xykhf xhcccababykkhf xh ykb k 要使三阶龙格-库塔方法具有三阶精度，必须使其局部截断误差为 O(h4)将 k1, k2, k3 代入 yn+1 的表达式中，在 (xn, yn) 处用二元泰勒公式展开，与 y(xn+1) 在 xn 处的泰勒展开式比较第42页/共48页123221331322233222233232211 21 31 6cccababbc ac ac ac ac b a 8 个未知参数，个未知参数，6 个个方程，有无穷多组方程，有无穷多组解解14166611211231213122(,)(,)(,)2nnnnnnnnyykkkkhf xykhf xh ykkhf xh ykk 三阶龙格库塔公式三阶龙格库塔公式第43页/共48页1123412132335RK1(22)6(,)11(,)2211(,)22(,),()nnnnnnnnnnyykkkkkhf xykhf xh ykkhf xh ykkhf xh ykO h类似以上推导过程，可求得经典四阶算式该公式结构整齐 局部截断误差为。第44页/共48页附注：附注： 二阶二阶Runge-Kutta方法的局部截断误差方法的局部截断误差 只能达到只能达到3();O h 五阶五阶Runge-Kutta方法的局部截断误差方法的局部截断误差 只能达到只能达到 四阶四阶Runge-Kutta方法的局部截断误差方法的局部截断误差 只能达到只能达到 三阶三阶Runge-Kutta方法的局部截断误差方法的局部截断误差 只能达到只能达到 4();O h5();O h6().O h第45页/共48页附注：附注： 龙格龙格-库塔法的主要运算在于计算库塔法的主要运算在于计算 的值，即计的值，即计算算 的值。的值。Butcher 于于1965年给出了计算量与可达年给出了计算量与可达到的最高精度阶数的关系：到的最高精度阶数的关系：iKf753可达到的最高精度可达到的最高精度642每步须算每步须算Ki 的个数的个数3()O h4()O h5()O h7()O h8()O h6()O h)(2nhO8n 由于龙格由于龙格-库塔法的导出基于泰勒展开，故精库塔法的导出基于泰勒展开，故精太好的解，最好采用太好的解，最好采用低阶算法低阶算法而将步长而将步长h 取小取小。度主要受解函数的光滑性影响。对于光滑性不度主要受解函数的光滑性影响。对于光滑性不第46页/共48页感谢您的听讲！感谢您的听讲！第47页/共48页感谢您的观看！感谢您的观看！第48页/共48页