oulahoucha.m

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

function [yt,th]=oulahoucha(f,y0,t0,td,h)
y=y0;
t=t0;

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

for i=2:fix((td-t0)/h)+1
fn=double(eval(f));                           %solve fn
yt(:,i)=yt(:,i-1)+fn'.*h;                     %solve yn+1
ytemp=yt(:,i);                                %ytemp=yn+1
y=yt(:,i);                                    %use current y(yn+1) to obtain fn+1

fn_1=double(eval(f));                         %solve fn+1
yt(:,i)=yt(:,i-1)+fn_1'.*h;                   %use fn+1,solve image yn+1

 while abs(ytemp-yt(:,i))>0.00001             %if        
        ytemp=yt(:,i);                        %ytemp=yn+1
        y=yt(:,i);                            %use current y(yn+1) to obtain fn+1
        fn_1= double(eval(f));                %solve a new fn+1
        yt(:,i)=yt(:,i-1)+fn_1'.*h;           %solve a new image yn+1
 end

th(i)=th(i-1)+h;                              %obtain new time point
y=yt(:,i);                                    %save yn+1 to y,use current y(yn+1) to obtain fn+1
t=th(i);           
end

%n=length(f);
%C=zeros(1,n);C(1)=1;
%out=C*yt;