matlab应用二步和四步显示的Adamas方法求解下列微分方程的初值问题（用四阶Runge-Kutta方法提供出发值）应用二步和四步显示的Adamas方法求解下列微分方程的初值问题（用四阶Runge-Kutta方法提供

matlab应用二步和四步显示的Adamas方法求解下列微分方程的初值问题（用四阶Runge-Kutta方法提供出发值）应用二步和四步显示的Adamas方法求解下列微分方程的初值问题（用四阶Runge-Kutta方法提供
matlab应用二步和四步显示的Adamas方法求解下列微分方程的初值问题（用四阶Runge-Kutta方法提供出发值）
应用二步和四步显示的Adamas方法求解下列微分方程的初值问题（用四阶Runge-Kutta方法提供出发值）
y'=0.1(x^3+y^3) (0

matlab应用二步和四步显示的Adamas方法求解下列微分方程的初值问题（用四阶Runge-Kutta方法提供出发值）应用二步和四步显示的Adamas方法求解下列微分方程的初值问题（用四阶Runge-Kutta方法提供
给你RK和Adamas一些参考程序：
例 9.2 用经典四阶 方法计算:
解:

%ex9_2.m
%四阶经典R-K公式作数值计算
clc;
F='y-2*x/y';
a=0;
b=1;
h=0.1;
n=(b-a)/h;
X=a:h:b;
Y=zeros(1,n+1);
Y(1)=1;
for i=1:n
x=X(i);
y=Y(i);
K1=h*eval(F);
x=x+h/2;
y=y+K1/2;
K2=h*eval(F);
x=x;
y=Y(i)+K2/2;
K3=h*eval(F);
x=X(i)+h;
y=Y(i)+K3;
K4=h*eval(F);
Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6;
end

%准确解
temp=[];
f=dsolve('Dy=y-2*x/y','y(0)=1','x');
df=zeros(1,n+1);
for i=1:n+1
temp=subs(f,'x',X(i));
df(i)=double(vpa(temp));
end
disp(' 步长 四阶经典R-K法 准确值');
disp([X',Y',df']);
%画图观察效果
figure;
plot(X,df,'k*',X,Y,'--r');
grid on;
title('四阶经典R-K法解常微分方程');
legend('准确值','四阶经典R-K法');

运行上述程序得到如下结果:

步长 四阶经典R-K法 准确值
0 1.000000000000000 1.000000000000000
0.100000000000000 1.095445531693094 1.095445115010332
0.200000000000000 1.183216745505993 1.183215956619923
0.300000000000000 1.264912228340392 1.264911064067352
0.400000000000000 1.341642353750372 1.341640786499874
0.500000000000000 1.414215577890085 1.414213562373095
0.600000000000000 1.483242222771993 1.483239697419133
0.700000000000000 1.549196452302143 1.549193338482967
0.800000000000000 1.612455349658987 1.612451549659710
0.900000000000000 1.673324659016256 1.673320053068151
1.000000000000000 1.732056365165566 1.732050807568877
作出的函数图形如下:
这个结果比上述两种方法精度高得多.

例 9.3 分别就和确定线性多步法的系数,使方法具有最高的截断误差阶.
解:
直接按公式求解相应系数即可,本题未知系数较多,且数目不确定,不适合编制自动求解程序.
先对情况讨论,按线性多步法的步骤得到如下方程组:
>> [a0,b0,b1]=solve('a0=1','b0+b1=1','2*b1=1')
a0 =
1
b0 =
1/2
b1 =
1/2

再由误差公式的表达式得到误差阶数:
>>C3=1/factorial(3)*(a0)-1/factorial(3-1)*(b1)
C3 =
-1/12

对情形类似讨论,本文不再给出解答.

例 9.4 利用四步四阶显式法计算:
解:
前三步用经典四阶R-K法启动计算.编制如下求解程序:

%ex9_4.m
%Adams四步四阶显式法作常微分方程数值计算
%[a,b]为求解区间,h为步长
clc;
F='y-2*x/y';
a=0;
b=1;
h=0.1;
n=(b-a)/h;
X=a:h:b;
Y=zeros(1,n+1);
Y(1)=1;
%以四阶R-K法启动
for i=1:3
x=X(i);
y=Y(i);
K1=h*eval(F);
x=x+h/2;
y=y+K1/2;
K2=h*eval(F);
x=x;
y=Y(i)+K2/2;
K3=h*eval(F);
x=X(i)+h;
y=Y(i)+K3;
K4=h*eval(F);
Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6;
end
%Adams四步四阶显式法
for i=4:n
x=X(i-3);
y=Y(i-3);
f1=eval(F);
x=X(i-2);
y=Y(i-2);
f2=eval(F);
x=X(i-1);
y=Y(i-1);
f3=eval(F);
x=X(i);
y=Y(i);
f4=eval(F);
Y(i+1)=Y(i)+h*(55*f4-59*f3+37*f2-9*f1)/24;
end

%准确解
temp=[];
f=dsolve('Dy=y-2*x/y','y(0)=1','x');
df=zeros(1,n+1);
for i=1:n+1
temp=subs(f,'x',X(i));
df(i)=double(vpa(temp));
end
disp(' 步长 Adams四步四阶显式法 准确值');
disp([X',Y',df']);
%画图观察效果
figure;
plot(X,df,'k*',X,Y,'--r');
grid on;
title('Adams四步四阶显式法解常微分方程');
legend('准确值','Adams四步四阶显式法');

程序运行结果:

步长 Adams四步四阶显式法 准确值
0 1.000000000000000 1.000000000000000
0.100000000000000 1.095445531693094 1.095445115010332
0.200000000000000 1.183216745505993 1.183215956619923
0.300000000000000 1.264912228340392 1.264911064067352
0.400000000000000 1.341551759049205 1.341640786499874
0.500000000000000 1.414046421479413 1.414213562373095
0.600000000000000 1.483018909732277 1.483239697419133
0.700000000000000 1.548918873971137 1.549193338482967
0.800000000000000 1.612116428793334 1.612451549659710
0.900000000000000 1.672917033446480 1.673320053068151
1.000000000000000 1.731569752635566 1.732050807568877

它与准确值的计算结果对比图如下:


由图上可看到线性多步法的精度还是很高的.它的优点在于每次计算量大大减小,缺点是不能自启动,需借助其他方法启动.

例 9.5 利用预测—校正公式法求解:
解:
前面三步还是用经典四阶R-K方法启动计算,求解程序如下:

%ex9_5.m
%Adams校正-预测法作常微分方程数值计算
%[a,b]为求解区间,h为步长
clc;
F='y-2*x/y';
a=0;
b=1;
h=0.1;
n=(b-a)/h;
X=a:h:b;
Y=zeros(1,n+1);%Adams预测值
Y(1)=1;
%以四阶R-K法启动
for i=1:3
x=X(i);
y=Y(i);
K1=h*eval(F);
x=x+h/2;
y=y+K1/2;
K2=h*eval(F);
x=x;
y=Y(i)+K2/2;
K3=h*eval(F);
x=X(i)+h;
y=Y(i)+K3;
K4=h*eval(F);
Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6;
end
%Adams校正-预测法
Y1=Y;%Adams校正值
for i=4:n
x=X(i-3);
y=Y(i-3);
f1=eval(F);
x=X(i-2);
y=Y(i-2);
f2=eval(F);
x=X(i-1);
y=Y(i-1);
f3=eval(F);
x=X(i);
y=Y(i);
f4=eval(F);
Y(i+1)=Y(i)+h*(55*f4-59*f3+37*f2-9*f1)/24;%Adams预测值
x=X(i+1);
y=Y(i+1);
f0=eval(F);
Y1(i+1)=Y(i)+h*(9*f0+19*f4-5*f3+f2)/24;%校正值
end

%准确解
temp=[];
f=dsolve('Dy=y-2*x/y','y(0)=1','x');
df=zeros(1,n+1);
for i=1:n+1
temp=subs(f,'x',X(i));
df(i)=double(vpa(temp));
end
disp(' 步长 Adams预测值 Adams校正值 准确值');
disp([X',Y',Y1',df']);
%画图观察效果
figure;
plot(X,df,'k*',X,Y,'-.r',X,Y1,'--b');
grid on;
title('Adams校正-预测法解常微分方程');
legend('准确值','Adams预测值','Adams校正值');

运行上述程序得到如下结果:

步长 Adams预测值 Adams校正值 准确值
0 1.000000000000000 1.000000000000000 1.000000000000000
0.100000000000000 1.095445531693094 1.095445531693094 1.095445115010332
0.200000000000000 1.183216745505993 1.183216745505993 1.183215956619923
0.300000000000000 1.264912228340392 1.264912228340392 1.264911064067352
0.400000000000000 1.341551759049205 1.341641357193254 1.341640786499874
0.500000000000000 1.414046421479413 1.414107280831795 1.414213562373095
0.600000000000000 1.483018909732277 1.483044033257615 1.483239697419133
0.700000000000000 1.548918873971137 1.548934845800237 1.549193338482967
0.800000000000000 1.612116428793334 1.612129054676922 1.612451549659710
0.900000000000000 1.672917033446480 1.672925295781879 1.673320053068151
1.000000000000000 1.731569752635566 1.731575065330948 1.732050807568877

它们的曲线如下:

上述曲线观察得不是很清楚,我们进行局部放大后得到如下效果图:


相对预测值,校正值的曲线更接近精确值的曲线.