hamming.m

1．编写欧拉前差、后差、梯形公式。 2．编写二阶、三阶龙格库塔法通用程序。 3．编写汉明积分法通用程序. 4．编写用状态转移法对连续系统状态方程进行离散化的通用程序。

function [yt,th]=hamming(f,y0,t0,td,h)          %our book page 38 hamming method
y=y0;                                           %begin with rungekutta_4
t=t0;

fn=zeros(length(f),fix((td-t0)/h)+1);           %define length of fn

yt(:,1)=double(y0);                             %initialize yt,th
th(1)=double(t0);

fn(:,1)=double(eval(f));                        %initialize fn

for i=2:4
fn(:,i-1)=double(eval(f));                      %solve fn
k1=fn(:,i-1);                                   %solve k1: k1=fn

y=yt(:,i-1)+k1'.*h/2;                           %obtain yn+1/2*hk1
t=th(i-1)+h/2;                                  %obtain tn+1/2*h
k2=double(eval(f));                             %solve k2=f(tn+1/2*h,yn+1/2*hk1)

y=yt(:,i-1)+k2'.*h/2;                           %obtain yn+1/2*hk2
t=th(i-1)+h/2;                                  %obtain tn+2/3*h
k3=double(eval(f));                             %solve k3=f(tn+1/2*h,yn+1/2*hk2)

y=yt(:,i-1)+k3'.*h;                             %obtain yn+hk3
t=th(i-1)+h;                                    %obtain tn+h
k4=double(eval(f));                             %solve k4=f(tn+1/2*h,yn+1/2*hk2)

yt(:,i)=yt(:,i-1)+(k1+2*k2+2*k3+k4)'.*h/6;      %use yn,k1,k2,k3,k4,solve yn+1

th(i)=th(i-1)+h;                                %next time point
y=yt(:,i);                                      %this two is used to solve the next fn
t=th(i);                                        %as above 
end

j=1;                                                                           %do yt5,由于yt1为y0，故yt5为y4
fn(:,4)=double(eval(f));                                                       %solve fn4为实际f3 ,fn1为实际f0
y_366(:,j)=yt(:,1)+4*(2*fn(:,4)-fn(:,3)+2*fn(:,2))'.*h/3;                      %yt(:,1)=y0, page 38,function 3-66
th(5)=th(4)+h;
y=y_366(:,j);                                   
t=th(5);
fn(:,5)=double(eval(f));                                                       %solve fn5为实际f4
yt(:,5)=(9*yt(:,4)-yt(:,2))/8+3*(fn(:,5)+2*fn(:,4)-fn(:,3))'.*h/8;             %function 3-68，yt5为实际y4

while abs(y_366(:,j)-yt(:,5))>0.0001                                           %迭代校正
    j=j+1;
    y_366(:,j)=yt(:,5);                                                        %y_366的j相当于（i），用来反复迭代，校正yn+1
    y=yt(:,5);
    fn(:,5)=double(eval(f));                                                   %公式3-68中的f+1（i）  
    yt(:,5)=(9*yt(:,4)-yt(:,2))/8+3*(fn(:,5)+2*fn(:,4)-fn(:,3))'.*h/8;         %function3-68
end

y=yt(:,5);

for p=6:fix((td-t0)/h)+1                                                       %do after yt6,实际是求y5之后的数
    fn(:,p-1)=double(eval(f));
    j=1;
    y_367(:,j)=yt(:,p-4)+4*(2*fn(:,p-1)-fn(:,p-2)+2*fn(:,p-3))'.*h/3;          %function 3-66
    y_367(:,1)=y_367(:,1)+112*(yt(:,p-1)-y_366(:,1))/121;                      %function 3-67
    th(p)=th(p-1)+h;
    y=y_367(:,j);
    t=th(p);
    fn(:,p)=double(eval(f));
    yt(:,p)=(9*yt(:,p-1)-yt(:,p-3))/8+3*(fn(:,p)+2*fn(:,p-1)-fn(:,p-2))'.*h/8; %function 3-68
    y_366=y_367;                                                               %经过3-67校正过的y赋值给3-68
    while abs(y_366(:,j)-yt(:,p))>0.0001                                       %迭代校正
        j=j+1;
        y_366(:,j)=yt(:,p);
        y=yt(:,p);
        fn(:,p)=double(eval(f));
        yt(:,p)=(9*yt(:,p-1)-yt(:,p-3))/8+3*(fn(:,p)+2*fn(:,p-1)-fn(:,p-2))'.*h/8;
    end
    y=yt(:,p);
end