【数值分析】龙格-库塔方法与亚当姆斯方法相结合求解偏微分方程

1.实现龙格-库塔方法，并验证它的正确性；
2.先用龙格-库塔方法计算前3点，然后实现亚当姆斯预报-校正系统，并验证算法的正确性。

一、实现龙格-库塔方法，并验证它的正确性

1.4阶经典龙格-库塔方法

2.使用4阶经典龙格-库塔方法求解初值问题，取步长h=0.1,并将结果与精确解进行比较。

求精确解，并计算在h=0.1上的每个等分点的函数值

精确解为
，在各个等分点上的函数值填入表1；
用四阶经典龙格-库塔方法求解该初值问题

将得到的结果填入表1；

比较表1中精确解与龙格库塔方法得到的数值解可见龙格-库塔方法的正确性，此外还可以看出龙格-库塔方法精度较高。

二、先用龙格-库塔方法计算前3点，然后实现亚当姆斯预报-校正系统，并验证算法的正确性。

1.亚当姆斯预报-校正系统

3.用亚当姆斯预报-校正系统求解（一）.2给出的初值问题，并将结果与精确解进行比较。

三、结论

1.四阶经典龙格-库塔方法具有精度较高、易于编制计算机程序与控制步长等特点，且有较好的稳定性；但它需要计算函数值的次数较多，故只有在函数较简单或精度要求不很高时才单独使用。当函数较复杂时，经典龙格一库塔方法更适用于为多步法提供必需的初值。

2.亚当姆斯预报-校正系统具有计算函数值次数少，结果精度高的特点，并且容易进行误差估计，在进行比较复杂的计算时，效果往往较好，但需要单步法启动。

四、体会

1.龙格一库塔方法的推导基于泰勒展开方法,因而它要求所求的解具有较好的光滑性。如果解的光滑性差，那么使用四阶龙格-库塔方法求得的数值解的精度反而不如单步法中的改进的欧拉方法。在实际计算时，应当针对问题的实际特点选择合适的算法，既要考虑方法的难易、计算量大小，又要考虑方法是截断误差与收敛性、稳定性。

图 1 龙格-库塔方法的稳定区域

2.截断误差(局部截断误差、整体截断误差)可以用来确定各种公式的计算精确度，并由此构造更高阶精度的公式；数值稳定性可以确定步长的选择范围，从而保证计算方法的适应性和可靠性。

五、matlab代码

1.使用新版dsolve方法求微分方程的解析解

syms y(x);
eqn=diff(y,x)==y-2*x/y;%方程组式一
cond=y(0)==1;%方程组式二，条件
s=dsolve(eqn,cond)

2.求解析解的数值代入

function GetDsolves(x0,h)
XA=x0:h:1;
YA=[];
YA=(2.*XA + 1).^(1/2)

3.4阶经典龙格-库塔方法求数值解

function RKode(x0,y0,h)
%x0,y0为边界，h是步长。
X=x0:h:1;
Y=[];
Y(1)=y0;
for i=1:(1-x0)/h		%1-10。迭代次数可随步长做出调整
    K1(i)=h*MyDy(X(i),Y(i));
    K2(i)=h*MyDy(X(i)+h/2,Y(i)+K1(i)/2);
    K3(i)=h*MyDy(X(i)+h/2,Y(i)+K2(i)/2);
    K4(i)=h*MyDy(X(i)+h,Y(i)+K3(i));
    Y(i+1)=Y(i)+(K1(i)+2*K2(i)+2*K3(i)+K4(i))/6;
end
Y

4.亚当姆斯预报-校正系统求数值解

function Adamsode(x0,y0,h)
%x0,y0为边界，h是步长。
X=x0:h:1;
Y=[];
Y(1)=y0;%校正值
DY=[];%校正一阶导数值
for i=1:4%龙格-库塔方法求y1,y2,y3
    DY(i)=MyDy(X(i),Y(i));
    K1(i)=h*MyDy(X(i),Y(i));
    K2(i)=h*MyDy(X(i)+h/2,Y(i)+K1(i)/2);
    K3(i)=h*MyDy(X(i)+h/2,Y(i)+K2(i)/2);
    K4(i)=h*MyDy(X(i)+h,Y(i)+K3(i));
    Y(i+1)=Y(i)+(K1(i)+2*K2(i)+2*K3(i)+K4(i))/6;
end
Yyb=[];%预报值
DYyb=[];%预报一阶导数值
%Y=[];校正值
%DY=[];校正一阶导数值
for i=4:(1-x0)/h
    %预报
    Yyb(i+1)=Y(i)+h/24*(55*DY(i)-59*DY(i-1)+37*DY(i-2)-9*DY(i-3));
    DYyb(i+1)=MyDy(X(i+1),Yyb(i+1));
    %校正
    Y(i+1)=Y(i)+(9*DYyb(i+1)+19*DY(i)-5*DY(i-1)+DY(i-2))*h/24;
    DY(i+1)=MyDy(X(i+1),Y(i+1));
end
Yyb
Y