微分方程數(shù)值解筆記之1到6階的Adams外插內(nèi)插實(shí)現(xiàn)

微分?程數(shù)值解筆記之1到6階的Adams外插內(nèi)插實(shí)現(xiàn)5.Adams?法，可實(shí)現(xiàn)1?6階的外插/內(nèi)插：5.1主函數(shù)，Adams外插法或內(nèi)插法下的逐點(diǎn)誤差，及誤差對數(shù)函數(shù)圖像：clc;clear;closeall;%%y'=y-2*x/y?Adams來求解驗(yàn)證%要求匿名函數(shù)?持向量化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內(nèi)插法ifstrcmp(method,'out')stra='外插法';elseifstrcmp(method,'in')stra='內(nèi)插法';endaccuracy=5;%%h=1/10下逐點(diǎn)誤差%參數(shù)：外插/內(nèi)插，導(dǎo)函數(shù)，因變量初值，?變量左端點(diǎn)，步長，?變量右端點(diǎn)，精度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,'下的數(shù)值解圖像與解析解圖像']);legend('數(shù)值解','解析解');disp('逐點(diǎn)誤差：');disp(err);

%%求誤差與步長的函數(shù)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);%求整體誤差的階，等于對數(shù)下的誤差與步長函數(shù)的斜率Deh=(lerr(2:end)-lerr(1:end-1))./(lh(2:end)-lh(1:end-1));deh=mean(Deh);%畫對數(shù)下誤差函數(shù)與步長的函數(shù)圖像figure(4);plot(log(h),log(err));title([num2str(round(deh)),'階精度Adams',stra,'下的誤差對數(shù)函數(shù)圖像']);disp(['此?法具有',num2str(round(deh)),'階精度']);5.2?函數(shù)1，實(shí)現(xiàn)具有1到6階精度的外插及內(nèi)插Adams?法：%求任意1到6階的數(shù)值解；%?先利?外置函數(shù)RungeKutta_auto來求啟動所需要的點(diǎn)，再?Adams外插法或者內(nèi)插法來求剩余的點(diǎn)%預(yù)備?作：根據(jù)所選擇?法求出所需啟動點(diǎn)%外插法流程：（1）求系數(shù)a(2)求系數(shù)b(3)求剩余點(diǎn)%內(nèi)插法流程:(1)求系數(shù)a(2)求系數(shù)b（3）求系數(shù)a*(4)求系數(shù)b*(5)求剩余點(diǎn)P(EC)...%注意事項(xiàng)：要求導(dǎo)函數(shù)?持向量化%完成時間：2022/3/20嫁佛幡function[u]=Adams_auto(method,func,y0,x0,h,xmax,accuracy)%[待求點(diǎn)]=[外插/內(nèi)插，導(dǎo)函數(shù)，因變量初值，?變量左端點(diǎn)，步長，?變量右端點(diǎn)，精度]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')%單步法求啟動點(diǎn)，ay階精度需要ay個點(diǎn)u(1:ay)=RungeKutta_auto(ay,func,y0,x0,h,x(ay));%求系數(shù)a,b[a,b]=Adams_arg('out',ay);%求剩余點(diǎn)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%單步法求啟動點(diǎn)，ay階精度需要ay-1個點(diǎn)，u(1:ay-1)=RungeKutta_auto(ay,func,y0,x0,h,x(ay-1));%求系數(shù)a,b[a,b]=Adams_arg('out',ay-1);%求系數(shù)a1,b1[a1,b1]=Adams_arg('in',ay);%求剩余點(diǎn)forj=ay-1:N-1%預(yù)估u_pre=u(j)+h*func(x(j:-1:j-ay+2),u(j:-1:j-ay+2))*b;%校準(zhǔn)tol=1e-7delta=1;u_part=h*func(x(j:-1:j-ay+2),u(j:-1:j-ay+2))*b1(2:end);%避免重復(fù)運(yùn)算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?函數(shù)2，單步法求具有1到6階精度的啟動點(diǎn)：function[u]=RungeKutta_auto(accuracy,func,y0,x0,h,xmax)%[待求點(diǎn)函數(shù)值]=[精度，導(dǎo)函數(shù)，因變量初值，?變量左端點(diǎn)，步長，因變量右端點(diǎn)]%參數(shù)初始化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%中點(diǎn)法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%經(jīng)典四階龍格庫塔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?函數(shù)3，求Adams外插法或內(nèi)插法的系數(shù)a,b：function[A,B]=Adams_arg(method,n)%[Adams系數(shù)a,b]=[?法，精度n]%輸出兩個列向量%系數(shù)對?法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