rungekutta_4.m

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

function [yt,th]=rungekutta_4(f,y0,t0,td,h)   %our book page 37 function 3_58
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/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

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