amendhamming.m

研一刚上完数值分析

function[xn,C]=AmendHamming(fun,x0,y0,h,N)
%修正Hamming公式，华中科技大学出版，李红写的《数值分析》书中5.6.1节
%x0,y0为初始条件；
%h为区间步长；
%N为区间的个数；
%x为Xn构成的向量；
%y为Yn构成的向量；
[xn,yn]=ode23(fun,[x0:h:x0+h*N],y0);
f=fun(xn,yn);
P=zeros(11,1);
C=zeros(11,1);
m=zeros(11,1);
md=zeros(11,1);

for n=1:4
    C(n)=yn(n);
    P(n)=C(n);
end

for n=4:N
    P(n+1)=C(n-3)+4/3*h*(2*f(n)-f(n-1)+2*f(n-2));
    m(n+1)=P(n+1)-112/121*(P(n+1)-C(n+1));
    md(n+1)=fun(xn(n+1),m(n+1));
    C(n+1)=1/8*(9*C(n)-C(n-2)+3*h*(md(n+1)+2*f(n)-2*f(n-1)));
    C(n+1)=C(n+1)+9/121*(P(n+1)-C(n+1));
    f(n+1)=fun(xn(n+1),C(n+1));
end