MATLAB运用(4)-常微分方程的性态和初值问题

一、实验一： 常微分方程形态和R-K法稳定性试验

考察下面微分方程右端项中函数y前面的参数a对方程性态的影响（它使得方程为好条件的或坏条件的）和研究计算步长h对4阶经典R-K法计算稳定性的影响。
$$\{\begin{array}{l}{y}^{\prime }=ay-ax+1,0<x<1\\ y(0)=1\end{array}$$
其中，-50≤a≤50，其精确解为$y(x)=e^{ax}+x$

1、步骤

（1）对参数a分别取4个不同的数值：一个大的正值，一个小的正值，一个绝对值小的负值，一个绝对值大的负值。取步长h=0.01，分别用经典R-K法计算，将四组计算结果画在同一张图上，进行比较并说明相应初值问题的形态。
（2）取参数a为一个绝对值不大的负值和两个计算步长，一个步长使参数ah在经典R-K法的稳定域内，另一个步长在经典R-K法的稳定域外，分别用经典R-K法计算并比较计算结果。取全域等距的10个点上的计算值，列表说明。
（备注：经典R-K法的绝对稳定域区间：$-2.785\le \lambda h\le 0,\lambda =f_y^{\prime },f=y^{\prime }$）

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的稳定区域为-50≤a≤0。当a处于稳定性范围内时，结果收敛，为好条件。反之则为换条件。
下面为分别选取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时间

对于常微分方程：
$$\{\begin{array}{l}{y}^{\prime }=-y+cos2x-2sin2x+2x{e}^{-x},0<x<2\\ y(0)=1\end{array}$$
有精确解$y(x)=x^2{e}^{-x}=cos2x$
1、比较4阶经典R-K算法和4阶Adams预测-校正修正法的计算精度，以及所花费的时间。
要求选择一个步长h使得4阶经典R-K算法和4阶Adams预测-校正修正法均稳定，分别用这两种方法计算初值问题，列出10个等距节点上的计算值和精确值并比较。其中多步法所需要的初值由R-K法提供。
2、再取h=0.001，比较这两种方法的时间。

1、思路

4阶Adams预测-校正修正法
经典R-K法的绝对稳定域区间：$-2.785\le \lambda h\le 0$
4阶Adams预测-校正法的绝对稳定域区间：$-0.75\le \lambda h\le 0$
因此本实验两算法的绝对稳定域区间：$0\le h\le -0.75$

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时：

从理论推到上：
显式Adams预测-校正法的局部截断误差：$T(n)=\frac{251}{720}{h}^5{y}^{(5)}(\xi )$，因此有4阶精度
R-K算法的局部截断误差：$T(n+1)=h^3[\frac{1}{6}{y}^{\prime \prime }(n)+(\frac{1}{6}-\frac{{\lambda }^2{c}_2^2}{2})(f_{xx}+2f_{xy}{y}^{\prime }+f_{yy}(y^{\prime }{)}^2)]+O(h^4)$，也具有4阶精度。
从实验结果上看，理论可得到验证。而仔细观察数据，可得出经典R-K算法比4阶Adams预测-校正法精度稍微高一点。
从计算时间上看，则是4阶Adams预测-校正法会比经典R-K法快。

当h=0.001

对应的R-K法计算时间： 0.0068151
对应的Adams预测-校正法计算时间：0.0048398