Mercurial > octave
comparison scripts/ode/ode23.m @ 20901:afe9c529760d
2015 Code Sprint: move ode23 and runge_kutta_23 from odepkg to core
* scripts/ode/ode23.m: new file
* scripts/ode/private/runge_kutta_23.m: new file
* scripts/ode/module.mk: list new files
* doc/interpreter/diffeq.txi: mention ode23 among available solvers
* scripts/help/__unimplemented__.m: remove ode23 from list of unimplemented functions
author | Stefan Miereis <stefan.miereis@gmx.de> |
---|---|
date | Tue, 15 Dec 2015 13:59:17 +0100 |
parents | |
children | 73cf3434e8c9 |
comparison
equal
deleted
inserted
replaced
20900:4d14b0a22b29 | 20901:afe9c529760d |
---|---|
1 ## Copyright (C) 2014-2015, Jacopo Corno <jacopo.corno@gmail.com> | |
2 ## Copyright (C) 2013-2015, Roberto Porcu' <roberto.porcu@polimi.it> | |
3 ## Copyright (C) 2006-2015, Thomas Treichl <treichl@users.sourceforge.net> | |
4 ## | |
5 ## This file is part of Octave. | |
6 ## | |
7 ## Octave is free software; you can redistribute it and/or modify it | |
8 ## under the terms of the GNU General Public License as published by | |
9 ## the Free Software Foundation; either version 3 of the License, or (at | |
10 ## your option) any later version. | |
11 ## | |
12 ## Octave is distributed in the hope that it will be useful, but | |
13 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
14 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
15 ## General Public License for more details. | |
16 ## | |
17 ## You should have received a copy of the GNU General Public License | |
18 ## along with Octave; see the file COPYING. If not, see | |
19 ## <http://www.gnu.org/licenses/>. | |
20 | |
21 ## -*- texinfo -*- | |
22 ## @deftypefn {Function File} {[@var{t}, @var{y}] =} ode23 (@var{fun}, @var{trange}, @var{init}) | |
23 ## @deftypefnx {Function File} {[@var{t}, @var{y}] =} ode23 (@var{fun}, @var{trange}, @var{init}, @var{ode_opt}) | |
24 ## @deftypefnx {Function File} {[@var{t}, @var{y}] =} ode23 (@dots{}, @var{par1}, @var{par2}, @dots{}) | |
25 ## @deftypefnx {Function File} {[@var{t}, @var{y}, @var{te}, @var{ye}, @var{ie}] =} ode23 (@dots{}) | |
26 ## @deftypefnx {Function File} {@var{solution} =} ode23 (@dots{}) | |
27 ## | |
28 ## Solve a set of non-stiff Ordinary Differential Equations (non-stiff ODEs) | |
29 ## with the well known explicit Bogacki-Shampine method of order 3. For the | |
30 ## definition of this method see | |
31 ## @url{http://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods}. | |
32 ## | |
33 ## @var{fun} is a function handle, inline function, or string containing the | |
34 ## name of the function that defines the ODE: @code{y' = f(t,y)}. The function | |
35 ## must accept two inputs where the first is time @var{t} and the second is a | |
36 ## column vector of unknowns @var{y}. | |
37 ## | |
38 ## @var{trange} specifies the time interval over which the ODE will be | |
39 ## evaluated. Typically, it is a two-element vector specifying the initial and | |
40 ## final times (@code{[tinit, tfinal]}). If there are more than two elements | |
41 ## then the solution will also be evaluated at these intermediate time | |
42 ## instances unless the integrate function specified is | |
43 ## @command{integrate_n_steps}. | |
44 ## | |
45 ## By default, @code{ode23} uses an adaptive timestep with the | |
46 ## @code{integrate_adaptive} algorithm. The tolerance for the timestep | |
47 ## computation may be changed by using the option @qcode{"Tau"}, that has a | |
48 ## default value of @math{1e-6}. If the ODE option @qcode{"TimeStepSize"} is | |
49 ## not empty, then the stepper called will be @code{integrate_const}. If, in | |
50 ## addition, the option @qcode{"TimeStepNumber"} is also specified then the | |
51 ## integrate function @code{integrate_n_steps} will be used. | |
52 ## | |
53 ## @var{init} contains the initial value for the unknowns. If it is a row | |
54 ## vector then the solution @var{y} will be a matrix in which each column is | |
55 ## the solution for the corresponding initial value in @var{init}. | |
56 ## | |
57 ## The optional fourth argument @var{ode_opt} specifies non-default options to | |
58 ## the ODE solver. It is a structure generated by @code{odeset}. | |
59 ## | |
60 ## The function typically returns two outputs. Variable @var{t} is a | |
61 ## column vector and contains the times where the solution was found. The | |
62 ## output @var{y} is a matrix in which each column refers to a different | |
63 ## unknown of the problem and each row corresponds to a time in @var{t}. | |
64 ## | |
65 ## The output can also be returned as a structure @var{solution} which | |
66 ## has field @var{x} containing the time where the solution was evaluated and | |
67 ## field @var{y} containing the solution matrix for the times in @var{x}. | |
68 ## Use @code{fieldnames (@var{solution})} to see the other fields and additional | |
69 ## information returned. | |
70 ## | |
71 ## If using the @qcode{"Events"} option then three additional outputs may | |
72 ## be returned. @var{te} holds the time when an Event function returned a | |
73 ## zero. @var{ye} holds the value of the solution at time @var{te}. @var{ie} | |
74 ## contains an index indicating which Event function was triggered in the case | |
75 ## of multiple Event functions. | |
76 ## | |
77 ## This function can be called with two output arguments: @var{t} and @var{y}. | |
78 ## Variable @var{t} is a column vector and contains the time stamps, instead | |
79 ## @var{y} is a matrix in which each column refers to a different unknown of the | |
80 ## problem and the rows number is the same of @var{t} rows number so that each | |
81 ## row of @var{y} contains the values of all unknowns at the time value | |
82 ## contained in the corresponding row in @var{t}. | |
83 ## | |
84 ## Example: Solve the Van der Pol equation | |
85 ## | |
86 ## @example | |
87 ## @group | |
88 ## fvdp = @@(@var{t},@var{y}) [@var{y}(2); (1 - @var{y}(1)^2) * @var{y}(2) - @var{y}(1)]; | |
89 ## [@var{t},@var{y}] = ode23 (fvdp, [0 20], [2 0]); | |
90 ## @end group | |
91 ## @end example | |
92 ## @seealso{odeset, odeget} | |
93 ## @end deftypefn | |
94 | |
95 ## ChangeLog: | |
96 ## 20010703 the function file "ode23.m" was written by Marc Compere | |
97 ## under the GPL for the use with this software. This function has been | |
98 ## taken as a base for the following implementation. | |
99 ## 20060810, Thomas Treichl | |
100 ## This function was adapted to the new syntax that is used by the | |
101 ## new OdePkg for Octave and is compatible to Matlab's ode23. | |
102 | |
103 function varargout = ode23 (fun, trange, init, varargin) | |
104 | |
105 if (nargin < 3) | |
106 print_usage (); | |
107 endif | |
108 | |
109 order = 3; | |
110 solver = "ode23"; | |
111 | |
112 if (nargin >= 4) | |
113 if (! isstruct (varargin{1})) | |
114 ## varargin{1:len} are parameters for fun | |
115 odeopts = odeset (); | |
116 odeopts.funarguments = varargin; | |
117 elseif (length (varargin) > 1) | |
118 ## varargin{1} is an ODE options structure opt | |
119 odeopts = ode_struct_value_check ("ode23", varargin{1}, "ode23"); | |
120 odeopts.funarguments = {varargin{2:length(varargin)}}; | |
121 else # if (isstruct (varargin{1})) | |
122 odeopts = ode_struct_value_check ("ode23", varargin{1}, "ode23"); | |
123 odeopts.funarguments = {}; | |
124 endif | |
125 else # nargin == 3 | |
126 odeopts = odeset (); | |
127 odeopts.funarguments = {}; | |
128 endif | |
129 | |
130 if (! isvector (trange) || ! isnumeric (trange)) | |
131 error ("Octave:invalid-input-arg", | |
132 "ode23: TRANGE must be a numeric vector"); | |
133 endif | |
134 | |
135 TimeStepNumber = odeget (odeopts, "TimeStepNumber", [], "fast"); | |
136 TimeStepSize = odeget (odeopts, "TimeStepSize", [], "fast"); | |
137 if (length (trange) < 2 | |
138 && (isempty (TimeStepSize) || isempty (TimeStepNumber))) | |
139 error ("Octave:invalid-input-arg", | |
140 "ode23: TRANGE must contain at least 2 elements"); | |
141 elseif (trange(2) == trange(1)) | |
142 error ("Octave:invalid-input-arg", | |
143 "ode23: invalid time span, TRANGE(1) == TRANGE(2)"); | |
144 else | |
145 odeopts.direction = sign (trange(2) - trange(1)); | |
146 endif | |
147 trange = trange(:); | |
148 | |
149 if (! isvector (init) || ! isnumeric (init)) | |
150 error ("Octave:invalid-input-arg", | |
151 "ode23: INIT must be a numeric vector"); | |
152 endif | |
153 init = init(:); | |
154 | |
155 if (ischar (fun)) | |
156 try | |
157 fun = str2func (fun); | |
158 catch | |
159 warning (lasterr); | |
160 end_try_catch | |
161 endif | |
162 if (! isa (fun, "function_handle")) | |
163 error ("Octave:invalid-input-arg", | |
164 "ode23: FUN must be a valid function handle"); | |
165 endif | |
166 | |
167 ## Start preprocessing, have a look which options are set in odeopts, | |
168 ## check if an invalid or unused option is set | |
169 if (isempty (TimeStepNumber) && isempty (TimeStepSize)) | |
170 integrate_func = "adaptive"; | |
171 odeopts.stepsizefixed = false; | |
172 elseif (! isempty (TimeStepNumber) && ! isempty (TimeStepSize)) | |
173 integrate_func = "n_steps"; | |
174 odeopts.stepsizefixed = true; | |
175 if (sign (TimeStepSize) != odeopts.direction) | |
176 warning ("Octave:invalid-input-arg", | |
177 ["ode23: option \"TimeStepSize\" has the wrong sign, ", ... | |
178 "but will be corrected automatically\n"]); | |
179 TimeStepSize = -TimeStepSize; | |
180 endif | |
181 elseif (isempty (TimeStepNumber) && ! isempty (TimeStepSize)) | |
182 integrate_func = "const"; | |
183 odeopts.stepsizefixed = true; | |
184 if (sign (TimeStepSize) != odeopts.direction) | |
185 warning ("Octave:invalid-input-arg", | |
186 ["ode23: option \"TimeStepSize\" has the wrong sign, ", | |
187 "but will be corrected automatically\n"]); | |
188 TimeStepSize = -TimeStepSize; | |
189 endif | |
190 else | |
191 warning ("Octave:invalid-input-arg", | |
192 "ode23: assuming an adaptive integrate function\n"); | |
193 integrate_func = "adaptive"; | |
194 endif | |
195 | |
196 if (isempty (odeopts.RelTol) && ! odeopts.stepsizefixed) | |
197 odeopts.RelTol = 1e-3; | |
198 warning ("Octave:invalid-input-arg", | |
199 "ode23: option \"RelTol\" not set, new value %f will be used\n", | |
200 odeopts.RelTol); | |
201 elseif (! isempty (odeopts.RelTol) && odeopts.stepsizefixed) | |
202 warning ("Octave:invalid-input-arg", | |
203 ["ode23: option \"RelTol\" is ignored", ... | |
204 " when fixed time stamps are given\n"]); | |
205 endif | |
206 | |
207 if (isempty (odeopts.AbsTol) && ! odeopts.stepsizefixed) | |
208 odeopts.AbsTol = 1e-6; | |
209 warning ("Octave:invalid-input-arg", | |
210 "ode23: option \"AbsTol\" not set, new value %f will be used\n", | |
211 odeopts.AbsTol); | |
212 elseif (! isempty (odeopts.AbsTol) && odeopts.stepsizefixed) | |
213 warning ("Octave:invalid-input-arg", | |
214 ["ode23: option \"AbsTol\" is ignored", ... | |
215 " when fixed time stamps are given\n"]); | |
216 else | |
217 odeopts.AbsTol = odeopts.AbsTol(:); # Create column vector | |
218 endif | |
219 | |
220 odeopts.normcontrol = strcmp (odeopts.NormControl, "on"); | |
221 | |
222 if (! isempty (odeopts.NonNegative)) | |
223 if (isempty (odeopts.Mass)) | |
224 odeopts.havenonnegative = true; | |
225 else | |
226 odeopts.havenonnegative = false; | |
227 warning ("Octave:invalid-input-arg", | |
228 ["ode23: option \"NonNegative\" is ignored", ... | |
229 " when mass matrix is set\n"]); | |
230 endif | |
231 else | |
232 odeopts.havenonnegative = false; | |
233 endif | |
234 | |
235 if (isempty (odeopts.OutputFcn) && nargout == 0) | |
236 odeopts.OutputFcn = @odeplot; | |
237 odeopts.haveoutputfunction = true; | |
238 else | |
239 odeopts.haveoutputfunction = ! isempty (odeopts.OutputFcn); | |
240 endif | |
241 | |
242 odeopts.haveoutputselection = ! isempty (odeopts.OutputSel); | |
243 | |
244 if (odeopts.Refine > 0) | |
245 odeopts.haverefine = true; | |
246 else | |
247 odeopts.haverefine = false; | |
248 endif | |
249 | |
250 if (isempty (odeopts.InitialStep) && strcmp (integrate_func, "adaptive")) | |
251 odeopts.InitialStep = odeopts.direction* ... | |
252 starting_stepsize (order, fun, trange(1), init, odeopts.AbsTol, | |
253 odeopts.RelTol, odeopts.normcontrol); | |
254 warning ("Octave:invalid-input-arg", | |
255 ["ode23: option \"InitialStep\" not set,", ... | |
256 " estimated value %f will be used\n"], | |
257 odeopts.InitialStep); | |
258 elseif (isempty (odeopts.InitialStep)) | |
259 odeopts.InitialStep = TimeStepSize; | |
260 endif | |
261 | |
262 if (isempty (odeopts.MaxStep) && ! odeopts.stepsizefixed) | |
263 odeopts.MaxStep = abs (trange(end) - trange(1)) / 10; | |
264 warning ("Octave:invalid-input-arg", | |
265 "ode23: option \"MaxStep\" not set, new value %f will be used\n", | |
266 odeopts.MaxStep); | |
267 endif | |
268 | |
269 odeopts.haveeventfunction = ! isempty (odeopts.Events); | |
270 | |
271 ## The options 'Jacobian', 'JPattern' and 'Vectorized' will be ignored | |
272 ## by this solver because this solver uses an explicit Runge-Kutta method | |
273 ## and therefore no Jacobian calculation is necessary | |
274 if (! isempty (odeopts.Jacobian)) | |
275 warning ("Octave:invalid-input-arg", | |
276 "ode23: option \"Jacobian\" is ignored by this solver\n"); | |
277 endif | |
278 | |
279 if (! isempty (odeopts.JPattern)) | |
280 warning ("Octave:invalid-input-arg", | |
281 "ode23: option \"JPattern\" is ignored by this solver\n"); | |
282 endif | |
283 | |
284 if (! isempty (odeopts.Vectorized)) | |
285 warning ("Octave:invalid-input-arg", | |
286 "ode23: option \"Vectorized\" is ignored by this solver\n"); | |
287 endif | |
288 | |
289 if (! isempty (odeopts.Mass) && isnumeric (odeopts.Mass)) | |
290 havemasshandle = false; | |
291 mass = odeopts.Mass; # constant mass | |
292 elseif (isa (odeopts.Mass, "function_handle")) | |
293 havemasshandle = true; # mass defined by a function handle | |
294 else # no mass matrix - creating a diag-matrix of ones for mass | |
295 havemasshandle = false; # mass = diag (ones (length (init), 1), 0); | |
296 endif | |
297 | |
298 massdependence = ! strcmp (odeopts.MStateDependence, "none"); | |
299 | |
300 ## Other options that are not used by this solver. | |
301 if (! isempty (odeopts.MvPattern)) | |
302 warning ("Octave:invalid-input-arg", | |
303 "ode23: option \"MvPattern\" is ignored by this solver\n"); | |
304 endif | |
305 | |
306 if (! isempty (odeopts.MassSingular)) | |
307 warning ("Octave:invalid-input-arg", | |
308 "ode23: option \"MassSingular\" is ignored by this solver\n"); | |
309 endif | |
310 | |
311 if (! isempty (odeopts.InitialSlope)) | |
312 warning ("Octave:invalid-input-arg", | |
313 "ode23: option \"InitialSlope\" is ignored by this solver\n"); | |
314 endif | |
315 | |
316 if (! isempty (odeopts.MaxOrder)) | |
317 warning ("Octave:invalid-input-arg", | |
318 "ode23: option \"MaxOrder\" is ignored by this solver\n"); | |
319 endif | |
320 | |
321 if (! isempty (odeopts.BDF)) | |
322 warning ("Octave:invalid-input-arg", | |
323 "ode23: option \"BDF\" is ignored by this solver\n"); | |
324 endif | |
325 | |
326 ## Starting the initialisation of the core solver ode23 | |
327 | |
328 if (havemasshandle) # Handle only the dynamic mass matrix, | |
329 if (massdependence) # constant mass matrices have already | |
330 mass = @(t,x) odeopts.Mass (t, x, odeopts.funarguments{:}); | |
331 fun = @(t,x) mass (t, x, odeopts.funarguments{:}) ... | |
332 \ fun (t, x, odeopts.funarguments{:}); | |
333 else # if (massdependence == false) | |
334 mass = @(t) odeopts.Mass (t, odeopts.funarguments{:}); | |
335 fun = @(t,x) mass (t, odeopts.funarguments{:}) ... | |
336 \ fun (t, x, odeopts.funarguments{:}); | |
337 endif | |
338 endif | |
339 | |
340 switch (integrate_func) | |
341 case "adaptive" | |
342 solution = integrate_adaptive (@runge_kutta_23, ... | |
343 order, fun, trange, init, odeopts); | |
344 case "n_steps" | |
345 solution = integrate_n_steps (@runge_kutta_23, ... | |
346 fun, trange(1), init, ... | |
347 TimeStepSize, TimeStepNumber, odeopts); | |
348 case "const" | |
349 solution = integrate_const (@runge_kutta_23, ... | |
350 fun, trange, init, ... | |
351 TimeStepSize, odeopts); | |
352 endswitch | |
353 | |
354 ## Postprocessing, do whatever when terminating integration algorithm | |
355 if (odeopts.haveoutputfunction) # Cleanup plotter | |
356 feval (odeopts.OutputFcn, solution.t(end), ... | |
357 solution.x(end,:)', "done", odeopts.funarguments{:}); | |
358 endif | |
359 if (odeopts.haveeventfunction) # Cleanup event function handling | |
360 ode_event_handler (odeopts.Events, solution.t(end), ... | |
361 solution.x(end,:)', "done", odeopts.funarguments{:}); | |
362 endif | |
363 | |
364 ## Print additional information if option Stats is set | |
365 if (strcmp (odeopts.Stats, "on")) | |
366 havestats = true; | |
367 nsteps = solution.cntloop-2; # cntloop from 2..end | |
368 nfailed = (solution.cntcycles-1)-nsteps+1; # cntcycl from 1..end | |
369 nfevals = 4 * (solution.cntcycles-1); # number of ode evaluations | |
370 ndecomps = 0; # number of LU decompositions | |
371 npds = 0; # number of partial derivatives | |
372 nlinsols = 0; # no. of solutions of linear systems | |
373 ## Print cost statistics if no output argument is given | |
374 if (nargout == 0) | |
375 printf ("Number of successful steps: %d\n", nsteps); | |
376 printf ("Number of failed attempts: %d\n", nfailed); | |
377 printf ("Number of function calls: %d\n", nfevals); | |
378 endif | |
379 else | |
380 havestats = false; | |
381 endif | |
382 | |
383 if (nargout == 2) | |
384 varargout{1} = solution.t; # Time stamps are first output argument | |
385 varargout{2} = solution.x; # Results are second output argument | |
386 elseif (nargout == 1) | |
387 varargout{1}.x = solution.t; # Time stamps are saved in field x | |
388 varargout{1}.y = solution.x; # Results are saved in field y | |
389 varargout{1}.solver = solver; # Solver name is saved in field solver | |
390 if (odeopts.haveeventfunction) | |
391 varargout{1}.ie = solution.event{2}; # Index info which event occurred | |
392 varargout{1}.xe = solution.event{3}; # Time info when an event occurred | |
393 varargout{1}.ye = solution.event{4}; # Results when an event occurred | |
394 endif | |
395 if (havestats) | |
396 varargout{1}.stats = struct (); | |
397 varargout{1}.stats.nsteps = nsteps; | |
398 varargout{1}.stats.nfailed = nfailed; | |
399 varargout{1}.stats.nfevals = nfevals; | |
400 varargout{1}.stats.npds = npds; | |
401 varargout{1}.stats.ndecomps = ndecomps; | |
402 varargout{1}.stats.nlinsols = nlinsols; | |
403 endif | |
404 elseif (nargout == 5) | |
405 varargout = cell (1,5); | |
406 varargout{1} = solution.t; | |
407 varargout{2} = solution.x; | |
408 if (odeopts.haveeventfunction) | |
409 varargout{3} = solution.event{3}; # Time info when an event occurred | |
410 varargout{4} = solution.event{4}; # Results when an event occurred | |
411 varargout{5} = solution.event{2}; # Index info which event occurred | |
412 endif | |
413 endif | |
414 | |
415 endfunction | |
416 | |
417 | |
418 ## We are using the "Van der Pol" implementation for all tests that are done | |
419 ## for this function. | |
420 ## For further tests we also define a reference solution (computed at high | |
421 ## accuracy) | |
422 %!function ydot = fpol (t, y) # The Van der Pol | |
423 %! ydot = [y(2); (1 - y(1)^2) * y(2) - y(1)]; | |
424 %!endfunction | |
425 %!function ref = fref () # The computed reference sol | |
426 %! ref = [0.32331666704577, -1.83297456798624]; | |
427 %!endfunction | |
428 %!function jac = fjac (t, y, varargin) # its Jacobian | |
429 %! jac = [0, 1; -1 - 2 * y(1) * y(2), 1 - y(1)^2]; | |
430 %!endfunction | |
431 %!function jac = fjcc (t, y, varargin) # sparse type | |
432 %! jac = sparse ([0, 1; -1 - 2 * y(1) * y(2), 1 - y(1)^2]) | |
433 %!endfunction | |
434 %!function [val, trm, dir] = feve (t, y, varargin) | |
435 %! val = fpol (t, y, varargin); # We use the derivatives | |
436 %! trm = zeros (2,1); # that's why component 2 | |
437 %! dir = ones (2,1); # seems to not be exact | |
438 %!endfunction | |
439 %!function [val, trm, dir] = fevn (t, y, varargin) | |
440 %! val = fpol (t, y, varargin); # We use the derivatives | |
441 %! trm = ones (2,1); # that's why component 2 | |
442 %! dir = ones (2,1); # seems to not be exact | |
443 %!endfunction | |
444 %!function mas = fmas (t, y, varargin) | |
445 %! mas = [1, 0; 0, 1]; # Dummy mass matrix for tests | |
446 %!endfunction | |
447 %!function mas = fmsa (t, y, varargin) | |
448 %! mas = sparse ([1, 0; 0, 1]); # A sparse dummy matrix | |
449 %!endfunction | |
450 %!function out = fout (t, y, flag, varargin) | |
451 %! if (regexp (char (flag), "init") == 1) | |
452 %! if (any (size (t) != [2, 1])) error ("\"fout\" step \"init\""); endif | |
453 %! elseif (isempty (flag)) | |
454 %! if (any (size (t) != [1, 1])) error ("\"fout\" step \"calc\""); endif | |
455 %! out = false; | |
456 %! elseif (regexp (char (flag), "done") == 1) | |
457 %! if (any (size (t) != [1, 1])) error ("\"fout\" step \"done\""); endif | |
458 %! else | |
459 %! error ("\"fout\" invalid flag"); | |
460 %! endif | |
461 %!endfunction | |
462 %! | |
463 %! ## Turn off output of warning messages for all tests, | |
464 %! ## turn them on again after the last test is called | |
465 %!test | |
466 %! warning ("off", "Octave:invalid-input-arg", "local"); | |
467 %!error ## ouput argument | |
468 %! B = ode23 (1, [0 25], [3 15 1]); | |
469 %!error ## input argument number one | |
470 %! [t, y] = ode23 (1, [0 25], [3 15 1]); | |
471 %!error ## input argument number two | |
472 %! [t, y] = ode23 (@fpol, 1, [3 15 1]); | |
473 %!test ## two output arguments | |
474 %! [t, y] = ode23 (@fpol, [0 2], [2 0]); | |
475 %! assert ([t(end), y(end,:)], [2, fref], 1e-3); | |
476 %!test ## anonymous function instead of real function | |
477 %! fvdb = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; | |
478 %! [t, y] = ode23 (fvdb, [0 2], [2 0]); | |
479 %! assert ([t(end), y(end,:)], [2, fref], 1e-3); | |
480 %!test ## extra input arguments passed through | |
481 %! [t, y] = ode23 (@fpol, [0 2], [2 0], 12, 13, "KL"); | |
482 %! assert ([t(end), y(end,:)], [2, fref], 1e-3); | |
483 %!test ## empty OdePkg structure *but* extra input arguments | |
484 %! opt = odeset; | |
485 %! [t, y] = ode23 (@fpol, [0 2], [2 0], opt, 12, 13, "KL"); | |
486 %! assert ([t(end), y(end,:)], [2, fref], 1e-2); | |
487 %!error ## strange OdePkg structure | |
488 %! opt = struct ("foo", 1); | |
489 %! [t, y] = ode23 (@fpol, [0 2], [2 0], opt); | |
490 %!test ## Solve vdp in fixed step sizes | |
491 %! opt = odeset("TimeStepSize",0.1); | |
492 %! [t, y] = ode23 (@fpol, [0,2], [2 0], opt); | |
493 %! assert (t(:), [0:0.1:2]', 1e-3); | |
494 %!test ## Solve another anonymous function below zero | |
495 %! ref = [0, 14.77810590694212]; | |
496 %! [t, y] = ode23 (@(t,y) y, [-2 0], 2); | |
497 %! assert ([t(end), y(end,:)], ref, 1e-2); | |
498 %!test ## InitialStep option | |
499 %! opt = odeset ("InitialStep", 1e-8); | |
500 %! [t, y] = ode23 (@fpol, [0 0.2], [2 0], opt); | |
501 %! assert ([t(2)-t(1)], [1e-8], 1e-9); | |
502 %!test ## MaxStep option | |
503 %! opt = odeset ("MaxStep", 1e-3); | |
504 %! sol = ode23 (@fpol, [0 0.2], [2 0], opt); | |
505 %! assert ([sol.x(5)-sol.x(4)], [1e-3], 1e-4); | |
506 %!test ## Solve in backward direction starting at t=0 | |
507 %! ref = [-1.205364552835178, 0.951542399860817]; | |
508 %! sol = ode23 (@fpol, [0 -2], [2 0]); | |
509 %! assert ([sol.x(end), sol.y(end,:)], [-2, ref], 5e-3); | |
510 %!test ## Solve in backward direction starting at t=2 | |
511 %! ref = [-1.205364552835178, 0.951542399860817]; | |
512 %! sol = ode23 (@fpol, [2 0 -2], fref); | |
513 %! assert ([sol.x(end), sol.y(end,:)], [-2, ref], 2e-2); | |
514 %!test ## Solve another anonymous function in backward direction | |
515 %! ref = [-1, 0.367879437558975]; | |
516 %! sol = ode23 (@(t,y) y, [0 -1], 1); | |
517 %! assert ([sol.x(end), sol.y(end,:)], ref, 1e-2); | |
518 %!test ## Solve another anonymous function below zero | |
519 %! ref = [0, 14.77810590694212]; | |
520 %! sol = ode23 (@(t,y) y, [-2 0], 2); | |
521 %! assert ([sol.x(end), sol.y(end,:)], ref, 1e-2); | |
522 %!test ## Solve in backward direction starting at t=0 with MaxStep option | |
523 %! ref = [-1.205364552835178, 0.951542399860817]; | |
524 %! opt = odeset ("MaxStep", 1e-3); | |
525 %! sol = ode23 (@fpol, [0 -2], [2 0], opt); | |
526 %! assert ([abs(sol.x(8)-sol.x(7))], [1e-3], 1e-3); | |
527 %! assert ([sol.x(end), sol.y(end,:)], [-2, ref], 1e-3); | |
528 %!test ## AbsTol option | |
529 %! opt = odeset ("AbsTol", 1e-5); | |
530 %! sol = ode23 (@fpol, [0 2], [2 0], opt); | |
531 %! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); | |
532 %!test ## AbsTol and RelTol option | |
533 %! opt = odeset ("AbsTol", 1e-8, "RelTol", 1e-8); | |
534 %! sol = ode23 (@fpol, [0 2], [2 0], opt); | |
535 %! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); | |
536 %!test ## RelTol and NormControl option -- higher accuracy | |
537 %! opt = odeset ("RelTol", 1e-8, "NormControl", "on"); | |
538 %! sol = ode23 (@fpol, [0 2], [2 0], opt); | |
539 %! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-4); | |
540 %!test ## Keeps initial values while integrating | |
541 %! opt = odeset ("NonNegative", 2); | |
542 %! sol = ode23 (@fpol, [0 2], [2 0], opt); | |
543 %! assert ([sol.x(end), sol.y(end,:)], [2, 2, 0], 1e-1); | |
544 %!test ## Details of OutputSel and Refine can't be tested | |
545 %! opt = odeset ("OutputFcn", @fout, "OutputSel", 1, "Refine", 5); | |
546 %! sol = ode23 (@fpol, [0 2], [2 0], opt); | |
547 %!test ## Stats must add further elements in sol | |
548 %! opt = odeset ("Stats", "on"); | |
549 %! sol = ode23 (@fpol, [0 2], [2 0], opt); | |
550 %! assert (isfield (sol, "stats")); | |
551 %! assert (isfield (sol.stats, "nsteps")); | |
552 %!test ## Events option add further elements in sol | |
553 %! opt = odeset ("Events", @feve); | |
554 %! sol = ode23 (@fpol, [0 10], [2 0], opt); | |
555 %! assert (isfield (sol, "ie")); | |
556 %! assert (sol.ie(1), 2); | |
557 %! assert (isfield (sol, "xe")); | |
558 %! assert (isfield (sol, "ye")); | |
559 %!test ## Events option, now stop integration | |
560 %! opt = odeset ("Events", @fevn, "NormControl", "on"); | |
561 %! sol = ode23 (@fpol, [0 10], [2 0], opt); | |
562 %! assert ([sol.ie, sol.xe, sol.ye], ... | |
563 %! [2.0, 2.496110, -0.830550, -2.677589], .5e-1); | |
564 %!test ## Events option, five output arguments | |
565 %! opt = odeset ("Events", @fevn, "NormControl", "on"); | |
566 %! [t, y, vxe, ye, vie] = ode23 (@fpol, [0 10], [2 0], opt); | |
567 %! assert ([vie, vxe, ye], ... | |
568 %! [2.0, 2.496110, -0.830550, -2.677589], 1e-1); | |
569 %!test ## Jacobian option | |
570 %! opt = odeset ("Jacobian", @fjac); | |
571 %! sol = ode23 (@fpol, [0 2], [2 0], opt); | |
572 %! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); | |
573 %!test ## Jacobian option and sparse return value | |
574 %! opt = odeset ("Jacobian", @fjcc); | |
575 %! sol = ode23 (@fpol, [0 2], [2 0], opt); | |
576 %! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); | |
577 %! | |
578 %! ## test for JPattern option is missing | |
579 %! ## test for Vectorized option is missing | |
580 %! | |
581 %!test ## Mass option as function | |
582 %! opt = odeset ("Mass", @fmas); | |
583 %! sol = ode23 (@fpol, [0 2], [2 0], opt); | |
584 %! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); | |
585 %!test ## Mass option as matrix | |
586 %! opt = odeset ("Mass", eye (2,2)); | |
587 %! sol = ode23 (@fpol, [0 2], [2 0], opt); | |
588 %! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); | |
589 %!test ## Mass option as sparse matrix | |
590 %! opt = odeset ("Mass", sparse (eye (2,2))); | |
591 %! sol = ode23 (@fpol, [0 2], [2 0], opt); | |
592 %! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); | |
593 %!test ## Mass option as function and sparse matrix | |
594 %! opt = odeset ("Mass", @fmsa); | |
595 %! sol = ode23 (@fpol, [0 2], [2 0], opt); | |
596 %! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); | |
597 %!test ## Mass option as function and MStateDependence | |
598 %! opt = odeset ("Mass", @fmas, "MStateDependence", "strong"); | |
599 %! sol = ode23 (@fpol, [0 2], [2 0], opt); | |
600 %! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); | |
601 %!test ## Set BDF option to something else than default | |
602 %! opt = odeset ("BDF", "on"); | |
603 %! [t, y] = ode23 (@fpol, [0 2], [2 0], opt); | |
604 %! assert ([t(end), y(end,:)], [2, fref], 1e-3); | |
605 %! | |
606 %! ## test for MvPattern option is missing | |
607 %! ## test for InitialSlope option is missing | |
608 %! ## test for MaxOrder option is missing | |
609 %! | |
610 %! warning ("on", "Octave:InvalidArgument"); | |
611 | |
612 ## Local Variables: *** | |
613 ## mode: octave *** | |
614 ## End: *** |