rk1.m

来自「(有源代码)数值分析作业,本文主要包括两个部分,第一部分是常微分方程(ODE)的」· M 代码 · 共 27 行

function [x,y] = RK(ODEfile,xi,xf,h,yi,varargin)

% if isempty(h) | h == 0
%    h = (xf - xi)/99;
% end

n = fix(n);

yi = (yi(:).')';	% Make sure it's a column vector

x = [xi:h:xf];		% Vector of x values
% if x(end) ~= xf
%    x(end+1) = xf;
% end
% d = diff(x);		% Vector of x-increments

y(:,1) = yi;		% Initial condition

   for i = 1:length(x)-1
      k1 = d(i) * feval(ODEfile,x(i),y(:,i),varargin{:});
      k2 = d(i) * feval(ODEfile,x(i)+d(i)/2,y(:,i)+k1/2, ...
         varargin{:});
      k3 = d(i) * feval(ODEfile,x(i)+d(i)/2,y(:,i)+k2/2, ...
         varargin{:});
      k4 = d(i) * feval(ODEfile,x(i+1),y(:,i)+k3,varargin{:});
      y(:,i+1) = y(:,i) + (k1+2*k2+2*k3+k4)/6;
   end