产品
Adams方法，可实现1至6阶的外插/内插：
1主函数，Adams外插法或内插法下的逐点误差，及误差对数函数图像：
clc;clear;close all;
%% y'=y-2*x/y用Adams来求解验证
%要求匿名函数支持向量化
func=@(x,y) y-2*x./y;
y_exa=@(x) sqrt(2.*x+1);
x0=0;
h=1/10;
xmax=1;
y0=1;
x=x0:h:xmax; %自变量
if x(end)~=xmax
x=[x,x(end)+h];
end
method='in'; %out 外插法 ；in 内插法
if strcmp(method,'out')
stra='外插法';
elseif strcmp(method,'in')
stra='内插法';
end
accuracy=5;
%% h=1/10下逐点误差
%参数：外插/内插，导函数，因变量初值，自变量左端点，步长，自变量右端点，精度
y=Adams_auto(method,func,y0,x0,h,xmax,accuracy);
y_rel=y_exa(x);
err=abs(y_rel-y);
figure(3);
plot(x,y,'r-',x,y_rel,'b:');
title(['Adams',stra,'下的数值解图像与解析解图像']);
legend('数值解','解析解');
disp('逐点误差：');
disp(err);
%% 求误差与步长的函数
h=1/10*2.^(linspace(-3,-7,5)); %步长过大过小精度都会降低
hl=length(h);
err=zeros(1,length(h));
for j=1:hl
x1=x0:h(j):xmax;
y1=Adams_auto(method,func,y0,x0,h(j),xmax,accuracy);
y_rel_1=y_exa(x1);
err(j)=mean(abs(y_rel_1-y1));
end
lerr=log(err);
lh=log(h);
%求整体误差的阶，等于对数下的误差与步长函数的斜率
Deh=(lerr(2:end)-lerr(1:end-1))./(lh(2:end)-lh(1:end-1));
deh=mean(Deh);
%画对数下误差函数与步长的函数图像
figure(4);
plot(log(h),log(err));
title([num2str(round(deh)),'阶精度Adams',stra,'下的误差对数函数图像']);
disp(['此方法具有',num2str(round(deh)),'阶精度']);
2子函数1，实现具有1到6阶精度的外插及内插Adams方法：
%求任意1到6阶的数值解；
%首先利用外置函数RungeKutta_auto来求启动所需要的点，再用Adams外插法或者内插法来求剩余的点
%预备工作：根据所选择方法求出所需启动点
%外插法流程：（1）求系数a (2)求系数b (3)求剩余点
%内插法流程: (1)求系数a (2)求系数b （3）求系数a* (4)求系数b* (5)求剩余点P(EC)...
%注意事项：要求导函数支持向量化
%完成时间：2022/3/20 嫁佛幡
function [ u ] = Adams_auto( method,func,y0,x0,h,xmax,accuracy)
%[待求点]=[外插/内插，导函数，因变量初值，自变量左端点，步长，自变量右端点，精度]
ay=accuracy; %精度
tol=1e-7; %可容误差限度
x=x0:h:xmax; %自变量
if x(end)~=xmax
x=[x,x(end)+h];
end
N=length(x); %自变量长度
u=zeros(1,N); %因变量初始化
u(1)=y0;
%外插法
if strcmp(method,'out')
%单步法求启动点，ay阶精度需要ay个点
u(1:ay)=RungeKutta_auto(ay,func,y0,x0,h,x(ay));
%求系数a,b
[a,b]=Adams_arg('out',ay);
%求剩余点
for j=ay:N-1
u(j+1)=u(j)+h*func(x(j:-1:j-ay+1),u(j:-1:j-ay+1))*b ;
end
elseif strcmp(method,'in') && ay>=2
%单步法求启动点，ay阶精度需要ay-1个点，
u(1:ay-1)=RungeKutta_auto(ay,func,y0,x0,h,x(ay-1));
%求系数a,b
[a,b]=Adams_arg('out',ay-1);
%求系数a1,b1
[a1,b1]=Adams_arg('in',ay);
%求剩余点
for j=ay-1:N-1
%预估
u_pre=u(j)+h*func(x(j:-1:j-ay+2),u(j:-1:j-ay+2))*b ;
%校准tol=1e-7
delta=1;
u_part=h*func(x(j:-1:j-ay+2),u(j:-1:j-ay+2))*b1(2:end); %避免重复运算
while delta>tol
u_new=u(j)+h*func(x(j+1),u_pre)*b1(1)+u_part;
delta=abs(u_new-u_pre);
u_pre=u_new;
end
u(j+1)=u_new;
end
else
disp('enter out or in')
end
end
3子函数2，单步法求具有1到6阶精度的启动点：
function [ u ] = RungeKutta_auto(accuracy,func,y0,x0,h,xmax )
%[待求点函数值]=[精度，导函数，因变量初值，自变量左端点，步长，因变量右端点]
%参数初始化
ay=accuracy;
x=x0:h:xmax; %自变量
if x(end)~=xmax
x=[x,x(end)+h];
end
N=length(x); %自变量长度
u=zeros(1,N); %因变量初始化
u(1)=y0;
acmax=7; %最高精度需要的维度
a=zeros(1,acmax);
b=zeros(acmax,acmax-1);
c=zeros(1,acmax);
if ay==1
%Euler法
c(1)=1;
elseif ay==2
%中点法
a(2)=1/2;
b(2,1)=1/2;
c(1:2)=[0,1];
elseif ay==3
%Heun三阶方法
a(2:3)=[1/3,2/3];
b(2,1)=1/3;
b(3,2)=2/3;
c(1:3)=[1/4,0,3/4];
elseif ay==4
%经典四阶龙格库塔
a(2:4)=[1/2,1/2,1];
b(2,1)=1/2;
b(3,2)=1/2;
b(4,3)=1;
c(1:4)=[1,2,2,1]/6;
elseif ay==5
%Runge_Kutta_Fehlberg
a(2:6)=[1/4,3/8,12/13,1,1/2];
b(2,1)=1/4;
b(3,1:2)=[3/32,9/32];
b(4,1:3)=[1932/2197,-7200/2197,7296/2197];
b(5,1:4)=[439/216,-8,3680/513,-845/4104];
b(6,1:5)=[-8/27,2,-3544/2565,1859/4104,-11/40];
c(1:6)=[16/135,0,6656/12825,28561/56430,-9/50,2/55];
% c(1:5)=[25/216,0,1408/2565,2197/4101,-1/5]; %
elseif ay==6
%by H.A.Luther from "An explicit sixth-order runge-kutta formula "
a(2:7)=[1,1/2,2/3,(7-sqrt(21))/14,(7+sqrt(21))/14,1];
b(2,1)=1;
b(3,1:2)=[3/8,1/8];
b(4,1:3)=[8,2,8]/27;
b(5,1:4)=[3*(3*sqrt(21)-7),-8*(7-sqrt(21)),48*(7-sqrt(21)),-3*(21-sqrt(21))]/392;
b(6,1:5)=[-5*(231+51*sqrt(21)),-40*(7+sqrt(21)),-320*sqrt(21),3*(21+121*sqrt(21)),392*(6+sqrt(21))]/1960;
b(7,1:6)=[15*(22+7*sqrt(21)),120,40*(7*sqrt(21)-5),-63*(3*sqrt(21)-2),-14*(49+9*sqrt(21)),70*(7-sqrt(21))]/180;
c(1:7)=[9,0,64,0,49,49,9]/180;
end
for j=1:N-1
k1=func(x(j),u(j)) ;
k2=func(x(j)+h*a(2),u(j)+h*k1*b(2,1));
k3=func(x(j)+h*a(3),u(j)+h*[k1,k2]*b(3,1:2)');
k4=func(x(j)+h*a(4),u(j)+h*[k1,k2,k3]*b(4,1:3)');
k5=func(x(j)+h*a(5),u(j)+h*[k1,k2,k3,k4]*b(5,1:4)');
k6=func(x(j)+h*a(6),u(j)+h*[k1,k2,k3,k4,k5]*b(6,1:5)');
k7=func(x(j)+h*a(7),u(j)+h*[k1,k2,k3,k4,k5,k6]*b(7,1:6)');
u(j+1)=u(j)+h*[k1,k2,k3,k4,k5,k6,k7]*c';
end
end
4子函数3，求Adams外插法或内插法的系数a,b：
function [A,B] = Adams_arg( method,n )
%[Adams系数a,b]=[方法，精度n]
%输出两个列向量
%系数对比法 PA=e
P=zeros(n);
A=zeros(n,1);
B=zeros(n,1);
for j=1:(n)
P(j,j:end)=1./(1:(n+1-j));
end
if strcmp(method,'out')
e=ones(n,1);
elseif strcmp(method,'in')
e=zeros(n,1);
e(end)=1;
else
disp('enter out or in')
return
end
A=flipud(P\e); %{a1,...,an}
%Bn
for L=0:n-1
for j=L:n-1
B(L+1)=B(L+1)+A(j+1)*factorial(j)/(factorial(L)*factorial(j-L));
end
B(L+1)=(-1)^L*B(L+1);
end
免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删
155-2731-8020
