苏州大学数值分析课件-第7章常微分方程的数值解法.pptx

苏州大学数值分析课件-第7章常微分方程的数值解法.pptx

《苏州大学数值分析课件-第7章常微分方程的数值解法.pptx》由会员分享，可在线阅读，更多相关《苏州大学数值分析课件-第7章常微分方程的数值解法.pptx（76页珍藏版）》请在冰豆网上搜索。

1引言2泰勒级数法（TaylorSeriesMethod）3基于数值微分公式的方法4龙格库塔法（Runge-KuttaMethod）5线性多步法（LinearMulti-stepMethod）,第七章常微分方程的数值解法,苏州大学数值分析,第七章常微分方程的数值解法,1引言只含有一个未知数的微分或导数的方程称为常微分方程（OrdinaryDifferentialEquation），形如一个常微分方程如果有解，则有无数个解，称为通解（GeneralSolution）；其中的某一个确定的解称为特解（ParticularSolution）.考虑初值问题,上有定义且连续，同,定理7.1（初值问题解的存在唯一性）设f（x,y）在时满足李普希兹条件,其中L为李普希兹常数，则初值问题（7.1）的解存在且唯一，并且解y（x）连续可微。

定理7.2如果f（x,y）满足李普希兹条件，则初值问题是稳定的。

仅有少数常微分方程可求解析解，大部分只能求数值解，即求出在区间a,b上若干离散点处函数值的近似值本章介绍求解常微分方程（组）的常用数值方法。

2泰勒级数法（TaylorSeriesMethod）,泰勒展开式,对于初值问题,考虑其解y=y（x）在初值点x0处的p阶泰勒多项式,其中,将上述各阶导数值代入*式，即得y1的值；再以（x1,y1）为起点，建立x1点的Taylor多项式，计算出y2,如此继续，依次通过可求出数值解,例1取步长h=0.1,用二阶泰勒多项式求解初值问题,解：

由于故有于是，相应的二阶泰勒多项式的迭代式为,由初值x0=0,y0=1开始迭代，计算结果为,（此初值问题的解析解为定义7.1假定yn=y（xn）,在此前提下用某种数值方法精确计算得到yn+1,称Rn+1=y（xn+1）-yn+1为该数值方法的局部截断误差。

定义7.2若数值方法的局部截断误差为Rn+1=O（hp+1）,则称该方法具有p阶精度或p阶方法，其中p为正整数。

结论展开至p阶的泰勒级数法具有p阶精度。

3基于数值微分公式的方法对于初值问题中的微分方程，使用向前差商数值微分公式,由于故有,，得,用近似值代替整理即得欧拉（Euler）公式：

局部截断误差,对于初值问题,中的微分方程，,使用向后差商数值微分公式得后退欧拉公式：

局部截断误差对上述两种公式进行平均，得梯形公式（改进的欧拉公式）：

局部截断误差,后退欧拉（Euler）公式：

梯形公式：

欧拉（Euler）公式：

显式法,隐式法,隐式法的求解方法（以梯形公式为例）,为止，将,反复迭代，直至满足精度要求作为xn+1点的数值解。

4龙格库塔法（Runge-KuttaMethod）,在泰勒级数法当中,考虑用f（x,y）在不同点的值的线性组合来代替其中的导数项。

以二阶龙格库塔公式为例。

二阶泰勒展式为,假定上式等价于其中,参数,待定，要求方法达到二阶精度。

按照二元函数的泰勒展式在,将（xn,yn）点展开为,将其代入,中整理得,由于要求方法是二阶的，故令上式与（*）式中关于h同次幂的系数相等，得四个未知数三个方程！

选取c2为自由参量，则有根据C2取值的不同，可得到不同的迭代公式，例如取c2=1/2,得到改进的欧拉法：

取c2=1,得到中点公式：

（见课本P211-P213）,同理可得，经典四阶龙格库塔公式,例7.2用四阶经典龙格库塔法解初值问题,解：

略。

（参考课本P214）,5线性多步法（LinearMulti-stepMethod）,在常微分方程的数值解法中，如果求解yn+1时只用到前一点yn的值，则称此数值方法为单步法，如前面所述各方法；如果求解yn+k时要用到前面k（k1）个点yn+k-1,yn+k-2,yn，则称此方法为多步法或k步法。

线性多步法的建立方法利用数值积分的方法基于泰勒展开式的待定系数法,形如,的计算公式称为线性多步法，通,若则称上式为显式k步法；否则为隐式k步法。

常取,5.1.1利用数值积分的方法,考虑方程在区间,上的积分,即如果通过等距节点上的已知值作Lagrange插值（例如k=2时），可得,代入（*）式可得,上式称为显式阿当姆斯（Adams）公式，其余式为,其中的系数值bi（i=1,2,k）和误差常数B为（课本P216）,上的已知值,作Lagrange,例如，通过两节点插值，可得,代入（*）式可得,即得系数,5.1.1利用数值积分的方法,考虑方程在区间,上的积分,即如果通过等距节点上的已知值作Lagrange插值，代入（*）式可得隐式阿当姆斯公式其余式为其中系数值与误差常数见课本P217表7.4。

系数绝对值较小，因而隐式公式比显式公式的舍入误差小。

5.1.2基于泰勒展式的待定系数法,得即得欧拉公式（k=1时的阿当姆斯公式，属线性单步法）：

对于微分方程考虑泰勒展式,确定线性多步法计算公式的一般方法待定系数法考虑线性多步法的局部截断误差公式,要求该计算公式是p阶的。

将上式中y（xn+j）以及y（xn+j）在xn处展开,其中,要使该计算公式是p阶的，只要令据此求出即可。

代入上述误差公式，合并同类项得,例3建立显式线性二步法计算公式。

解：

由题意知而根据前面得到的关系式得,设,待定。

解得,从而得显式线性二步法的计算公式为,例4建立具有最高阶数的隐式线性二步法计算公式。

解：

由题意知而待定。

根据前面得到的关系式，最多可以建立5个方程,解得而从而得具有最高阶数的隐式线性二步法的计算公式为,吉尔方法（GearMethod）,形如,的线性k步法称为吉尔方法。

即普通线性k步法中取的情形。

按照前述待定系数法可求得，k阶k步吉尔方法的计算公式，其系数见课本P220表7.5,5.2线性多步法的起始值的求法,采用p阶线性k步法求解常微分方程的初值问题时，除了已知的y0之外，还必须再附加y1,y2,yk-1这k-1个起始值才能进行计算。

线性多步法的所有起始值也称为“表头”。

这些起始值的精度至少不应低于该线性k步法的局部截断误差O（hp+1）.通常的作法：

先采用单步法（如泰勒级数法）逐次获得满足精度要求的y1,y2,yk-1然后再用k步法公式进行计算。

5.3线性多步法的数值求解方法,显式公式的求解方法计算表头按照多步法公式计算例5取步长h=0.1,用四阶显式阿当姆斯公式求解初值问题,解：

四阶显式阿当姆斯公式为,整理得,要用此公式计算，除y0之外，还需知道y1,y2,y3的值，采用泰勒公式计算之。

由于相同阶数的泰勒级数法与阿当姆斯方法的局部截断误差是同阶的，故应采用不低于四阶的泰勒多项式进行计算。

四阶泰勒展开式为,其中,整理得,将x0=0,y0=1,xi=0.1h,i=1,2,3代入上式计算得,然后将上式“起始值”代入阿当姆斯公式计算得,5.3线性多步法的数值求解方法,5.3.2隐式公式的求解方法将隐式公式简记为,其中为不含及常用的隐式公式求解方法：

迭代解法预估校正法,的项。

1.迭代解法,2.预估校正法,先用显式方法,作为,的值，再用,对于隐式公式求得再按照进行迭代，把满足精度要求的,上述方法求例7见课本225例7.7,作业,1.P249第1、3、4题,2.预习chap81,SolutionusingMATLAB:

FirstOrderEquations,TheODEp（x,y）+q（x,y）y=0maybesolvedusingeitherofthefunctions:

dsolve（p（x,y）+q（x,y）*Dy=0,x）forthegeneralsolution;dsolve（p（x,y）+q（x,y）*Dy=0,y（x0）=y0,x）foraparticularsolution.,Example1,SolutionusingMATLAB:

SecondOrderEquations,TheODEy+p（x）y+q（x）y=r（x）maybesolvedusingeitherofthefunctions:

dsolve（D2y+p（x）*Dy+q（x）*y=r（x）,x）forthegeneralsolution;dsolve（D2y+p（x）*Dy+q（x）*y=r（x）,y（x0）=y0,Dy（x0）=v0,x）（initialvalue）ordsolve（D2y+p（x）*Dy+q（x）*y=r（x）,y（x0）=y0,y（x1）=y1,x）（boundaryvalue）foraparticularsolution.,Attention1:

WhentheindependentvariableofODEist,wecanomitthelastargumenttinthedsolvefunction,becausetisadefaultinMATLAB.,Example2,Attention2:

Ifnoclosed-form（explicit）solutionisfound,animplicitsolutionisattempted.Whenanimplicitsolutionisreturned,awarningisgiven.Suchas,inputY=dsolve（Dy=y2*（1-y）willproduceWarning:

Explicitsolutioncouldnotbefound;implicitsolutionreturned.Y=t+1/y-log（y）+log（-1+y）+C1=0,Attention3:

Ifneitheranexplicitnorimplicitsolutioncanbecomputed,thenawarningisgivenandtheemptysymisreturned.Suchas,inputY=dsolve（Dy=1+y2+t3）willgetWarning:

Explicitsolutioncouldnotbefound.ans=emptysymWecanonlygetthenumericalsolutioninthatcase.,1.2ODEBasicNumericalMethod,1.2.1EulersMethod：

（BF5.2，P256）GivenafirstorderODEwithinitialvaluecondition:

wetrytogetanapproximationofy（x）.Approach:

Accordingtotheequationandtheinitialvalue,trytogettheapproximatevaluesofy（x）onlyatafinitenumberofpointsornodes.Ifthenodesarethickenough,theobtaineddiscretepointscanapproximatethefunctiony（x）.,EulersMethod:

1.Dividetheintervala,bintoNpiecesequally,2.Taylorseries,3.Weget,KeyAlgorithmofEulersMethod:

Example3SolvetheODEbyEulersMethod.,Solution:

1.createan“M-file”,pileyourprograminthem-filewindow,3.produceafigureoftheapproximatecurve.,4.produceanExcelspreadsheet:

a.Doubleclickthevectorintheworkspacewindow.,b.CheckalldataintheforegoingArrayEditorwindow,copyandpastethemintothecorrespondingplaceofaExcelspreadsheet.,c.copyandpastethefigureintothecorrespondingplaceofaExcelspreadsheet.,1.2.2HigherOrderTaylorMethod（BF5.3,P266）,Toapproximatethesolutionoftheinitial-valueproblem:

at（N+1）equallyspacednotesintheintervala,b.TaylorsmethodisanimprovementonEulersmethodinfact.InEulersmethodthefirsttwotermsoftheTaylorseriesareused.Theapproximationmaybeimprovedbytakinghigherorderterms.,TaylorsMethod:

Dividetheintervala,bintoNpiecesequallyTaylorseries,3.Weget,KeyAlgorithmofTaylorsMethod:

Example4,Solution:

1.createan“M-file”,pileyourprograminthem-filewindow,3.produceafigureoftheapproximatecurve.,4.produceanExcelspreadsheet:

1.3ODEStandardNumericalMethodRunge-KuttaMethods（BF5.4，P272）,InRunge-KuttamethodsthehigherderivativesrequiredintheaboveTaylorseriesarereplacedbyaseriesofapproximategradientsintheintervalxn,xn+1.,InEulersmethodyn+1=yn+h*f（xn,yn）ienewy=oldyplushtimesLHSgradient.Abetterapproximationistotaketheaverageofgradientsateachsideoftheintervalxn,xn+1,sothatyn+1=yn+h*f（xn,yn）+f（xn+1,yn+1）andwithinthisyn+1isapproximatedbyyn+h*f（xn,yn）.Thisisnormallyre-writtenasyn+1=yn+k1+k2wherek1=h*f（xn,yn）andk2=h*f（xn+1,yn+k1）.,1.3.1SecondOrderRunge-KuttaMethod（RK2）,KeyAlgorithmofRK2:

Example5,Solution:

1.createan“M-file”,pileyourprograminthem-filewindow,3.produceafigureoftheapproximatecurve.,Example6SolvethefollowingODEusingEulers,TaylorsandRK2Methods.,Solution:

1.createan“M-file”,pileyourprograminthem-filewindow,3.produceafigureoftheapproximatecurve.,1.3.2ForthOrderRunge-KuttaMethod（RK4）,InEulersmethodyn+1=yn+h*f（xn,yn）,i.e.newy=oldyplushtimesLHSgradient.InRK2methodyn+1=yn+h*f（xn,yn）+f（xn+1,yn+1）,twoapproximationstothegradientareused:

attheLHS&RHSoftheintervalxn,xn+1.Inthefourthordermethod,4approximationstothegradientareused:

oneattheLHS,twoatthecentreandoneattheRHS.,KeyAlgorithmOfRK4Method:

Example7SolvethefollowingODEbytheForthOrderRunge-KuttaMethod.,Solution:

1.createan“M-file”,pileyourprograminthem-filewindow,3.produceafigureoftheapproximatecurve.,1.4HigherOrderODES:

（BF5.9，P313）,ManyODEshavehigherderivatives:

withy（x0）andy（x0）given.Byintroducinganewvariablez=y,thenz=ytsecondequationbecomeapairoffirstorderequationswithy（x0）andz（x0）given.Thenwecansolvetheequationsbytheprevioustheory.,Example8Solution:

rewritetheequationasapairoffirstorderequationsImplementanumericalmethodegEuler,forsimplicity,compileyourprograminthem-filewindow,produceafigureoftheapproximatecurve.,Approximatethesolutiontotheinitial-valueproblemwhereh=0.2bytheEulersmethod;theTaylorsmethodofordertwo;（3）theRunge-Kuttamethodoforderfour.,Exercise: