Mercurial > octave

comparison scripts/ode/ode15i.m @ 30893:e1788b1a315f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Use "fcn" as preferred abbreviation for "function" in m-files. * accumarray.m, accumdim.m, quadl.m, quadv.m, randi.m, structfun.m, __is_function__.m, uigetfile.m, uimenu.m, uiputfile.m, doc_cache_create.m, colorspace_conversion_input_check.m, imageIO.m, argnames.m, vectorize.m, vectorize.m, normest1.m, inputname.m, nthargout.m, display_info_file.m, decic.m, ode15i.m, ode15s.m, ode23.m, ode23s.m, ode45.m, odeset.m, check_default_input.m, integrate_adaptive.m, ode_event_handler.m, runge_kutta_23.m, runge_kutta_23s.m, runge_kutta_45_dorpri.m, runge_kutta_interpolate.m, starting_stepsize.m, __all_opts__.m, fminbnd.m, fminsearch.m, fminunc.m, fsolve.m, fzero.m, sqp.m, fplot.m, plotyy.m, __bar__.m, __ezplot__.m, flat_entry.html, profexport.m, movfun.m, bicg.m, bicgstab.m, cgs.m, eigs.m, gmres.m, pcg.m, __alltohandles__.m, __sprand__.m, qmr.m, tfqmr.m, dump_demos.m: Replace "func", "fun", "fn" in documentation and variable names with "fcn".
author Rik <rik@octave.org>
date Mon, 04 Apr 2022 18:14:56 -0700
parents 796f54d4ddbf
children 360d330cc30e
comparison
equal deleted inserted replaced
30892:1a3cc2811090 30893:e1788b1a315f
22 ## <https://www.gnu.org/licenses/>. 22 ## <https://www.gnu.org/licenses/>.
23 ## 23 ##
24 ######################################################################## 24 ########################################################################
25 25
26 ## -*- texinfo -*- 26 ## -*- texinfo -*-
27 ## @deftypefn {} {[@var{t}, @var{y}] =} ode15i (@var{fun}, @var{trange}, @var{y0}, @var{yp0}) 27 ## @deftypefn {} {[@var{t}, @var{y}] =} ode15i (@var{fcn}, @var{trange}, @var{y0}, @var{yp0})
28 ## @deftypefnx {} {[@var{t}, @var{y}] =} ode15i (@var{fun}, @var{trange}, @var{y0}, @var{yp0}, @var{ode_opt}) 28 ## @deftypefnx {} {[@var{t}, @var{y}] =} ode15i (@var{fcn}, @var{trange}, @var{y0}, @var{yp0}, @var{ode_opt})
29 ## @deftypefnx {} {[@var{t}, @var{y}, @var{te}, @var{ye}, @var{ie}] =} ode15i (@dots{}) 29 ## @deftypefnx {} {[@var{t}, @var{y}, @var{te}, @var{ye}, @var{ie}] =} ode15i (@dots{})
30 ## @deftypefnx {} {@var{solution} =} ode15i (@dots{}) 30 ## @deftypefnx {} {@var{solution} =} ode15i (@dots{})
31 ## @deftypefnx {} {} ode15i (@dots{}) 31 ## @deftypefnx {} {} ode15i (@dots{})
32 ## Solve a set of fully-implicit Ordinary Differential Equations (ODEs) or 32 ## Solve a set of fully-implicit Ordinary Differential Equations (ODEs) or
33 ## index 1 Differential Algebraic Equations (DAEs). 33 ## index 1 Differential Algebraic Equations (DAEs).
34 ## 34 ##
35 ## @code{ode15i} uses a variable step, variable order BDF (Backward 35 ## @code{ode15i} uses a variable step, variable order BDF (Backward
36 ## Differentiation Formula) method that ranges from order 1 to 5. 36 ## Differentiation Formula) method that ranges from order 1 to 5.
37 ## 37 ##
38 ## @var{fun} is a function handle, inline function, or string containing the 38 ## @var{fcn} is a function handle, inline function, or string containing the
39 ## name of the function that defines the ODE: @code{0 = f(t,y,yp)}. The 39 ## name of the function that defines the ODE: @code{0 = f(t,y,yp)}. The
40 ## function must accept three inputs where the first is time @var{t}, the 40 ## function must accept three inputs where the first is time @var{t}, the
41 ## second is the function value @var{y} (a column vector), and the third 41 ## second is the function value @var{y} (a column vector), and the third
42 ## is the derivative value @var{yp} (a column vector). 42 ## is the derivative value @var{yp} (a column vector).
43 ## 43 ##
94 ## @end group 94 ## @end group
95 ## @end smallexample 95 ## @end smallexample
96 ## @seealso{decic, odeset, odeget} 96 ## @seealso{decic, odeset, odeget}
97 ## @end deftypefn 97 ## @end deftypefn
98 98
99 function varargout = ode15i (fun, trange, y0, yp0, varargin) 99 function varargout = ode15i (fcn, trange, y0, yp0, varargin)
100 100
101 if (nargin < 4) 101 if (nargin < 4)
102 print_usage (); 102 print_usage ();
103 endif 103 endif
104 104
110 options = varargin{1}; 110 options = varargin{1};
111 else 111 else
112 options = odeset (); 112 options = odeset ();
113 endif 113 endif
114 114
115 ## Check fun, trange, y0, yp0 115 ## Check fcn, trange, y0, yp0
116 fun = check_default_input (fun, trange, solver, y0, yp0); 116 fcn = check_default_input (fcn, trange, solver, y0, yp0);
117 117
118 if (! isempty (options.Jacobian)) 118 if (! isempty (options.Jacobian))
119 if (ischar (options.Jacobian)) 119 if (ischar (options.Jacobian))
120 if (! exist (options.Jacobian)) 120 if (! exist (options.Jacobian))
121 error ("Octave:invalid-input-arg", 121 error ("Octave:invalid-input-arg",
170 classes, attributes, solver); 170 classes, attributes, solver);
171 171
172 ## Jacobian 172 ## Jacobian
173 options.havejac = false; 173 options.havejac = false;
174 options.havejacsparse = false; 174 options.havejacsparse = false;
175 options.havejacfun = false; 175 options.havejacfcn = false;
176 176
177 if (! isempty (options.Jacobian)) 177 if (! isempty (options.Jacobian))
178 options.havejac = true; 178 options.havejac = true;
179 if (iscell (options.Jacobian)) 179 if (iscell (options.Jacobian))
180 if (numel (options.Jacobian) == 2) 180 if (numel (options.Jacobian) == 2)
194 error ("Octave:invalid-input-arg", 194 error ("Octave:invalid-input-arg",
195 'ode15i: invalid value assigned to field "Jacobian"'); 195 'ode15i: invalid value assigned to field "Jacobian"');
196 endif 196 endif
197 197
198 elseif (is_function_handle (options.Jacobian)) 198 elseif (is_function_handle (options.Jacobian))
199 options.havejacfun = true; 199 options.havejacfcn = true;
200 if (nargin (options.Jacobian) == 3) 200 if (nargin (options.Jacobian) == 3)
201 [J1, J2] = options.Jacobian (trange(1), y0, yp0); 201 [J1, J2] = options.Jacobian (trange(1), y0, yp0);
202 202
203 if ( ! issquare (J1) || rows (J1) != n 203 if ( ! issquare (J1) || rows (J1) != n
204 || ! isnumeric (J1) || ! isreal (J1) 204 || ! isnumeric (J1) || ! isreal (J1)
206 || ! isnumeric (J2) || ! isreal (J2)) 206 || ! isnumeric (J2) || ! isreal (J2))
207 error ("Octave:invalid-input-arg", 207 error ("Octave:invalid-input-arg",
208 'ode15i: "Jacobian" function must evaluate to a real square matrix'); 208 'ode15i: "Jacobian" function must evaluate to a real square matrix');
209 endif 209 endif
210 if (issparse (J1) && issparse (J2)) 210 if (issparse (J1) && issparse (J2))
211 options.havejacsparse = true; # Jac is sparse fun 211 options.havejacsparse = true; # Jac is sparse fcn
212 endif 212 endif
213 else 213 else
214 error ("Octave:invalid-input-arg", 214 error ("Octave:invalid-input-arg",
215 'ode15i: invalid value assigned to field "Jacobian"'); 215 'ode15i: invalid value assigned to field "Jacobian"');
216 endif 216 endif
254 254
255 ## Events 255 ## Events
256 options.haveeventfunction = ! isempty (options.Events); 256 options.haveeventfunction = ! isempty (options.Events);
257 257
258 ## 3 arguments in the event callback of ode15i 258 ## 3 arguments in the event callback of ode15i
259 [t, y, te, ye, ie] = __ode15__ (fun, trange, y0, yp0, options, 3); 259 [t, y, te, ye, ie] = __ode15__ (fcn, trange, y0, yp0, options, 3);
260 260
261 if (nargout == 2) 261 if (nargout == 2)
262 varargout{1} = t; 262 varargout{1} = t;
263 varargout{2} = y; 263 varargout{2} = y;
264 elseif (nargout == 1) 264 elseif (nargout == 1)
284 endfunction 284 endfunction
285 285
286 286
287 %!demo 287 %!demo
288 %! ## Solve Robertson's equations with ode15i 288 %! ## Solve Robertson's equations with ode15i
289 %! fun = @(t, y, yp) [-(yp(1) + 0.04*y(1) - 1e4*y(2)*y(3)); 289 %! fcn = @(t, y, yp) [-(yp(1) + 0.04*y(1) - 1e4*y(2)*y(3));
290 %! -(yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2); 290 %! -(yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2);
291 %! y(1) + y(2) + y(3) - 1]; 291 %! y(1) + y(2) + y(3) - 1];
292 %! 292 %!
293 %! opt = odeset ("RelTol", 1e-4, "AbsTol", [1e-8, 1e-14, 1e-6]); 293 %! opt = odeset ("RelTol", 1e-4, "AbsTol", [1e-8, 1e-14, 1e-6]);
294 %! y0 = [1; 0; 0]; 294 %! y0 = [1; 0; 0];
295 %! yp0 = [-1e-4; 1e-4; 0]; 295 %! yp0 = [-1e-4; 1e-4; 0];
296 %! tspan = [0 4*logspace(-6, 6)]; 296 %! tspan = [0 4*logspace(-6, 6)];
297 %! 297 %!
298 %! [t, y] = ode15i (fun, tspan, y0, yp0, opt); 298 %! [t, y] = ode15i (fcn, tspan, y0, yp0, opt);
299 %! 299 %!
300 %! y(:,2) = 1e4 * y(:, 2); 300 %! y(:,2) = 1e4 * y(:, 2);
301 %! figure (2); 301 %! figure (2);
302 %! semilogx (t, y, "o"); 302 %! semilogx (t, y, "o");
303 %! xlabel ("time"); 303 %! xlabel ("time");
440 %!testif HAVE_SUNDIALS 440 %!testif HAVE_SUNDIALS
441 %! opt = odeset ("RelTol", 1e-6); 441 %! opt = odeset ("RelTol", 1e-6);
442 %! [t, y] = ode15i (@rob, [0, 100], [1; 0; 0], [-1e-4; 1e-4; 0], opt); 442 %! [t, y] = ode15i (@rob, [0, 100], [1; 0; 0], [-1e-4; 1e-4; 0], opt);
443 %! assert ([t(end), y(end,:)], fref, 1e-3); 443 %! assert ([t(end), y(end,:)], fref, 1e-3);
444 444
445 ## Jacobian fun dense 445 ## Jacobian fcn dense
446 %!testif HAVE_SUNDIALS 446 %!testif HAVE_SUNDIALS
447 %! opt = odeset ("Jacobian", @jacfundense); 447 %! opt = odeset ("Jacobian", @jacfundense);
448 %! [t, y] = ode15i (@rob, [0, 100], [1; 0; 0], [-1e-4; 1e-4; 0], opt); 448 %! [t, y] = ode15i (@rob, [0, 100], [1; 0; 0], [-1e-4; 1e-4; 0], opt);
449 %! assert ([t(end), y(end,:)], fref, 1e-3); 449 %! assert ([t(end), y(end,:)], fref, 1e-3);
450 450
451 ## Jacobian fun dense as string 451 ## Jacobian fcn dense as string
452 %!testif HAVE_SUNDIALS 452 %!testif HAVE_SUNDIALS
453 %! opt = odeset ("Jacobian", "jacfundense"); 453 %! opt = odeset ("Jacobian", "jacfundense");
454 %! [t, y] = ode15i (@rob, [0, 100], [1; 0; 0], [-1e-4; 1e-4; 0], opt); 454 %! [t, y] = ode15i (@rob, [0, 100], [1; 0; 0], [-1e-4; 1e-4; 0], opt);
455 %! assert ([t(end), y(end,:)], fref, 1e-3); 455 %! assert ([t(end), y(end,:)], fref, 1e-3);
456 456
457 ## Jacobian fun sparse 457 ## Jacobian fcn sparse
458 %!testif HAVE_SUNDIALS_SUNLINSOL_KLU 458 %!testif HAVE_SUNDIALS_SUNLINSOL_KLU
459 %! opt = odeset ("Jacobian", @jacfunsparse, "AbsTol", 1e-7, "RelTol", 1e-7); 459 %! opt = odeset ("Jacobian", @jacfunsparse, "AbsTol", 1e-7, "RelTol", 1e-7);
460 %! [t, y] = ode15i (@rob, [0, 100], [1; 0; 0], [-1e-4; 1e-4; 0], opt); 460 %! [t, y] = ode15i (@rob, [0, 100], [1; 0; 0], [-1e-4; 1e-4; 0], opt);
461 %! assert ([t(end), y(end,:)], fref, 1e-3); 461 %! assert ([t(end), y(end,:)], fref, 1e-3);
462 462
530 %! A = eye (2); 530 %! A = eye (2);
531 %! [tout, yout] = ode15i (@(t, y, yp) A * y - A * yp, ... 531 %! [tout, yout] = ode15i (@(t, y, yp) A * y - A * yp, ...
532 %! [0, 1], [1, 1], [1; 1]); 532 %! [0, 1], [1, 1], [1; 1]);
533 %! assert (size (yout), [20, 2]); 533 %! assert (size (yout), [20, 2]);
534 534
535 ## Jacobian fun wrong dimension 535 ## Jacobian fcn wrong dimension
536 %!testif HAVE_SUNDIALS 536 %!testif HAVE_SUNDIALS
537 %! opt = odeset ("Jacobian", @jacwrong); 537 %! opt = odeset ("Jacobian", @jacwrong);
538 %! fail ("[t, y] = ode15i (@rob, [0, 4e6], [1; 0; 0], [-1e-4; 1e-4; 0], opt)", 538 %! fail ("[t, y] = ode15i (@rob, [0, 4e6], [1; 0; 0], [-1e-4; 1e-4; 0], opt)",
539 %! '"Jacobian" function must evaluate to a real square matrix'); 539 %! '"Jacobian" function must evaluate to a real square matrix');
540 540
583 %!testif HAVE_SUNDIALS 583 %!testif HAVE_SUNDIALS
584 %! opt = odeset ("Jacobian", "_5yVNhWVJWJn47RKnzxPsyb_"); 584 %! opt = odeset ("Jacobian", "_5yVNhWVJWJn47RKnzxPsyb_");
585 %! fail ("[t, y] = ode15i (@rob, [0, 4e6], [1; 0; 0], [-1e-4; 1e-4; 0], opt)", 585 %! fail ("[t, y] = ode15i (@rob, [0, 4e6], [1; 0; 0], [-1e-4; 1e-4; 0], opt)",
586 %! '"Jacobian" function "_5yVNhWVJWJn47RKnzxPsyb_" not found'); 586 %! '"Jacobian" function "_5yVNhWVJWJn47RKnzxPsyb_" not found');
587 587
588 %!function ydot = fun (t, y, yp) 588 %!function ydot = fcn (t, y, yp)
589 %! ydot = [y - yp]; 589 %! ydot = [y - yp];
590 %!endfunction 590 %!endfunction
591 591
592 %!testif HAVE_SUNDIALS 592 %!testif HAVE_SUNDIALS
593 %! fail ("ode15i ()", "Invalid call to ode15i"); 593 %! fail ("ode15i ()", "Invalid call to ode15i");
600 600
601 %!testif HAVE_SUNDIALS 601 %!testif HAVE_SUNDIALS
602 %! fail ("ode15i (1, 1, 1)", "Invalid call to ode15i"); 602 %! fail ("ode15i (1, 1, 1)", "Invalid call to ode15i");
603 603
604 %!testif HAVE_SUNDIALS 604 %!testif HAVE_SUNDIALS
605 %! fail ("ode15i (1, 1, 1, 1)", "ode15i: fun must be of class:"); 605 %! fail ("ode15i (1, 1, 1, 1)", "ode15i: fcn must be of class:");
606 606
607 %!testif HAVE_SUNDIALS 607 %!testif HAVE_SUNDIALS
608 %! fail ("ode15i (1, 1, 1, 1, 1)", "ode15i: fun must be of class:"); 608 %! fail ("ode15i (1, 1, 1, 1, 1)", "ode15i: fcn must be of class:");
609 609
610 %!testif HAVE_SUNDIALS 610 %!testif HAVE_SUNDIALS
611 %! fail ("ode15i (1, 1, 1, 1, 1, 1)", "ode15i: fun must be of class:"); 611 %! fail ("ode15i (1, 1, 1, 1, 1, 1)", "ode15i: fcn must be of class:");
612 612
613 %!testif HAVE_SUNDIALS 613 %!testif HAVE_SUNDIALS
614 %! fail ("ode15i (@fun, 1, 1, 1)", 614 %! fail ("ode15i (@fcn, 1, 1, 1)",
615 %! "ode15i: invalid value assigned to field 'trange'"); 615 %! "ode15i: invalid value assigned to field 'trange'");
616 616
617 %!testif HAVE_SUNDIALS 617 %!testif HAVE_SUNDIALS
618 %! fail ("ode15i (@fun, [1, 1], 1, 1)", 618 %! fail ("ode15i (@fcn, [1, 1], 1, 1)",
619 %! "ode15i: invalid value assigned to field 'trange'"); 619 %! "ode15i: invalid value assigned to field 'trange'");
620 620
621 %!testif HAVE_SUNDIALS 621 %!testif HAVE_SUNDIALS
622 %! fail ("ode15i (@fun, [1, 2], 1, [1, 2])", 622 %! fail ("ode15i (@fcn, [1, 2], 1, [1, 2])",
623 %! "ode15i: y0 must have 2 elements"); 623 %! "ode15i: y0 must have 2 elements");
624 624
625 %!testif HAVE_SUNDIALS 625 %!testif HAVE_SUNDIALS
626 %! opt = odeset ("RelTol", "_5yVNhWVJWJn47RKnzxPsyb_"); 626 %! opt = odeset ("RelTol", "_5yVNhWVJWJn47RKnzxPsyb_");
627 %! fail ("[t, y] = ode15i (@fun, [0, 2], 2, 2, opt)", 627 %! fail ("[t, y] = ode15i (@fcn, [0, 2], 2, 2, opt)",
628 %! "ode15i: RelTol must be of class:"); 628 %! "ode15i: RelTol must be of class:");
629 629
630 %!testif HAVE_SUNDIALS 630 %!testif HAVE_SUNDIALS
631 %! opt = odeset ("RelTol", [1, 2]); 631 %! opt = odeset ("RelTol", [1, 2]);
632 %! fail ("[t, y] = ode15i (@fun, [0, 2], 2, 2, opt)", 632 %! fail ("[t, y] = ode15i (@fcn, [0, 2], 2, 2, opt)",
633 %! "ode15i: RelTol must be scalar"); 633 %! "ode15i: RelTol must be scalar");
634 634
635 %!testif HAVE_SUNDIALS 635 %!testif HAVE_SUNDIALS
636 %! opt = odeset ("RelTol", -2); 636 %! opt = odeset ("RelTol", -2);
637 %! fail ("[t, y] = ode15i (@fun, [0, 2], 2, 2, opt)", 637 %! fail ("[t, y] = ode15i (@fcn, [0, 2], 2, 2, opt)",
638 %! "ode15i: RelTol must be positive"); 638 %! "ode15i: RelTol must be positive");
639 639
640 %!testif HAVE_SUNDIALS 640 %!testif HAVE_SUNDIALS
641 %! opt = odeset ("AbsTol", "_5yVNhWVJWJn47RKnzxPsyb_"); 641 %! opt = odeset ("AbsTol", "_5yVNhWVJWJn47RKnzxPsyb_");
642 %! fail ("[t, y] = ode15i (@fun, [0, 2], 2, 2, opt)", 642 %! fail ("[t, y] = ode15i (@fcn, [0, 2], 2, 2, opt)",
643 %! "ode15i: AbsTol must be of class:"); 643 %! "ode15i: AbsTol must be of class:");
644 644
645 %!testif HAVE_SUNDIALS 645 %!testif HAVE_SUNDIALS
646 %! opt = odeset ("AbsTol", -1); 646 %! opt = odeset ("AbsTol", -1);
647 %! fail ("[t, y] = ode15i (@fun, [0, 2], 2, 2, opt)", 647 %! fail ("[t, y] = ode15i (@fcn, [0, 2], 2, 2, opt)",
648 %! "ode15i: AbsTol must be positive"); 648 %! "ode15i: AbsTol must be positive");
649 649
650 %!testif HAVE_SUNDIALS 650 %!testif HAVE_SUNDIALS
651 %! opt = odeset ("AbsTol", [1, 1, 1]); 651 %! opt = odeset ("AbsTol", [1, 1, 1]);
652 %! fail ("[t, y] = ode15i (@fun, [0, 2], 2, 2, opt)", 652 %! fail ("[t, y] = ode15i (@fcn, [0, 2], 2, 2, opt)",
653 %! 'ode15i: invalid value assigned to field "AbsTol"'); 653 %! 'ode15i: invalid value assigned to field "AbsTol"');
654 654
655 %!testif HAVE_SUNDIALS 655 %!testif HAVE_SUNDIALS
656 %! A = zeros (2); 656 %! A = zeros (2);
657 %! fail ("ode15i (@(t, y, yp) A * y - A * yp, [0, 1], eye (2), [1, 1])", 657 %! fail ("ode15i (@(t, y, yp) A * y - A * yp, [0, 1], eye (2), [1, 1])",