Hi Nils,
yes, I've further questions :) I nearly made all of your changes, but if
I try to solve the equation
f = ode45(@(t,y)y, [0 -1], 1)
then the endpoint still is not hit exactly. So I wonder what I am still
doing wrong? Can you take another look at my modifiations (cf. attachment).
The second thing is thet with the change of line 435 from
vstepsize = max (vodeoptions.MaxStep, ...
to
vstepsize = max (-vodeoptions.MaxStep, ...
I cannot solve the equation from before at all. All the other changes
should be in the file...
Thanks for help again
best regards
Thomas
Am 05.06.11 18:10, schrieb Nils Strunk:
Hi Thomas,
ad (2): You're right. I already noticed the same and tried to fix it.
I think, the best way to avoid this problem is to change line 318 from
vstepsize = vtimestop - vdirection * vtimestamp;
to
vstepsize = vdirection * abs(abs(vtimestop) - abs(vtimestamp));
I didn't carefully check if it solves the problem in general, but I
think it will.
I changed a few more things, that you may carry over to the official
version as well:
(1) If you solve an ODE backwards, and want to define InitialStep and
MaxStep, then you'll notice that InitialStep should be negativ,
whereas MaxStep couldn't be defined negative (so is supposed to be
positive). On the one hand it's confusing and on the other hand this
leads to problems in the algorithm. If you want to fix wrong signs in
InitialStep (so, if you choose the correct sign in the definition of
InitialStep doesn't matters any more), you have to replace line 272
vstepsize = vodeoptions.InitialStep;
by
if (sign(vodeoptions.InitialStep) == vdirection)
vstepsize = vodeoptions.InitialStep;
else %# fix wrong direction of InitialStep
vstepsize = -vodeoptions.InitialStep;
end
Furthermore (and this is really a bug!) you should change line 435
vstepsize = max (vodeoptions.MaxStep, ...
to
vstepsize = max (-vodeoptions.MaxStep, ...
(2) In line 349 there is the following variable defined:
vSaveVUForRefine = vu;
which isn't used in common. So the following if-clause will prevent to
save that variable if it's not needed:
if (vhaveoutputfunction && vhaverefine)
vSaveVUForRefine = vu;
end
(3) Also a minor change: Add "\n" at the end of "fprintf" in lines
510-512.
If you encounter any problems with these changes, or if you have any
questions regarding to this, please don't hesitate to write me.
Kind regards,
Nils
Am 04.06.2011 10:23, schrieb Thomas Treichl:
Hi Nils,
I've taken a look at the solver problem that you've reported:
(1) Looks good and I think that makes sense.
(2) Does not work in general, further, if I take your example
f = ode45(@(t,y)y, [0 -1], 1)
and your modified line 316 of the form
if (vdirection * (vtimestamp + vstepsize)> vdirection * vtimestop)
then the output on my machine is
error: Solving has not been successful. The iterative integration
loop exited at time t = -1.000000 before endpoint at tend = -1.000000
was reached. This happened because the iterative integration loop
does not find a valid solution at this time stamp. Try to reduce the
value of "InitialStep" and/or "MaxStep" with the command "odeset".
Did you use any options to make that work? You're right if you say
that the endpoint is not hidden correctly! Looks like we still need
another condition. Suggestions?
Best regards
Thomas
%# Copyright (C) 2006-2009, Thomas Treichl <trei...@users.sourceforge.net>
%# OdePkg - A package for solving ordinary differential equations and more
%#
%# This program is free software; you can redistribute it and/or modify
%# it under the terms of the GNU General Public License as published by
%# the Free Software Foundation; either version 2 of the License, or
%# (at your option) any later version.
%#
%# This program is distributed in the hope that it will be useful,
%# but WITHOUT ANY WARRANTY; without even the implied warranty of
%# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%# GNU General Public License for more details.
%#
%# You should have received a copy of the GNU General Public License
%# along with this program; If not, see <http://www.gnu.org/licenses/>.
%# -*- texinfo -*-
%# @deftypefn {Function File} {[@var{}] =} ode45 (@var{@@fun}, @var{slot},
@var{init}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
%# @deftypefnx {Command} {[@var{sol}] =} ode45 (@var{@@fun}, @var{slot},
@var{init}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
%# @deftypefnx {Command} {[@var{t}, @var{y}, [@var{xe}, @var{ye}, @var{ie}]] =}
ode45 (@var{@@fun}, @var{slot}, @var{init}, [@var{opt}], [@var{par1},
@var{par2}, @dots{}])
%#
%# This function file can be used to solve a set of non--stiff ordinary
differential equations (non--stiff ODEs) or non--stiff differential algebraic
equations (non--stiff DAEs) with the well known explicit Runge--Kutta method of
order (4,5).
%#
%# If this function is called with no return argument then plot the solution
over time in a figure window while solving the set of ODEs that are defined in
a function and specified by the function handle @var{@@fun}. The second input
argument @var{slot} is a double vector that defines the time slot, @var{init}
is a double vector that defines the initial values of the states, @var{opt} can
optionally be a structure array that keeps the options created with the command
@command{odeset} and @var{par1}, @var{par2}, @dots{} can optionally be other
input arguments of any type that have to be passed to the function defined by
@var{@@fun}.
%#
%# If this function is called with one return argument then return the solution
@var{sol} of type structure array after solving the set of ODEs. The solution
@var{sol} has the fields @var{x} of type double column vector for the steps
chosen by the solver, @var{y} of type double column vector for the solutions at
each time step of @var{x}, @var{solver} of type string for the solver name and
optionally the extended time stamp information @var{xe}, the extended solution
information @var{ye} and the extended index information @var{ie} all of type
double column vector that keep the informations of the event function if an
event function handle is set in the option argument @var{opt}.
%#
%# If this function is called with more than one return argument then return
the time stamps @var{t}, the solution values @var{y} and optionally the
extended time stamp information @var{xe}, the extended solution information
@var{ye} and the extended index information @var{ie} all of type double column
vector.
%#
%# For example, solve an anonymous implementation of the Van der Pol equation
%#
%# @example
%# fvdb = @@(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
%#
%# vopt = odeset ("RelTol", 1e-3, "AbsTol", 1e-3, \
%# "NormControl", "on", "OutputFcn", @@odeplot);
%# ode45 (fvdb, [0 20], [2 0], vopt);
%# @end example
%# @end deftypefn
%#
%# @seealso{odepkg}
%# ChangeLog:
%# 20010703 the function file "ode45.m" was written by Marc Compere
%# under the GPL for the use with this software. This function has been
%# taken as a base for the following implementation.
%# 20060810, Thomas Treichl
%# This function was adapted to the new syntax that is used by the
%# new OdePkg for Octave and is compatible to Matlab's ode45.
function [varargout] = ode45 (vfun, vslot, vinit, varargin)
if (nargin == 0) %# Check number and types of all input arguments
help ('ode45');
error ('OdePkg:InvalidArgument', ...
'Number of input arguments must be greater than zero');
elseif (nargin < 3)
print_usage;
elseif ~(isa (vfun, 'function_handle') || isa (vfun, 'inline'))
error ('OdePkg:InvalidArgument', ...
'First input argument must be a valid function handle');
elseif (~isvector (vslot) || length (vslot) < 2)
error ('OdePkg:InvalidArgument', ...
'Second input argument must be a valid vector');
elseif (~isvector (vinit) || ~isnumeric (vinit))
error ('OdePkg:InvalidArgument', ...
'Third input argument must be a valid numerical value');
elseif (nargin >= 4)
if (~isstruct (varargin{1}))
%# varargin{1:len} are parameters for vfun
vodeoptions = odeset;
vfunarguments = varargin;
elseif (length (varargin) > 1)
%# varargin{1} is an OdePkg options structure vopt
vodeoptions = odepkg_structure_check (varargin{1}, 'ode45');
vfunarguments = {varargin{2:length(varargin)}};
else %# if (isstruct (varargin{1}))
vodeoptions = odepkg_structure_check (varargin{1}, 'ode45');
vfunarguments = {};
end
else %# if (nargin == 3)
vodeoptions = odeset;
vfunarguments = {};
end
%# Start preprocessing, have a look which options are set in
%# vodeoptions, check if an invalid or unused option is set
vslot = vslot(:).'; %# Create a row vector
vinit = vinit(:).'; %# Create a row vector
if (length (vslot) > 2) %# Step size checking
vstepsizefixed = true;
else
vstepsizefixed = false;
end
%# Get the default options that can be set with 'odeset' temporarily
vodetemp = odeset;
%# Implementation of the option RelTol has been finished. This option
%# can be set by the user to another value than default value.
if (isempty (vodeoptions.RelTol) && ~vstepsizefixed)
vodeoptions.RelTol = 1e-6;
warning ('OdePkg:InvalidArgument', ...
'Option "RelTol" not set, new value %f is used', vodeoptions.RelTol);
elseif (~isempty (vodeoptions.RelTol) && vstepsizefixed)
warning ('OdePkg:InvalidArgument', ...
'Option "RelTol" will be ignored if fixed time stamps are given');
end
%# Implementation of the option AbsTol has been finished. This option
%# can be set by the user to another value than default value.
if (isempty (vodeoptions.AbsTol) && ~vstepsizefixed)
vodeoptions.AbsTol = 1e-6;
warning ('OdePkg:InvalidArgument', ...
'Option "AbsTol" not set, new value %f is used', vodeoptions.AbsTol);
elseif (~isempty (vodeoptions.AbsTol) && vstepsizefixed)
warning ('OdePkg:InvalidArgument', ...
'Option "AbsTol" will be ignored if fixed time stamps are given');
else
vodeoptions.AbsTol = vodeoptions.AbsTol(:); %# Create column vector
end
%# Implementation of the option NormControl has been finished. This
%# option can be set by the user to another value than default value.
if (strcmp (vodeoptions.NormControl, 'on')) vnormcontrol = true;
else vnormcontrol = false; end
%# Implementation of the option NonNegative has been finished. This
%# option can be set by the user to another value than default value.
if (~isempty (vodeoptions.NonNegative))
if (isempty (vodeoptions.Mass)), vhavenonnegative = true;
else
vhavenonnegative = false;
warning ('OdePkg:InvalidArgument', ...
'Option "NonNegative" will be ignored if mass matrix is set');
end
else vhavenonnegative = false;
end
%# Implementation of the option OutputFcn has been finished. This
%# option can be set by the user to another value than default value.
if (isempty (vodeoptions.OutputFcn) && nargout == 0)
vodeoptions.OutputFcn = @odeplot;
vhaveoutputfunction = true;
elseif (isempty (vodeoptions.OutputFcn)), vhaveoutputfunction = false;
else vhaveoutputfunction = true;
end
%# Implementation of the option OutputSel has been finished. This
%# option can be set by the user to another value than default value.
if (~isempty (vodeoptions.OutputSel)), vhaveoutputselection = true;
else vhaveoutputselection = false; end
%# Implementation of the option OutputSave has been finished. This
%# option can be set by the user to another value than default value.
if (isempty (vodeoptions.OutputSave)), vodeoptions.OutputSave = 1;
end
%# Implementation of the option Refine has been finished. This option
%# can be set by the user to another value than default value.
if (vodeoptions.Refine > 0), vhaverefine = true;
else vhaverefine = false; end
%# Implementation of the option Stats has been finished. This option
%# can be set by the user to another value than default value.
%# Implementation of the option InitialStep has been finished. This
%# option can be set by the user to another value than default value.
if (isempty (vodeoptions.InitialStep) && ~vstepsizefixed)
vodeoptions.InitialStep = (vslot(1,2) - vslot(1,1)) / 10;
vodeoptions.InitialStep = vodeoptions.InitialStep / 10^vodeoptions.Refine;
warning ('OdePkg:InvalidArgument', ...
'Option "InitialStep" not set, new value %f is used',
vodeoptions.InitialStep);
end
%# Implementation of the option MaxStep has been finished. This option
%# can be set by the user to another value than default value.
if (isempty (vodeoptions.MaxStep) && ~vstepsizefixed)
vodeoptions.MaxStep = (vslot(1,2) - vslot(1,1)) / 10;
warning ('OdePkg:InvalidArgument', ...
'Option "MaxStep" not set, new value %f is used', vodeoptions.MaxStep);
end
%# Implementation of the option Events has been finished. This option
%# can be set by the user to another value than default value.
if (~isempty (vodeoptions.Events)), vhaveeventfunction = true;
else vhaveeventfunction = false; end
%# The options 'Jacobian', 'JPattern' and 'Vectorized' will be ignored
%# by this solver because this solver uses an explicit Runge-Kutta
%# method and therefore no Jacobian calculation is necessary
if (~isequal (vodeoptions.Jacobian, vodetemp.Jacobian))
warning ('OdePkg:InvalidArgument', ...
'Option "Jacobian" will be ignored by this solver');
end
if (~isequal (vodeoptions.JPattern, vodetemp.JPattern))
warning ('OdePkg:InvalidArgument', ...
'Option "JPattern" will be ignored by this solver');
end
if (~isequal (vodeoptions.Vectorized, vodetemp.Vectorized))
warning ('OdePkg:InvalidArgument', ...
'Option "Vectorized" will be ignored by this solver');
end
if (~isequal (vodeoptions.NewtonTol, vodetemp.NewtonTol))
warning ('OdePkg:InvalidArgument', ...
'Option "NewtonTol" will be ignored by this solver');
end
if (~isequal (vodeoptions.MaxNewtonIterations,...
vodetemp.MaxNewtonIterations))
warning ('OdePkg:InvalidArgument', ...
'Option "MaxNewtonIterations" will be ignored by this solver');
end
%# Implementation of the option Mass has been finished. This option
%# can be set by the user to another value than default value.
if (~isempty (vodeoptions.Mass) && ismatrix (vodeoptions.Mass))
vhavemasshandle = false; vmass = vodeoptions.Mass; %# constant mass
elseif (isa (vodeoptions.Mass, 'function_handle'))
vhavemasshandle = true; %# mass defined by a function handle
else %# no mass matrix - creating a diag-matrix of ones for mass
vhavemasshandle = false; %# vmass = diag (ones (length (vinit), 1), 0);
end
%# Implementation of the option MStateDependence has been finished.
%# This option can be set by the user to another value than default
%# value.
if (strcmp (vodeoptions.MStateDependence, 'none'))
vmassdependence = false;
else vmassdependence = true;
end
%# Other options that are not used by this solver. Print a warning
%# message to tell the user that the option(s) is/are ignored.
if (~isequal (vodeoptions.MvPattern, vodetemp.MvPattern))
warning ('OdePkg:InvalidArgument', ...
'Option "MvPattern" will be ignored by this solver');
end
if (~isequal (vodeoptions.MassSingular, vodetemp.MassSingular))
warning ('OdePkg:InvalidArgument', ...
'Option "MassSingular" will be ignored by this solver');
end
if (~isequal (vodeoptions.InitialSlope, vodetemp.InitialSlope))
warning ('OdePkg:InvalidArgument', ...
'Option "InitialSlope" will be ignored by this solver');
end
if (~isequal (vodeoptions.MaxOrder, vodetemp.MaxOrder))
warning ('OdePkg:InvalidArgument', ...
'Option "MaxOrder" will be ignored by this solver');
end
if (~isequal (vodeoptions.BDF, vodetemp.BDF))
warning ('OdePkg:InvalidArgument', ...
'Option "BDF" will be ignored by this solver');
end
%# Starting the initialisation of the core solver ode45
vtimestamp = vslot(1,1); %# timestamp = start time
vtimelength = length (vslot); %# length needed if fixed steps
vtimestop = vslot(1,vtimelength); %# stop time = last value
%# 20110611, reported by Nils Strunk
%# Make it possible to solve equations from negativ to zero,
%# eg. vres = ode45 (@(t,y) y, [-2 0], 2)
vdirection = sign (vtimestop - vtimestamp); %# Direction flag
if (~vstepsizefixed)
if (sign (vodeoptions.InitialStep) == vdirection)
vstepsize = vodeoptions.InitialStep;
else %# Fix wrong direction of InitialStep.
vstepsize = - vodeoptions.InitialStep;
end
vminstepsize = (vtimestop - vtimestamp) / (1/eps);
else %# If step size is given then use the fixed time steps
vstepsize = vslot(1,2) - vslot(1,1);
vminstepsize = sign (vstepsize) * eps;
end
vretvaltime = vtimestamp; %# first timestamp output
vretvalresult = vinit; %# first solution output
%# Initialize the OutputFcn
if (vhaveoutputfunction)
if (vhaveoutputselection) vretout = vretvalresult(vodeoptions.OutputSel);
else vretout = vretvalresult; end
feval (vodeoptions.OutputFcn, vslot.', ...
vretout.', 'init', vfunarguments{:});
end
%# Initialize the EventFcn
if (vhaveeventfunction)
odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
vretvalresult.', 'init', vfunarguments{:});
end
vpow = 1/5; %# 20071016, reported by Luis Randez
va = [0, 0, 0, 0, 0; %# The Runge-Kutta-Fehlberg 4(5) coefficients
1/4, 0, 0, 0, 0; %# Coefficients proved on 20060827
3/32, 9/32, 0, 0, 0; %# See p.91 in Ascher & Petzold
1932/2197, -7200/2197, 7296/2197, 0, 0;
439/216, -8, 3680/513, -845/4104, 0;
-8/27, 2, -3544/2565, 1859/4104, -11/40];
%# 4th and 5th order b-coefficients
vb4 = [25/216; 0; 1408/2565; 2197/4104; -1/5; 0];
vb5 = [16/135; 0; 6656/12825; 28561/56430; -9/50; 2/55];
vc = sum (va, 2);
%# The solver main loop - stop if the endpoint has been reached
vcntloop = 2; vcntcycles = 1; vu = vinit; vk = vu.' * zeros(1,6);
vcntiter = 0; vunhandledtermination = true; vcntsave = 2;
while ((vdirection * (vtimestamp) < vdirection * (vtimestop)) && ...
(vdirection * (vstepsize) >= vdirection * (vminstepsize)))
%# Hit the endpoint of the time slot exactely
if ((vtimestamp + vstepsize) > vdirection * vtimestop)
%# vstepsize = vtimestop - vdirection * vtimestamp;
%# 20110611, reported by Nils Strunk
%# The endpoint of the time slot must be hit exactly,
%# vres = ode45 (@(t,y) y, [0 -1], 1);
vstepsize = vdirection * abs (abs (vtimestop) - abs (vtimestamp));
end
%# Estimate the six results when using this solver
for j = 1:6
vthetime = vtimestamp + vc(j,1) * vstepsize;
vtheinput = vu.' + vstepsize * vk(:,1:j-1) * va(j,1:j-1).';
if (vhavemasshandle) %# Handle only the dynamic mass matrix,
if (vmassdependence) %# constant mass matrices have already
vmass = feval ... %# been set before (if any)
(vodeoptions.Mass, vthetime, vtheinput, vfunarguments{:});
else %# if (vmassdependence == false)
vmass = feval ... %# then we only have the time argument
(vodeoptions.Mass, vthetime, vfunarguments{:});
end
vk(:,j) = vmass \ feval ...
(vfun, vthetime, vtheinput, vfunarguments{:});
else
vk(:,j) = feval ...
(vfun, vthetime, vtheinput, vfunarguments{:});
end
end
%# Compute the 4th and the 5th order estimation
y4 = vu.' + vstepsize * (vk * vb4);
y5 = vu.' + vstepsize * (vk * vb5);
if (vhavenonnegative)
vu(vodeoptions.NonNegative) = abs (vu(vodeoptions.NonNegative));
y4(vodeoptions.NonNegative) = abs (y4(vodeoptions.NonNegative));
y5(vodeoptions.NonNegative) = abs (y5(vodeoptions.NonNegative));
end
if (vhaveoutputfunction && vhaverefine)
vSaveVUForRefine = vu;
end
%# Calculate the absolute local truncation error and the acceptable error
if (~vstepsizefixed)
if (~vnormcontrol)
vdelta = abs (y5 - y4);
vtau = max (vodeoptions.RelTol * abs (vu.'), vodeoptions.AbsTol);
else
vdelta = norm (y5 - y4, Inf);
vtau = max (vodeoptions.RelTol * max (norm (vu.', Inf), 1.0), ...
vodeoptions.AbsTol);
end
else %# if (vstepsizefixed == true)
vdelta = 1; vtau = 2;
end
%# If the error is acceptable then update the vretval variables
if (all (vdelta <= vtau))
vtimestamp = vtimestamp + vstepsize;
vu = y5.'; %# MC2001: the higher order estimation as "local extrapolation"
%# Save the solution every vodeoptions.OutputSave steps
if (mod (vcntloop-1,vodeoptions.OutputSave) == 0)
vretvaltime(vcntsave,:) = vtimestamp;
vretvalresult(vcntsave,:) = vu;
vcntsave = vcntsave + 1;
end
vcntloop = vcntloop + 1; vcntiter = 0;
%# Call plot only if a valid result has been found, therefore this
%# code fragment has moved here. Stop integration if plot function
%# returns false
if (vhaveoutputfunction)
for vcnt = 0:vodeoptions.Refine %# Approximation between told and t
if (vhaverefine) %# Do interpolation
vapproxtime = (vcnt + 1) * vstepsize / (vodeoptions.Refine + 2);
vapproxvals = vSaveVUForRefine.' + vapproxtime * (vk * vb5);
vapproxtime = (vtimestamp - vstepsize) + vapproxtime;
else
vapproxvals = vu.';
vapproxtime = vtimestamp;
end
if (vhaveoutputselection)
vapproxvals = vapproxvals(vodeoptions.OutputSel);
end
vpltret = feval (vodeoptions.OutputFcn, vapproxtime, ...
vapproxvals, [], vfunarguments{:});
if vpltret %# Leave refinement loop
break;
end
end
if (vpltret) %# Leave main loop
vunhandledtermination = false;
break;
end
end
%# Call event only if a valid result has been found, therefore this
%# code fragment has moved here. Stop integration if veventbreak is
%# true
if (vhaveeventfunction)
vevent = ...
odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
vu(:), [], vfunarguments{:});
if (~isempty (vevent{1}) && vevent{1} == 1)
vretvaltime(vcntloop-1,:) = vevent{3}(end,:);
vretvalresult(vcntloop-1,:) = vevent{4}(end,:);
vunhandledtermination = false; break;
end
end
end %# If the error is acceptable ...
%# Update the step size for the next integration step
if (~vstepsizefixed)
%# 20080425, reported by Marco Caliari
%# vdelta cannot be negative (because of the absolute value that
%# has been introduced) but it could be 0, then replace the zeros
%# with the maximum value of vdelta
vdelta(find (vdelta == 0)) = max (vdelta);
%# It could happen that max (vdelta) == 0 (ie. that the original
%# vdelta was 0), in that case we double the previous vstepsize
vdelta(find (vdelta == 0)) = max (vtau) .* (0.4 ^ (1 / vpow));
if (vdirection == 1)
vstepsize = min (vodeoptions.MaxStep, ...
min (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow));
else
vstepsize = max (vodeoptions.MaxStep, ...
max (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow));
end
else %# if (vstepsizefixed)
if (vcntloop <= vtimelength)
vstepsize = vslot(vcntloop) - vslot(vcntloop-1);
else %# Get out of the main integration loop
break;
end
end
%# Update counters that count the number of iteration cycles
vcntcycles = vcntcycles + 1; %# Needed for cost statistics
vcntiter = vcntiter + 1; %# Needed to find iteration problems
%# Stop solving because the last 1000 steps no successful valid
%# value has been found
if (vcntiter >= 5000)
error (['Solving has not been successful. The iterative', ...
' integration loop exited at time t = %f before endpoint at', ...
' tend = %f was reached. This happened because the iterative', ...
' integration loop does not find a valid solution at this time', ...
' stamp. Try to reduce the value of "InitialStep" and/or', ...
' "MaxStep" with the command "odeset".\n'], vtimestamp, vtimestop);
end
end %# The main loop
%# Check if integration of the ode has been successful
if (vdirection * vtimestamp < vdirection * vtimestop)
if (vunhandledtermination == true)
error ('OdePkg:InvalidArgument', ...
['Solving has not been successful. The iterative', ...
' integration loop exited at time t = %f', ...
' before endpoint at tend = %f was reached. This may', ...
' happen if the stepsize grows smaller than defined in', ...
' vminstepsize. Try to reduce the value of "InitialStep" and/or', ...
' "MaxStep" with the command "odeset".\n'], vtimestamp, vtimestop);
else
warning ('OdePkg:InvalidArgument', ...
['Solver has been stopped by a call of "break" in', ...
' the main iteration loop at time t = %f before endpoint at', ...
' tend = %f was reached. This may happen because the @odeplot', ...
' function returned "true" or the @event function returned "true".'],
...
vtimestamp, vtimestop);
end
end
%# Postprocessing, do whatever when terminating integration algorithm
if (vhaveoutputfunction) %# Cleanup plotter
feval (vodeoptions.OutputFcn, vtimestamp, ...
vu.', 'done', vfunarguments{:});
end
if (vhaveeventfunction) %# Cleanup event function handling
odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
vu.', 'done', vfunarguments{:});
end
%# Save the last step, if not already saved
if (mod (vcntloop-2,vodeoptions.OutputSave) ~= 0)
vretvaltime(vcntsave,:) = vtimestamp;
vretvalresult(vcntsave,:) = vu;
end
%# Print additional information if option Stats is set
if (strcmp (vodeoptions.Stats, 'on'))
vhavestats = true;
vnsteps = vcntloop-2; %# vcntloop from 2..end
vnfailed = (vcntcycles-1)-(vcntloop-2)+1; %# vcntcycl from 1..end
vnfevals = 6*(vcntcycles-1); %# number of ode evaluations
vndecomps = 0; %# number of LU decompositions
vnpds = 0; %# number of partial derivatives
vnlinsols = 0; %# no. of solutions of linear
systems
%# Print cost statistics if no output argument is given
if (nargout == 0)
vmsg = fprintf (1, 'Number of successful steps: %d\n', vnsteps);
vmsg = fprintf (1, 'Number of failed attempts: %d\n', vnfailed);
vmsg = fprintf (1, 'Number of function calls: %d\n', vnfevals);
end
else
vhavestats = false;
end
if (nargout == 1) %# Sort output variables, depends on nargout
varargout{1}.x = vretvaltime; %# Time stamps are saved in field x
varargout{1}.y = vretvalresult; %# Results are saved in field y
varargout{1}.solver = 'ode45'; %# Solver name is saved in field solver
if (vhaveeventfunction)
varargout{1}.ie = vevent{2}; %# Index info which event occured
varargout{1}.xe = vevent{3}; %# Time info when an event occured
varargout{1}.ye = vevent{4}; %# Results when an event occured
end
if (vhavestats)
varargout{1}.stats = struct;
varargout{1}.stats.nsteps = vnsteps;
varargout{1}.stats.nfailed = vnfailed;
varargout{1}.stats.nfevals = vnfevals;
varargout{1}.stats.npds = vnpds;
varargout{1}.stats.ndecomps = vndecomps;
varargout{1}.stats.nlinsols = vnlinsols;
end
elseif (nargout == 2)
varargout{1} = vretvaltime; %# Time stamps are first output argument
varargout{2} = vretvalresult; %# Results are second output argument
elseif (nargout == 5)
varargout{1} = vretvaltime; %# Same as (nargout == 2)
varargout{2} = vretvalresult; %# Same as (nargout == 2)
varargout{3} = []; %# LabMat doesn't accept lines like
varargout{4} = []; %# varargout{3} = varargout{4} = [];
varargout{5} = [];
if (vhaveeventfunction)
varargout{3} = vevent{3}; %# Time info when an event occured
varargout{4} = vevent{4}; %# Results when an event occured
varargout{5} = vevent{2}; %# Index info which event occured
end
end
end
%! # We are using the "Van der Pol" implementation for all tests that
%! # are done for this function. We also define a Jacobian, Events,
%! # pseudo-Mass implementation. For further tests we also define a
%! # reference solution (computed at high accuracy) and an OutputFcn
%!function [ydot] = fpol (vt, vy, varargin) %# The Van der Pol
%! ydot = [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
%!function [vjac] = fjac (vt, vy, varargin) %# its Jacobian
%! vjac = [0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2];
%!function [vjac] = fjcc (vt, vy, varargin) %# sparse type
%! vjac = sparse ([0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2]);
%!function [vval, vtrm, vdir] = feve (vt, vy, varargin)
%! vval = fpol (vt, vy, varargin); %# We use the derivatives
%! vtrm = zeros (2,1); %# that's why component 2
%! vdir = ones (2,1); %# seems to not be exact
%!function [vval, vtrm, vdir] = fevn (vt, vy, varargin)
%! vval = fpol (vt, vy, varargin); %# We use the derivatives
%! vtrm = ones (2,1); %# that's why component 2
%! vdir = ones (2,1); %# seems to not be exact
%!function [vmas] = fmas (vt, vy)
%! vmas = [1, 0; 0, 1]; %# Dummy mass matrix for tests
%!function [vmas] = fmsa (vt, vy)
%! vmas = sparse ([1, 0; 0, 1]); %# A sparse dummy matrix
%!function [vref] = fref () %# The computed reference sol
%! vref = [0.32331666704577, -1.83297456798624];
%!function [vout] = fout (vt, vy, vflag, varargin)
%! if (regexp (char (vflag), 'init') == 1)
%! if (any (size (vt) ~= [2, 1])) error ('"fout" step "init"'); end
%! elseif (isempty (vflag))
%! if (any (size (vt) ~= [1, 1])) error ('"fout" step "calc"'); end
%! vout = false;
%! elseif (regexp (char (vflag), 'done') == 1)
%! if (any (size (vt) ~= [1, 1])) error ('"fout" step "done"'); end
%! else error ('"fout" invalid vflag');
%! end
%!
%! %# Turn off output of warning messages for all tests, turn them on
%! %# again if the last test is called
%!error %# input argument number one
%! warning ('off', 'OdePkg:InvalidArgument');
%! B = ode45 (1, [0 25], [3 15 1]);
%!error %# input argument number two
%! B = ode45 (@fpol, 1, [3 15 1]);
%!error %# input argument number three
%! B = ode45 (@flor, [0 25], 1);
%!test %# one output argument
%! vsol = ode45 (@fpol, [0 2], [2 0]);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%! assert (isfield (vsol, 'solver'));
%! assert (vsol.solver, 'ode45');
%!test %# two output arguments
%! [vt, vy] = ode45 (@fpol, [0 2], [2 0]);
%! assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
%!test %# five output arguments and no Events
%! [vt, vy, vxe, vye, vie] = ode45 (@fpol, [0 2], [2 0]);
%! assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
%! assert ([vie, vxe, vye], []);
%!test %# anonymous function instead of real function
%! fvdb = @(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
%! vsol = ode45 (fvdb, [0 2], [2 0]);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# extra input arguments passed trhough
%! vsol = ode45 (@fpol, [0 2], [2 0], 12, 13, 'KL');
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# empty OdePkg structure *but* extra input arguments
%! vopt = odeset;
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt, 12, 13, 'KL');
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!error %# strange OdePkg structure
%! vopt = struct ('foo', 1);
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%!test %# Solve vdp in fixed step sizes
%! vsol = ode45 (@fpol, [0:0.1:2], [2 0]);
%! assert (vsol.x(:), [0:0.1:2]');
%! assert (vsol.y(end,:), fref, 1e-3);
%!test %# Solve in backward direction starting at t=0
%! %# vref = [-1.2054034414, 0.9514292694];
%! vsol = ode45 (@fpol, [0 -2], [2 0]);
%! %# assert ([vsol.x(end), vsol.y(end,:)], [-2, fref], 1e-3);
%!test %# Solve in backward direction starting at t=2
%! %# vref = [-1.2154183302, 0.9433018000];
%! vsol = ode45 (@fpol, [2 -2], [0.3233166627 -1.8329746843]);
%! %# assert ([vsol.x(end), vsol.y(end,:)], [-2, fref], 1e-3);
%!test %# AbsTol option
%! vopt = odeset ('AbsTol', 1e-5);
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# AbsTol and RelTol option
%! vopt = odeset ('AbsTol', 1e-8, 'RelTol', 1e-8);
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# RelTol and NormControl option -- higher accuracy
%! vopt = odeset ('RelTol', 1e-8, 'NormControl', 'on');
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-5);
%!test %# Keeps initial values while integrating
%! vopt = odeset ('NonNegative', 2);
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, 2, 0], 0.5);
%!test %# Details of OutputSel and Refine can't be tested
%! vopt = odeset ('OutputFcn', @fout, 'OutputSel', 1, 'Refine', 5);
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%!test %# Details of OutputSave can't be tested
%! vopt = odeset ('OutputSave', 1, 'OutputSel', 1);
%! vsla = ode45 (@fpol, [0 2], [2 0], vopt);
%! vopt = odeset ('OutputSave', 2);
%! vslb = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert (length (vsla.x) > length (vslb.x))
%!test %# Stats must add further elements in vsol
%! vopt = odeset ('Stats', 'on');
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert (isfield (vsol, 'stats'));
%! assert (isfield (vsol.stats, 'nsteps'));
%!test %# InitialStep option
%! vopt = odeset ('InitialStep', 1e-8);
%! vsol = ode45 (@fpol, [0 0.2], [2 0], vopt);
%! assert ([vsol.x(2)-vsol.x(1)], [1e-8], 1e-9);
%!test %# MaxStep option
%! vopt = odeset ('MaxStep', 1e-2);
%! vsol = ode45 (@fpol, [0 0.2], [2 0], vopt);
%! assert ([vsol.x(5)-vsol.x(4)], [1e-2], 1e-3);
%!test %# Events option add further elements in vsol
%! vopt = odeset ('Events', @feve);
%! vsol = ode45 (@fpol, [0 10], [2 0], vopt);
%! assert (isfield (vsol, 'ie'));
%! assert (vsol.ie(1), 2);
%! assert (isfield (vsol, 'xe'));
%! assert (isfield (vsol, 'ye'));
%!test %# Events option, now stop integration
%! vopt = odeset ('Events', @fevn, 'NormControl', 'on');
%! vsol = ode45 (@fpol, [0 10], [2 0], vopt);
%! assert ([vsol.ie, vsol.xe, vsol.ye], ...
%! [2.0, 2.496110, -0.830550, -2.677589], .5e-1);
%!test %# Events option, five output arguments
%! vopt = odeset ('Events', @fevn, 'NormControl', 'on');
%! [vt, vy, vxe, vye, vie] = ode45 (@fpol, [0 10], [2 0], vopt);
%! assert ([vie, vxe, vye], ...
%! [2.0, 2.496110, -0.830550, -2.677589], 1e-1);
%!test %# Jacobian option
%! vopt = odeset ('Jacobian', @fjac);
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# Jacobian option and sparse return value
%! vopt = odeset ('Jacobian', @fjcc);
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!
%! %# test for JPattern option is missing
%! %# test for Vectorized option is missing
%! %# test for NewtonTol option is missing
%! %# test for MaxNewtonIterations option is missing
%!
%!test %# Mass option as function
%! vopt = odeset ('Mass', @fmas);
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# Mass option as matrix
%! vopt = odeset ('Mass', eye (2,2));
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# Mass option as sparse matrix
%! vopt = odeset ('Mass', sparse (eye (2,2)));
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# Mass option as function and sparse matrix
%! vopt = odeset ('Mass', @fmsa);
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# Mass option as function and MStateDependence
%! vopt = odeset ('Mass', @fmas, 'MStateDependence', 'strong');
%! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
%!test %# Set BDF option to something else than default
%! vopt = odeset ('BDF', 'on');
%! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt);
%! assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
%!
%! %# test for MvPattern option is missing
%! %# test for InitialSlope option is missing
%! %# test for MaxOrder option is missing
%!
%! warning ('on', 'OdePkg:InvalidArgument');
%# Local Variables: ***
%# mode: octave ***
%# End: ***
------------------------------------------------------------------------------
EditLive Enterprise is the world's most technically advanced content
authoring tool. Experience the power of Track Changes, Inline Image
Editing and ensure content is compliant with Accessibility Checking.
http://p.sf.net/sfu/ephox-dev2dev
_______________________________________________
Octave-dev mailing list
Octave-dev@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/octave-dev
