rk_4.m

内涵大量基于MATLAB的仿真算法

function [tout,yout] = rk_4(odefile,t0,h,th,y0)
%RK_4 实现了 4 阶 Runge-Kutta 定步长微分方程求解。
%
%  [tout,yout] = rk_4(odefile,t0,h,th,y0)
%
%  odefile 为描述微分方程的 M 函数， t0, h 和 th 分别为
%  初始时间、步长和终止时间， y0 为状态变量初值向量，返回
%  的 tout 和 yout 分别为时间向量和状态变量解矩阵。

%Designed by Prof D Xue (c) 2000
tout=[t0:h:th]'; yout=[];
for t=tout'
   k1=h*eval([odefile '(t,y0)']);
   k2=h*eval([odefile '(t+h/2,y0+0.5*k1)']);
   k3=h*eval([odefile '(t+h/2,y0+0.5*k2)']);
   k4=h*eval([odefile '(t+h,y0+k3)']);
   y0=y0+(k1+2*k2+2*k3+k4)/6;
   yout=[yout; y0'];
end