Adams与Milline-Hamming预测校正公式例题

例：y(0)=1(0x1)本题采用adams与Miline-Haming预测校正公式求解的一种解法。f=sym('y-2*x/y');%sym为定义符号变量函数。例f=sym(‘sin2’)输出结果为f=sin2x=0;y=1;%初始条件h=0.1;%步长a=0;b=1;%上下线n=(b-a)/h+1;%步数sjs=[x,y];%先定义一个数值解；sjs符号自定义即可，没有要求adx=0;%adams方法的x初始值；%下面的是龙格库塔方法求解，这里求得的是[0.1]中的所有解，为和后面的方法进行比较fori=1:n-1RK(i)=y;L1=subs(f,{'x','y'},{x,y});L2=subs(f,{'x','y'},{x+h/2,y+h*L1/2});%subs是替代函数，也就是用后面的代替前面的%例如：L2中，用(x+h/2,y+h*L1/2)代替原(x,y);此时函数应为f(x+h/2,y+h*L1/2);%下面同理L3=subs(f,{'x','y'},{x+h/2,y+h*L2/2});L4=subs(f,{'x','y'},{x+h,y+h*L3});y=y+h*(L1+2*L2+2*L3+L4)/6;x=x+h;adx(i+1)=x;sjs=[sjs;x,y];endady(1)=RK(1);mh(1)=RK(1);%ady为adams方法的y值，而mh为另一种方法的y值ady(2)=RK(2);mh(2)=RK(2);%因为这两种方法均要用龙格库塔方法求出前几项ady(3)=RK(3);mh(3)=RK(3);ady(4)=RK(4);mh(4)=RK(4);%下面是adams预测-校正方法forj=5:npre1=ady(j-1)+h*(55*subs(f,{'x','y'},{adx(j-1),ady(j-1)})-59*subs(f,{'x','y'},{adx(j-2),ady(j-2)})+37*subs(f,{'x','y'},{adx(j-3),ady(j-3)})-9*subs(f,{'x','y'},{adx(j-4),ady(j-4)}))/24;ady(j)=ady(j-1)+h*(9*subs(f,{'x','y'},{adx(j),pre1})+19*subs(f,{'x','y'},{adx(j-1),ady(j-1)})-5*subs(f,{'x','y'},{adx(j-2),ady(j-2)})+subs(f,{'x','y'},{adx(j-3),ady(j-3)}))/24;end%下面是Miline-Hamming预测校正方法fork=5:11pre2=mh(k-4)+4*h*(2*subs(f,{'x','y'},{adx(k-1),mh(k-1)})-subs(f,{'x','y'},{adx(k-2),mh(k-2)})+2*subs(f,{'x','y'},{adx(k-3),mh(k-3)}))/3;mh(k)=(9*mh(k-1)-mh(k-3))/8+3*h/8*(subs(f,{'x','y'},{adx(k),pre2})+2*subs(f,{'x','y'},{adx(k-1),mh(k-1)})-subs(f,{'x','y'},{adx(k-2),mh(k-2)}));end%输出结果并打印图像sjs%RK公式结果plot(sjs(:,1),sjs(:,2),'r')%RK公式图像holdonady%adams公式结果plot(sjs(:,1),ady,'ob')%adams公式图像holdonexact=sqrt(1+2*sjs(:,1))%精确解plot(sjs(:,1),exact,'*g');%精确解图像mh%mh公式结果holdonplot(sjs(:,1),mh,'k+');%mh公式图像