计算方法与实习-常微分方程数值解法.

第七章常微分方程数值解法第一节问题的提出含有一个或多个导数及其函数的方程式称为微分方程，在工程中常遇到求解微分方程的问题。很多微分方程的解不能用初等函数来表示，有时即使能够用解析式表示其解，但计算量太大而不实用(表达式过于复杂)。需要用数值方法来求解，一般只要求得到若干个点上的近似值或者解的简单的近似表达式(精度要求满足即可)。常微分方程（一个自变量）微分方程偏微分方程（一个以上自变量）一.定解问题实际中求解常微分方程的所谓定解问题有两类：初值问题和边值问题定解指已知因变量和/或其导数在某些点上是已知的(约束条件)。1.边解问题约束条件为已知，在自变量的任一非初值上，已知函数值和/或其导数值，如y”=f(x,y,y’)求解yy(a)=，y(b)=常常可以将边解问题转化为初值问题求解。2.初值问题约束条件为在自变量的初值上已知函数值，如：y’=f(x,y)dy/dx=f(x,y)xx0y(x0)=y0y(x0)=y0求解y(x)，以满足上述两式，即在ax0x1……xnb上的y(xi)的近似值yi(i=0,1,2,…,n)。通常取等距节点，即h=xi+1-xi，有xi=x0+ih(i=0,1,2,…,n)初值问题的数值解法特点：按节点顺序依次推进，由已知的y0,y1,…,yi，求出yi+1，这可以通过递推公式得到。单步法：利用前一个单步的信息(一个点)，在y=f(x)上找下一点yi，有欧拉法，龙格－库格法。初值问题的常见解法预测校正法：多步法，利用一个以上的前点信息求f(x)上的下一个yi，常用迭代法，如改进欧拉法，阿当姆斯法。第二节欧拉法y’=f(x,y)求解计算y(xi)(i=1,2,…,n)的近似值yiy(x0)=y0一、欧拉折线法1.解一阶方程初值问题的几何意义y’=f(x,y)(x,y)R，R={axb，-y+}也即解y=y(x)所表示的积分曲线y=∫f(x,y)dx+c上，每一点(x，y)的切线斜率等于函数f(x,y)在该点的值，即：y’(xk)=f(xk，yk)y(x0)=y0，表示积分曲线从P0（x0，y0）点出发且在P0（x0，y0）的切线斜率为f（x0，y0），故设想积分曲线y=y(x)在x=x0附近可用切线近似代替。2.欧拉折线法在（x0，y0）点上的切线方程为y=y0+f(x0,y0)(x-x0)设方程与x=x1交点P1（x1，y1）纵坐标y1取为y(x1)的近似值，则有y(x1)y1=y0+f(x0,y0)(x1-x0)=y0+hf(x0,y0)同理：在(x1，y1)上的切线方程y=y1+f(x1,y1)(x-x1)与x=x2交点P2(x2,y2)的纵坐标y2取为y(x2)的近似值有y(x2)y2=y1+f(x1,y1)(x2-x1)=y1+hf(x1,y1)上述的h=x2－x1=x1－x0过Pi(xi,yi)作斜率为f(xi,yi)的切线方程与x=xi+1的交点，有欧拉公式y(xi+1)yi+1=yi+hf(xi，yi)(i=0,1,…,n-1)xi=x0+ih(i=1,2,…,n)例1：求解初值问题y’=-2xy0x1.8y(0)=1取h=0.1解：f(x,y)=-2xy，x0=0，y(x0)=1，h=0.1所以y(xi+1)yi+1=yi+hf(xi,yi)=yi+0.1(-2xiyi)=(1-0.2xi)yi计算见P142表7-2-1。3.误差若y(xi)=yi，则y(xi+1)－yi+10由泰勒展开式在xi上展开有y(xi+1)=y(xi+h)=y(xi)+y‘(xi)h+(h2/2!)y”()而yi+1=yi+hf(xi,yi)=y(xi)+hf(xi,y(xi))=y(xi)+hy’(xi)所以y(xi+1)－yi+1=(h2/2)y”()=O(h2)这是局部截断误差(因为认为yi=y(xi))实际上可能yiy(xi)，从而有误差的积累，所以欧拉方法虽然简单，但精度不够，很少采用，只是说明这个方法，有助理解其他方法。二、改进欧拉法1.梯形公式y’=f(x,y)在[xi,xi+1]上的积分有即现被积函数中y(x)未知，用数值积分(矩形公式)11),('iiiixxxxdxyxfdxy1))(,()()(1iixxiidxxyxfxyxy11),())(,()()(1iiiixxixxiiidxyxfydxxyxfxyxyhxyxfdxxyxfiixxii))(,())(,(1用y(xi)yi，则y(xi+1)yi+1=yi+hf(xi,yi)若用梯形公式求积分，则所以实际上是取（xi，yi）及（xi+1，yi+1）两点的平均斜率作为在（xi，yi）点的斜率代入欧拉公式中得到上式的。现有yi+1在公式右边，如用欧拉公式求出y(xi+1)的一个粗糙近似值yi+1（预测值），代入梯形公式中得yi+1（校正值）。))(,())(,(2))(,(111iiiixxxyxfxyxfhdxxyxfii))(,())(,(2)(1111iiiiiiixyxfxyxfhyyxyyi+1（校正值）较yi+1（预测值）准确，即有改进的欧拉公式（预测-校正公式）yi+1=yi+hf(xi，yi)yi+1=yi+h/2[f(xi,y(xi)),f(xi+1,y(xi+1))]为便于编程，改为yp=yi+hf（xi，yi）yi+1=yi+h/2（k1+k2）yc=yi+hf（xi+1，yp）或k1=f（xi，yi）yi+1=1/2（yp+yc）k2=f（xi+1，yi+hk1）此外为提高精度，可迭代求yi+1，即有yi+1（0）=yi+hf（xi，yi）yi+1（k+1）=yi+h/2[f（xi，yi），f（xi+1，yi+1（k））]直到yi+1（k+1）-yi+1（k）时，取y（xi+1）yi+1（k+1）例2用改进欧拉法求例1中的初值问题解：yp=yi+hf（xi，yi）=（1－0.2xi）yiyc=yi+hf（xi+1，yp）=yi－0.2（xi+0.1）ypyi+1=1/2（yp+yc）计算见P144表7-2-22.误差y(xi+1)=y(xi+h)=y(xi)+hy’(xi)+(h2/2!)y”(xi)+……k1=f（xi，yi）=f（xi，y（xi））=y’（xi）k2=f（xi+1，yi+hk1）=f（xi+h，y（xi）+hk1）=f（xi，y（xi））+h[f（xi，y（xi））/x]+hk1[f（xi，y（xi））/y]+……=y’（xi）+h[f（xi，y（xi））/x+y’（xi）f（xi，y（xi））/y]+O(h2)=y’（xi）+hy”（xi）+O(h2)所以：yi+1=yi+h/2(k1+k2)=y(xi)+h/2[y’(xi)+y’(xi)+hy”(xi)+O(h3)]=y(xi)+hy’(xi)+h/2y”(xi)+O(h3)∴y(xi+1)－yi+1=O(h3)局部截断误差O(hp+1)时，方法具有p阶精度。因此欧拉公式有一阶精度，而改进欧拉公式有二阶精度。3.改进欧拉法N-S图及程序读入x,y读入数据h,xn输出x,yx≤xnyl=y+hf(x,y)x=x+h输出x,y结束定义函数f(x,y)=….y=y+0.5*h*[f(x,y)+f(x+h,yl)]第三节龙格－库塔方法一、基本思想根据微分中值定理有：y’(xi+h)=[y(xi+1)－y’(xi)]/h(01，h=xi+1-xi)即f(xi+h,y(xi+h)=[y(xi+1)－yi]/h是平均斜率∴y(xi+1)=yi+hf(xi+h,y(xi+h))=yi+hk*记k*=f(xi+h,y(xi+h))而计算k*的方法是欧拉公式中k*=f(xi,yi)精度低改进欧拉公式中k*=0.5*[f(xi,yi)+f(xi+1,yi+1)]=0.5*[k1+f(xi+1,yi+hk1)]=0.5*[k1+k2]多取了一个点使精度得以提高。如果设法在[xi，xi+1]上多预报几个点的斜率值，然后将它们的加权平均值作为k*=f(xi+h,y(xi+h))的近似值，则有可能构造出更高精度的计算公式，这就是龙格－库塔方法的基本思想。二、二阶龙格－库塔公式在[xi，xi+1]内有xi+l=xi+lh（0l1）取k*=1k1+2k2(是xi,xi+1两个点的斜率值k1,k2加权平均值）则yi+1=yi+hk*=yi+h（1k1+2k2）用欧拉公式，取k1=f（xi，yi）而k2=f（xi+l，yi+l）=f（xi+l，yi+lhk1）yi+1=yi+hf（1k1+2k2）∴k1=f（xi，yi）k2=f（xi+l，yi+lhk1）现计算1，2及l。若要使上述公式具有二阶精度，即y(xi+1)－yi+1=O(h3)设y(xi)=yi，展开y(xi+1)=y’(xi)+hy’(xi)+h2/2y”(xi)+O(h3)∵k1=f(xi,yi)=y’(xi)k2=f(xi+l,yi+lhk1)=f(xi,yi)+lh[fx(xi,yi)+fy(xi,yi)+f(xi,yi)]+O(h2)=y’(xi)+lhy”(xi)+O(h2)∴yi+1=y(xi)+h(1+2)y’(xi)+h22ly”(xi)+O(h3)具有二阶精度，只需1+2=1三个未知量，两个方程2l=1/2（1）当l=1时，即xi+l=xi+1，解得1=2=1/2yi+1=yi+h/2（k1+k2）∴k1=f（xi，yi）k2=f（xi+l，yi+hk1）即为改进欧拉公式（2）取l=1/2，则1=0，2=1，有变形的欧拉公式yi+1=yi+hk2∴k1=f（xi，yi）k2=f（xi+h/2，yi+h/2k1）三、高阶龙格－库塔公式在[xi，xi+1]上除xi，xi+1外再增加一点xi+m=xi+mh（lm1）设xi+m处的斜率为k3，则：k*=1k1+2k2+3k3yi+1=yi+hk*=yi+h(1k1+2k2+3k3)k1=f（xi，yi）k2=f（xi+lh，yi+lhk1）k3=f（xi+m，yi+m）=f（xi+mh，yi+mh（μ1k1+μ2k2））得yi+1=yi+h（1k1+2k2+3k3）k1=f（xi，yi）k2=f（xi+lh，yi+lhk1）k3=f（xi+mh，yi+mh（μ1k1+μ2k2））利用泰勒展开方法，选取1,2,3,l,m及μ1,μ2使上述公式具有三阶精度O(h4)。得μ1+μ2=11+2+3=1有7个未知数，5个方程，2l+3m=1/2可得一族解，所得公式2l2+3m2=1/3为三阶龙格－库塔公式3lmμ2=1/61.库塔公式取l=1/2,m=1,则1=1/6,2=2/3,3=1/6,μ1=-1,μ2=2得库塔公式yi+1=yi+h/6（k1+4k2+k3）k1=f（xi，yi）k2=f（xi+h/2，yi+h/2k1）k3=f（xi+h，yi-hk1+2hk2）2.四阶龙格－库塔公式在[xi，xi+1]上用四个点的ki(i=1,2,3,4)做k*的加权平均值，可得一族经典四阶龙格－库塔公式。yi+1=yi+h/6（k1+2k2+2k3+k4）k1=f（xi，yi）k2=f（xi+h/2，yi+h/2k1）k3=f（xi+h/2，yi+h/2k2）k4=f（xi+h，yi+hk3）例3：用四阶龙格－库塔方法求解y’=－2xy0x1.8y(0)=1取h=0.2局部截断误差O(h5)，精度提高。解：yi+1=yi+0.2/6（k1+2k2+2k3+k4）k1=-2xiyik2=f（xi+0