Mercurial > octave-nkf

comparison scripts/ode/ode45.m @ 20581:e368ce72a844

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Use Octave coding conventions for ode* functions. Use '!' instead of '~' for logical not. Use '##' only for comments which standalone on a line. Use double quotes instead of single quotes where possible. Wrap long lines to < 80 characters. * ode45.m, odeget.m, odeset.m: Use Octave coding conventions.
author Rik <rik@octave.org>
date Sat, 03 Oct 2015 22:10:45 -0700
parents 25623ef2ff4f
children 45151de7423f
comparison
equal deleted inserted replaced
20580:25623ef2ff4f 20581:e368ce72a844
70 ## @end group 70 ## @end group
71 ## @end example 71 ## @end example
72 ## @seealso{odeset, odeget} 72 ## @seealso{odeset, odeget}
73 ## @end deftypefn 73 ## @end deftypefn
74 74
75 function [varargout] = ode45 (vfun, vslot, vinit, varargin) 75 function varargout = ode45 (vfun, vslot, vinit, varargin)
76 76
77 vorder = 5; % runge_kutta_45_dorpri uses local extrapolation 77 vorder = 5; # runge_kutta_45_dorpri uses local extrapolation
78 vsolver = "ode45"; 78 vsolver = "ode45";
79 79
80 if (nargin == 0) ## Check number and types of all input arguments 80 if (nargin == 0) # Check number and types of all input arguments
81 help (vsolver); 81 help (vsolver);
82 error ("OdePkg:InvalidArgument", ... 82 error ("OdePkg:InvalidArgument", ...
83 "number of input arguments must be greater than zero"); 83 "number of input arguments must be greater than zero");
84 endif 84 endif
85 85
86 if (nargin < 3) 86 if (nargin < 3)
87 print_usage; 87 print_usage;
88 endif 88 endif
89 89
90 if (nargin >= 4) 90 if (nargin >= 4)
91 if (~isstruct (varargin{1})) 91 if (! isstruct (varargin{1}))
92 ## varargin{1:len} are parameters for vfun 92 ## varargin{1:len} are parameters for vfun
93 vodeoptions = odeset; 93 vodeoptions = odeset;
94 vodeoptions.vfunarguments = varargin; 94 vodeoptions.vfunarguments = varargin;
95 elseif (length (varargin) > 1) 95 elseif (length (varargin) > 1)
96 ## varargin{1} is an OdePkg options structure vopt 96 ## varargin{1} is an OdePkg options structure vopt
97 vodeoptions = odepkg_structure_check (varargin{1}, "ode45"); 97 vodeoptions = odepkg_structure_check (varargin{1}, "ode45");
98 vodeoptions.vfunarguments = {varargin{2:length(varargin)}}; 98 vodeoptions.vfunarguments = {varargin{2:length(varargin)}};
99 else ## if (isstruct (varargin{1})) 99 else # if (isstruct (varargin{1}))
100 vodeoptions = odepkg_structure_check (varargin{1}, "ode45"); 100 vodeoptions = odepkg_structure_check (varargin{1}, "ode45");
101 vodeoptions.vfunarguments = {}; 101 vodeoptions.vfunarguments = {};
102 endif 102 endif
103 else ## if (nargin == 3) 103 else # nargin == 3
104 vodeoptions = odeset; 104 vodeoptions = odeset;
105 vodeoptions.vfunarguments = {}; 105 vodeoptions.vfunarguments = {};
106 endif 106 endif
107 107
108 if (~isvector (vslot) || ~isnumeric (vslot)) 108 if (! isvector (vslot) || ! isnumeric (vslot))
109 error ("OdePkg:InvalidArgument", ... 109 error ("OdePkg:InvalidArgument", ...
110 "second input argument must be a valid vector"); 110 "second input argument must be a valid vector");
111 endif 111 endif
112 112
113 if (length (vslot) < 2 && ... 113 if (length (vslot) < 2 && ...
121 else 121 else
122 vodeoptions.vdirection = sign (vslot(2) - vslot(1)); 122 vodeoptions.vdirection = sign (vslot(2) - vslot(1));
123 endif 123 endif
124 vslot = vslot(:); 124 vslot = vslot(:);
125 125
126 if (~isvector (vinit) || ~isnumeric (vinit)) 126 if (! isvector (vinit) || ! isnumeric (vinit))
127 error ("OdePkg:InvalidArgument", ... 127 error ("OdePkg:InvalidArgument", ...
128 "third input argument must be a valid numerical value"); 128 "third input argument must be a valid numerical value");
129 endif 129 endif
130 vinit = vinit(:); 130 vinit = vinit(:);
131 131
132 if ~(isa (vfun, "function_handle") || isa (vfun, "inline")) 132 if (! (isa (vfun, "function_handle") || isa (vfun, "inline")))
133 error ("OdePkg:InvalidArgument", ... 133 error ("OdePkg:InvalidArgument", ...
134 "first input argument must be a valid function handle"); 134 "first input argument must be a valid function handle");
135 endif 135 endif
136 136
137 ## Start preprocessing, have a look which options are set in 137 ## Start preprocessing, have a look which options are set in vodeoptions,
138 ## vodeoptions, check if an invalid or unused option is set 138 ## check if an invalid or unused option is set
139 if (isempty (vodeoptions.TimeStepNumber) ... 139 if (isempty (vodeoptions.TimeStepNumber) ...
140 && isempty (vodeoptions.TimeStepSize)) 140 && isempty (vodeoptions.TimeStepSize))
141 integrate_func = "adaptive"; 141 integrate_func = "adaptive";
142 vodeoptions.vstepsizefixed = false; 142 vodeoptions.vstepsizefixed = false;
143 elseif (~isempty (vodeoptions.TimeStepNumber) ... 143 elseif (! isempty (vodeoptions.TimeStepNumber) ...
144 && ~isempty (vodeoptions.TimeStepSize)) 144 && ! isempty (vodeoptions.TimeStepSize))
145 integrate_func = "n_steps"; 145 integrate_func = "n_steps";
146 vodeoptions.vstepsizefixed = true; 146 vodeoptions.vstepsizefixed = true;
147 if (sign (vodeoptions.TimeStepSize) != vodeoptions.vdirection) 147 if (sign (vodeoptions.TimeStepSize) != vodeoptions.vdirection)
148 warning ("OdePkg:InvalidArgument", ... 148 warning ("OdePkg:InvalidArgument", ...
149 "option ''TimeStepSize'' has a wrong sign", ... 149 "option 'TimeStepSize' has a wrong sign", ...
150 "it will be corrected automatically"); 150 "it will be corrected automatically");
151 vodeoptions.TimeStepSize = (-1)*vodeoptions.TimeStepSize; 151 vodeoptions.TimeStepSize = (-1)*vodeoptions.TimeStepSize;
152 endif 152 endif
153 elseif (isempty (vodeoptions.TimeStepNumber) && ~isempty (vodeoptions.TimeStepSize)) 153 elseif (isempty (vodeoptions.TimeStepNumber) ...
154 && ! isempty (vodeoptions.TimeStepSize))
154 integrate_func = "const"; 155 integrate_func = "const";
155 vodeoptions.vstepsizefixed = true; 156 vodeoptions.vstepsizefixed = true;
156 if (sign (vodeoptions.TimeStepSize) != vodeoptions.vdirection) 157 if (sign (vodeoptions.TimeStepSize) != vodeoptions.vdirection)
157 warning ("OdePkg:InvalidArgument", ... 158 warning ("OdePkg:InvalidArgument", ...
158 "option ''TimeStepSize'' has a wrong sign", ... 159 "option 'TimeStepSize' has a wrong sign", ...
159 "it will be corrected automatically"); 160 "it will be corrected automatically");
160 vodeoptions.TimeStepSize = (-1)*vodeoptions.TimeStepSize; 161 vodeoptions.TimeStepSize = (-1)*vodeoptions.TimeStepSize;
161 endif 162 endif
162 else 163 else
163 warning ("OdePkg:InvalidArgument", ... 164 warning ("OdePkg:InvalidArgument", ...
166 endif 167 endif
167 168
168 ## Get the default options that can be set with "odeset" temporarily 169 ## Get the default options that can be set with "odeset" temporarily
169 vodetemp = odeset; 170 vodetemp = odeset;
170 171
171 ## Implementation of the option RelTol has been finished. This option 172 ## Implementation of the option RelTol has been finished.
172 ## can be set by the user to another value than default value. 173 ## This option can be set by the user to another value than default value.
173 if (isempty (vodeoptions.RelTol) && ~vodeoptions.vstepsizefixed) 174 if (isempty (vodeoptions.RelTol) && ! vodeoptions.vstepsizefixed)
174 vodeoptions.RelTol = 1e-3; 175 vodeoptions.RelTol = 1e-3;
175 warning ("OdePkg:InvalidArgument", ... 176 warning ("OdePkg:InvalidArgument", ...
176 "Option ''RelTol'' not set, new value %f is used", ... 177 "Option 'RelTol' not set, new value %f is used", ...
177 vodeoptions.RelTol); 178 vodeoptions.RelTol);
178 elseif (~isempty (vodeoptions.RelTol) && vodeoptions.vstepsizefixed) 179 elseif (! isempty (vodeoptions.RelTol) && vodeoptions.vstepsizefixed)
179 warning ("OdePkg:InvalidArgument", ... 180 warning ("OdePkg:InvalidArgument", ...
180 "Option ''RelTol'' will be ignored if fixed time stamps are given"); 181 "Option 'RelTol' will be ignored if fixed time stamps are given");
181 endif 182 endif
182 183
183 ## Implementation of the option AbsTol has been finished. This option 184 ## Implementation of the option AbsTol has been finished.
184 ## can be set by the user to another value than default value. 185 ## This option can be set by the user to another value than default value.
185 if (isempty (vodeoptions.AbsTol) && ~vodeoptions.vstepsizefixed) 186 if (isempty (vodeoptions.AbsTol) && ! vodeoptions.vstepsizefixed)
186 vodeoptions.AbsTol = 1e-6; 187 vodeoptions.AbsTol = 1e-6;
187 warning ("OdePkg:InvalidArgument", ... 188 warning ("OdePkg:InvalidArgument", ...
188 "Option ''AbsTol'' not set, new value %f is used", ... 189 "Option 'AbsTol' not set, new value %f is used", ...
189 vodeoptions.AbsTol); 190 vodeoptions.AbsTol);
190 elseif (~isempty (vodeoptions.AbsTol) && vodeoptions.vstepsizefixed) 191 elseif (! isempty (vodeoptions.AbsTol) && vodeoptions.vstepsizefixed)
191 warning ("OdePkg:InvalidArgument", ... 192 warning ("OdePkg:InvalidArgument", ...
192 "Option ''AbsTol'' will be ignored if fixed time stamps are given"); 193 "Option 'AbsTol' will be ignored if fixed time stamps are given");
193 else 194 else
194 vodeoptions.AbsTol = vodeoptions.AbsTol(:); ## Create column vector 195 vodeoptions.AbsTol = vodeoptions.AbsTol(:); # Create column vector
195 endif 196 endif
196 197
197 ## Implementation of the option NormControl has been finished. This 198 ## Implementation of the option NormControl has been finished.
198 ## option can be set by the user to another value than default value. 199 ## This option can be set by the user to another value than default value.
199 if (strcmp (vodeoptions.NormControl, "on")) 200 if (strcmp (vodeoptions.NormControl, "on"))
200 vodeoptions.vnormcontrol = true; 201 vodeoptions.vnormcontrol = true;
201 else 202 else
202 vodeoptions.vnormcontrol = false; 203 vodeoptions.vnormcontrol = false;
203 endif 204 endif
204 205
205 ## Implementation of the option NonNegative has been finished. This 206 ## Implementation of the option NonNegative has been finished.
206 ## option can be set by the user to another value than default value. 207 ## This option can be set by the user to another value than default value.
207 if (~isempty (vodeoptions.NonNegative)) 208 if (! isempty (vodeoptions.NonNegative))
208 if (isempty (vodeoptions.Mass)) 209 if (isempty (vodeoptions.Mass))
209 vodeoptions.vhavenonnegative = true; 210 vodeoptions.vhavenonnegative = true;
210 else 211 else
211 vodeoptions.vhavenonnegative = false; 212 vodeoptions.vhavenonnegative = false;
212 warning ("OdePkg:InvalidArgument", ... 213 warning ("OdePkg:InvalidArgument", ...
214 endif 215 endif
215 else 216 else
216 vodeoptions.vhavenonnegative = false; 217 vodeoptions.vhavenonnegative = false;
217 endif 218 endif
218 219
219 ## Implementation of the option OutputFcn has been finished. This 220 ## Implementation of the option OutputFcn has been finished.
220 ## option can be set by the user to another value than default value. 221 ## This option can be set by the user to another value than default value.
221 if (isempty (vodeoptions.OutputFcn) && nargout == 0) 222 if (isempty (vodeoptions.OutputFcn) && nargout == 0)
222 vodeoptions.OutputFcn = @odeplot; 223 vodeoptions.OutputFcn = @odeplot;
223 vodeoptions.vhaveoutputfunction = true; 224 vodeoptions.vhaveoutputfunction = true;
224 elseif (isempty (vodeoptions.OutputFcn)) 225 elseif (isempty (vodeoptions.OutputFcn))
225 vodeoptions.vhaveoutputfunction = false; 226 vodeoptions.vhaveoutputfunction = false;
226 else 227 else
227 vodeoptions.vhaveoutputfunction = true; 228 vodeoptions.vhaveoutputfunction = true;
228 endif 229 endif
229 230
230 ## Implementation of the option OutputSel has been finished. This 231 ## Implementation of the option OutputSel has been finished.
231 ## option can be set by the user to another value than default value. 232 ## This option can be set by the user to another value than default value.
232 if (~isempty (vodeoptions.OutputSel)) 233 if (! isempty (vodeoptions.OutputSel))
233 vodeoptions.vhaveoutputselection = true; 234 vodeoptions.vhaveoutputselection = true;
234 else 235 else
235 vodeoptions.vhaveoutputselection = false; 236 vodeoptions.vhaveoutputselection = false;
236 endif 237 endif
237 238
238 ## Implementation of the option OutputSave has been finished. This 239 ## Implementation of the option OutputSave has been finished.
239 ## option can be set by the user to another value than default value. 240 ## This option can be set by the user to another value than default value.
240 if (isempty (vodeoptions.OutputSave)) 241 if (isempty (vodeoptions.OutputSave))
241 vodeoptions.OutputSave = 1; 242 vodeoptions.OutputSave = 1;
242 endif 243 endif
243 244
244 ## Implementation of the option Refine has been finished. This option 245 ## Implementation of the option Refine has been finished.
245 ## can be set by the user to another value than default value. 246 ## This option can be set by the user to another value than default value.
246 if (vodeoptions.Refine > 0) 247 if (vodeoptions.Refine > 0)
247 vodeoptions.vhaverefine = true; 248 vodeoptions.vhaverefine = true;
248 else 249 else
249 vodeoptions.vhaverefine = false; 250 vodeoptions.vhaverefine = false;
250 endif 251 endif
251 252
252 ## Implementation of the option Stats has been finished. This option 253 ## Implementation of the option Stats has been finished.
253 ## can be set by the user to another value than default value. 254 ## This option can be set by the user to another value than default value.
254 255
255 ## Implementation of the option InitialStep has been finished. This 256 ## Implementation of the option InitialStep has been finished.
256 ## option can be set by the user to another value than default value. 257 ## This option can be set by the user to another value than default value.
257 if (isempty (vodeoptions.InitialStep) && strcmp (integrate_func, "adaptive")) 258 if (isempty (vodeoptions.InitialStep) && strcmp (integrate_func, "adaptive"))
258 vodeoptions.InitialStep = vodeoptions.vdirection* ... 259 vodeoptions.InitialStep = ...
259 starting_stepsize (vorder, vfun, vslot(1), vinit, vodeoptions.AbsTol, ... 260 vodeoptions.vdirection * starting_stepsize (vorder, vfun, vslot(1), vinit,
260 vodeoptions.RelTol, vodeoptions.vnormcontrol); 261 vodeoptions.AbsTol,
261 warning ("OdePkg:InvalidArgument", ... 262 vodeoptions.RelTol,
262 "option ''InitialStep'' not set, estimated value %f is used", ... 263 vodeoptions.vnormcontrol);
264 warning ("OdePkg:InvalidArgument", ...
265 "option 'InitialStep' not set, estimated value %f is used", ...
263 vodeoptions.InitialStep); 266 vodeoptions.InitialStep);
264 elseif(isempty (vodeoptions.InitialStep)) 267 elseif(isempty (vodeoptions.InitialStep))
265 vodeoptions.InitialStep = odeget (vodeoptions, "TimeStepSize"); 268 vodeoptions.InitialStep = odeget (vodeoptions, "TimeStepSize");
266 endif 269 endif
267 270
268 ## Implementation of the option MaxStep has been finished. This option 271 ## Implementation of the option MaxStep has been finished. This option
269 ## can be set by the user to another value than default value. 272 ## can be set by the user to another value than default value.
270 if (isempty (vodeoptions.MaxStep) && ~vodeoptions.vstepsizefixed) 273 if (isempty (vodeoptions.MaxStep) && ! vodeoptions.vstepsizefixed)
271 vodeoptions.MaxStep = abs (vslot(end) - vslot(1)) / 10; 274 vodeoptions.MaxStep = abs (vslot(end) - vslot(1)) / 10;
272 warning ("OdePkg:InvalidArgument", ... 275 warning ("OdePkg:InvalidArgument", ...
273 "Option ''MaxStep'' not set, new value %f is used", ... 276 "Option 'MaxStep' not set, new value %f is used", ...
274 vodeoptions.MaxStep); 277 vodeoptions.MaxStep);
275 endif 278 endif
276 279
277 ## Implementation of the option Events has been finished. This option 280 ## Implementation of the option Events has been finished.
278 ## can be set by the user to another value than default value. 281 ## This option can be set by the user to another value than default value.
279 if (~isempty (vodeoptions.Events)) 282 if (! isempty (vodeoptions.Events))
280 vodeoptions.vhaveeventfunction = true; 283 vodeoptions.vhaveeventfunction = true;
281 else 284 else
282 vodeoptions.vhaveeventfunction = false; 285 vodeoptions.vhaveeventfunction = false;
283 endif 286 endif
284 287
285 ## The options 'Jacobian', 'JPattern' and 'Vectorized' will be ignored 288 ## The options "Jacobian", "JPattern" and "Vectorized" will be ignored
286 ## by this solver because this solver uses an explicit Runge-Kutta 289 ## by this solver because this solver uses an explicit Runge-Kutta method
287 ## method and therefore no Jacobian calculation is necessary 290 ## and therefore no Jacobian calculation is necessary.
288 if (~isequal (vodeoptions.Jacobian, vodetemp.Jacobian)) 291 if (! isequal (vodeoptions.Jacobian, vodetemp.Jacobian))
289 warning ("OdePkg:InvalidArgument", ... 292 warning ("OdePkg:InvalidArgument", ...
290 "option ''Jacobian'' will be ignored by this solver"); 293 "option 'Jacobian' will be ignored by this solver");
291 endif 294 endif
292 295
293 if (~isequal (vodeoptions.JPattern, vodetemp.JPattern)) 296 if (! isequal (vodeoptions.JPattern, vodetemp.JPattern))
294 warning ("OdePkg:InvalidArgument", ... 297 warning ("OdePkg:InvalidArgument", ...
295 "option ''JPattern'' will be ignored by this solver"); 298 "option 'JPattern' will be ignored by this solver");
296 endif 299 endif
297 300
298 if (~isequal (vodeoptions.Vectorized, vodetemp.Vectorized)) 301 if (! isequal (vodeoptions.Vectorized, vodetemp.Vectorized))
299 warning ("OdePkg:InvalidArgument", ... 302 warning ("OdePkg:InvalidArgument", ...
300 "option ''Vectorized'' will be ignored by this solver"); 303 "option 'Vectorized' will be ignored by this solver");
301 endif 304 endif
302 305
303 if (~isequal (vodeoptions.NewtonTol, vodetemp.NewtonTol)) 306 if (! isequal (vodeoptions.NewtonTol, vodetemp.NewtonTol))
304 warning ("OdePkg:InvalidArgument", ... 307 warning ("OdePkg:InvalidArgument", ...
305 "option ''NewtonTol'' will be ignored by this solver"); 308 "option 'NewtonTol' will be ignored by this solver");
306 endif 309 endif
307 310
308 if (~isequal (vodeoptions.MaxNewtonIterations,vodetemp.MaxNewtonIterations)) 311 if (! isequal (vodeoptions.MaxNewtonIterations,vodetemp.MaxNewtonIterations))
309 warning ("OdePkg:InvalidArgument", ... 312 warning ("OdePkg:InvalidArgument", ...
310 "option ''MaxNewtonIterations'' will be ignored by this solver"); 313 "option 'MaxNewtonIterations' will be ignored by this solver");
311 endif 314 endif
312 315
313 ## Implementation of the option Mass has been finished. This option 316 ## Implementation of the option Mass has been finished.
314 ## can be set by the user to another value than default value. 317 ## This option can be set by the user to another value than default value.
315 if (~isempty (vodeoptions.Mass) && isnumeric (vodeoptions.Mass)) 318 if (! isempty (vodeoptions.Mass) && isnumeric (vodeoptions.Mass))
316 vhavemasshandle = false; 319 vhavemasshandle = false;
317 vmass = vodeoptions.Mass; ## constant mass 320 vmass = vodeoptions.Mass; # constant mass
318 elseif (isa (vodeoptions.Mass, "function_handle")) 321 elseif (isa (vodeoptions.Mass, "function_handle"))
319 vhavemasshandle = true; ## mass defined by a function handle 322 vhavemasshandle = true; # mass defined by a function handle
320 else ## no mass matrix - creating a diag-matrix of ones for mass 323 else # no mass matrix - creating a diag-matrix of ones for mass
321 vhavemasshandle = false; ## vmass = diag (ones (length (vinit), 1), 0); 324 vhavemasshandle = false; # vmass = diag (ones (length (vinit), 1), 0);
322 endif 325 endif
323 326
324 ## Implementation of the option MStateDependence has been finished. 327 ## Implementation of the option MStateDependence has been finished.
325 ## This option can be set by the user to another value than default 328 ## This option can be set by the user to another value than default value.
326 ## value.
327 if (strcmp (vodeoptions.MStateDependence, "none")) 329 if (strcmp (vodeoptions.MStateDependence, "none"))
328 vmassdependence = false; 330 vmassdependence = false;
329 else 331 else
330 vmassdependence = true; 332 vmassdependence = true;
331 endif 333 endif
332 334
333 ## Other options that are not used by this solver. Print a warning 335 ## Other options that are not used by this solver.
334 ## message to tell the user that the option(s) is/are ignored. 336 ## Print a warning message to tell the user that the option is ignored.
335 if (~isequal (vodeoptions.MvPattern, vodetemp.MvPattern)) 337 if (! isequal (vodeoptions.MvPattern, vodetemp.MvPattern))
336 warning ("OdePkg:InvalidArgument", ... 338 warning ("OdePkg:InvalidArgument", ...
337 "option ''MvPattern'' will be ignored by this solver"); 339 "option 'MvPattern' will be ignored by this solver");
338 endif 340 endif
339 341
340 if (~isequal (vodeoptions.MassSingular, vodetemp.MassSingular)) 342 if (! isequal (vodeoptions.MassSingular, vodetemp.MassSingular))
341 warning ("OdePkg:InvalidArgument", ... 343 warning ("OdePkg:InvalidArgument", ...
342 "option ''MassSingular'' will be ignored by this solver"); 344 "option 'MassSingular' will be ignored by this solver");
343 endif 345 endif
344 346
345 if (~isequal (vodeoptions.InitialSlope, vodetemp.InitialSlope)) 347 if (! isequal (vodeoptions.InitialSlope, vodetemp.InitialSlope))
346 warning ("OdePkg:InvalidArgument", ... 348 warning ("OdePkg:InvalidArgument", ...
347 "option ''InitialSlope'' will be ignored by this solver"); 349 "option 'InitialSlope' will be ignored by this solver");
348 endif 350 endif
349 351
350 if (~isequal (vodeoptions.MaxOrder, vodetemp.MaxOrder)) 352 if (! isequal (vodeoptions.MaxOrder, vodetemp.MaxOrder))
351 warning ("OdePkg:InvalidArgument", ... 353 warning ("OdePkg:InvalidArgument", ...
352 "option ''MaxOrder'' will be ignored by this solver"); 354 "option 'MaxOrder' will be ignored by this solver");
353 endif 355 endif
354 356
355 if (~isequal (vodeoptions.BDF, vodetemp.BDF)) 357 if (! isequal (vodeoptions.BDF, vodetemp.BDF))
356 warning ("OdePkg:InvalidArgument", ... 358 warning ("OdePkg:InvalidArgument", ...
357 "option ''BDF'' will be ignored by this solver"); 359 "option 'BDF' will be ignored by this solver");
358 endif 360 endif
359 361
360 ## Starting the initialisation of the core solver ode45 362 ## Starting the initialisation of the core solver ode45
361 SubOpts = vodeoptions; 363 SubOpts = vodeoptions;
362 364
363 if (vhavemasshandle) ## Handle only the dynamic mass matrix, 365 if (vhavemasshandle) # Handle only the dynamic mass matrix,
364 if (vmassdependence) ## constant mass matrices have already 366 if (vmassdependence) # constant mass matrices have already
365 vmass = @(t,x) vodeoptions.Mass (t, x, vodeoptions.vfunarguments{:}); 367 vmass = @(t,x) vodeoptions.Mass (t, x, vodeoptions.vfunarguments{:});
366 vfun = @(t,x) vmass (t, x, vodeoptions.vfunarguments{:}) ... 368 vfun = @(t,x) vmass (t, x, vodeoptions.vfunarguments{:}) ...
367 \ vfun (t, x, vodeoptions.vfunarguments{:}); 369 \ vfun (t, x, vodeoptions.vfunarguments{:});
368 else ## if (vmassdependence == false) 370 else # if (vmassdependence == false)
369 vmass = @(t) vodeoptions.Mass (t, vodeoptions.vfunarguments{:}); 371 vmass = @(t) vodeoptions.Mass (t, vodeoptions.vfunarguments{:});
370 vfun = @(t,x) vmass (t, vodeoptions.vfunarguments{:}) ... 372 vfun = @(t,x) vmass (t, vodeoptions.vfunarguments{:}) ...
371 \ vfun (t, x, vodeoptions.vfunarguments{:}); 373 \ vfun (t, x, vodeoptions.vfunarguments{:});
372 endif 374 endif
373 endif 375 endif
386 vfun, vslot, vinit, ... 388 vfun, vslot, vinit, ...
387 vodeoptions.TimeStepSize, SubOpts); 389 vodeoptions.TimeStepSize, SubOpts);
388 endswitch 390 endswitch
389 391
390 ## Postprocessing, do whatever when terminating integration algorithm 392 ## Postprocessing, do whatever when terminating integration algorithm
391 if (vodeoptions.vhaveoutputfunction) ## Cleanup plotter 393 if (vodeoptions.vhaveoutputfunction) # Cleanup plotter
392 feval (vodeoptions.OutputFcn, solution.t(end), ... 394 feval (vodeoptions.OutputFcn, solution.t(end), ...
393 solution.x(end,:)', "done", vodeoptions.vfunarguments{:}); 395 solution.x(end,:)', "done", vodeoptions.vfunarguments{:});
394 endif 396 endif
395 if (vodeoptions.vhaveeventfunction) ## Cleanup event function handling 397 if (vodeoptions.vhaveeventfunction) # Cleanup event function handling
396 odepkg_event_handle (vodeoptions.Events, solution.t(end), ... 398 odepkg_event_handle (vodeoptions.Events, solution.t(end), ...
397 solution.x(end,:)', "done", vodeoptions.vfunarguments{:}); 399 solution.x(end,:)', "done",
400 vodeoptions.vfunarguments{:});
398 endif 401 endif
399 402
400 ## Print additional information if option Stats is set 403 ## Print additional information if option Stats is set
401 if (strcmp (vodeoptions.Stats, "on")) 404 if (strcmp (vodeoptions.Stats, "on"))
402 vhavestats = true; 405 vhavestats = true;
403 vnsteps = solution.vcntloop-2; ## vcntloop from 2..end 406 vnsteps = solution.vcntloop-2; # vcntloop from 2..end
404 vnfailed = (solution.vcntcycles-1)-(solution.vcntloop-2)+1; ## vcntcycl from 1..end 407 vnfailed = (solution.vcntcycles-1)-(solution.vcntloop-2)+1; # vcntcycl from 1..end
405 vnfevals = 7*(solution.vcntcycles-1); ## number of ode evaluations 408 vnfevals = 7*(solution.vcntcycles-1); # number of ode evaluations
406 vndecomps = 0; ## number of LU decompositions 409 vndecomps = 0; # number of LU decompositions
407 vnpds = 0; ## number of partial derivatives 410 vnpds = 0; # number of partial derivatives
408 vnlinsols = 0; ## no. of solutions of linear systems 411 vnlinsols = 0; # no. of linear systems solutions
409 ## Print cost statistics if no output argument is given 412 ## Print cost statistics if no output argument is given
410 if (nargout == 0) 413 if (nargout == 0)
411 vmsg = fprintf (1, "Number of successful steps: %d\n", vnsteps); 414 vmsg = fprintf (1, "Number of successful steps: %d\n", vnsteps);
412 vmsg = fprintf (1, "Number of failed attempts: %d\n", vnfailed); 415 vmsg = fprintf (1, "Number of failed attempts: %d\n", vnfailed);
413 vmsg = fprintf (1, "Number of function calls: %d\n", vnfevals); 416 vmsg = fprintf (1, "Number of function calls: %d\n", vnfevals);
414 endif 417 endif
415 else 418 else
416 vhavestats = false; 419 vhavestats = false;
417 endif 420 endif
418 421
419 if (nargout == 1) ## Sort output variables, depends on nargout 422 if (nargout == 1) # Sort output variables, depends on nargout
420 varargout{1}.x = solution.t; ## Time stamps are saved in field x 423 varargout{1}.x = solution.t; # Time stamps are saved in field x
421 varargout{1}.y = solution.x; ## Results are saved in field y 424 varargout{1}.y = solution.x; # Results are saved in field y
422 varargout{1}.solver = vsolver; ## Solver name is saved in field solver 425 varargout{1}.solver = vsolver; # Solver name is saved in field solver
423 if (vodeoptions.vhaveeventfunction) 426 if (vodeoptions.vhaveeventfunction)
424 varargout{1}.ie = solution.vevent{2}; ## Index info which event occured 427 varargout{1}.ie = solution.vevent{2}; # Index info which event occurred
425 varargout{1}.xe = solution.vevent{3}; ## Time info when an event occured 428 varargout{1}.xe = solution.vevent{3}; # Time info when an event occurred
426 varargout{1}.ye = solution.vevent{4}; ## Results when an event occured 429 varargout{1}.ye = solution.vevent{4}; # Results when an event occurred
427 endif 430 endif
428 if (vhavestats) 431 if (vhavestats)
429 varargout{1}.stats = struct; 432 varargout{1}.stats = struct;
430 varargout{1}.stats.nsteps = vnsteps; 433 varargout{1}.stats.nsteps = vnsteps;
431 varargout{1}.stats.nfailed = vnfailed; 434 varargout{1}.stats.nfailed = vnfailed;
433 varargout{1}.stats.npds = vnpds; 436 varargout{1}.stats.npds = vnpds;
434 varargout{1}.stats.ndecomps = vndecomps; 437 varargout{1}.stats.ndecomps = vndecomps;
435 varargout{1}.stats.nlinsols = vnlinsols; 438 varargout{1}.stats.nlinsols = vnlinsols;
436 endif 439 endif
437 elseif (nargout == 2) 440 elseif (nargout == 2)
438 varargout{1} = solution.t; ## Time stamps are first output argument 441 varargout{1} = solution.t; # Time stamps are first output argument
439 varargout{2} = solution.x; ## Results are second output argument 442 varargout{2} = solution.x; # Results are second output argument
440 elseif (nargout == 5) 443 elseif (nargout == 5)
441 varargout{1} = solution.t; ## Same as (nargout == 2) 444 varargout{1} = solution.t; # Same as (nargout == 2)
442 varargout{2} = solution.x; ## Same as (nargout == 2) 445 varargout{2} = solution.x; # Same as (nargout == 2)
443 varargout{3} = []; ## LabMat doesn't accept lines like 446 varargout{3} = []; # LabMat doesn't accept lines like
444 varargout{4} = []; ## varargout{3} = varargout{4} = []; 447 varargout{4} = []; # varargout{3} = varargout{4} = [];
445 varargout{5} = []; 448 varargout{5} = [];
446 if (vodeoptions.vhaveeventfunction) 449 if (vodeoptions.vhaveeventfunction)
447 varargout{3} = solution.vevent{3}; ## Time info when an event occured 450 varargout{3} = solution.vevent{3}; # Time info when an event occurred
448 varargout{4} = solution.vevent{4}; ## Results when an event occured 451 varargout{4} = solution.vevent{4}; # Results when an event occurred
449 varargout{5} = solution.vevent{2}; ## Index info which event occured 452 varargout{5} = solution.vevent{2}; # Index info which event occurred
450 endif 453 endif
451 endif 454 endif
452 455
453 endfunction 456 endfunction
454 457
455 %! # We are using the "Van der Pol" implementation for all tests that 458
456 %! # are done for this function. 459 ## We are using the "Van der Pol" implementation for all tests that are done
457 %! # For further tests we also define a reference solution (computed at high accuracy) 460 ## for this function.
458 %!function [ydot] = fpol (vt, vy) ## The Van der Pol 461 ## For further tests we also define a reference solution (computed at high
459 %! ydot = [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)]; 462 ## accuracy)
460 %!endfunction 463 %!function [ydot] = fpol (vt, vy) # The Van der Pol
461 %!function [vref] = fref () ## The computed reference sol 464 %! ydot = [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
462 %! vref = [0.32331666704577, -1.83297456798624]; 465 %!endfunction
463 %!endfunction 466 %!function [vref] = fref () # The computed reference solution
464 %!function [vjac] = fjac (vt, vy, varargin) ## its Jacobian 467 %! vref = [0.32331666704577, -1.83297456798624];
465 %! vjac = [0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2]; 468 %!endfunction
466 %!function [vjac] = fjcc (vt, vy, varargin) ## sparse type 469 %!function [vjac] = fjac (vt, vy, varargin) # its Jacobian
467 %! vjac = sparse ([0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2]) 470 %! vjac = [0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2];
471 %!endfunction
472 %!function [vjac] = fjcc (vt, vy, varargin) # sparse type
473 %! vjac = sparse ([0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2])
474 %!endfunction
468 %!function [vval, vtrm, vdir] = feve (vt, vy, varargin) 475 %!function [vval, vtrm, vdir] = feve (vt, vy, varargin)
469 %! vval = fpol (vt, vy, varargin); ## We use the derivatives 476 %! vval = fpol (vt, vy, varargin); # We use the derivatives
470 %! vtrm = zeros (2,1); ## that's why component 2 477 %! vtrm = zeros (2,1); # that's why component 2
471 %! vdir = ones (2,1); ## seems to not be exact 478 %! vdir = ones (2,1); # seems to not be exact
479 %!endfunction
472 %!function [vval, vtrm, vdir] = fevn (vt, vy, varargin) 480 %!function [vval, vtrm, vdir] = fevn (vt, vy, varargin)
473 %! vval = fpol (vt, vy, varargin); ## We use the derivatives 481 %! vval = fpol (vt, vy, varargin); # We use the derivatives
474 %! vtrm = ones (2,1); ## that's why component 2 482 %! vtrm = ones (2,1); # that's why component 2
475 %! vdir = ones (2,1); ## seems to not be exact 483 %! vdir = ones (2,1); # seems to not be exact
484 %!endfunction
476 %!function [vmas] = fmas (vt, vy, varargin) 485 %!function [vmas] = fmas (vt, vy, varargin)
477 %! vmas = [1, 0; 0, 1]; ## Dummy mass matrix for tests 486 %! vmas = [1, 0; 0, 1]; # Dummy mass matrix for tests
487 %!endfunction
478 %!function [vmas] = fmsa (vt, vy, varargin) 488 %!function [vmas] = fmsa (vt, vy, varargin)
479 %! vmas = sparse ([1, 0; 0, 1]); ## A sparse dummy matrix 489 %! vmas = sparse ([1, 0; 0, 1]); # A sparse dummy matrix
490 %!endfunction
480 %!function [vout] = fout (vt, vy, vflag, varargin) 491 %!function [vout] = fout (vt, vy, vflag, varargin)
481 %! if (regexp (char (vflag), 'init') == 1) 492 %! if (regexp (char (vflag), 'init') == 1)
482 %! if (any (size (vt) ~= [2, 1])) error ('"fout" step "init"'); end 493 %! if (any (size (vt) != [2, 1])) error ('"fout" step "init"'); endif
483 %! elseif (isempty (vflag)) 494 %! elseif (isempty (vflag))
484 %! if (any (size (vt) ~= [1, 1])) error ('"fout" step "calc"'); end 495 %! if (any (size (vt) != [1, 1])) error ('"fout" step "calc"'); endif
485 %! vout = false; 496 %! vout = false;
486 %! elseif (regexp (char (vflag), 'done') == 1) 497 %! elseif (regexp (char (vflag), 'done') == 1)
487 %! if (any (size (vt) ~= [1, 1])) error ('"fout" step "done"'); end 498 %! if (any (size (vt) != [1, 1])) error ('"fout" step "done"'); endif
488 %! else error ('"fout" invalid vflag'); 499 %! else
489 %! end 500 %! error ('"fout" invalid vflag');
501 %! endif
502 %!endfunction
490 %! 503 %!
491 %! ## Turn off output of warning messages for all tests, turn them on 504 %! ## Turn off output of warning messages for all tests, turn them on
492 %! ## again if the last test is called 505 %! ## again if the last test is called
493 %!error ## ouput argument 506 %!error # ouput argument
494 %! warning ('off', 'OdePkg:InvalidArgument'); 507 %! warning ("off", "OdePkg:InvalidArgument");
495 %! B = ode45 (1, [0 25], [3 15 1]); 508 %! B = ode45 (1, [0 25], [3 15 1]);
496 %!error ## input argument number one 509 %!error # input argument number one
497 %! [vt, vy] = ode45 (1, [0 25], [3 15 1]); 510 %! [vt, vy] = ode45 (1, [0 25], [3 15 1]);
498 %!error ## input argument number two 511 %!error # input argument number two
499 %! [vt, vy] = ode45 (@fpol, 1, [3 15 1]); 512 %! [vt, vy] = ode45 (@fpol, 1, [3 15 1]);
500 %!test ## two output arguments 513 %!test # two output arguments
501 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0]); 514 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0]);
502 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-2); 515 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-2);
503 %!test ## not too many steps 516 %!test # not too many steps
504 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0]); 517 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0]);
505 %! assert (size (vt) < 20); 518 %! assert (size (vt) < 20);
506 %!test ## anonymous function instead of real function 519 %!test # anonymous function instead of real function
507 %! fvdb = @(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)]; 520 %! fvdb = @(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
508 %! [vt, vy] = ode45 (fvdb, [0 2], [2 0]); 521 %! [vt, vy] = ode45 (fvdb, [0 2], [2 0]);
509 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-2); 522 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-2);
510 %!test ## extra input arguments passed through 523 %!test # extra input arguments passed through
511 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0], 12, 13, 'KL'); 524 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0], 12, 13, "KL");
512 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-2); 525 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-2);
513 %!test ## empty OdePkg structure *but* extra input arguments 526 %!test # empty OdePkg structure *but* extra input arguments
514 %! vopt = odeset; 527 %! vopt = odeset;
515 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt, 12, 13, 'KL'); 528 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt, 12, 13, "KL");
516 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-2); 529 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-2);
517 %!error ## strange OdePkg structure 530 %!error # strange OdePkg structure
518 %! vopt = struct ('foo', 1); 531 %! vopt = struct ("foo", 1);
519 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt); 532 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt);
520 %!test ## Solve vdp in fixed step sizes 533 %!test # Solve vdp in fixed step sizes
521 %! vopt = odeset('TimeStepSize', 0.1); 534 %! vopt = odeset("TimeStepSize", 0.1);
522 %! [vt, vy] = ode45 (@fpol, [0,2], [2 0], vopt); 535 %! [vt, vy] = ode45 (@fpol, [0,2], [2 0], vopt);
523 %! assert (vt(:), [0:0.1:2]', 1e-2); 536 %! assert (vt(:), [0:0.1:2]', 1e-2);
524 %!test ## Solve another anonymous function below zero 537 %!test # Solve another anonymous function below zero
525 %! vref = [0, 14.77810590694212]; 538 %! vref = [0, 14.77810590694212];
526 %! [vt, vy] = ode45 (@(t,y) y, [-2 0], 2); 539 %! [vt, vy] = ode45 (@(t,y) y, [-2 0], 2);
527 %! assert ([vt(end), vy(end,:)], vref, 1e-1); 540 %! assert ([vt(end), vy(end,:)], vref, 1e-1);
528 %!test ## InitialStep option 541 %!test # InitialStep option
529 %! vopt = odeset ('InitialStep', 1e-8); 542 %! vopt = odeset ("InitialStep", 1e-8);
530 %! [vt, vy] = ode45 (@fpol, [0 0.2], [2 0], vopt); 543 %! [vt, vy] = ode45 (@fpol, [0 0.2], [2 0], vopt);
531 %! assert ([vt(2)-vt(1)], [1e-8], 1e-9); 544 %! assert ([vt(2)-vt(1)], [1e-8], 1e-9);
532 %!test ## MaxStep option 545 %!test # MaxStep option
533 %! vopt = odeset ('MaxStep', 1e-3); 546 %! vopt = odeset ("MaxStep", 1e-3);
534 %! vsol = ode45 (@fpol, [0 0.2], [2 0], vopt); 547 %! vsol = ode45 (@fpol, [0 0.2], [2 0], vopt);
535 %! assert ([vsol.x(5)-vsol.x(4)], [1e-3], 1e-3); 548 %! assert ([vsol.x(5)-vsol.x(4)], [1e-3], 1e-3);
536 %!test ## Solve with intermidiate step 549 %!test # Solve with intermidiate step
537 %! vsol = ode45 (@fpol, [0 1 2], [2 0]); 550 %! vsol = ode45 (@fpol, [0 1 2], [2 0]);
538 %! assert (any((vsol.x-1) == 0)); 551 %! assert (any((vsol.x-1) == 0));
539 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); 552 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
540 %!test ## Solve in backward direction starting at t=0 553 %!test # Solve in backward direction starting at t=0
541 %! vref = [-1.205364552835178, 0.951542399860817]; 554 %! vref = [-1.205364552835178, 0.951542399860817];
542 %! vsol = ode45 (@fpol, [0 -2], [2 0]); 555 %! vsol = ode45 (@fpol, [0 -2], [2 0]);
543 %! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-2); 556 %! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-2);
544 %!test ## Solve in backward direction starting at t=2 557 %!test # Solve in backward direction starting at t=2
545 %! vref = [-1.205364552835178, 0.951542399860817]; 558 %! vref = [-1.205364552835178, 0.951542399860817];
546 %! vsol = ode45 (@fpol, [2 -2], fref); 559 %! vsol = ode45 (@fpol, [2 -2], fref);
547 %! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-2); 560 %! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-2);
548 %!test ## Solve in backward direction starting at t=2, with intermidiate step 561 %!test # Solve in backward direction starting at t=2, with intermidiate step
549 %! vref = [-1.205364552835178, 0.951542399860817]; 562 %! vref = [-1.205364552835178, 0.951542399860817];
550 %! vsol = ode45 (@fpol, [2 0 -2], fref); 563 %! vsol = ode45 (@fpol, [2 0 -2], fref);
551 %! idx = find(vsol.x < 0, 1, 'first') - 1; 564 %! idx = find(vsol.x < 0, 1, "first") - 1;
552 %! assert ([vsol.x(idx), vsol.y(idx,:)], [0 2 0], 1e-2); 565 %! assert ([vsol.x(idx), vsol.y(idx,:)], [0 2 0], 1e-2);
553 %! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-2); 566 %! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-2);
554 %!test ## Solve another anonymous function in backward direction 567 %!test # Solve another anonymous function in backward direction
555 %! vref = [-1, 0.367879437558975]; 568 %! vref = [-1, 0.367879437558975];
556 %! vsol = ode45 (@(t,y) y, [0 -1], 1); 569 %! vsol = ode45 (@(t,y) y, [0 -1], 1);
557 %! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3); 570 %! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3);
558 %!test ## Solve another anonymous function below zero 571 %!test # Solve another anonymous function below zero
559 %! vref = [0, 14.77810590694212]; 572 %! vref = [0, 14.77810590694212];
560 %! vsol = ode45 (@(t,y) y, [-2 0], 2); 573 %! vsol = ode45 (@(t,y) y, [-2 0], 2);
561 %! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3); 574 %! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3);
562 %!test ## Solve in backward direction starting at t=0 with MaxStep option 575 %!test # Solve in backward direction starting at t=0 with MaxStep option
563 %! vref = [-1.205364552835178, 0.951542399860817]; 576 %! vref = [-1.205364552835178, 0.951542399860817];
564 %! vopt = odeset ('MaxStep', 1e-3); 577 %! vopt = odeset ("MaxStep", 1e-3);
565 %! vsol = ode45 (@fpol, [0 -2], [2 0], vopt); 578 %! vsol = ode45 (@fpol, [0 -2], [2 0], vopt);
566 %! assert ([abs(vsol.x(8)-vsol.x(7))], [1e-3], 1e-3); 579 %! assert ([abs(vsol.x(8)-vsol.x(7))], [1e-3], 1e-3);
567 %! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3); 580 %! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3);
568 %!test ## AbsTol option 581 %!test # AbsTol option
569 %! vopt = odeset ('AbsTol', 1e-5); 582 %! vopt = odeset ("AbsTol", 1e-5);
570 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt); 583 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
571 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); 584 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
572 %!test ## AbsTol and RelTol option 585 %!test # AbsTol and RelTol option
573 %! vopt = odeset ('AbsTol', 1e-8, 'RelTol', 1e-8); 586 %! vopt = odeset ("AbsTol", 1e-8, "RelTol", 1e-8);
574 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt); 587 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
575 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); 588 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
576 %!test ## RelTol and NormControl option -- higher accuracy 589 %!test # RelTol and NormControl option -- higher accuracy
577 %! vopt = odeset ('RelTol', 1e-8, 'NormControl', 'on'); 590 %! vopt = odeset ("RelTol", 1e-8, "NormControl", "on");
578 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt); 591 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
579 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-5); 592 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-5);
580 %!test ## Keeps initial values while integrating 593 %!test # Keeps initial values while integrating
581 %! vopt = odeset ('NonNegative', 2); 594 %! vopt = odeset ("NonNegative", 2);
582 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt); 595 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
583 %! assert ([vsol.x(end), vsol.y(end,:)], [2, 2, 0], 0.5); 596 %! assert ([vsol.x(end), vsol.y(end,:)], [2, 2, 0], 0.5);
584 %!test ## Details of OutputSel and Refine can't be tested 597 %!test # Details of OutputSel and Refine can't be tested
585 %! vopt = odeset ('OutputFcn', @fout, 'OutputSel', 1, 'Refine', 5); 598 %! vopt = odeset ("OutputFcn", @fout, "OutputSel", 1, "Refine", 5);
586 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt); 599 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
587 %!test ## Details of OutputSave can't be tested 600 %!test # Details of OutputSave can't be tested
588 %! vopt = odeset ('OutputSave', 1, 'OutputSel', 1); 601 %! vopt = odeset ("OutputSave", 1, "OutputSel", 1);
589 %! vsla = ode45 (@fpol, [0 2], [2 0], vopt); 602 %! vsla = ode45 (@fpol, [0 2], [2 0], vopt);
590 %! vopt = odeset ('OutputSave', 2); 603 %! vopt = odeset ("OutputSave", 2);
591 %! vslb = ode45 (@fpol, [0 2], [2 0], vopt); 604 %! vslb = ode45 (@fpol, [0 2], [2 0], vopt);
592 %! assert (length (vsla.x) + 1 >= 2 * length (vslb.x)) 605 %! assert (length (vsla.x) + 1 >= 2 * length (vslb.x))
593 %!test ## Stats must add further elements in vsol 606 %!test # Stats must add further elements in vsol
594 %! vopt = odeset ('Stats', 'on'); 607 %! vopt = odeset ("Stats", "on");
595 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt); 608 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
596 %! assert (isfield (vsol, 'stats')); 609 %! assert (isfield (vsol, "stats"));
597 %! assert (isfield (vsol.stats, 'nsteps')); 610 %! assert (isfield (vsol.stats, "nsteps"));
598 %!test ## Events option add further elements in vsol 611 %!test # Events option add further elements in vsol
599 %! vopt = odeset ('Events', @feve); 612 %! vopt = odeset ("Events", @feve);
600 %! vsol = ode45 (@fpol, [0 10], [2 0], vopt); 613 %! vsol = ode45 (@fpol, [0 10], [2 0], vopt);
601 %! assert (isfield (vsol, 'ie')); 614 %! assert (isfield (vsol, "ie"));
602 %! assert (vsol.ie(1), 2); 615 %! assert (vsol.ie(1), 2);
603 %! assert (isfield (vsol, 'xe')); 616 %! assert (isfield (vsol, "xe"));
604 %! assert (isfield (vsol, 'ye')); 617 %! assert (isfield (vsol, "ye"));
605 %!test ## Events option, now stop integration 618 %!test # Events option, now stop integration
606 %! vopt = odeset ('Events', @fevn, 'NormControl', 'on'); 619 %! vopt = odeset ("Events", @fevn, "NormControl", "on");
607 %! vsol = ode45 (@fpol, [0 10], [2 0], vopt); 620 %! vsol = ode45 (@fpol, [0 10], [2 0], vopt);
608 %! assert ([vsol.ie, vsol.xe, vsol.ye], ... 621 %! assert ([vsol.ie, vsol.xe, vsol.ye], ...
609 %! [2.0, 2.496110, -0.830550, -2.677589], 6e-1); 622 %! [2.0, 2.496110, -0.830550, -2.677589], 6e-1);
610 %!test ## Events option, five output arguments 623 %!test # Events option, five output arguments
611 %! vopt = odeset ('Events', @fevn, 'NormControl', 'on'); 624 %! vopt = odeset ("Events", @fevn, "NormControl", "on");
612 %! [vt, vy, vxe, vye, vie] = ode45 (@fpol, [0 10], [2 0], vopt); 625 %! [vt, vy, vxe, vye, vie] = ode45 (@fpol, [0 10], [2 0], vopt);
613 %! assert ([vie, vxe, vye], ... 626 %! assert ([vie, vxe, vye], ...
614 %! [2.0, 2.496110, -0.830550, -2.677589], 6e-1); 627 %! [2.0, 2.496110, -0.830550, -2.677589], 6e-1);
615 %!test ## Jacobian option 628 %!test # Jacobian option
616 %! vopt = odeset ('Jacobian', @fjac); 629 %! vopt = odeset ("Jacobian", @fjac);
617 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt); 630 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
618 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); 631 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
619 %!test ## Jacobian option and sparse return value 632 %!test # Jacobian option and sparse return value
620 %! vopt = odeset ('Jacobian', @fjcc); 633 %! vopt = odeset ("Jacobian", @fjcc);
621 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt); 634 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
622 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); 635 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
623 %! 636 %!
624 %! ## test for JPattern option is missing 637 %!## test for JPattern option is missing
625 %! ## test for Vectorized option is missing 638 %!## test for Vectorized option is missing
626 %! ## test for NewtonTol option is missing 639 %!## test for NewtonTol option is missing
627 %! ## test for MaxNewtonIterations option is missing 640 %!## test for MaxNewtonIterations option is missing
628 %! 641 %!
629 %!test ## Mass option as function 642 %!test # Mass option as function
630 %! vopt = odeset ('Mass', @fmas); 643 %! vopt = odeset ("Mass", @fmas);
631 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt); 644 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
632 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); 645 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
633 %!test ## Mass option as matrix 646 %!test # Mass option as matrix
634 %! vopt = odeset ('Mass', eye (2,2)); 647 %! vopt = odeset ("Mass", eye (2,2));
635 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt); 648 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
636 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); 649 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
637 %!test ## Mass option as sparse matrix 650 %!test # Mass option as sparse matrix
638 %! vopt = odeset ('Mass', sparse (eye (2,2))); 651 %! vopt = odeset ("Mass", sparse (eye (2,2)));
639 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt); 652 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
640 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); 653 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
641 %!test ## Mass option as function and sparse matrix 654 %!test # Mass option as function and sparse matrix
642 %! vopt = odeset ('Mass', @fmsa); 655 %! vopt = odeset ("Mass", @fmsa);
643 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt); 656 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
644 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); 657 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
645 %!test ## Mass option as function and MStateDependence 658 %!test # Mass option as function and MStateDependence
646 %! vopt = odeset ('Mass', @fmas, 'MStateDependence', 'strong'); 659 %! vopt = odeset ("Mass", @fmas, "MStateDependence", "strong");
647 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt); 660 %! vsol = ode45 (@fpol, [0 2], [2 0], vopt);
648 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); 661 %! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
649 %!test ## Set BDF option to something else than default 662 %!test # Set BDF option to something else than default
650 %! vopt = odeset ('BDF', 'on'); 663 %! vopt = odeset ("BDF", "on");
651 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt); 664 %! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt);
652 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-3); 665 %! assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
653 %! 666 %!
654 %! ## test for MvPattern option is missing 667 %!## test for MvPattern option is missing
655 %! ## test for InitialSlope option is missing 668 %!## test for InitialSlope option is missing
656 %! ## test for MaxOrder option is missing 669 %!## test for MaxOrder option is missing
657 %! 670 %!
658 %! warning ('on', 'OdePkg:InvalidArgument'); 671 %! warning ("on", "OdePkg:InvalidArgument");
659 672
660 ## Local Variables: ***
661 ## mode: octave ***
662 ## End: ***