微分方程数值解笔记之1到6阶的Adams外插内插实现

微分⽅程数值解笔记之1到6阶的Adams外插内插实现5.Adams⽅法，可实现1⾄6阶的外插/内插：5.1主函数，Adams外插法或内插法下的逐点误差，及误差对数函数图像：clc;clear;closeall;%%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;ifx(end)~=xmaxx=[x,x(end)+h];end%⾃变量method='in';%out外插法；in内插法ifstrcmp(method,'out')stra='外插法';elseifstrcmp(method,'in')stra='内插法';endaccuracy=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));forj=1:hlx1=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));endlerr=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)),'阶精度']);5.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;ifx(end)~=xmax

x=[x,x(end)+h];endN=length(x);u=zeros(1,N);u(1)=y0;%⾃变量长度%因变量初始化%外插法ifstrcmp(method,'out')%单步法求启动点，ay阶精度需要ay个点u(1:ay)=RungeKutta_auto(ay,func,y0,x0,h,x(ay));%求系数a,b[a,b]=Adams_arg('out',ay);%求剩余点forj=ay:N-1u(j+1)=u(j)+h*func(x(j:-1:j-ay+1),u(j:-1:j-ay+1))*b;endelseifstrcmp(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);%求剩余点forj=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-7delta=1;u_part=h*func(x(j:-1:j-ay+2),u(j:-1:j-ay+2))*b1(2:end);%避免重复运算whiledelta>tolu_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;endelsedisp('enteroutorin')endend5.3⼦函数2，单步法求具有1到6阶精度的启动点：function[u]=RungeKutta_auto(accuracy,func,y0,x0,h,xmax)%[待求点函数值]=[精度，导函数，因变量初值，⾃变量左端点，步长，因变量右端点]%参数初始化ay=accuracy;x=x0:h:xmax;ifx(end)~=xmaxx=[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);ifay==1%Euler法c(1)=1;elseifay==2%中点法a(2)=1/2;b(2,1)=1/2;c(1:2)=[0,1];elseifay==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];elseifay==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;elseifay==5%Runge_Kutta_Fehlberga(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];%elseifay==6%byH.A.Lutherfrom"Anexplicitsixth-orderrunge-kuttaformula"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;endforj=1:N-1k1=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';endend5.4⼦函数3，求Adams外插法或内插法的系数a,b：function[A,B]=Adams_arg(method,n)%[Adams系数a,b]=[⽅法，精度n]%输出两个列向量%系数对⽐法PA=eP=zeros(n);A=zeros(n,1);B=zeros(n,1);forj=1:(n)P(j,j:end)=1./(1:(n+1-j));endifstrcmp(method,'out')e=ones(n,1);elseif