93 lines
2.6 KiB
Mathematica
93 lines
2.6 KiB
Mathematica
|
function Y = ode5(odefun,tspan,y0,varargin)
|
||
|
|
%ODE5 Solve differential equations with a non-adaptive method of order 5.
|
||
|
|
% Y = ODE5(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates
|
||
|
|
% the system of differential equations y' = f(t,y) by stepping from T0 to
|
||
|
|
% T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.
|
||
|
|
% The vector Y0 is the initial conditions at T0. Each row in the solution
|
||
|
|
% array Y corresponds to a time specified in TSPAN.
|
||
|
|
%
|
||
|
|
% Y = ODE5(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters
|
||
|
|
% P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).
|
||
|
|
%
|
||
|
|
% This is a non-adaptive solver. The step sequence is determined by TSPAN
|
||
|
|
% but the derivative function ODEFUN is evaluated multiple times per step.
|
||
|
|
% The solver implements the Dormand-Prince method of order 5 in a general
|
||
|
|
% framework of explicit Runge-Kutta methods.
|
||
|
|
%
|
||
|
|
% Example
|
||
|
|
% tspan = 0:0.1:20;
|
||
|
|
% y = ode5(@vdp1,tspan,[2 0]);
|
||
|
|
% plot(tspan,y(:,1));
|
||
|
|
% solves the system y' = vdp1(t,y) with a constant step size of 0.1,
|
||
|
|
% and plots the first component of the solution.
|
||
|
|
|
||
|
|
if ~isnumeric(tspan)
|
||
|
|
error('TSPAN should be a vector of integration steps.');
|
||
|
|
end
|
||
|
|
|
||
|
|
if ~isnumeric(y0)
|
||
|
|
error('Y0 should be a vector of initial conditions.');
|
||
|
|
end
|
||
|
|
|
||
|
|
h = diff(tspan);
|
||
|
|
if any(sign(h(1))*h <= 0)
|
||
|
|
error('Entries of TSPAN are not in order.')
|
||
|
|
end
|
||
|
|
|
||
|
|
try
|
||
|
|
f0 = feval(odefun,tspan(1),y0,varargin{:});
|
||
|
|
catch
|
||
|
|
msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];
|
||
|
|
error(msg);
|
||
|
|
end
|
||
|
|
|
||
|
|
y0 = y0(:); % Make a column vector.
|
||
|
|
if ~isequal(size(y0),size(f0))
|
||
|
|
error('Inconsistent sizes of Y0 and f(t0,y0).');
|
||
|
|
end
|
||
|
|
|
||
|
|
neq = length(y0);
|
||
|
|
N = length(tspan);
|
||
|
|
Y = zeros(neq,N);
|
||
|
|
|
||
|
|
% Method coefficients -- Butcher's tableau
|
||
|
|
%
|
||
|
|
% C | A
|
||
|
|
% --+---
|
||
|
|
% | B
|
||
|
|
|
||
|
|
C = [1/5; 3/10; 4/5; 8/9; 1];
|
||
|
|
|
||
|
|
A = [ 1/5, 0, 0, 0, 0
|
||
|
|
3/40, 9/40, 0, 0, 0
|
||
|
|
44/45 -56/15, 32/9, 0, 0
|
||
|
|
19372/6561, -25360/2187, 64448/6561, -212/729, 0
|
||
|
|
9017/3168, -355/33, 46732/5247, 49/176, -5103/18656];
|
||
|
|
|
||
|
|
B = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84];
|
||
|
|
|
||
|
|
% More convenient storage
|
||
|
|
A = A.';
|
||
|
|
B = B(:);
|
||
|
|
|
||
|
|
nstages = length(B);
|
||
|
|
F = zeros(neq,nstages);
|
||
|
|
|
||
|
|
Y(:,1) = y0;
|
||
|
|
for i = 2:N
|
||
|
|
ti = tspan(i-1);
|
||
|
|
hi = h(i-1);
|
||
|
|
yi = Y(:,i-1);
|
||
|
|
|
||
|
|
% General explicit Runge-Kutta framework
|
||
|
|
F(:,1) = feval(odefun,ti,yi,varargin{:});
|
||
|
|
for stage = 2:nstages
|
||
|
|
tstage = ti + C(stage-1)*hi;
|
||
|
|
ystage = yi + F(:,1:stage-1)*(hi*A(1:stage-1,stage-1));
|
||
|
|
F(:,stage) = feval(odefun,tstage,ystage,varargin{:});
|
||
|
|
end
|
||
|
|
Y(:,i) = yi + F*(hi*B);
|
||
|
|
|
||
|
|
end
|
||
|
|
Y = Y.';
|