PageRenderTime 52ms CodeModel.GetById 19ms RepoModel.GetById 1ms app.codeStats 0ms

/tags/R2009-06-07/main/odepkg/inst/ode23d.m

#
MATLAB | 716 lines | 424 code | 45 blank | 247 comment | 50 complexity | 8874bc7f2038f77a9d575c07c4af8399 MD5 | raw file
Possible License(s): GPL-2.0, BSD-3-Clause, LGPL-2.1, GPL-3.0, LGPL-3.0
  1. %# Copyright (C) 2008, Thomas Treichl <treichl@users.sourceforge.net>
  2. %# OdePkg - A package for solving ordinary differential equations and more
  3. %#
  4. %# This program is free software; you can redistribute it and/or modify
  5. %# it under the terms of the GNU General Public License as published by
  6. %# the Free Software Foundation; either version 2 of the License, or
  7. %# (at your option) any later version.
  8. %#
  9. %# This program is distributed in the hope that it will be useful,
  10. %# but WITHOUT ANY WARRANTY; without even the implied warranty of
  11. %# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  12. %# GNU General Public License for more details.
  13. %#
  14. %# You should have received a copy of the GNU General Public License
  15. %# along with this program; If not, see <http://www.gnu.org/licenses/>.
  16. %# -*- texinfo -*-
  17. %# @deftypefn {Function File} {[@var{}] =} ode23d (@var{@@fun}, @var{slot}, @var{init}, @var{lags}, @var{hist}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
  18. %# @deftypefnx {Command} {[@var{sol}] =} ode23d (@var{@@fun}, @var{slot}, @var{init}, @var{lags}, @var{hist}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
  19. %# @deftypefnx {Command} {[@var{t}, @var{y}, [@var{xe}, @var{ye}, @var{ie}]] =} ode23d (@var{@@fun}, @var{slot}, @var{init}, @var{lags}, @var{hist}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
  20. %#
  21. %# This function file can be used to solve a set of non--stiff delay differential equations (non--stiff DDEs) with a modified version of the well known explicit Runge--Kutta method of order (2,3).
  22. %#
  23. %# If this function is called with no return argument then plot the solution over time in a figure window while solving the set of DDEs 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{lags} is a double vector that describes the lags of time, @var{hist} is a double matrix and describes the history of the DDEs, @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}.
  24. %#
  25. %# If this function is called with one return argument then return the solution @var{sol} of type structure array after solving the set of DDEs. 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}.
  26. %#
  27. %# 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.
  28. %#
  29. %# For example, solve an anonymous implementation of a chaotic behavior
  30. %# @example
  31. %# fcao = @@(vt, vy, vz) [2 * vz / (1 + vz^9.65) - vy];
  32. %#
  33. %# vopt = odeset ("NormControl", "on", "RelTol", 1e-3);
  34. %# vsol = ode23d (fcao, [0, 100], 0.5, 2, 0.5, vopt);
  35. %#
  36. %# vlag = interp1 (vsol.x, vsol.y, vsol.x - 2);
  37. %# plot (vsol.y, vlag); legend ("fcao (t,y,z)");
  38. %# @end example
  39. %# @end deftypefn
  40. %#
  41. %# @seealso{odepkg}
  42. function [varargout] = ode23d (vfun, vslot, vinit, vlags, vhist, varargin)
  43. if (nargin == 0) %# Check number and types of all input arguments
  44. help ('ode23d');
  45. error ('OdePkg:InvalidArgument', ...
  46. 'Number of input arguments must be greater than zero');
  47. elseif (nargin < 5)
  48. print_usage;
  49. elseif (~isa (vfun, 'function_handle'))
  50. error ('OdePkg:InvalidArgument', ...
  51. 'First input argument must be a valid function handle');
  52. elseif (~isvector (vslot) || length (vslot) < 2)
  53. error ('OdePkg:InvalidArgument', ...
  54. 'Second input argument must be a valid vector');
  55. elseif (~isvector (vinit) || ~isnumeric (vinit))
  56. error ('OdePkg:InvalidArgument', ...
  57. 'Third input argument must be a valid numerical value');
  58. elseif (~isvector (vlags) || ~isnumeric (vlags))
  59. error ('OdePkg:InvalidArgument', ...
  60. 'Fourth input argument must be a valid numerical value');
  61. elseif ~(isnumeric (vhist) || isa (vhist, 'function_handle'))
  62. error ('OdePkg:InvalidArgument', ...
  63. 'Fifth input argument must either be numeric or a function handle');
  64. elseif (nargin >= 6)
  65. if (~isstruct (varargin{1}))
  66. %# varargin{1:len} are parameters for vfun
  67. vodeoptions = odeset;
  68. vfunarguments = varargin;
  69. elseif (length (varargin) > 1)
  70. %# varargin{1} is an OdePkg options structure vopt
  71. vodeoptions = odepkg_structure_check (varargin{1}, 'ode23d');
  72. vfunarguments = {varargin{2:length(varargin)}};
  73. else %# if (isstruct (varargin{1}))
  74. vodeoptions = odepkg_structure_check (varargin{1}, 'ode23d');
  75. vfunarguments = {};
  76. end
  77. else %# if (nargin == 5)
  78. vodeoptions = odeset;
  79. vfunarguments = {};
  80. end
  81. %# Start preprocessing, have a look which options have been set in
  82. %# vodeoptions. Check if an invalid or unused option has been set and
  83. %# print warnings.
  84. vslot = vslot(:)'; %# Create a row vector
  85. vinit = vinit(:)'; %# Create a row vector
  86. vlags = vlags(:)'; %# Create a row vector
  87. %# Check if the user has given fixed points of time
  88. if (length (vslot) > 2), vstepsizegiven = true; %# Step size checking
  89. else vstepsizegiven = false; end
  90. %# Get the default options that can be set with 'odeset' temporarily
  91. vodetemp = odeset;
  92. %# Implementation of the option RelTol has been finished. This option
  93. %# can be set by the user to another value than default value.
  94. if (isempty (vodeoptions.RelTol) && ~vstepsizegiven)
  95. vodeoptions.RelTol = 1e-6;
  96. warning ('OdePkg:InvalidOption', ...
  97. 'Option "RelTol" not set, new value %f is used', vodeoptions.RelTol);
  98. elseif (~isempty (vodeoptions.RelTol) && vstepsizegiven)
  99. warning ('OdePkg:InvalidOption', ...
  100. 'Option "RelTol" will be ignored if fixed time stamps are given');
  101. %# This implementation has been added to odepkg_structure_check.m
  102. %# elseif (~isscalar (vodeoptions.RelTol) && ~vstepsizegiven)
  103. %# error ('OdePkg:InvalidOption', ...
  104. %# 'Option "RelTol" must be set to a scalar value for this solver');
  105. end
  106. %# Implementation of the option AbsTol has been finished. This option
  107. %# can be set by the user to another value than default value.
  108. if (isempty (vodeoptions.AbsTol) && ~vstepsizegiven)
  109. vodeoptions.AbsTol = 1e-6;
  110. warning ('OdePkg:InvalidOption', ...
  111. 'Option "AbsTol" not set, new value %f is used', vodeoptions.AbsTol);
  112. elseif (~isempty (vodeoptions.AbsTol) && vstepsizegiven)
  113. warning ('OdePkg:InvalidOption', ...
  114. 'Option "AbsTol" will be ignored if fixed time stamps are given');
  115. else %# create column vector
  116. vodeoptions.AbsTol = vodeoptions.AbsTol(:);
  117. end
  118. %# Implementation of the option NormControl has been finished. This
  119. %# option can be set by the user to another value than default value.
  120. if (strcmp (vodeoptions.NormControl, 'on')), vnormcontrol = true;
  121. else vnormcontrol = false;
  122. end
  123. %# Implementation of the option NonNegative has been finished. This
  124. %# option can be set by the user to another value than default value.
  125. if (~isempty (vodeoptions.NonNegative))
  126. if (isempty (vodeoptions.Mass)), vhavenonnegative = true;
  127. else
  128. vhavenonnegative = false;
  129. warning ('OdePkg:InvalidOption', ...
  130. 'Option "NonNegative" will be ignored if mass matrix is set');
  131. end
  132. else vhavenonnegative = false;
  133. end
  134. %# Implementation of the option OutputFcn has been finished. This
  135. %# option can be set by the user to another value than default value.
  136. if (isempty (vodeoptions.OutputFcn) && nargout == 0)
  137. vodeoptions.OutputFcn = @odeplot;
  138. vhaveoutputfunction = true;
  139. elseif (isempty (vodeoptions.OutputFcn)), vhaveoutputfunction = false;
  140. else vhaveoutputfunction = true;
  141. end
  142. %# Implementation of the option OutputSel has been finished. This
  143. %# option can be set by the user to another value than default value.
  144. if (~isempty (vodeoptions.OutputSel)), vhaveoutputselection = true;
  145. else vhaveoutputselection = false; end
  146. %# Implementation of the option Refine has been finished. This option
  147. %# can be set by the user to another value than default value.
  148. if (isequal (vodeoptions.Refine, vodetemp.Refine)), vhaverefine = true;
  149. else vhaverefine = false; end
  150. %# Implementation of the option Stats has been finished. This option
  151. %# can be set by the user to another value than default value.
  152. %# Implementation of the option InitialStep has been finished. This
  153. %# option can be set by the user to another value than default value.
  154. if (isempty (vodeoptions.InitialStep) && ~vstepsizegiven)
  155. vodeoptions.InitialStep = abs (vslot(1,1) - vslot(1,2)) / 10;
  156. vodeoptions.InitialStep = vodeoptions.InitialStep / 10^vodeoptions.Refine;
  157. warning ('OdePkg:InvalidOption', ...
  158. 'Option "InitialStep" not set, new value %f is used', vodeoptions.InitialStep);
  159. end
  160. %# Implementation of the option MaxStep has been finished. This option
  161. %# can be set by the user to another value than default value.
  162. if (isempty (vodeoptions.MaxStep) && ~vstepsizegiven)
  163. vodeoptions.MaxStep = abs (vslot(1,1) - vslot(1,length (vslot))) / 10;
  164. %# vodeoptions.MaxStep = vodeoptions.MaxStep / 10^vodeoptions.Refine;
  165. warning ('OdePkg:InvalidOption', ...
  166. 'Option "MaxStep" not set, new value %f is used', vodeoptions.MaxStep);
  167. end
  168. %# Implementation of the option Events has been finished. This option
  169. %# can be set by the user to another value than default value.
  170. if (~isempty (vodeoptions.Events)), vhaveeventfunction = true;
  171. else vhaveeventfunction = false; end
  172. %# The options 'Jacobian', 'JPattern' and 'Vectorized' will be ignored
  173. %# by this solver because this solver uses an explicit Runge-Kutta
  174. %# method and therefore no Jacobian calculation is necessary
  175. if (~isequal (vodeoptions.Jacobian, vodetemp.Jacobian))
  176. warning ('OdePkg:InvalidOption', ...
  177. 'Option "Jacobian" will be ignored by this solver');
  178. end
  179. if (~isequal (vodeoptions.JPattern, vodetemp.JPattern))
  180. warning ('OdePkg:InvalidOption', ...
  181. 'Option "JPattern" will be ignored by this solver');
  182. end
  183. if (~isequal (vodeoptions.Vectorized, vodetemp.Vectorized))
  184. warning ('OdePkg:InvalidOption', ...
  185. 'Option "Vectorized" will be ignored by this solver');
  186. end
  187. %# Implementation of the option Mass has been finished. This option
  188. %# can be set by the user to another value than default value.
  189. if (~isempty (vodeoptions.Mass) && ismatrix (vodeoptions.Mass))
  190. vhavemasshandle = false; vmass = vodeoptions.Mass; %# constant mass
  191. elseif (isa (vodeoptions.Mass, 'function_handle'))
  192. vhavemasshandle = true; %# mass defined by a function handle
  193. else %# no mass matrix - creating a diag-matrix of ones for mass
  194. vhavemasshandle = false; %# vmass = diag (ones (length (vinit), 1), 0);
  195. end
  196. %# Implementation of the option MStateDependence has been finished.
  197. %# This option can be set by the user to another value than default
  198. %# value.
  199. if (strcmp (vodeoptions.MStateDependence, 'none'))
  200. vmassdependence = false;
  201. else vmassdependence = true;
  202. end
  203. %# Other options that are not used by this solver. Print a warning
  204. %# message to tell the user that the option(s) is/are ignored.
  205. if (~isequal (vodeoptions.MvPattern, vodetemp.MvPattern))
  206. warning ('OdePkg:InvalidOption', ...
  207. 'Option "MvPattern" will be ignored by this solver');
  208. end
  209. if (~isequal (vodeoptions.MassSingular, vodetemp.MassSingular))
  210. warning ('OdePkg:InvalidOption', ...
  211. 'Option "MassSingular" will be ignored by this solver');
  212. end
  213. if (~isequal (vodeoptions.InitialSlope, vodetemp.InitialSlope))
  214. warning ('OdePkg:InvalidOption', ...
  215. 'Option "InitialSlope" will be ignored by this solver');
  216. end
  217. if (~isequal (vodeoptions.MaxOrder, vodetemp.MaxOrder))
  218. warning ('OdePkg:InvalidOption', ...
  219. 'Option "MaxOrder" will be ignored by this solver');
  220. end
  221. if (~isequal (vodeoptions.BDF, vodetemp.BDF))
  222. warning ('OdePkg:InvalidOption', ...
  223. 'Option "BDF" will be ignored by this solver');
  224. end
  225. %# Starting the initialisation of the core solver ode23d
  226. vtimestamp = vslot(1,1); %# timestamp = start time
  227. vtimelength = length (vslot); %# length needed if fixed steps
  228. vtimestop = vslot(1,vtimelength); %# stop time = last value
  229. if (~vstepsizegiven)
  230. vstepsize = vodeoptions.InitialStep;
  231. vminstepsize = (vtimestop - vtimestamp) / (1/eps);
  232. else %# If step size is given then use the fixed time steps
  233. vstepsize = abs (vslot(1,1) - vslot(1,2));
  234. vminstepsize = eps; %# vslot(1,2) - vslot(1,1) - eps;
  235. end
  236. vretvaltime = vtimestamp; %# first timestamp output
  237. if (vhaveoutputselection) %# first solution output
  238. vretvalresult = vinit(vodeoptions.OutputSel);
  239. else vretvalresult = vinit;
  240. end
  241. %# Initialize the OutputFcn
  242. if (vhaveoutputfunction)
  243. feval (vodeoptions.OutputFcn, vslot', ...
  244. vretvalresult', 'init', vfunarguments{:});
  245. end
  246. %# Initialize the History
  247. if (isnumeric (vhist))
  248. vhmat = vhist;
  249. vhavehistnumeric = true;
  250. else %# it must be a function handle
  251. for vcnt = 1:length (vlags);
  252. vhmat(:,vcnt) = feval (vhist, (vslot(1)-vlags(vcnt)), vfunarguments{:});
  253. end
  254. vhavehistnumeric = false;
  255. end
  256. %# Initialize DDE variables for history calculation
  257. vsaveddetime = [vtimestamp - vlags, vtimestamp]';
  258. vsaveddeinput = [vhmat, vinit']';
  259. vsavedderesult = [vhmat, vinit']';
  260. %# Initialize the EventFcn
  261. if (vhaveeventfunction)
  262. odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
  263. {vretvalresult', vhmat}, 'init', vfunarguments{:});
  264. end
  265. vpow = 1/3; %# 20071016, reported by Luis Randez
  266. va = [ 0, 0, 0; %# The Runge-Kutta-Fehlberg 2(3) coefficients
  267. 1/2, 0, 0; %# Coefficients proved on 20060827
  268. -1, 2, 0]; %# See p.91 in Ascher & Petzold
  269. vb2 = [0; 1; 0]; %# 2nd and 3rd order
  270. vb3 = [1/6; 2/3; 1/6]; %# b-coefficients
  271. vc = sum (va, 2);
  272. %# The solver main loop - stop if the endpoint has been reached
  273. vcntloop = 2; vcntcycles = 1; vu = vinit; vk = vu' * zeros(1,3);
  274. vcntiter = 0; vunhandledtermination = true;
  275. while ((vtimestamp < vtimestop && vstepsize >= vminstepsize))
  276. %# Hit the endpoint of the time slot exactely
  277. if ((vtimestamp + vstepsize) > vtimestop)
  278. vstepsize = vtimestop - vtimestamp; end
  279. %# Estimate the three results when using this solver
  280. for j = 1:3
  281. vthetime = vtimestamp + vc(j,1) * vstepsize;
  282. vtheinput = vu' + vstepsize * vk(:,1:j-1) * va(j,1:j-1)';
  283. %# Claculate the history values (or get them from an external
  284. %# function) that are needed for the next step of solving
  285. if (vhavehistnumeric)
  286. for vcnt = 1:length (vlags)
  287. %# Direct implementation of a 'quadrature cubic Hermite interpolation'
  288. %# found at the Faculty for Mathematics of the University of Stuttgart
  289. %# http://mo.mathematik.uni-stuttgart.de/inhalt/aussage/aussage1269
  290. vnumb = find (vthetime - vlags(vcnt) >= vsaveddetime);
  291. velem = min (vnumb(end), length (vsaveddetime) - 1);
  292. vstep = vsaveddetime(velem+1) - vsaveddetime(velem);
  293. vdiff = (vthetime - vlags(vcnt) - vsaveddetime(velem)) / vstep;
  294. vsubs = 1 - vdiff;
  295. %# Calculation of the coefficients for the interpolation algorithm
  296. vua = (1 + 2 * vdiff) * vsubs^2;
  297. vub = (3 - 2 * vdiff) * vdiff^2;
  298. vva = vstep * vdiff * vsubs^2;
  299. vvb = -vstep * vsubs * vdiff^2;
  300. vhmat(:,vcnt) = vua * vsaveddeinput(velem,:)' + ...
  301. vub * vsaveddeinput(velem+1,:)' + ...
  302. vva * vsavedderesult(velem,:)' + ...
  303. vvb * vsavedderesult(velem+1,:)';
  304. end
  305. else %# the history must be a function handle
  306. for vcnt = 1:length (vlags)
  307. vhmat(:,vcnt) = feval ...
  308. (vhist, vthetime - vlags(vcnt), vfunarguments{:});
  309. end
  310. end
  311. if (vhavemasshandle) %# Handle only the dynamic mass matrix,
  312. if (vmassdependence) %# constant mass matrices have already
  313. vmass = feval ... %# been set before (if any)
  314. (vodeoptions.Mass, vthetime, vtheinput, vfunarguments{:});
  315. else %# if (vmassdependence == false)
  316. vmass = feval ... %# then we only have the time argument
  317. (vodeoptions.Mass, vthetime, vfunarguments{:});
  318. end
  319. vk(:,j) = vmass \ feval ...
  320. (vfun, vthetime, vtheinput, vhmat, vfunarguments{:});
  321. else
  322. vk(:,j) = feval ...
  323. (vfun, vthetime, vtheinput, vhmat, vfunarguments{:});
  324. end
  325. end
  326. %# Compute the 2nd and the 3rd order estimation
  327. y2 = vu' + vstepsize * (vk * vb2);
  328. y3 = vu' + vstepsize * (vk * vb3);
  329. if (vhavenonnegative)
  330. vu(vodeoptions.NonNegative) = abs (vu(vodeoptions.NonNegative));
  331. y2(vodeoptions.NonNegative) = abs (y2(vodeoptions.NonNegative));
  332. y3(vodeoptions.NonNegative) = abs (y3(vodeoptions.NonNegative));
  333. end
  334. vSaveVUForRefine = vu;
  335. %# Calculate the absolute local truncation error and the acceptable error
  336. if (~vstepsizegiven)
  337. if (~vnormcontrol)
  338. vdelta = y3 - y2;
  339. vtau = max (vodeoptions.RelTol * vu', vodeoptions.AbsTol);
  340. else
  341. vdelta = norm (y3 - y2, Inf);
  342. vtau = max (vodeoptions.RelTol * max (norm (vu', Inf), 1.0), ...
  343. vodeoptions.AbsTol);
  344. end
  345. else %# if (vstepsizegiven == true)
  346. vdelta = 1; vtau = 2;
  347. end
  348. %# If the error is acceptable then update the vretval variables
  349. if (all (vdelta <= vtau))
  350. vtimestamp = vtimestamp + vstepsize;
  351. vu = y3'; %# MC2001: the higher order estimation as "local extrapolation"
  352. vretvaltime(vcntloop,:) = vtimestamp;
  353. if (vhaveoutputselection)
  354. vretvalresult(vcntloop,:) = vu(vodeoptions.OutputSel);
  355. else
  356. vretvalresult(vcntloop,:) = vu;
  357. end
  358. vcntloop = vcntloop + 1; vcntiter = 0;
  359. %# Update DDE values for next history calculation
  360. vsaveddetime(end+1) = vtimestamp;
  361. vsaveddeinput(end+1,:) = vtheinput';
  362. vsavedderesult(end+1,:) = vu;
  363. %# Call plot only if a valid result has been found, therefore this
  364. %# code fragment has moved here. Stop integration if plot function
  365. %# returns false
  366. if (vhaveoutputfunction)
  367. if (vhaverefine) %# Do interpolation
  368. for vcnt = 0:vodeoptions.Refine %# Approximation between told and t
  369. vapproxtime = (vcnt + 1) * vstepsize / (vodeoptions.Refine + 2);
  370. vapproxvals = vSaveVUForRefine' + vapproxtime * (vk * vb3);
  371. if (vhaveoutputselection)
  372. vapproxvals = vapproxvals(vodeoptions.OutputSel);
  373. end
  374. feval (vodeoptions.OutputFcn, (vtimestamp - vstepsize) + vapproxtime, ...
  375. vapproxvals, [], vfunarguments{:});
  376. end
  377. end
  378. vpltret = feval (vodeoptions.OutputFcn, vtimestamp, ...
  379. vretvalresult(vcntloop-1,:)', [], vfunarguments{:});
  380. if (vpltret), vunhandledtermination = false; break; end
  381. end
  382. %# Call event only if a valid result has been found, therefore this
  383. %# code fragment has moved here. Stop integration if veventbreak is
  384. %# true
  385. if (vhaveeventfunction)
  386. vevent = ...
  387. odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
  388. {vu(:), vhmat}, [], vfunarguments{:});
  389. if (~isempty (vevent{1}) && vevent{1} == 1)
  390. vretvaltime(vcntloop-1,:) = vevent{3}(end,:);
  391. vretvalresult(vcntloop-1,:) = vevent{4}(end,:);
  392. vunhandledtermination = false; break;
  393. end
  394. end
  395. end %# If the error is acceptable ...
  396. %# Update the step size for the next integration step
  397. if (~vstepsizegiven)
  398. %# vdelta may be 0 or even negative - could be an iteration problem
  399. vdelta = max (vdelta, eps);
  400. vstepsize = min (vodeoptions.MaxStep, ...
  401. min (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow));
  402. elseif (vstepsizegiven)
  403. if (vcntloop < vtimelength)
  404. vstepsize = vslot(1,vcntloop-1) - vslot(1,vcntloop-2);
  405. end
  406. end
  407. %# Update counters that count the number of iteration cycles
  408. vcntcycles = vcntcycles + 1; %# Needed for postprocessing
  409. vcntiter = vcntiter + 1; %# Needed to find iteration problems
  410. %# Stop solving because the last 1000 steps no successful valid
  411. %# value has been found
  412. if (vcntiter >= 5000)
  413. error (['Solving has not been successful. The iterative', ...
  414. ' integration loop exited at time t = %f before endpoint at', ...
  415. ' tend = %f was reached. This happened because the iterative', ...
  416. ' integration loop does not find a valid solution at this time', ...
  417. ' stamp. Try to reduce the value of "InitialStep" and/or', ...
  418. ' "MaxStep" with the command "odeset".\n'], vtimestamp, vtimestop);
  419. end
  420. end %# The main loop
  421. %# Check if integration of the ode has been successful
  422. if (vtimestamp < vtimestop)
  423. if (vunhandledtermination == true)
  424. error (['Solving has not been successful. The iterative', ...
  425. ' integration loop exited at time t = %f', ...
  426. ' before endpoint at tend = %f was reached. This may', ...
  427. ' happen if the stepsize grows smaller than defined in', ...
  428. ' vminstepsize. Try to reduce the value of "InitialStep" and/or', ...
  429. ' "MaxStep" with the command "odeset".\n'], vtimestamp, vtimestop);
  430. else
  431. warning ('OdePkg:HideWarning', ...
  432. ['Solver has been stopped by a call of "break" in', ...
  433. ' the main iteration loop at time t = %f before endpoint at', ...
  434. ' tend = %f was reached. This may happen because the @odeplot', ...
  435. ' function returned "true" or the @event function returned "true".'], ...
  436. vtimestamp, vtimestop);
  437. end
  438. end
  439. %# Postprocessing, do whatever when terminating integration algorithm
  440. if (vhaveoutputfunction) %# Cleanup plotter
  441. feval (vodeoptions.OutputFcn, vtimestamp, ...
  442. vretvalresult(vcntloop-1,:)', 'done', vfunarguments{:});
  443. end
  444. if (vhaveeventfunction) %# Cleanup event function handling
  445. odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
  446. {vretvalresult(vcntloop-1,:), vhmat}, 'done', vfunarguments{:});
  447. end
  448. %# Print additional information if option Stats is set
  449. if (strcmp (vodeoptions.Stats, 'on'))
  450. vhavestats = true;
  451. vnsteps = vcntloop-2; %# vcntloop from 2..end
  452. vnfailed = (vcntcycles-1)-(vcntloop-2)+1; %# vcntcycl from 1..end
  453. vnfevals = 3*(vcntcycles-1); %# number of ode evaluations
  454. vndecomps = 0; %# number of LU decompositions
  455. vnpds = 0; %# number of partial derivatives
  456. vnlinsols = 0; %# no. of solutions of linear systems
  457. %# Print cost statistics if no output argument is given
  458. if (nargout == 0)
  459. vmsg = fprintf (1, 'Number of successful steps: %d', vnsteps);
  460. vmsg = fprintf (1, 'Number of failed attempts: %d', vnfailed);
  461. vmsg = fprintf (1, 'Number of function calls: %d', vnfevals);
  462. end
  463. else vhavestats = false;
  464. end
  465. if (nargout == 1) %# Sort output variables, depends on nargout
  466. varargout{1}.x = vretvaltime; %# Time stamps are saved in field x
  467. varargout{1}.y = vretvalresult; %# Results are saved in field y
  468. varargout{1}.solver = 'ode23d'; %# Solver name is saved in field solver
  469. if (vhaveeventfunction)
  470. varargout{1}.ie = vevent{2}; %# Index info which event occured
  471. varargout{1}.xe = vevent{3}; %# Time info when an event occured
  472. varargout{1}.ye = vevent{4}; %# Results when an event occured
  473. end
  474. if (vhavestats)
  475. varargout{1}.stats = struct;
  476. varargout{1}.stats.nsteps = vnsteps;
  477. varargout{1}.stats.nfailed = vnfailed;
  478. varargout{1}.stats.nfevals = vnfevals;
  479. varargout{1}.stats.npds = vnpds;
  480. varargout{1}.stats.ndecomps = vndecomps;
  481. varargout{1}.stats.nlinsols = vnlinsols;
  482. end
  483. elseif (nargout == 2)
  484. varargout{1} = vretvaltime; %# Time stamps are first output argument
  485. varargout{2} = vretvalresult; %# Results are second output argument
  486. elseif (nargout == 5)
  487. varargout{1} = vretvaltime; %# Same as (nargout == 2)
  488. varargout{2} = vretvalresult; %# Same as (nargout == 2)
  489. varargout{3} = []; %# LabMat doesn't accept lines like
  490. varargout{4} = []; %# varargout{3} = varargout{4} = [];
  491. varargout{5} = [];
  492. if (vhaveeventfunction)
  493. varargout{3} = vevent{3}; %# Time info when an event occured
  494. varargout{4} = vevent{4}; %# Results when an event occured
  495. varargout{5} = vevent{2}; %# Index info which event occured
  496. end
  497. %# else nothing will be returned, varargout{1} undefined
  498. end
  499. %! # We are using a "pseudo-DDE" implementation for all tests that
  500. %! # are done for this function. We also define an Events and a
  501. %! # pseudo-Mass implementation. For further tests we also define a
  502. %! # reference solution (computed at high accuracy) and an OutputFcn.
  503. %!function [vyd] = fexp (vt, vy, vz, varargin)
  504. %! vyd(1,1) = exp (- vt) - vz(1); %# The DDEs that are
  505. %! vyd(2,1) = vy(1) - vz(2); %# used for all examples
  506. %!function [vval, vtrm, vdir] = feve (vt, vy, vz, varargin)
  507. %! vval = fexp (vt, vy, vz); %# We use the derivatives
  508. %! vtrm = zeros (2,1); %# don't stop solving here
  509. %! vdir = ones (2,1); %# in positive direction
  510. %!function [vval, vtrm, vdir] = fevn (vt, vy, vz, varargin)
  511. %! vval = fexp (vt, vy, vz); %# We use the derivatives
  512. %! vtrm = ones (2,1); %# stop solving here
  513. %! vdir = ones (2,1); %# in positive direction
  514. %!function [vmas] = fmas (vt, vy, vz, varargin)
  515. %! vmas = [1, 0; 0, 1]; %# Dummy mass matrix for tests
  516. %!function [vmas] = fmsa (vt, vy, vz, varargin)
  517. %! vmas = sparse ([1, 0; 0, 1]); %# A dummy sparse matrix
  518. %!function [vref] = fref () %# The reference solution
  519. %! vref = [0.12194462133618, 0.01652432423938];
  520. %!function [vout] = fout (vt, vy, vflag, varargin)
  521. %! if (regexp (char (vflag), 'init') == 1)
  522. %! if (any (size (vt) ~= [2, 1])) error ('"fout" step "init"'); end
  523. %! elseif (isempty (vflag))
  524. %! if (any (size (vt) ~= [1, 1])) error ('"fout" step "calc"'); end
  525. %! vout = false;
  526. %! elseif (regexp (char (vflag), 'done') == 1)
  527. %! if (any (size (vt) ~= [1, 1])) error ('"fout" step "done"'); end
  528. %! else error ('"fout" invalid vflag');
  529. %! end
  530. %!
  531. %! %# Turn off output of warning messages for all tests, turn them on
  532. %! %# again if the last test is called
  533. %!error %# input argument number one
  534. %! warning ('off', 'OdePkg:InvalidOption');
  535. %! B = ode23d (1, [0 5], [1; 0], 1, [1; 0]);
  536. %!error %# input argument number two
  537. %! B = ode23d (@fexp, 1, [1; 0], 1, [1; 0]);
  538. %!error %# input argument number three
  539. %! B = ode23d (@fexp, [0 5], 1, 1, [1; 0]);
  540. %!error %# input argument number four
  541. %! B = ode23d (@fexp, [0 5], [1; 0], [1; 1], [1; 0]);
  542. %!error %# input argument number five
  543. %! B = ode23d (@fexp, [0 5], [1; 0], 1, 1);
  544. %!test %# one output argument
  545. %! vsol = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0]);
  546. %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 1e-1);
  547. %! assert (isfield (vsol, 'solver'));
  548. %! assert (vsol.solver, 'ode23d');
  549. %!test %# two output arguments
  550. %! [vt, vy] = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0]);
  551. %! assert ([vt(end), vy(end,:)], [5, fref], 1e-1);
  552. %!test %# five output arguments and no Events
  553. %! [vt, vy, vxe, vye, vie] = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0]);
  554. %! assert ([vt(end), vy(end,:)], [5, fref], 1e-1);
  555. %! assert ([vie, vxe, vye], []);
  556. %!test %# anonymous function instead of real function
  557. %! faym = @(vt, vy, vz) [exp(-vt) - vz(1); vy(1) - vz(2)];
  558. %! vsol = ode23d (faym, [0 5], [1; 0], 1, [1; 0]);
  559. %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 1e-1);
  560. %!test %# extra input arguments passed trhough
  561. %! vsol = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0], 'KL');
  562. %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 1e-1);
  563. %!test %# empty OdePkg structure *but* extra input arguments
  564. %! vopt = odeset;
  565. %! vsol = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt, 12, 13, 'KL');
  566. %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 1e-1);
  567. %!error %# strange OdePkg structure
  568. %! vopt = struct ('foo', 1);
  569. %! vsol = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
  570. %!test %# AbsTol option
  571. %! vopt = odeset ('AbsTol', 1e-5);
  572. %! vsol = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
  573. %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 1e-1);
  574. %!test %# AbsTol and RelTol option
  575. %! vopt = odeset ('AbsTol', 1e-7, 'RelTol', 1e-7);
  576. %! vsol = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
  577. %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 1e-1);
  578. %!test %# RelTol and NormControl option
  579. %! vopt = odeset ('AbsTol', 1e-7, 'NormControl', 'on');
  580. %! vsol = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
  581. %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], .5e-1);
  582. %!test %# NonNegative for second component
  583. %! vopt = odeset ('NonNegative', 1);
  584. %! vsol = ode23d (@fexp, [0 2.5], [1; 0], 1, [1; 0], vopt);
  585. %! assert ([vsol.x(end), vsol.y(end,:)], [2.5, 0.001, 0.237], 1e-1);
  586. %!test %# Details of OutputSel and Refine can't be tested
  587. %! vopt = odeset ('OutputFcn', @fout, 'OutputSel', 1, 'Refine', 5);
  588. %! vsol = ode23d (@fexp, [0 2.5], [1; 0], 1, [1; 0], vopt);
  589. %!test %# Stats must add further elements in vsol
  590. %! vopt = odeset ('Stats', 'on');
  591. %! vsol = ode23d (@fexp, [0 2.5], [1; 0], 1, [1; 0], vopt);
  592. %! assert (isfield (vsol, 'stats'));
  593. %! assert (isfield (vsol.stats, 'nsteps'));
  594. %!test %# InitialStep option
  595. %! vopt = odeset ('InitialStep', 1e-8);
  596. %! vsol = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
  597. %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 1e-1);
  598. %!test %# MaxStep option
  599. %! vopt = odeset ('MaxStep', 1e-2);
  600. %! vsol = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
  601. %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 1e-1);
  602. %!test %# Events option add further elements in vsol
  603. %! vopt = odeset ('Events', @feve);
  604. %! vsol = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
  605. %! assert (isfield (vsol, 'ie'));
  606. %! assert (vsol.ie, [1; 1]);
  607. %! assert (isfield (vsol, 'xe'));
  608. %! assert (isfield (vsol, 'ye'));
  609. %!test %# Events option, now stop integration
  610. %! warning ('off', 'OdePkg:HideWarning');
  611. %! vopt = odeset ('Events', @fevn, 'NormControl', 'on');
  612. %! vsol = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
  613. %! assert ([vsol.ie, vsol.xe, vsol.ye], ...
  614. %! [1.0000, 2.9219, -0.2127, -0.2671], 1e-1);
  615. %!test %# Events option, five output arguments
  616. %! vopt = odeset ('Events', @fevn, 'NormControl', 'on');
  617. %! [vt, vy, vxe, vye, vie] = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
  618. %! assert ([vie, vxe, vye], ...
  619. %! [1.0000, 2.9219, -0.2127, -0.2671], 1e-1);
  620. %!
  621. %! %# test for Jacobian option is missing
  622. %! %# test for Jacobian (being a sparse matrix) is missing
  623. %! %# test for JPattern option is missing
  624. %! %# test for Vectorized option is missing
  625. %!
  626. %!test %# Mass option as function
  627. %! vopt = odeset ('Mass', eye (2,2));
  628. %! vsol = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
  629. %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 1e-1);
  630. %!test %# Mass option as matrix
  631. %! vopt = odeset ('Mass', eye (2,2));
  632. %! vsol = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
  633. %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 1e-1);
  634. %!test %# Mass option as sparse matrix
  635. %! vopt = odeset ('Mass', sparse (eye (2,2)));
  636. %! vsol = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
  637. %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 1e-1);
  638. %!test %# Mass option as function and sparse matrix
  639. %! vopt = odeset ('Mass', @fmsa);
  640. %! vsol = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
  641. %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 1e-1);
  642. %!test %# Mass option as function and MStateDependence
  643. %! vopt = odeset ('Mass', @fmas, 'MStateDependence', 'strong');
  644. %! vsol = ode23d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
  645. %! assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 1e-1);
  646. %!
  647. %! %# test for MvPattern option is missing
  648. %! %# test for InitialSlope option is missing
  649. %! %# test for MaxOrder option is missing
  650. %! %# test for BDF option is missing
  651. %!
  652. %! warning ('on', 'OdePkg:InvalidOption');
  653. %# Local Variables: ***
  654. %# mode: octave ***
  655. %# End: ***