rungekutta_3.m

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

function [yt,th]=rungekutta_3(f,y0,t0,td,h)   %our book page 36 function 3_57
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
k1=fn;                                          %solve k1=fn

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

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

yt(:,i)=yt(:,i-1)+(k1+3*k3)'.*h/4;              %use yn,k1,k2,k3,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

%n=length(f);                                    %state space function's output
%C=zeros(1,n);C(1)=1;
%out=C*yt;