oulapingjun.m

来自「1．编写欧拉前差、后差、梯形公式。 2．编写二阶、三阶龙格库塔法通用程序。 」· M 代码 · 共 31 行

function [yt,th]=oulapingjun(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+fn)'.*h/2;            %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+fn)'.*h/2;           %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;