一、实验一： 常微分方程形态和R-K法稳定性试验
考察下面微分方程右端项中函数y前面的参数a对方程性态的影响（它使得方程为好条件的或坏条件的）和研究计算步长h对4阶经典R-K法计算稳定性的影响。
其中，-50≤a≤50，其精确解为y ( x ) =
+ x
1、步骤
（1）对参数a分别取4个不同的数值：一个大的正值，一个小的正值，一个绝对值小的负值，一个绝对值大的负值。取步长h=0.01，分别用经典R-K法计算，将四组计算结果画在同一张图上，进行比较并说明相应初值问题的形态。
（2）取参数a为一个绝对值不大的负值和两个计算步长，一个步长使参数ah在经典R-K法的稳定域内，另一个步长在经典R-K法的稳定域外，分别用经典R-K法计算并比较计算结果。取全域等距的10个点上的计算值，列表说明。
2、程序
function test_5_1 %输入参数a和步长h %利用经典4阶R-K法计算，并做图 result=inputdlg({'请输入[-50,50]之间的整数'},'text_5_1',1,{'-50'}); a=str2num(char(result)); if (a<-50)|| (a>50) errordlg({'超出范围，请重新输入参数a'}); return; end result=inputdlg({'请输入(0,1)之间的步长h'},'text_5_1',1,{'0.01'}); h=str2num(char(result)); if (h<=0)|| (h>=1) errordlg({'超出范围，请重新输入参数a'}); return; end x=0:h:1; y=1; N=length(x); fun=@(x,y)a*y-a*x+1; for i=1:N-1 k1=fun(x(i),y(i)); k2=fun(x(i)+h/2,y(i)+h*k1/2); k3=fun(x(i)+h/2,y(i)+h*k2/2); k4=fun(x(i)+h,y(i)+h*k3); y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4); end plot(x,y,'r-')%将plot改为semilogy，y就是对数坐标轴 hold on y0=exp(a*x)+x; plot(x,y0,'b*') xlabel('x'); ylabel('y'); legend('R-K近似值','精确值') r=abs((y(101)-y0(101))/y0(101));%误差分析 sprintf('当h是%e，a是%d时的绝对误差是%e',h,a,r)
3、结果
第一问：
由于步长选取为h=0.01，而4阶经典R-K法的绝对稳定域是-2.785≤λh≤0。因此，当h=0.01时，-278.5≤a≤0，该方法才稳定。又要求-50≤a≤50，因此a的稳定区域为-5
围内时，结果收敛，为好条件。反之则为换条件。
下面为分别选取a=-45，-1，1，45下的图形。
从上图上并不能直接看出产生多少误差（不明显）。因此程序中加入一段绝对误差的计算。我们假设一开始的初值没有误差，最后迭代到101次产生了的绝对误差如下：
a=-45,r=0
a=-1,r=2.2600e-11
a=1,r=6.0416e-11
a=45,r=1.0535e-02
可见在稳定区域内或靠近稳定区域内，误差很小。但是远离稳定区域会带来较大的误差。
第二问
取a=-10，则h的稳定区域为0<h≤0.2785。取步长h分别为0.1和0.3，作图如下
由表中等距的10个点的计算值和图像可以看出，h值超过稳定区域时结果发散，与精确值相差很大。而h在稳定区域内时，计算结果和精确值几乎一样。
二、实验二：比较同阶精度的单步法和预测-校正多步法的计算精度和所花费计算机CPU时间
对于常微分方程：
1、比较4阶经典R-K算法和4阶Adams预测-校正修正法的计算精度，以及所花费的时间。
要求选择一个步长h使得4阶经典R-K算法和4阶Adams预测-校正修正法均稳定，分别用这两种方法计算初值问题，列出10个等距节点上的计算值和精确值并比较。其中多步法所需要的初值由R-K法提供。
2、再取h=0.001，比较这两种方法的时间。
1、思路
4阶Adams预测-校正修正法
2、程序
function test_5_2 result=inputdlg({'请输入步长h'},'test_5_2',1,{'0.1'}); h=str2num(char(result)); fun=@(x,y)-y+cos(2*x)-2*sin(2*x)+2*x.*exp(-x); %经典R-K法 tic; [x,y1]=R_K(fun,h,0,2,1); time_RK=toc; %4阶Adams预测-校正法 tic; N=length(x); y=x;%x与y同维,将x赋值给y，随后运算中逐个替换为相应的元素。这样比重新构造一个向量再逐次插入元素快。 f=x; y(1)=y1(1);y(2)=y1(2);y(3)=y1(3);y(4)=y1(4);%以R-K法获得的初值 ycn=0;ypn=0; f(1)=fun(x(1),y(1));f(2)=fun(x(2),y(2));f(3)=fun(x(3),y(3)); for n=4:N-1 f(n)=fun(x(n),y(n)); ypn1=y(n)+h/24*(55*f(n)-59*f(n-1)+37*f(n-2)-9*f(n-3)); ypmn1=ypn1+251/270*(ycn-ypn); ycn1=y(n)+h/24*(9*fun(x(n+1),ypmn1)+19*f(n)-5*f(n-1)+f(n-2)); y(n+1)=ycn1-19/270*(ycn1-ypn1); ycn=ycn1; ypn=ypn1; end time_Ad=toc; y0=x.^2.*exp(-x)+cos(2*x); subplot(1,2,1); plot(x,y1,'b--') hold on plot(x,y0,'g*') xlabel('x');ylabel('y');legend('R-K法','精确值'); subplot(1,2,2); plot(x,y,'b--') hold on plot(x,y0,'g*') xlabel('x');ylabel('y');legend('Adams预测-校正法','精确值'); d=floor(N/10); disp(['十个等距节点： ',num2str(x(1:d:N))]); disp(['对应的R-K法计算值： ',num2str(y1(1:d:N))]); disp(['对应的Adams预测-校正法计算值：',num2str(y(1:d:N))]); disp(['精确值： ',num2str(y0(1:d:N))]); disp(['对应的R-K法计算时间： ',num2str(time_RK)]); disp(['对应的Adams预测-校正法计算时间：',num2str(time_Ad)]); function [x,y]=R_K(fun,h,a,b,y0) x=a:h:b; y=y0; N=length(x); for i=1:N-1%R-K法得到初值y(n) k1=fun(x(i),y(i)); k2=fun(x(i)+h/2,y(i)+h*k1/2); k3=fun(x(i)+h/2,y(i)+h*k2/2); k4=fun(x(i)+h,y(i)+h*k3); y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4); end
3、结果
当h=0.1时：
从理论推到上：
从实验结果上看，理论可得到验证。而仔细观察数据，可得出经典R-K算法比4阶Adams预测-校正法精度稍微高一点。
从计算时间上看，则是4阶Adams预测-校正法会比经典R-K法快。
当h=0.001
对应的R-K法计算时间： 0.0068151
对应的Adams预测-校正法计算时间：0.0048398
