Article directory
-
- summary
- References
- Code interpretation
-
-
- ode45 function header
- Initial parameter settings
- DP54 format parameter settings
- Initial step size setting
- DP54 format calculation
- error estimate
- Update step size
- Accept step size
- End single step loop
- End of program
-
- Detailed review
- summary
- Replenish
Summary
This article introduces the variable step strategy used in matlab’s built-in function ode45 and the interpretation of the corresponding code.
The matlab version used by the author is MTALAB R2018b
Reference materials
matlab ordinary differential equation solver technical documentation
The MATLAB ODE Suite, L. F. Shampine and M. W. Reichelt, SIAM Journal on Scientific Computing, 18-1, 1997
Original paper in DP54 format
J. R. DORMAND AND P. J. PRINCE, A family of embedded Runge-Kutta formulae, J. Comp. Appl. Math., 6 (1980), pp. 19-26.
Code interpretation
ode45 function header
% ode is the right-hand term of the equation, tspan is the solution interval, y0 is the initial value, and options is the parameter structure returned by odeset function varargout = ode45(ode,tspan,y0,options,varargin)
Initial parameter settings
solver_name = 'ode45'; % solver name % Stats nsteps = 0; % number of solution steps nfailed = 0; % number of re-iterations nfevals = 0; % number of function evaluations %Output FcnHandlesUsed = isa(ode,'function_handle'); % Whether to use function handles output_sol = (FcnHandlesUsed & amp; & amp; (nargout==1));%sol=odeXX(...) returns the function handle sol=[]; output_ty = (~output_sol & amp; & amp; (nargout > 0));%[t,y,...]=odeXX(...) returns parameter matrix % There might be no output requested... % Handle solver arguments % parameter settings [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ... options, threshold, rtol, normcontrol, normy, hmax, htry, htspan, dataType] = ... odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Go to odarguments function:
function [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, args, odeFcn, ...
options, threshold, rtol, normcontrol, normy, hmax, htry, htspan, ...
dataType ] = ...
odarguments(FcnHandlesUsed, solver, ode, tspan, y0, options, extras)
tspan = tspan(:);
ntspan = length(tspan);
if ntspan == 1 % Integrate from 0 to tspan
t0 = 0;
next = 1; % Next entry in tspan.
else
t0 = tspan(1);
next = 2; % next entry in tspan
end
htspan = abs(tspan(next) - t0);
tfinal = tspan(end);
y0 = y0(:);
neq = length(y0);
tdir = sign(tfinal - t0);
f0 = feval(ode,t0,y0,args{<!-- -->:});
dataType = superiorfloat(t0,y0,f0); % highest precision type
% Get the error control options, and set defaults.
rtol = odeget(options,'RelTol',1e-3,'fast'); % relative error setting
if rtol < 100 * eps(dataType)
rtol = 100 * eps(dataType); % The error cannot be less than 100 times the minimum accuracy that can be recognized by the current type
end
atol = odeget(options,'AbsTol',1e-6,'fast'); % absolute error setting
atol = atol(:);
threshold = atol / rtol; % absolute error and relative error ratio, convenient for calculation and use
% By default, hmax is 1/10 of the interval.
% Maximum step size setting, the default is 1/10 of the interval
safehmax = 16.0*eps(dataType)*max(abs(t0),abs(tfinal)); % numerical safety step size
defaulthmax = max(0.1*(abs(tfinal-t0)), safehmax);
hmax = min(abs(tfinal-t0), abs(odeget(options,'MaxStep',defaulthmax,'fast')));
htry = abs(odeget(options,'InitialStep',[],'fast')); % initial step setting
odeFcn = ode;
end
Return to ode45 function:
refine = max(1,odeget(options,'Refine',4,'fast')); % number of interpolations = refine-1
if ntspan > 2
outputAt = 1; % output only at tspan points Interpolation is required only at specified points
elseif refine <= 1
outputAt = 2; % computed points, no refinement no interpolation required
else
outputAt = 3; % computed points, with refinement requires interpolation
S = (1:refine-1) / refine;
end
t = t0;
y = y0;
noout = 0;
tout = []; yout = []; % initialization of output parameters
if nargout > 0
if ntspan > 2 % output only at tspan points Output values only at specified points
tout = zeros(1,ntspan,dataType);
yout = zeros(neq,ntspan,dataType);
else % alloc in chunks
chunk = min(max(100,50*refine), refine + floor((2^13)/neq));
tout = zeros(1,chunk,dataType);
yout = zeros(neq,chunk,dataType);
end
noout = 1;
tout(nout) = t;
yout(:,nout) = y;
end
Pay attention to the refine parameter in the above code. Not all of the output of matlab is obtained by forwarding the corresponding format. Some solutions (default is 3/4) are obtained by interpolation formulas
DP54 format parameter settings
% Initialize method parameters. pow = 1/5; A = [1/5, 3/10, 4/5, 8/9, 1, 1]; % Still used by restarting criteria % B = [ % 1/5 3/40 44/45 19372/6561 9017/3168 35/384 % 0 9/40 -56/15 -25360/2187 -355/33 0 % 0 0 32/9 64448/6561 46732/5247 500/1113 % 0 0 0 -212/729 49/176 125/192 % 0 0 0 0 -5103/18656 -2187/6784 % 0 0 0 0 0 11/84 % 0 0 0 0 0 0 % ]; % E = [71/57600; 0; -71/16695; 71/1920; -17253/339200; 22/525; -1/40]; % Same values as above extracted as scalars (1 and 0 values committed) a2=cast(1/5,dataType); % Convert data type a3=cast(3/10,dataType); a4=cast(4/5,dataType); a5=cast(8/9,dataType); b11=cast(1/5,dataType); b21=cast(3/40,dataType); b31=cast(44/45,dataType); b41=cast(19372/6561,dataType); b51=cast(9017/3168,dataType); b61=cast(35/384,dataType); b22=cast(9/40,dataType); b32=cast(-56/15,dataType); b42=cast(-25360/2187,dataType); b52=cast(-355/33,dataType); b33=cast(32/9,dataType); b43=cast(64448/6561,dataType); b53=cast(46732/5247,dataType); b63=cast(500/1113,dataType); b44=cast(-212/729,dataType); b54=cast(49/176,dataType); b64=cast(125/192,dataType); b55=cast(-5103/18656,dataType); b65=cast(-2187/6784,dataType); b66=cast(11/84,dataType); e1=cast(71/57600,dataType); e3=cast(-71/16695,dataType); e4=cast(71/1920,dataType); e5=cast(-17253/339200,dataType); e6=cast(22/525,dataType); e7=cast(-1/40,dataType);
Initial step size setting
hmin = 16*eps(t); % minimum step size
if isempty(htry)
% Compute an initial step size h using y'(t).
% Calculate the initial step size based on the value of the right-hand term of the initial point
absh = min(hmax, htspan);
rh = norm(f0 ./ max(abs(y),threshold),inf) / (0.8 * rtol^pow);
if absh * rh > 1
absh = 1 / rh;
end
absh = max(absh, hmin);
else
absh = min(hmax, max(hmin, htry));
end
f1 = f0;
The above initial step size is generally
0.8
?
m
a
x
(
absolute error
,
Relative error
?
y
)
f
?
Relative error
4
5
0.8*{{max(absolute error, relative error*y)}\over{f*{relative error}^{{4}\over{5} }}}
0.8?f? Relative error 54?max (absolute error, relative error?y)?
DP54 format calculation
% Cleanup the main ode function call
FcnUsed = isa(odeFcn,'function_handle');
odeFcn_main = odefcncleanup(FcnUsed,odeFcn,{<!-- -->});
% THE MAIN LOOP
done = false;
while~done
% By default, hmin is a small number such that t + hmin is only slightly
% different than t. It might be 0 if t is 0.
hmin = 16*eps(t);
absh = min(hmax, max(hmin, absh)); % couldn't limit absh until new hmin
h = tdir * absh; % current step size
% Stretch the step if within 10% of tfinal-t.
% Whether it reaches the right endpoint of the solution interval
if 1.1*absh >= abs(tfinal - t)
h = tfinal - t;
absh = abs(h);
done = true;
end
% LOOP FOR ADVANCING ONE STEP.
%Single step loop
nofailed = true; % no failed attempts
while true
y2 = y + h .* (b11.*f1 );
t2 = t + h .* a2;
f2 = odeFcn_main(t2, y2);
y3 = y + h .* (b21.*f1 + b22.*f2 );
t3 = t + h .* a3;
f3 = odeFcn_main(t3, y3);
y4 = y + h .* (b31.*f1 + b32.*f2 + b33.*f3 );
t4 = t + h .* a4;
f4 = odeFcn_main(t4, y4);
y5 = y + h .* (b41.*f1 + b42.*f2 + b43.*f3 + b44.*f4 );
t5 = t + h .* a5;
f5 = odeFcn_main(t5, y5);
y6 = y + h .* (b51.*f1 + b52.*f2 + b53.*f3 + b54.*f4 + b55.*f5 );
t6 = t + h;
f6 = odeFcn_main(t6, y6);
tnew = t + h;
if done
tnew = tfinal; % Hit end point exactly.
end
h = tnew - t; % Purify h.
ynew = y + h.* ( b61.*f1 + b63.*f3 + b64.*f4 + b65.*f5 + b66.*f6 ); % Calculated numerical solution
f7 = odeFcn_main(tnew,ynew);
nfevals = nfevals + 6; % performed 6 valuations
Error estimate
% Estimate the error.
NNrejectStep = false;
fE = f1*e1 + f3*e3 + f4*e4 + f5*e5 + f6*e6 + f7*e7; % error estimate
err = absh * norm((fE) ./ max(max(abs(y),abs(ynew)),threshold),inf);
% Accept the solution only if the weighted error is no more than the
% tolerance rtol. Estimate an h that will yield an error of rtol on
% the next step or the next try at taking this step, as the case may be,
% and use 0.8 of this value to avoid failures.
if err > rtol % Failed step
nfailed = nfailed + 1;
The condition for not accepting this step size is
h
?
f
E
>
m
a
x
(
absolute error
,
Relative error
?
m
a
x
(
y
n
,
y
n
+
1
)
)
h*f_E>max(absolute error, relative error*max(y_n,y_{n + 1}))
h?fE?>max(absolute error, relative error?max(yn?,yn + 1?))
Update step size
if nofailed
nofailed = false;
absh = max(hmin, absh * max(0.1, 0.8*(rtol/err)^pow));
else
absh = max(hmin, 0.5 * absh);
end
h = tdir * absh;
done = false;
The step size update formula is
h
n
e
w
=
h
?
0.8
?
(
m
a
x
(
absolute error
,
Relative error
?
m
a
x
(
y
n
,
y
n
+
1
)
)
h
?
f
E
)
1
5
h_{new}=h*0.8*({{max(absolute error, relative error*max(y_n,y_{n + 1}))}\over{h*f_E}}) ^{{1}\over{5}}
hnew?=h?0.8?(h?fE?max(absolute error, relative error?max(yn?,yn + 1?))?)51?
The minimum update is
h
n
e
w
=
h
?
0.1
h_{new}=h*0.1
hnew?=h?0.1
If the step size is still not accepted, then
h
n
e
w
=
h
?
0.5
h_{new}=h*0.5
hnew?=h?0.5
Accept step size
else % Successful step
break; % jump out of while true loop
end
end
nsteps = nsteps + 1;
% generate output
if output_ty || haveOutputFcn
switch outputAt
case 2 % computed points, no refinement, no interpolation, output directly
nout_new = 1;
tout_new = tnew;
yout_new = ynew;
case 3 % computed points, with refinement for interpolation
tref = t + (tnew-t)*S;
nout_new = refine;
tout_new = [tref, tnew];
yntrp45 = ntrp45split(tref,t,y,h,f1,f3,f4,f5,f6,f7,[]); % ode45 interpolation function
yout_new = [yntrp45, ynew];
case 1 % output only at tspan points interpolation output only at specific points
nout_new = 0;
tout_new = [];
yout_new = [];
while next <= ntspan
if tdir * (tnew - tspan(next)) < 0
break; % Interpolate when the evaluation point crosses the output point
end
nout_new = nout_new + 1;
tout_new = [tout_new, tspan(next)];
if tspan(next) == tnew
yout_new = [yout_new, ynew];
else
yntrp45 = ntrp45split(tspan(next),t,y,h,f1,f3,f4,f5,f6,f7,[]);
yout_new = [yout_new, yntrp45];
end
next = next + 1;
end
end
% write output
if nout_new > 0
if output_ty
oldnout = nout;
nout = nout + nout_new;
if nout > length(tout)
tout = [tout, zeros(1,chunk,dataType)]; % requires chunk >= refine
yout = [yout, zeros(neq,chunk,dataType)];
end
idx = (oldnout + 1):nout;
tout(idx) = tout_new;
yout(:,idx) = yout_new;
end
end
end
% If there were no failures compute a new h.
% The step size still needs to be changed when this step size is accepted
if nofailed
% Note that absh may shrink by 0.8, and that err may be 0.
temp = 1.25*(err/rtol)^pow;
if temp > 0.2
absh = absh / temp;
else
absh = 5.0*absh;
end
end
The step size update formula here is the same as above, but the step size is required to be increased by up to five times.
End single step loop
% Advance the integration one step. t = tnew; y = ynew; f1 = f7; % Already have f(tnew,ynew) end
End of program
solver_output = {<!-- -->};
if (nout > 0) % produce output
if isempty(sol) % output [t,y,...]
solver_output{<!-- -->1} = tout(1:nout).';
solver_output{<!-- -->2} = yout(:,1:nout).';
solver_output{<!-- -->end + 1} = [nsteps, nfailed, nfevals]; % Column vector
end
end
if nargout > 0
varargout = solver_output; % ode45 output
end
Details review
- The necessary parameters include the right-hand term of the equation (ode), the initial value (y0), and the solution interval (tspan)
- Optional parameters include absolute error (atol, default 1e-6), relative error (rtol, default 1e-3), maximum step size (hmax, default 1/10 of the interval), initial step size (htry), interpolation flag (refine, default is 4)
- When the initial step is not set, an initial step will be calculated based on the value of the initial right-hand term.
- Whether to perform interpolation will be selected based on the interpolation flag.
- The step update rate is limited to between 0.1 and 5.0
- Step size update formula
h
n
e
w
=
h
?
0.8
?
(
m
a
x
(
absolute error
,
Relative error
?
m
a
x
(
y
n
,
y
n
+
1
)
)
h
?
f
E
)
1
5
h_{new}=h*0.8*({{max(absolute error, relative error*max(y_n,y_{n + 1}))}\over{h*f_E}}) ^{{1}\over{5}}
hnew?=h?0.8?(h?fE?max(absolute error, relative error?max(yn?,yn + 1?))?)51?
Summary
This article only retains the variable step size and general ordinary differential equation solution parts in the ode45 source code, and removes other irrelevant variable step content (such as the settings and functions of other optional parameters; in the form
M
(
t
)
y
‘
=
f
(
x
,
y
)
M(t)y’=f(x,y)
M(t)y′=f(x,y) solution of ordinary differential equations, etc.), but this does not affect the operation of the variable step part of the above code. The variable step strategy of ode45 is not complicated. The technical document in the reference material was issued in 1997. Most of the space is about matlab’s efforts in solving rigid problems. ode45 only accounts for a small part of it. Even the variable step size strategy is not mentioned. The original paper on DP54 format was published in 1980. The variable step size strategy mentioned in the article is very similar to the one used now. When solving equations, ode45 will improve the solving efficiency and the smoothness of the drawing through interpolation, and does not have too strict requirements on the solving accuracy.
Supplement
Two private functions of MATLAB (odefcncleanup and ntrp45split) appear in the above code. Users cannot call them directly. They can be copied to the user directory and called.