常微分方程数值解法 5

主要内容1.Euler方法2.改进Euler法3.龙格库塔(Runge-Kutta)法4、一阶微分方程组与高阶微分方程的数值解法5、相容性，收敛性和稳定性6、线性多步法章节练习,第九章常微分方程数值解法,目的：用数值解逼近解析解,由常微分方程理论知，初值问题（9-1）有解且唯一。,处得近似值的方法，称为问题（9-1）的数值解。为步长。,常用几种方法:用差商近似导数；用数值积分方法；用Taylor多项式近似。,数值解法-求问题（9-1）的解y(x)在若干点,现将微分方程离散化-建立数值解法,(1)用差商近似导数,可由初值y0逐次算出.称离散化问题（9-3）为差分方程初值问题。,（2）用数值积分方法,（3）用Taylor多项式近似,函数f(x)在xn处展开，取一次Taylor多项式近似，得,Euler方法-用差分方程初值问题的解近似微分方程初值问题（9-1）。,可由初值y0逐次算出.,用Euler方法求下面初值问题的数值解,9-1、Euler方法,例1、,解：,若取h=0.1,N=10,可由初值y0=0逐次算出.见下表,注：,（1）Euler方法-Euler折线法,（2）,1-2Euler方法的误差估计,（1）截断误差,（2）截断整体误差（不含舍入误差）,结合前面有,若用梯形公式代替微分方程（9-4）右端积分，有,梯形公式的局部截断误差,2-1改进的Euler方法,2-1-1梯形公式,由于梯形公式仍然是隐式方法，一般用迭代法求解，若每一步只代一次，相当将Euler公式与梯形公式结合使用。,为上机编程方便，改写为,算法9-1（见教材209）,2-2改进Euler法(梯形法),改进Euler法求解下面初值问题的数值解。（取h=0.1,N=10）,用公式（9-7）,若取h=0.1,N=10,可由初值y0=0逐次算出.见下表,例2、,解：,改进Euler法优于Euler法.(二阶方法),思考与练习,近似值与准确值比较，Euler法有2位有效数字，改进Euler法有3位有效数字,3.龙格库塔(Runge-Kutta)法,由Euler公式,改进Euler公式,3-2(Runge-Kutta)法的构造,一般有,确定相应参数的原则：近似公式在（Xn,Yn)处的泰勒展式与y(x)在Xn处的泰勒展式的前面的项尽可能重合，以提高精度。,当P=2时的模型,上式在（Xn,Yn)处展开,比较系数有：,改进Euler公式,（二阶）中点公式,（3）三阶(Runge-Kutta)公式,（4）四阶经典(Runge-Kutta)公式,算法9-2（见教材212）,用四阶经典(Runge-Kutta)法求解下面初值问题的数值解。（取h=0.2,N=5）,依题意四阶经典(Runge-Kutta)公式应为,解：,例3、,由初值y0=0逐次算出下表,特点：用四阶经典(Runge-Kutta)法计算的结果，其精度高于改进Euler法.,例1.,解:,利用经典的四阶龙格-库塔公式有,计算结果如下表。,该结果具有四位有效数字，和本章第一节例1比，可见四阶龙格-库塔方法的精度高于Euler公式和改进Euler（预估-校正）公式。,%a=0;y0=1;h=0.2;n=0;x_rec(1)=a;y_rec(1)=y0;fprintf(x%2.0f=%10.6fy%2.0f=%10.6fn,n,x_rec(1),n,y_rec(1)x=a;y=y0;whilex1n=n+1k1=y-2*x/y;k2=(y+h*k1/2)-(2*(x+h/2)/(y+h*k1/2);k3=(y+h*k2/2)-(2*(x+h/2)/(y+h*k2/2);k4=(y+h*k3)-(2*(x+h)/(y+h*k3);y=y+h*(k1+2*k2+2*k3+k4)/6;x=a+n*h;y_rec(n+1)=y;x_rec(n+1)=x;fprintf(x%2.0f=%10.6fy%2.0f=%10.6fn,n,x_rec(n+1),n,y_rec(n+1)endplot(x_rec,y_rec),x0=0.000000y0=1.000000n=1x1=0.200000y1=1.183229n=2x2=0.400000y2=1.341667n=3x3=0.600000y3=1.483281n=4x4=0.800000y4=1.612514n=5x5=1.000000y5=1.732142,%a=0;y0=2;h=0.25;n=0;x_rec(1)=a;y_rec(1)=y0;fprintf(x%2.0f=%10.6fy%2.0f=%10.6fn,n,x_rec(1),n,y_rec(1)x=a;y=y0;whilex5n=n+1k1=-x*y2;k2=-(x+h/2)*(y+h*k1/2)2;k3=-(x+h/2)*(y+h*k2/2)2;k4=-(x+h)*(y+h*k3)2;y=y+h*(k1+2*k2+2*k3+k4)/6;x=a+n*h;y_rec(n+1)=y;x_rec(n+1)=x;fprintf(x%2.0f=%10.6y%2.0f=%10.6fn,n,x_rec(n+1),n,y_rec(n+1)end,例1.,plot(x_rec,y_rec)x0=0.000000y0=2.000000n=1x1=0.250000y1=1.882308n=2x2=0.500000y2=1.599896n=3x3=0.750000y3=1.279948n=4x4=1.000000y4=1.000027n=5x5=1.250000y5=0.780556n=6x6=1.500000y6=0.615459n=7x7=1.750000y7=0.492374n=8x8=2.000000y8=0.400054n=9x9=2.250000y9=0.329940n=10 x10=2.500000y10=0.275895,n=11x11=2.750000y11=0.233602n=12x12=3.000000y12=0.200020n=13x13=3.250000y13=0.172989n=14x14=3.500000y14=0.150956n=15x15=3.750000y15=0.132790n=16x16=4.000000y16=0.117655n=17x17=4.250000y17=0.104924n=18x18=4.500000y18=0.094123n=19x19=4.750000y19=0.084885n=20 x20=5.000000y20=0.076927,4.一阶微分方程组与高阶微分方程的数值解法,4-1一阶微分方程组的数值解法,设一阶微分方程组的初值问题,向量形式,(1)改进Euler公式,特别有：,(2)四阶经典(Runge-Kutta)公式,一般地，用于单步法公式,当m=2时，有,四阶经典(Runge-Kutta)公式,4-2高阶微分方程的数值解法,具体方法：通过变量代换把其高阶微分方程化为一阶微分方程组的初值问题。（降阶法）,设m阶常微分方程初值问题,方法：用对一阶微分方程组的初值问题求解即可。,用四阶经典(Runge-Kutta)方法求解（取h=0.1）,解：,例5、,用公式（*）、（*）有:,从而,逐次算出下表,原问题的精确解为它们之间的误差见上表的最后一列。,注：刚性问题：若对一阶微分方程组,则称（*）为刚性方程组或Stiff方程组.,特别：对有刚性的一阶微分方程组，用上述方法求解会遇到根本的困难，这是由数值方法自身的稳定性限制所决定的。（由h确定）,5-1.问题:离散化是否合理？,5、相容性，收敛性和稳定性,2结论,5-2.稳定性,绝对稳定区域-,绝对稳定区域包含左半平面-该方法称为A稳定。,（1）Euler方法稳定性,由Euler方法知，计算公式为,则Euler方法的绝对稳定区域为,特点：,-1,-2,（2）梯形公式的稳定性,由梯形公式知,类似于前面梯形公式的方法,绝对稳定区域为左半平面。,（3）RK方法的稳定性,仿前面梯形公式的方法，将RK方法用于讨论的模型，可得二、三和四阶RK方法的绝对稳定区域。,特别有：,-1,-2,i,2i,3i,-i,-2i,-3i,N=1,N=2,N=3,N=4,由四阶RK方法的计算公式,解：,例5、,由初值y0=1逐次算出下表,由表中可看出，当h=0.1,计算结果的相对误差很大，但计算过程是稳定的。当h=0.1,计算过程是不稳定的。,事实上，,在考虑稳定性对步长的限制时，常用,6、线性多步法,为提高精确度和计算速度，由单步法一般公式,可推广跟一般的线性多步法公式,方法：通过相应的条件，代入公式比较系数，可求得参数。,（2）阿达姆斯外插公式;,介绍三个问题,（1）线性多步方法的基本思想；,（3）阿达姆斯内插公式.,（1）线性多步方法的基本思想,基本思想。,（2）阿达姆斯(Adams)外插公式及误差,上式公式称为线性四步阿达姆斯显式公式，也称为阿达姆斯外插公式。,其局部截断误差就是数值积分的误差：,注意：阿达姆斯外插公式要进行计算，必须提供初值。在实际计算中，常用四阶龙格-库塔方法计算出这些初值，然后再用阿达姆斯外插公式计算。,例1,解先由四阶龙格-库塔公式求出初值，结果为,然后上式计算其他节点处的值，结果为,（3）阿达姆斯内插公式,公式上式(*1)称为线性三步阿达姆斯公式，或称为阿达姆斯内插公式，是四阶方法。,公式(*1)是隐式公式，不便于直接