南京工业大学-计算机在材料中的应用常微分方程数值解法

7常微分方程数值解法,7.1问题的提出含有一个或多个导数及其函数的方程式称为微分方程，在工程中常遇到求解微分方程的问题。很多微分方程的解不能用初等函数来表示，有时即使能够用解析式表示其解，但计算量太大而不实用(表达式过于复杂)。需要用数值方法来求解，一般只要求得到若干个点上的近似值或者解的简单的近似表达式(精度要求满足即可)。常微分方程（一个自变量）微分方程偏微分方程（一个以上自变量）,7.1.1定解问题实际中求解常微分方程的所谓定解问题有两类：初值问题和边值问题定解指已知因变量和/或其导数在某些点上是已知的(约束条件)。(1)边值问题约束条件为已知，在自变量的任一非初值上，已知函数值和/或其导数值，如y”=f(x,y,y)求解yy(a)=，y(b)=常常可以将边值问题转化为初值问题求解。,(2)初值问题约束条件为在自变量的初值上已知函数值，如：y=f(x,y)dy/dx=f(x,y)xx0y(x0)=y0y(x0)=y0求解y(x)，以满足上述两式，即在ax0、x1、xnb上的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，常用迭代法如改进欧拉法、阿当姆斯法。,7.2欧拉方法y=f(x,y)求解计算y(xi)(i=1,2,n)的近似值yiy(x0)=y07.2.1欧拉公式(1)解一阶方程初值问题的几何意义y=f(x,y)(x,y)R，R=axb，-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附近可用切线近似代替。P163图7-2-1(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=x2x1=x1x0,过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)P164例7.1：求解初值问题y=-2xy0x1.8y(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计算见P164表7-2-1。,(3)误差若y(xi)=yi，则y(xi+1)yi+10由泰勒展开式在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”(i)=O(h2)这是局部截断误差(因为认为yi=y(xi)实际上可能yiy(xi)，从而有误差的积累，所以欧拉方法虽然简单，但精度不够，很少采用，只是说明这个方法，有助理解其他方法。,7.2.2梯形公式(改进欧拉法)(1)梯形公式y=f(x,y)在xi,xi+1上的积分有即现被积函数中y(x)未知，用数值积分(矩形公式),用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(校正值)。,yi+1(校正值）较yi+1(预测值）准确，即有改进欧拉公式(预测-校正公式）yi+1=yi+hf(xi,yi)yi+1=yi+h/2f(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/2f(xi,yi)，f(xi+1,yi+1(k)直到yi+1(k+1)-yi+1(k)时，取y(xi+1)yi+1(k+1),P166例7.2用改进欧拉法求例7.1中的初值问题解：yp=yi+hf(xi,yi)=(10.2xi)yiyc=yi+hf(xi+1,yp)=yi0.2(xi+0.1)ypyi+1=1/2(yp+yc)计算见P167表7-2-2(2)误差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)+hf(xi,y(xi)/x+hk1f(xi,y(xi)/y+,=y(xi)+hf(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/2y(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,y,xxn,yt=y+hf(x,y),x=x+h,输出x,y,结束,定义函数f(x,y)=.,y=y+0.5*h*f(x,y)+f(x+h,yt),7.3龙格库塔方法7.3.1基本思想根据微分中值定理有：y(xi+h)=y(xi+1)y(xi)/h(01,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)的近似值，则有可能构造出更高精度的计算公式，这就是龙格库塔方法的基本思想。7.3.2二阶龙格库塔公式在xi,xi+1内有xi+w=xi+wh（0w1)取k*=1k1+2k2(是xi,xi+1两个点的斜率值k1,k2加权平均值)则yi+1=yi+hk*=yi+h(1k1+2k2)用欧拉公式，取k1=f(xi,yi)而k2=f(xi+w,yi+w)=f(xi+w,yi+whk1)yi+1=yi+hf(1k1+2k2)k1=f(xi,yi)k2=f(xi+w,yi+whk1),现计算1,2及w。若要使上述公式具有二阶精度，即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+w,yi+whk1)=f(xi,yi)+whfx(xi,yi)+fy(xi,yi)+f(xi,yi)+O(h2)=y(xi)+why”(xi)+O(h2)yi+1=y(xi)+h(1+2)y(xi)+h22wy”(xi)+O(h3),具有二阶精度，只需1+2=1三个未知量，两个方程w2=1/2(1)当w=1时，即xi+w=xi+1，解得1=2=1/2yi+1=yi+h/2(k1+k2)k1=f(xi,yi)k2=f(xi+l,yi+hk1)即为改进欧拉公式(2)取w=1/2，则1=0，2=1，有变形的欧拉公式yi+1=yi+hk2k1=f(xi,yi)k2=f(xi+1/2,yi+h/2k1),7.3.3高阶龙格库塔公式在xi,xi+1上除xi,xi+w外再增加一点xi+m=xi+mh(wm1)设xi+m处的斜率为k3，则：k*=1k1+2k2+3k3yi+1=yi+hk*=yi+h(1k1+2k2+3k3)k1=f(xi,yi)k2=f(xi+wh,yi+whk1)k3=f(xi+mh,yi+m)=f(xi+mh,yi+mh(1k1+2k2)得yi+1=yi+h(1k1+2k2+3k3)k1=f(xi,yi)k2=f(xi+wh,yi+whk1)k3=f(xi+mh,yi+mh(1k1+2k2)利用泰勒展开方法，选取1,2,3,w,m及1,2使上述公式具有三阶精度O(h4)。,得1+2=11+2+3=1有7个未知数，5个方程，2w+3m=1/2可得一族解，所得公式2w2+3m2=1/3为三阶龙格库塔公式3wm2=1/6(1)库塔公式取w=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)P173例7.3用四阶龙格库塔方法求解例7.1y=2xy0x1.8y(0)=1取h=0.2,局部截断误差O(h5)，精度提高。,解：yi+1=yi+0.2/6(k1+2k2+2k3+k4)k1=-2xiyik2=f(xi+0.1,yi+0.1k1)=-2(xi+0.1)(yi+0.1k1)k3=f(xi+0.1,yi+0.1k2)=-2(xi+0.1)(yi+0.1k2)k4=f(xi+0.2,yi+0.2k3)=-2(xi+0.2)(yi+0.2k3)计算结果见P174表7-3-1，在表7-3-2中对比了几种方法的误差。表7-3-3说明再提高阶数没有多大意义。四阶是兼顾了精度和计算量的理想公式。,(3)四阶龙格-库塔公式的几何意义,xi,xi+h,A,B,E,C,D,R,S,Pi(xi，yi),hk1,斜率为f(xi，yi)的PA，AM=hki，PA中点R(xi+h/2,yi+h/2ki)的斜率作PB,BM=hk2，有PB中点S(xi+h/2,yi+h/2k2)的斜率作PC，CM=hk3，C点(xi+h,yi+hk3)的斜率作PD，DM=hk4，取EM=hk*=h/6*(k1+2k2+2k3+k4)y(xi+h)=yi+h=yi+1/6*(k1+2k2+2k3+k4),M,hk2,hk3,hk4,(4)自动定步长四阶龙格库塔设从xi出发，先以h为步长，用四阶龙格库塔公式求得y(xi+1)的近似值yi+1(h)，则y(xi+1)yi+1(h)ch5再以h/2为步长，两次计算后(先算yi+h/2，再算yi+1)可得yi+1y(xi+1)yi+1(h/2)c(h/2)5=c/16*h5有y(xi+1)yi+1(h/2)yi+1(h/2)yi+1(h)/15即利用事后误差估计法|yi+1(h/2)yi+1(h)|，判断是否停止计算。,(5)四阶龙格库塔法N-S图及程序,读入初值x,y,读入数据h,读入xn,xxn,k1,k2,k3,k4,y=y+h/6(k1+2k2+2k3+k4),x=x+h,结束,定义函数f(x,y),输出k1,k2,k3,k4,x,y,7.4线形多步法,单步法中，求解yi+1时，实际已求出了y0,y1yi。利用多个y值，期望得到较高精度的yi+1y(xi+1)=y(xi)+xi+1f(x，y(x)dxxi用更精确的求积方法，可用更高次的插值多次来代f(x,y)，选取的插值节点不同，则解法不同。7.4.1阿当姆斯内插公式除取xi,xi+1两个节点外，其他节点取在xi,xi+1时，f的值未知，故选取xi,xi+1以外的节点为xi-1,xi-2等。现取xi-1,xi-2,xi,xi+1四个插值节点，由于计算yi+1时，yi-2,yi-1,yi,已求解出来，插值多项式,L3(x)=(x-xi)(x-xi-1)(x-xi-2)/(xi+1-xi)(xi+1-xi-1)(xi+1-xi-2)f(xi+1,y(xi+1)+(x-xi+1)(x-xi-1)(x-xi-2)/(xi-xi+1i)(xi-xi-1)(xi-xi-2)f(xi,y(xi)+(x-xi+1)(x-xi)(x-xi-2)/(xi-1-xi+1)(xi-1-xi)(xi-1-xi-2)f(xi-1,y(xi-1)+(x-xi+1)(x-xi)(x-xi-1)/(xi-2-xi)(xi-2-xi-1)(xi-2-xi+1)f(xi-2,y(xi-2)f(x,y(x)L3(x)所以yi+1y(xi)+xi+1L3(x)dxxi等距节点时，xi+1-xi=xi-xi-1=xi+1-xi-2=h，设yi=f(xi,y(xi)令x=xi+th0t1,yi+1=yi+h/6f(xi+1,y(xi+1)1t(t-1)(t+2)dt0-h/2f(xi,y(xi)1(t-1)(t+1)(t+2)dt0+h/2f(xi-1,y(xi-1)1(t-1)t(t+2)dt0-h/6f(xi-2,y(xi-2)1(t-1)t(t+1)dt0=y(xi)+h/249yi+1+19yi-5yi-1+yi-2阿当姆斯内插公式是四阶公式y(xi-1)-yi+1=(-19/720)h5yi(5)+=0(h5),7.4.2阿当姆斯外推公式,阿当姆斯内插公式中由于选取了xi+1点为节点，从而使公式成为隐式，若取xi-3，xi-2，xi-1，xi作为插值节点，即：iiL3(x)=(x-xj)/(xk-xj)f(xj-y(xj)k=i-3j=I-3jk所以y(xj+1)=yi+1=xi+1L3(x)dxxi=yi+h/2455yi+59yi-1+37yi-2-9yi-3阿当姆斯外推公式y(xi+1)-yi+1=(25/720)h5yi(5)+=0(h5),7.4.3阿当姆斯预测校正式计算f(x,y)的次数，算yi+1时yi-1，yi-2，yi-3已在计算yi时得到，只要算yi就可，无需算四次，但需用其它方法先行计算y1，y2，y3后才能用该法。依次计算y4，y5，所以不具有自启动性，此外h要改变也很麻烦。若将外推公式作预测式，内插公式作校正，可得阿当姆斯预测校正式:yi+1(0)=yi+h/2455yi+59yi-1+37yi-2-9yi-3yi+1=yi+h/249f(xi+1，yi+1(0)+19yi+5yi-1-yi-2,P179例7.4用阿当姆斯法求解例7.1y=-2xy0x1.2y(0)=1取h=0.1解：四阶龙格-库塔求出，y1，y2，y3，由y0，y1，y2，y3阿当姆斯预测校正公式。计算见P179表7-4-1，表7-4-2预测校正值比外推法精度高的多。,7.5一阶方程组与高阶方程,7.5.1一阶方程组u=(x,u,v),u(x0)=u0v=(x,u,v),v(x0)=v0初值解记y=(u,v)T,y0=(u0,v0)TF=(,)Ty=F(x,y)y(x0)=y0用龙格库塔解,yi+1=yi+h/6(k1+2k2+2