9 常微分方程数值解

/*NumericalMethodsforOrdinaryDifferentialEquations*/,考虑一阶常微分方程的初值问题/*Initial-ValueProblem*/:,只要f(x,y)在a,bR1上连续，且关于y满足Lipschitz条件，即存在与x,y无关的常数L使对任意定义在a,b上的y1(x)和y2(x)都成立，则上述IVP存在唯一解。,要计算出解函数y(x)在一系列节点a=x0x1xn=b处的近似值,节点间距为步长，通常采用等距节点，即取hi=h(常数)。,第9章常微分方程数值解,1引言,2欧拉方法,/*EulersMethod*/,欧拉公式：,欧拉法的局部截断误差：,欧拉法具有1阶精度。,Ri的主项/*leadingterm*/,亦称为欧拉折线法/*Eulerspolygonalarcmethod*/,例1：,解：h=0.2,xi=1+ih，精确解为：y=x2-2xy1=y0+hf(x0,y0)=-1y2=y1+hf(x1,y1)=-0.9333y3=y2+hf(x2,y2)=-0.8y4=y3+hf(x3,y3)=-0.6y5=y4+hf(x4,y4)=-0.3333y6=y5+hf(x5,y5)=0,可以看出误差随着计算在积累。,步长h=0.2,1EulersMethod,隐式欧拉法/*implicitEulermethod*/,由于未知数yi+1同时出现在等式的两边，不能直接得到，故称为隐式/*implicit*/欧拉公式，而前者称为显式/*explicit*/欧拉公式。,一般先用显式计算一个初值，再迭代求解。,隐式欧拉法的局部截断误差：,即隐式欧拉公式具有1阶精度。,欧拉公式的改进：,1EulersMethod,/*trapezoidformula*/,显/隐式两种算法的平均,注：的确有局部截断误差，即梯形公式具有2阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式，计算时不得不用到迭代法。,中点欧拉公式/*midpointformula*/,假设，则可以导出即中点公式具有2阶精度。,需要2个初值y0和y1来启动递推过程，这样的算法称为双步法/*double-stepmethod*/，而前面的三种算法都是单步法/*single-stepmethod*/。,梯形公式,1EulersMethod,简单,精度低,稳定性最好,精度低,计算量大,精度提高,计算量大,精度提高,显式,多一个初值,可能影响精度,Cantyougivemeaformulawithalltheadvantagesyetwithoutanyofthedisadvantages?,Doyouthinkitpossible?,Well,callmegreedy,OK,letsmakeitpossible.,/*modifiedEulersmethod*/,注：此法亦称为预测-校正法/*predictor-correctormethod*/。可以证明该算法具有2阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。后面将看到，它的稳定性高于显式欧拉法。,1EulersMethod,改进欧拉法,例2：,解：梯形公式为：,可以预测校正方法来求：,解yk+1出来比较困难，遇到的是一个二次方程，,即：,/*Runge-KuttaMethod*/,建立高精度的单步递推格式。,单步递推法的基本思想是从(xi,yi)点出发，以某一斜率沿直线达到(xi+1,yi+1)点。欧拉法及其各种变形所能达到的最高精度为2阶。,斜率一定取K1K2的平均值吗？,步长一定是一个h吗？,3龙格-库塔法,2Runge-KuttaMethod,首先希望能确定系数1、2、p，使得到的算法格式有2阶精度，即在的前提假设下，使得,Step1:将K2在(xi,yi)点作Taylor展开,Step2:将K2代入第1式，得到,2Runge-KuttaMethod,Step3:将yi+1与y(xi+1)在xi点的泰勒展开作比较,要求，则必须有：,这里有个未知数，个方程。,3,2,存在无穷多个解。所有满足上式的格式统称为2阶龙格-库塔格式。,注意到，就是改进的欧拉法。,Q:为获得更高的精度，应该如何进一步推广？,其中i(i=1,m)，i(i=2,m)和ij(i=2,m;j=1,i1)均为待定系数，确定这些系数的步骤与前面相似。,2Runge-KuttaMethod,最常用为四级4阶经典龙格-库塔法/*ClassicalRunge-KuttaMethod*/：,例3：,其精确解为：,用经典龙格库塔方法，h=0.1,(1)求y1，此时,解：,以下计算用表格列出：,附：三阶龙格库塔方法:,共有八个参数：1,2,3,1,2,21,31,32，,1、三阶龙格库塔公式格式如下：,将f在xk点泰勒展开：,将K1,K2按展开式代入K3得：(高于h2的不计),将y(xk+1)在xk点泰勒展开：,其中:,欲使