常微分方程数值解法

第九章常微分方程的数值解法§1、引言§2、初值问题的数值解法--单步法§3、龙格-库塔方法§4、收敛性与稳定性§5、初值问题的数值解法―多步法§6、方程组和刚性方程§7、习题和总结主要内容主要内容研究的问题数值解法的意义§1、引言现实世界中大多数事物内部联系非常复杂找出其状态和状态变化规律之间的相互联系，也即一个或一些函数与他们的导数之间的关系此种关系的数学表达就为微分方程1.什么是微分方程?其状态随着时间、地点、条件的不同而不同如何利用数值方法求解微分方程（组）的问题。2.数值求解微分方程的意义如何建立数学模型已在建模课程中得到讨论，各类微分方程本身和他们的解所具有的特性已在常微分方程及数学物理方程中得以解释，本章专门讨论寻找解析解的过程称为求解微分方程组。一个或一组具有所要求阶连续导数的解析函数，将它代入微分方程（组），恰使其所有条件都得到满足的解称为解析解(或古典解),称为真解或解。3.什么是微分方程(组)的解析解?3.什么是微分方程(组)的解析解?4.什么是微分方程的数值解?虽然求解微分方程有许多解析方法,但解析方法只能够求解一些特殊类型的方程,从实际意义上来讲我们更关心的是某些特定的自变量在某一个定义范围内的一系列离散点上的近似值.寻找数值解的过程称为数值求解微分方程。把这样一组近似解称为微分方程在该范围内的数值解在大量的实际方程中出现的函数起码的连续性都无法保证，更何况要求阶的导数求解数值解很多微分方程根本求不到问题的解析解！重要手段。常微分方程的数值解法常用来求近似解根据提供的算法通过计算机便捷地实现5.常微分方程数值解法的特点数值解法得到的近似解(含误差)是一个离散的函数表.00(,)(1.1)()(1.2)yfxyaxbyxy本章主要讨论一阶常微方程的初值问题6.基本知识其中f(x,y)是已知函数，(1.2)是定解条件也称为初值条件。各种数值解法则称f(x,y)对y满足李普希兹条件，L称为Lipschitz常数.常微分方程的理论指出：当f(x,y)定义在区域G=(a≤x≤b,｜y｜＜∞)若存在正的常数L使：1212|(,)(,)|||(1.3)fxyfxyLyy12[,],xabyy使得对任意的及都成立就可保证方程解的存在唯一性(Lipschitz)条件我们以下的讨论，都在满足上述条件下进行.若f(x,y)在区域G连续，关于y一阶常微分方程的初值问题的解存在且唯一.满足李普希兹条件一阶常微分方程组常表述为：1111(,,,)()(,,,)mmmmyfxyyaxbyfxyy1010()()mmyxyx方程组初值条件写成向量形式为(1)()00000(,)(),(,,)mTyfxyaxbyxyxxx高阶常微分方程定解问题如二阶定解问题：(,,)()()yfxyyaxbyaya也就解决了高阶方程的定解问题.这些解法都可以写成向量形式用于一阶常微分方程组的初值问题.简单的数值方法与基本概念1.简单欧拉法（Euler）2．后退的欧拉法3．梯形法4．改进Euler法§2、初值问题的数值解法―单步法1.简单的欧拉(Euler)方法考虑模型：00(,)(1.1)()(1.2)yfxyaxbyxy在精度要求不高时通过欧拉方法的讨论弄清常微方程初值问题数值解法的一些基本概念和构造方法的思路.欧拉方法最简单而直观实用方法2.欧拉方法的导出把区间［a,b］分为n个小区间步长为要计算出解函数y(x)在一系列节点()iiyyx)/()(iiixaihhhban,一般取即等距节点处的近似值01naxxxb1(-)iiihxxN等分00(,)(1.1)()(1.2)yfxyaxbyxy对微分方程(1.1)两端从1nnxx到进行积分11(,())nnnnxxxxydxfxyxdx11()()(,())nnxnnxyxyxfxyxdx右端积分用左矩形数值求积公式2()()()()()2baggxdxbagaba()(,())gxfxyx令1)1(,())((,())nnnnxxnnfxyxnnyyfxyxh得x0x11()()()(,)nnnnnnyxyxhyxyhfxy)1,...,0(),(1niyxfhyyiiii亦称为欧拉折线法/*Euler’spolygonalarcmethod*/每步计算只用到1ny或用向前差商近似导数1()()()nnnyxyxyxh100021111(,)(,)(,)nnnnyyhfxyyyhfxyyyhfxy依上述公式逐次计算可得：ny例题3.欧拉公式有明显的几何意义依此类推得到一折线故也称Euler为单步法。ny公式右端只含有已知项所以又称为显格式的单步法。000011122(,)()(,)()(,)(,)xyyxxyyxxxxyxy过点的曲线是解在作的切线与直线交于再作切线交于也称欧拉折线法.就是用这条折线近似地代替曲线欧拉方法()yx()yyxxynx1nxnp1np1npx从上述几何意义上得知，由Euler法所得的折线明显偏离了积分曲线，可见此方法非常粗糙。4.欧拉法的局部截断误差：在假设第i步计算是精确的前提下，考虑截断误差称为局部截断误差/*localtruncationerror*/。1()iiiRyxy定义若某算法的局部截断误差为O(hp+1)，则称该算法有p阶精度。Ri的主项/*leadingterm*/欧拉法的局部截断误差：232()()hiyxOh欧拉法具有1阶精度。)],([)]()()()([)(32112iiiihiiiiiyxhfyhOxyxyhxyyxyR定义在假设yi=y(xi)，即第i步计算是精确的前提下，考虑的截断误差Ri=y(xi+1)yi+1称为局部截断误差/*localtruncationerror*/。如果单步差分公式的局部截断误差为O(hp+1),则称该公式为p阶方法.这里p为非负整数.显然,阶数越高,方法的精度越高.Taylor展开式,一元函数的Taylor展开式为:321!3)(!2)()()()()(hxyhxyhxyxyhxyxynnnnnn23()[()()()()][(,)]112hRyxyyxhyxyxOhyhfxyiiiiiiiii232()()hiyxOh若某算法的局部截断误差为O(hp+1)，则称该算法有p阶精度。Ri的主项/*leadingterm*/5.欧拉公式的改进：隐式欧拉法/*implicitEulermethod*/向后差商近似导数hxyxyxy)()()(011x0x1))(,()(1101xyxfhyxy)1,...,0(),(111niyxfhyyiiii由于未知数yi+1同时出现在等式的两边，不能直接得到，故称为隐式/*implicit*/欧拉公式，而前者称为显式/*explicit*/欧拉公式。一般先用显式计算一个初值，再迭代求解。隐式欧拉法的局部截断误差：11)(iiiyxyR)()(322hOxyih即隐式欧拉公式具有1阶精度。6.梯形公式/*trapezoidformula*/—显、隐式两种算法的平均)1,...,0()],(),([2111niyxfyxfhyyiiiiii注：的确有局部截断误差，即梯形公式具有2阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。)()(311hOyxyRiii中点欧拉公式/*midpointformula*/中心差商近似导数hxyxyxy2)()()(021x0x2x1))(,(2)()(1102xyxfhxyxy1,...,1),(211niyxfhyyiiii假设，则可以导出即中点公式具有2阶精度。)(),(11iiiixyyxyy)()(311hOyxyRiii需要2个初值y0和y1来启动递推过程，这样的算法称为双步法/*double-stepmethod*/，而前面的三种算法都是单步法/*single-stepmethod*/。方法显式欧拉隐式欧拉梯形公式中点公式简单精度低稳定性最好精度低,计算量大精度提高计算量大精度提高,显式多一个初值,可能影响精度改进欧拉法/*modifiedEuler’smethod*/Step1:先用显式欧拉公式作预测，算出),(1iiiiyxfhyyStep2:再将代入隐式梯形公式的右边作校正，得到1iy)],(),([2111iiiiiiyxfyxfhyy注：此法亦称为预测-校正法/*predictor-correctormethod*/。可以证明该算法具有2阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。后面将看到，它的稳定性高于显式欧拉法。)1,...,0(),(,),(211niyxfhyxfyxfhyyiiiiiiii§3龙格-库塔法/*Runge-KuttaMethod*/考察改进的欧拉法，可以将其改写为：),(),(2121121211hKyhxfKyxfKKKhyyiiiiii斜率一定取K1K2的平均值吗？步长一定是一个h吗？单步递推法的基本思想是从(xi,yi)点出发，以某一斜率沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形所能达到的最高精度为2阶。建立高精度的单步递推格式。首先希望能确定系数1、2、p，使得到的算法格式有2阶精度，即在的前提假设下，使得)(iixyy)()(311hOyxyRiiiStep1:将K2在(xi,yi)点作Taylor展开)(),(),(),(),(2112hOyxfphKyxphfyxfphKyphxfKiiyiixiiii)()()(2hOxyphxyii将改进欧拉法推广为：),(),(][12122111phKyphxfKyxfKKKhyyiiiiii),(),(),(),(),(),()(yxfyxfyxfdxdyyxfyxfyxfdxdxyyxyxStep2:将K2代入第1式，得到)()()()()]()()([)(322212211hOxyphxyhyhOxyphxyxyhyyiiiiiiiiStep3:将yi+1与y(xi+1)在xi点的泰勒展开作比较)()()()(322211hOxyphxyhyyiiii)()(2)()()(321hOxyhxyhxyxyiiii要求，则必须有：)()(311hOyxyRiii21,1221p这里有个未知数，个方程。32存在无穷多个解。所有满足上式的格式统称为2阶龙格-库塔格式。21,121p注意到，就是改进的欧拉法。Q:为获得更高的精度，应该如何进一步推广？其中i(i=1,…,m)，i(i=2,…,m)和ij(i=2,…,m;j=1,…,i1)均为待定系数，确定这些系数的步骤与前面相似。)...,(......),(),(),(]...[1122112321313312122122111mmmmmmimiiiiiimmiihKhKhKyhxfKhKhKyhxfKhKyhxfKyxfKKKKhyy最常用为四级4阶经典龙格-库塔法/*ClassicalRunge-KuttaMethod*/：