1@c DO NOT EDIT! Generated automatically by munge-texi.pl. 2 3@c Copyright (C) 1996-2019 John W. Eaton 4@c 5@c This file is part of Octave. 6@c 7@c Octave is free software: you can redistribute it and/or modify it 8@c under the terms of the GNU General Public License as published by 9@c the Free Software Foundation, either version 3 of the License, or 10@c (at your option) any later version. 11@c 12@c Octave is distributed in the hope that it will be useful, but 13@c WITHOUT ANY WARRANTY; without even the implied warranty of 14@c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15@c GNU General Public License for more details. 16@c 17@c You should have received a copy of the GNU General Public License 18@c along with Octave; see the file COPYING. If not, see 19@c <https://www.gnu.org/licenses/>. 20 21@node Differential Equations 22@chapter Differential Equations 23 24Octave has built-in functions for solving ordinary differential equations 25(ODEs), and differential-algebraic equations (DAEs). 26 27@menu 28* Ordinary Differential Equations:: 29* Differential-Algebraic Equations:: 30* Matlab-compatible solvers:: 31@end menu 32 33@cindex differential equations 34@cindex ODE 35@cindex DAE 36 37@node Ordinary Differential Equations 38@section Ordinary Differential Equations 39 40The function @code{lsode} can be used to solve ODEs of the form 41@tex 42$$ 43 {dx\over dt} = f (x, t) 44$$ 45@end tex 46@ifnottex 47 48@example 49@group 50dx 51-- = f (x, t) 52dt 53@end group 54@end example 55 56@end ifnottex 57 58@noindent 59using @nospell{Hindmarsh's} ODE solver @sc{lsode}. 60 61@c lsode libinterp/corefcn/lsode.cc 62@anchor{XREFlsode} 63@deftypefn {} {[@var{x}, @var{istate}, @var{msg}] =} lsode (@var{fcn}, @var{x_0}, @var{t}) 64@deftypefnx {} {[@var{x}, @var{istate}, @var{msg}] =} lsode (@var{fcn}, @var{x_0}, @var{t}, @var{t_crit}) 65Ordinary Differential Equation (ODE) solver. 66 67The set of differential equations to solve is 68@tex 69$$ {dx \over dt} = f (x, t) $$ 70with 71$$ x(t_0) = x_0 $$ 72@end tex 73@ifnottex 74 75@example 76@group 77dx 78-- = f (x, t) 79dt 80@end group 81@end example 82 83@noindent 84with 85 86@example 87x(t_0) = x_0 88@end example 89 90@end ifnottex 91The solution is returned in the matrix @var{x}, with each row 92corresponding to an element of the vector @var{t}. The first element 93of @var{t} should be @math{t_0} and should correspond to the initial 94state of the system @var{x_0}, so that the first row of the output 95is @var{x_0}. 96 97The first argument, @var{fcn}, is a string, inline, or function handle 98that names the function @math{f} to call to compute the vector of right 99hand sides for the set of equations. The function must have the form 100 101@example 102@var{xdot} = f (@var{x}, @var{t}) 103@end example 104 105@noindent 106in which @var{xdot} and @var{x} are vectors and @var{t} is a scalar. 107 108If @var{fcn} is a two-element string array or a two-element cell array 109of strings, inline functions, or function handles, the first element names 110the function @math{f} described above, and the second element names a 111function to compute the Jacobian of @math{f}. The Jacobian function 112must have the form 113 114@example 115@var{jac} = j (@var{x}, @var{t}) 116@end example 117 118@noindent 119in which @var{jac} is the matrix of partial derivatives 120@tex 121$$ J = {\partial f_i \over \partial x_j} = \left[\matrix{ 122{\partial f_1 \over \partial x_1} 123 & {\partial f_1 \over \partial x_2} 124 & \cdots 125 & {\partial f_1 \over \partial x_N} \cr 126{\partial f_2 \over \partial x_1} 127 & {\partial f_2 \over \partial x_2} 128 & \cdots 129 & {\partial f_2 \over \partial x_N} \cr 130 \vdots & \vdots & \ddots & \vdots \cr 131{\partial f_3 \over \partial x_1} 132 & {\partial f_3 \over \partial x_2} 133 & \cdots 134 & {\partial f_3 \over \partial x_N} \cr}\right]$$ 135@end tex 136@ifnottex 137 138@example 139@group 140 | df_1 df_1 df_1 | 141 | ---- ---- ... ---- | 142 | dx_1 dx_2 dx_N | 143 | | 144 | df_2 df_2 df_2 | 145 | ---- ---- ... ---- | 146 df_i | dx_1 dx_2 dx_N | 147jac = ---- = | | 148 dx_j | . . . . | 149 | . . . . | 150 | . . . . | 151 | | 152 | df_N df_N df_N | 153 | ---- ---- ... ---- | 154 | dx_1 dx_2 dx_N | 155@end group 156@end example 157 158@end ifnottex 159 160The second argument specifies the initial state of the system @math{x_0}. The 161third argument is a vector, @var{t}, specifying the time values for which a 162solution is sought. 163 164The fourth argument is optional, and may be used to specify a set of 165times that the ODE solver should not integrate past. It is useful for 166avoiding difficulties with singularities and points where there is a 167discontinuity in the derivative. 168 169After a successful computation, the value of @var{istate} will be 2 170(consistent with the Fortran version of @sc{lsode}). 171 172If the computation is not successful, @var{istate} will be something 173other than 2 and @var{msg} will contain additional information. 174 175You can use the function @code{lsode_options} to set optional 176parameters for @code{lsode}. 177 178See @nospell{Alan C. Hindmarsh}, 179@cite{ODEPACK, A Systematized Collection of ODE Solvers}, 180in Scientific Computing, @nospell{R. S. Stepleman}, editor, (1983) 181or @url{https://computing.llnl.gov/projects/odepack} 182for more information about the inner workings of @code{lsode}. 183 184Example: Solve the @nospell{Van der Pol} equation 185 186@example 187@group 188fvdp = @@(@var{y},@var{t}) [@var{y}(2); (1 - @var{y}(1)^2) * @var{y}(2) - @var{y}(1)]; 189@var{t} = linspace (0, 20, 100); 190@var{y} = lsode (fvdp, [2; 0], @var{t}); 191@end group 192@end example 193@xseealso{@ref{XREFdaspk,,daspk}, @ref{XREFdassl,,dassl}, @ref{XREFdasrt,,dasrt}} 194@end deftypefn 195 196 197@c lsode_options libinterp/corefcn/LSODE-opts.cc 198@anchor{XREFlsode_options} 199@deftypefn {} {} lsode_options () 200@deftypefnx {} {val =} lsode_options (@var{opt}) 201@deftypefnx {} {} lsode_options (@var{opt}, @var{val}) 202Query or set options for the function @code{lsode}. 203 204When called with no arguments, the names of all available options and 205their current values are displayed. 206 207Given one argument, return the value of the option @var{opt}. 208 209When called with two arguments, @code{lsode_options} sets the option 210@var{opt} to value @var{val}. 211 212Options include 213 214@table @asis 215@item @qcode{"absolute tolerance"} 216Absolute tolerance. May be either vector or scalar. If a vector, it 217must match the dimension of the state vector. 218 219@item @qcode{"relative tolerance"} 220Relative tolerance parameter. Unlike the absolute tolerance, this 221parameter may only be a scalar. 222 223The local error test applied at each integration step is 224 225@example 226@group 227 abs (local error in x(i)) <= ... 228 rtol * abs (y(i)) + atol(i) 229@end group 230@end example 231 232@item @qcode{"integration method"} 233A string specifying the method of integration to use to solve the ODE 234system. Valid values are 235 236@table @asis 237@item @qcode{"adams"} 238@itemx @qcode{"non-stiff"} 239No Jacobian used (even if it is available). 240 241@item @qcode{"bdf"} 242@itemx @qcode{"stiff"} 243Use stiff backward differentiation formula (BDF) method. If a 244function to compute the Jacobian is not supplied, @code{lsode} will 245compute a finite difference approximation of the Jacobian matrix. 246@end table 247 248@item @qcode{"initial step size"} 249The step size to be attempted on the first step (default is determined 250automatically). 251 252@item @qcode{"maximum order"} 253Restrict the maximum order of the solution method. If using the Adams 254method, this option must be between 1 and 12. Otherwise, it must be 255between 1 and 5, inclusive. 256 257@item @qcode{"maximum step size"} 258Setting the maximum stepsize will avoid passing over very large 259regions (default is not specified). 260 261@item @qcode{"minimum step size"} 262The minimum absolute step size allowed (default is 0). 263 264@item @qcode{"step limit"} 265Maximum number of steps allowed (default is 100000). 266@end table 267@end deftypefn 268 269 270Here is an example of solving a set of three differential equations using 271@code{lsode}. Given the function 272 273@cindex oregonator 274 275@example 276@group 277## oregonator differential equation 278function xdot = f (x, t) 279 280 xdot = zeros (3,1); 281 282 xdot(1) = 77.27 * (x(2) - x(1)*x(2) + x(1) ... 283 - 8.375e-06*x(1)^2); 284 xdot(2) = (x(3) - x(1)*x(2) - x(2)) / 77.27; 285 xdot(3) = 0.161*(x(1) - x(3)); 286 287endfunction 288@end group 289@end example 290 291@noindent 292and the initial condition @code{x0 = [ 4; 1.1; 4 ]}, the set of 293equations can be integrated using the command 294 295@example 296@group 297t = linspace (0, 500, 1000); 298 299y = lsode ("f", x0, t); 300@end group 301@end example 302 303If you try this, you will see that the value of the result changes 304dramatically between @var{t} = 0 and 5, and again around @var{t} = 305. 305A more efficient set of output points might be 306 307@example 308@group 309t = [0, logspace(-1, log10(303), 150), ... 310 logspace(log10(304), log10(500), 150)]; 311@end group 312@end example 313 314An m-file for the differential equation used above is included with the 315Octave distribution in the examples directory under the name 316@file{oregonator.m}. 317 318 319@node Differential-Algebraic Equations 320@section Differential-Algebraic Equations 321 322The function @code{daspk} can be used to solve DAEs of the form 323@tex 324$$ 325 0 = f (\dot{x}, x, t), \qquad x(t=0) = x_0, \dot{x}(t=0) = \dot{x}_0 326$$ 327@end tex 328@ifnottex 329 330@example 3310 = f (x-dot, x, t), x(t=0) = x_0, x-dot(t=0) = x-dot_0 332@end example 333 334@end ifnottex 335 336@noindent 337where 338@tex 339$\dot{x} = {dx \over dt}$ 340@end tex 341@ifnottex 342@math{x-dot} 343@end ifnottex 344is the derivative of @math{x}. The equation is solved using 345@nospell{Petzold's} DAE solver @sc{daspk}. 346 347@c daspk libinterp/corefcn/daspk.cc 348@anchor{XREFdaspk} 349@deftypefn {} {[@var{x}, @var{xdot}, @var{istate}, @var{msg}] =} daspk (@var{fcn}, @var{x_0}, @var{xdot_0}, @var{t}, @var{t_crit}) 350Solve a set of differential-algebraic equations. 351 352@code{daspk} solves the set of equations 353@tex 354$$ 0 = f (x, \dot{x}, t) $$ 355with 356$$ x(t_0) = x_0, \dot{x}(t_0) = \dot{x}_0 $$ 357@end tex 358@ifnottex 359 360@example 3610 = f (x, xdot, t) 362@end example 363 364@noindent 365with 366 367@example 368x(t_0) = x_0, xdot(t_0) = xdot_0 369@end example 370 371@end ifnottex 372The solution is returned in the matrices @var{x} and @var{xdot}, 373with each row in the result matrices corresponding to one of the 374elements in the vector @var{t}. The first element of @var{t} 375should be @math{t_0} and correspond to the initial state of the 376system @var{x_0} and its derivative @var{xdot_0}, so that the first 377row of the output @var{x} is @var{x_0} and the first row 378of the output @var{xdot} is @var{xdot_0}. 379 380The first argument, @var{fcn}, is a string, inline, or function handle 381that names the function @math{f} to call to compute the vector of 382residuals for the set of equations. It must have the form 383 384@example 385@var{res} = f (@var{x}, @var{xdot}, @var{t}) 386@end example 387 388@noindent 389in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a 390scalar. 391 392If @var{fcn} is a two-element string array or a two-element cell array 393of strings, inline functions, or function handles, the first element names 394the function @math{f} described above, and the second element names a 395function to compute the modified Jacobian 396@tex 397$$ 398J = {\partial f \over \partial x} 399 + c {\partial f \over \partial \dot{x}} 400$$ 401@end tex 402@ifnottex 403 404@example 405@group 406 df df 407jac = -- + c ------ 408 dx d xdot 409@end group 410@end example 411 412@end ifnottex 413 414The modified Jacobian function must have the form 415 416@example 417@group 418 419@var{jac} = j (@var{x}, @var{xdot}, @var{t}, @var{c}) 420 421@end group 422@end example 423 424The second and third arguments to @code{daspk} specify the initial 425condition of the states and their derivatives, and the fourth argument 426specifies a vector of output times at which the solution is desired, 427including the time corresponding to the initial condition. 428 429The set of initial states and derivatives are not strictly required to 430be consistent. If they are not consistent, you must use the 431@code{daspk_options} function to provide additional information so 432that @code{daspk} can compute a consistent starting point. 433 434The fifth argument is optional, and may be used to specify a set of 435times that the DAE solver should not integrate past. It is useful for 436avoiding difficulties with singularities and points where there is a 437discontinuity in the derivative. 438 439After a successful computation, the value of @var{istate} will be 440greater than zero (consistent with the Fortran version of @sc{daspk}). 441 442If the computation is not successful, the value of @var{istate} will be 443less than zero and @var{msg} will contain additional information. 444 445You can use the function @code{daspk_options} to set optional 446parameters for @code{daspk}. 447@xseealso{@ref{XREFdassl,,dassl}} 448@end deftypefn 449 450 451@c daspk_options libinterp/corefcn/DASPK-opts.cc 452@anchor{XREFdaspk_options} 453@deftypefn {} {} daspk_options () 454@deftypefnx {} {val =} daspk_options (@var{opt}) 455@deftypefnx {} {} daspk_options (@var{opt}, @var{val}) 456Query or set options for the function @code{daspk}. 457 458When called with no arguments, the names of all available options and 459their current values are displayed. 460 461Given one argument, return the value of the option @var{opt}. 462 463When called with two arguments, @code{daspk_options} sets the option 464@var{opt} to value @var{val}. 465 466Options include 467 468@table @asis 469@item @qcode{"absolute tolerance"} 470Absolute tolerance. May be either vector or scalar. If a vector, it 471must match the dimension of the state vector, and the relative 472tolerance must also be a vector of the same length. 473 474@item @qcode{"relative tolerance"} 475Relative tolerance. May be either vector or scalar. If a vector, it 476must match the dimension of the state vector, and the absolute 477tolerance must also be a vector of the same length. 478 479The local error test applied at each integration step is 480 481@example 482@group 483 abs (local error in x(i)) 484 <= rtol(i) * abs (Y(i)) + atol(i) 485@end group 486@end example 487 488@item @qcode{"compute consistent initial condition"} 489Denoting the differential variables in the state vector by @samp{Y_d} 490and the algebraic variables by @samp{Y_a}, @code{ddaspk} can solve 491one of two initialization problems: 492 493@enumerate 494@item Given Y_d, calculate Y_a and Y'_d 495 496@item Given Y', calculate Y. 497@end enumerate 498 499In either case, initial values for the given components are input, and 500initial guesses for the unknown components must also be provided as 501input. Set this option to 1 to solve the first problem, or 2 to solve 502the second (the default is 0, so you must provide a set of 503initial conditions that are consistent). 504 505If this option is set to a nonzero value, you must also set the 506@qcode{"algebraic variables"} option to declare which variables in the 507problem are algebraic. 508 509@item @qcode{"use initial condition heuristics"} 510Set to a nonzero value to use the initial condition heuristics options 511described below. 512 513@item @qcode{"initial condition heuristics"} 514A vector of the following parameters that can be used to control the 515initial condition calculation. 516 517@table @code 518@item MXNIT 519Maximum number of Newton iterations (default is 5). 520 521@item MXNJ 522Maximum number of Jacobian evaluations (default is 6). 523 524@item MXNH 525Maximum number of values of the artificial stepsize parameter to be 526tried if the @qcode{"compute consistent initial condition"} option has 527been set to 1 (default is 5). 528 529Note that the maximum total number of Newton iterations allowed is 530@code{MXNIT*MXNJ*MXNH} if the @qcode{"compute consistent initial 531condition"} option has been set to 1 and @code{MXNIT*MXNJ} if it is 532set to 2. 533 534@item LSOFF 535Set to a nonzero value to disable the linesearch algorithm (default is 5360). 537 538@item STPTOL 539Minimum scaled step in linesearch algorithm (default is eps^(2/3)). 540 541@item EPINIT 542Swing factor in the Newton iteration convergence test. The test is 543applied to the residual vector, premultiplied by the approximate 544Jacobian. For convergence, the weighted RMS norm of this vector 545(scaled by the error weights) must be less than @code{EPINIT*EPCON}, 546where @code{EPCON} = 0.33 is the analogous test constant used in the 547time steps. The default is @code{EPINIT} = 0.01. 548@end table 549 550@item @qcode{"print initial condition info"} 551Set this option to a nonzero value to display detailed information 552about the initial condition calculation (default is 0). 553 554@item @qcode{"exclude algebraic variables from error test"} 555Set to a nonzero value to exclude algebraic variables from the error 556test. You must also set the @qcode{"algebraic variables"} option to 557declare which variables in the problem are algebraic (default is 0). 558 559@item @qcode{"algebraic variables"} 560A vector of the same length as the state vector. A nonzero element 561indicates that the corresponding element of the state vector is an 562algebraic variable (i.e., its derivative does not appear explicitly 563in the equation set). 564 565This option is required by the 566@qcode{"compute consistent initial condition"} and 567@qcode{"exclude algebraic variables from error test"} options. 568 569@item @qcode{"enforce inequality constraints"} 570Set to one of the following values to enforce the inequality 571constraints specified by the @qcode{"inequality constraint types"} 572option (default is 0). 573 574@enumerate 575@item To have constraint checking only in the initial condition calculation. 576 577@item To enforce constraint checking during the integration. 578 579@item To enforce both options 1 and 2. 580@end enumerate 581 582@item @qcode{"inequality constraint types"} 583A vector of the same length as the state specifying the type of 584inequality constraint. Each element of the vector corresponds to an 585element of the state and should be assigned one of the following 586codes 587 588@table @asis 589@item -2 590Less than zero. 591 592@item -1 593Less than or equal to zero. 594 595@item 0 596Not constrained. 597 598@item 1 599Greater than or equal to zero. 600 601@item 2 602Greater than zero. 603@end table 604 605This option only has an effect if the 606@qcode{"enforce inequality constraints"} option is nonzero. 607 608@item @qcode{"initial step size"} 609Differential-algebraic problems may occasionally suffer from severe 610scaling difficulties on the first step. If you know a great deal 611about the scaling of your problem, you can help to alleviate this 612problem by specifying an initial stepsize (default is computed 613automatically). 614 615@item @qcode{"maximum order"} 616Restrict the maximum order of the solution method. This option must 617be between 1 and 5, inclusive (default is 5). 618 619@item @qcode{"maximum step size"} 620Setting the maximum stepsize will avoid passing over very large 621regions (default is not specified). 622@end table 623@end deftypefn 624 625 626Octave also includes @sc{dassl}, an earlier version of @sc{daspk}, 627and @sc{dasrt}, which can be used to solve DAEs with constraints 628(stopping conditions). 629 630@c dassl libinterp/corefcn/dassl.cc 631@anchor{XREFdassl} 632@deftypefn {} {[@var{x}, @var{xdot}, @var{istate}, @var{msg}] =} dassl (@var{fcn}, @var{x_0}, @var{xdot_0}, @var{t}, @var{t_crit}) 633Solve a set of differential-algebraic equations. 634 635@code{dassl} solves the set of equations 636@tex 637$$ 0 = f (x, \dot{x}, t) $$ 638with 639$$ x(t_0) = x_0, \dot{x}(t_0) = \dot{x}_0 $$ 640@end tex 641@ifnottex 642 643@example 6440 = f (x, xdot, t) 645@end example 646 647@noindent 648with 649 650@example 651x(t_0) = x_0, xdot(t_0) = xdot_0 652@end example 653 654@end ifnottex 655The solution is returned in the matrices @var{x} and @var{xdot}, 656with each row in the result matrices corresponding to one of the 657elements in the vector @var{t}. The first element of @var{t} 658should be @math{t_0} and correspond to the initial state of the 659system @var{x_0} and its derivative @var{xdot_0}, so that the first 660row of the output @var{x} is @var{x_0} and the first row 661of the output @var{xdot} is @var{xdot_0}. 662 663The first argument, @var{fcn}, is a string, inline, or function handle 664that names the function @math{f} to call to compute the vector of 665residuals for the set of equations. It must have the form 666 667@example 668@var{res} = f (@var{x}, @var{xdot}, @var{t}) 669@end example 670 671@noindent 672in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a 673scalar. 674 675If @var{fcn} is a two-element string array or a two-element cell array 676of strings, inline functions, or function handles, the first element names 677the function @math{f} described above, and the second element names a 678function to compute the modified Jacobian 679 680@tex 681$$ 682J = {\partial f \over \partial x} 683 + c {\partial f \over \partial \dot{x}} 684$$ 685@end tex 686@ifnottex 687 688@example 689@group 690 df df 691jac = -- + c ------ 692 dx d xdot 693@end group 694@end example 695 696@end ifnottex 697 698The modified Jacobian function must have the form 699 700@example 701@group 702 703@var{jac} = j (@var{x}, @var{xdot}, @var{t}, @var{c}) 704 705@end group 706@end example 707 708The second and third arguments to @code{dassl} specify the initial 709condition of the states and their derivatives, and the fourth argument 710specifies a vector of output times at which the solution is desired, 711including the time corresponding to the initial condition. 712 713The set of initial states and derivatives are not strictly required to 714be consistent. In practice, however, @sc{dassl} is not very good at 715determining a consistent set for you, so it is best if you ensure that 716the initial values result in the function evaluating to zero. 717 718The fifth argument is optional, and may be used to specify a set of 719times that the DAE solver should not integrate past. It is useful for 720avoiding difficulties with singularities and points where there is a 721discontinuity in the derivative. 722 723After a successful computation, the value of @var{istate} will be 724greater than zero (consistent with the Fortran version of @sc{dassl}). 725 726If the computation is not successful, the value of @var{istate} will be 727less than zero and @var{msg} will contain additional information. 728 729You can use the function @code{dassl_options} to set optional 730parameters for @code{dassl}. 731@xseealso{@ref{XREFdaspk,,daspk}, @ref{XREFdasrt,,dasrt}, @ref{XREFlsode,,lsode}} 732@end deftypefn 733 734 735@c dassl_options libinterp/corefcn/DASSL-opts.cc 736@anchor{XREFdassl_options} 737@deftypefn {} {} dassl_options () 738@deftypefnx {} {val =} dassl_options (@var{opt}) 739@deftypefnx {} {} dassl_options (@var{opt}, @var{val}) 740Query or set options for the function @code{dassl}. 741 742When called with no arguments, the names of all available options and 743their current values are displayed. 744 745Given one argument, return the value of the option @var{opt}. 746 747When called with two arguments, @code{dassl_options} sets the option 748@var{opt} to value @var{val}. 749 750Options include 751 752@table @asis 753@item @qcode{"absolute tolerance"} 754Absolute tolerance. May be either vector or scalar. If a vector, it 755must match the dimension of the state vector, and the relative 756tolerance must also be a vector of the same length. 757 758@item @qcode{"relative tolerance"} 759Relative tolerance. May be either vector or scalar. If a vector, it 760must match the dimension of the state vector, and the absolute 761tolerance must also be a vector of the same length. 762 763The local error test applied at each integration step is 764 765@example 766@group 767 abs (local error in x(i)) 768 <= rtol(i) * abs (Y(i)) + atol(i) 769@end group 770@end example 771 772@item @qcode{"compute consistent initial condition"} 773If nonzero, @code{dassl} will attempt to compute a consistent set of initial 774conditions. This is generally not reliable, so it is best to provide 775a consistent set and leave this option set to zero. 776 777@item @qcode{"enforce nonnegativity constraints"} 778If you know that the solutions to your equations will always be 779non-negative, it may help to set this parameter to a nonzero 780value. However, it is probably best to try leaving this option set to 781zero first, and only setting it to a nonzero value if that doesn't 782work very well. 783 784@item @qcode{"initial step size"} 785Differential-algebraic problems may occasionally suffer from severe 786scaling difficulties on the first step. If you know a great deal 787about the scaling of your problem, you can help to alleviate this 788problem by specifying an initial stepsize. 789 790@item @qcode{"maximum order"} 791Restrict the maximum order of the solution method. This option must 792be between 1 and 5, inclusive. 793 794@item @qcode{"maximum step size"} 795Setting the maximum stepsize will avoid passing over very large 796regions (default is not specified). 797 798@item @qcode{"step limit"} 799Maximum number of integration steps to attempt on a single call to the 800underlying Fortran code. 801@end table 802@end deftypefn 803 804 805@c dasrt libinterp/corefcn/dasrt.cc 806@anchor{XREFdasrt} 807@deftypefn {} {[@var{x}, @var{xdot}, @var{t_out}, @var{istat}, @var{msg}] =} dasrt (@var{fcn}, @var{g}, @var{x_0}, @var{xdot_0}, @var{t}) 808@deftypefnx {} {@dots{} =} dasrt (@var{fcn}, @var{g}, @var{x_0}, @var{xdot_0}, @var{t}, @var{t_crit}) 809@deftypefnx {} {@dots{} =} dasrt (@var{fcn}, @var{x_0}, @var{xdot_0}, @var{t}) 810@deftypefnx {} {@dots{} =} dasrt (@var{fcn}, @var{x_0}, @var{xdot_0}, @var{t}, @var{t_crit}) 811Solve a set of differential-algebraic equations. 812 813@code{dasrt} solves the set of equations 814@tex 815$$ 0 = f (x, \dot{x}, t) $$ 816with 817$$ x(t_0) = x_0, \dot{x}(t_0) = \dot{x}_0 $$ 818@end tex 819@ifnottex 820 821@example 8220 = f (x, xdot, t) 823@end example 824 825@noindent 826with 827 828@example 829x(t_0) = x_0, xdot(t_0) = xdot_0 830@end example 831 832@end ifnottex 833with functional stopping criteria (root solving). 834 835The solution is returned in the matrices @var{x} and @var{xdot}, 836with each row in the result matrices corresponding to one of the 837elements in the vector @var{t_out}. The first element of @var{t} 838should be @math{t_0} and correspond to the initial state of the 839system @var{x_0} and its derivative @var{xdot_0}, so that the first 840row of the output @var{x} is @var{x_0} and the first row 841of the output @var{xdot} is @var{xdot_0}. 842 843The vector @var{t} provides an upper limit on the length of the 844integration. If the stopping condition is met, the vector 845@var{t_out} will be shorter than @var{t}, and the final element of 846@var{t_out} will be the point at which the stopping condition was met, 847and may not correspond to any element of the vector @var{t}. 848 849The first argument, @var{fcn}, is a string, inline, or function handle 850that names the function @math{f} to call to compute the vector of 851residuals for the set of equations. It must have the form 852 853@example 854@var{res} = f (@var{x}, @var{xdot}, @var{t}) 855@end example 856 857@noindent 858in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a 859scalar. 860 861If @var{fcn} is a two-element string array or a two-element cell array 862of strings, inline functions, or function handles, the first element names 863the function @math{f} described above, and the second element names a 864function to compute the modified Jacobian 865 866@tex 867$$ 868J = {\partial f \over \partial x} 869 + c {\partial f \over \partial \dot{x}} 870$$ 871@end tex 872@ifnottex 873 874@example 875@group 876 df df 877jac = -- + c ------ 878 dx d xdot 879@end group 880@end example 881 882@end ifnottex 883 884The modified Jacobian function must have the form 885 886@example 887@group 888 889@var{jac} = j (@var{x}, @var{xdot}, @var{t}, @var{c}) 890 891@end group 892@end example 893 894The optional second argument names a function that defines the 895constraint functions whose roots are desired during the integration. 896This function must have the form 897 898@example 899@var{g_out} = g (@var{x}, @var{t}) 900@end example 901 902@noindent 903and return a vector of the constraint function values. 904If the value of any of the constraint functions changes sign, @sc{dasrt} 905will attempt to stop the integration at the point of the sign change. 906 907If the name of the constraint function is omitted, @code{dasrt} solves 908the same problem as @code{daspk} or @code{dassl}. 909 910Note that because of numerical errors in the constraint functions 911due to round-off and integration error, @sc{dasrt} may return false 912roots, or return the same root at two or more nearly equal values of 913@var{T}. If such false roots are suspected, the user should consider 914smaller error tolerances or higher precision in the evaluation of the 915constraint functions. 916 917If a root of some constraint function defines the end of the problem, 918the input to @sc{dasrt} should nevertheless allow integration to a 919point slightly past that root, so that @sc{dasrt} can locate the root 920by interpolation. 921 922The third and fourth arguments to @code{dasrt} specify the initial 923condition of the states and their derivatives, and the fourth argument 924specifies a vector of output times at which the solution is desired, 925including the time corresponding to the initial condition. 926 927The set of initial states and derivatives are not strictly required to 928be consistent. In practice, however, @sc{dassl} is not very good at 929determining a consistent set for you, so it is best if you ensure that 930the initial values result in the function evaluating to zero. 931 932The sixth argument is optional, and may be used to specify a set of 933times that the DAE solver should not integrate past. It is useful for 934avoiding difficulties with singularities and points where there is a 935discontinuity in the derivative. 936 937After a successful computation, the value of @var{istate} will be 938greater than zero (consistent with the Fortran version of @sc{dassl}). 939 940If the computation is not successful, the value of @var{istate} will be 941less than zero and @var{msg} will contain additional information. 942 943You can use the function @code{dasrt_options} to set optional 944parameters for @code{dasrt}. 945@xseealso{@ref{XREFdasrt_options,,dasrt_options}, @ref{XREFdaspk,,daspk}, @ref{XREFdasrt,,dasrt}, @ref{XREFlsode,,lsode}} 946@end deftypefn 947 948 949@c dasrt_options libinterp/corefcn/DASRT-opts.cc 950@anchor{XREFdasrt_options} 951@deftypefn {} {} dasrt_options () 952@deftypefnx {} {val =} dasrt_options (@var{opt}) 953@deftypefnx {} {} dasrt_options (@var{opt}, @var{val}) 954Query or set options for the function @code{dasrt}. 955 956When called with no arguments, the names of all available options and 957their current values are displayed. 958 959Given one argument, return the value of the option @var{opt}. 960 961When called with two arguments, @code{dasrt_options} sets the option 962@var{opt} to value @var{val}. 963 964Options include 965 966@table @asis 967@item @qcode{"absolute tolerance"} 968Absolute tolerance. May be either vector or scalar. If a vector, it 969must match the dimension of the state vector, and the relative 970tolerance must also be a vector of the same length. 971 972@item @qcode{"relative tolerance"} 973Relative tolerance. May be either vector or scalar. If a vector, it 974must match the dimension of the state vector, and the absolute 975tolerance must also be a vector of the same length. 976 977The local error test applied at each integration step is 978 979@example 980@group 981 abs (local error in x(i)) <= ... 982 rtol(i) * abs (Y(i)) + atol(i) 983@end group 984@end example 985 986@item @qcode{"initial step size"} 987Differential-algebraic problems may occasionally suffer from severe 988scaling difficulties on the first step. If you know a great deal 989about the scaling of your problem, you can help to alleviate this 990problem by specifying an initial stepsize. 991 992@item @qcode{"maximum order"} 993Restrict the maximum order of the solution method. This option must 994be between 1 and 5, inclusive. 995 996@item @qcode{"maximum step size"} 997Setting the maximum stepsize will avoid passing over very large 998regions. 999 1000@item @qcode{"step limit"} 1001Maximum number of integration steps to attempt on a single call to the 1002underlying Fortran code. 1003@end table 1004@end deftypefn 1005 1006 1007See @nospell{K. E. Brenan}, et al., @cite{Numerical Solution of Initial-Value 1008Problems in Differential-Algebraic Equations}, North-Holland (1989), 1009DOI: @url{https://doi.org/10.1137/1.9781611971224}, 1010for more information about the implementation of @sc{dassl}. 1011 1012 1013@node Matlab-compatible solvers 1014@section Matlab-compatible solvers 1015 1016Octave also provides a set of solvers for initial value problems for ordinary 1017differential equations (ODEs) that have a @sc{matlab}-compatible interface. 1018The options for this class of methods are set using the functions. 1019 1020@itemize 1021 @item @ref{XREFodeset,,odeset} 1022 1023 @item @ref{XREFodeget,,odeget} 1024@end itemize 1025 1026Currently implemented solvers are: 1027 1028@itemize 1029 @item @nospell{Runge-Kutta} methods 1030 1031 @itemize 1032 @item @ref{XREFode45,,ode45} integrates a system of non-stiff ODEs or 1033 index-1 differential-algebraic equations (DAEs) using the high-order, 1034 variable-step @nospell{Dormand-Prince} method. It requires six function 1035 evaluations per integration step, but may take larger steps on smooth 1036 problems than @code{ode23}: potentially offering improved efficiency at 1037 smaller tolerances. 1038 1039 @item @ref{XREFode23,,ode23} integrates a system of non-stiff ODEs or (or 1040 index-1 DAEs). It uses the third-order @nospell{Bogacki-Shampine} method 1041 and adapts the local step size in order to satisfy a user-specified 1042 tolerance. The solver requires three function evaluations per integration 1043 step. 1044 1045 @item @ref{XREFode23s,,ode23s} integrates a system of stiff ODEs (or 1046 index-1 DAEs) using a modified second-order @nospell{Rosenbrock} method. 1047 @end itemize 1048 1049 @item Linear multistep methods 1050 1051 @itemize 1052 @item @ref{XREFode15s,,ode15s} integrates a system of stiff ODEs (or 1053 index-1 DAEs) using a variable step, variable order method based on 1054 Backward Difference Formulas (BDF). 1055 1056 @item @ref{XREFode15i,,ode15i} integrates a system of fully-implicit ODEs 1057 (or index-1 DAEs) using the same variable step, variable order method as 1058 @code{ode15s}. @ref{XREFdecic,,decic} can be used to compute consistent 1059 initial conditions for @code{ode15i}. 1060 @end itemize 1061@end itemize 1062 1063Detailed information on the solvers are given in @nospell{L. F. Shampine} and 1064@nospell{M. W. Reichelt}, @cite{The MATLAB ODE Suite}, SIAM Journal on 1065Scientific Computing, Vol. 18, 1997, pp. 1–22, 1066DOI: @url{https://doi.org/10.1137/S1064827594276424}. 1067 1068@c ode45 scripts/ode/ode45.m 1069@anchor{XREFode45} 1070@deftypefn {} {[@var{t}, @var{y}] =} ode45 (@var{fun}, @var{trange}, @var{init}) 1071@deftypefnx {} {[@var{t}, @var{y}] =} ode45 (@var{fun}, @var{trange}, @var{init}, @var{ode_opt}) 1072@deftypefnx {} {[@var{t}, @var{y}, @var{te}, @var{ye}, @var{ie}] =} ode45 (@dots{}) 1073@deftypefnx {} {@var{solution} =} ode45 (@dots{}) 1074@deftypefnx {} {} ode45 (@dots{}) 1075 1076Solve a set of non-stiff Ordinary Differential Equations (non-stiff ODEs) 1077with the well known explicit @nospell{Dormand-Prince} method of order 4. 1078 1079@var{fun} is a function handle, inline function, or string containing the 1080name of the function that defines the ODE: @code{y' = f(t,y)}. The function 1081must accept two inputs where the first is time @var{t} and the second is a 1082column vector of unknowns @var{y}. 1083 1084@var{trange} specifies the time interval over which the ODE will be 1085evaluated. Typically, it is a two-element vector specifying the initial and 1086final times (@code{[tinit, tfinal]}). If there are more than two elements 1087then the solution will also be evaluated at these intermediate time 1088instances. 1089 1090By default, @code{ode45} uses an adaptive timestep with the 1091@code{integrate_adaptive} algorithm. The tolerance for the timestep 1092computation may be changed by using the options @qcode{"RelTol"} and 1093@qcode{"AbsTol"}. 1094 1095@var{init} contains the initial value for the unknowns. If it is a row 1096vector then the solution @var{y} will be a matrix in which each column is 1097the solution for the corresponding initial value in @var{init}. 1098 1099The optional fourth argument @var{ode_opt} specifies non-default options to 1100the ODE solver. It is a structure generated by @code{odeset}. 1101 1102The function typically returns two outputs. Variable @var{t} is a 1103column vector and contains the times where the solution was found. The 1104output @var{y} is a matrix in which each column refers to a different 1105unknown of the problem and each row corresponds to a time in @var{t}. 1106 1107The output can also be returned as a structure @var{solution} which has a 1108field @var{x} containing a row vector of times where the solution was 1109evaluated and a field @var{y} containing the solution matrix such that each 1110column corresponds to a time in @var{x}. Use 1111@w{@code{fieldnames (@var{solution})}} to see the other fields and 1112additional information returned. 1113 1114If no output arguments are requested, and no @qcode{"OutputFcn"} is 1115specified in @var{ode_opt}, then the @qcode{"OutputFcn"} is set to 1116@code{odeplot} and the results of the solver are plotted immediately. 1117 1118If using the @qcode{"Events"} option then three additional outputs may be 1119returned. @var{te} holds the time when an Event function returned a zero. 1120@var{ye} holds the value of the solution at time @var{te}. @var{ie} 1121contains an index indicating which Event function was triggered in the case 1122of multiple Event functions. 1123 1124Example: Solve the @nospell{Van der Pol} equation 1125 1126@example 1127@group 1128fvdp = @@(@var{t},@var{y}) [@var{y}(2); (1 - @var{y}(1)^2) * @var{y}(2) - @var{y}(1)]; 1129[@var{t},@var{y}] = ode45 (fvdp, [0, 20], [2, 0]); 1130@end group 1131@end example 1132@xseealso{@ref{XREFodeset,,odeset}, @ref{XREFodeget,,odeget}, @ref{XREFode23,,ode23}, @ref{XREFode15s,,ode15s}} 1133@end deftypefn 1134 1135 1136@c ode23 scripts/ode/ode23.m 1137@anchor{XREFode23} 1138@deftypefn {} {[@var{t}, @var{y}] =} ode23 (@var{fun}, @var{trange}, @var{init}) 1139@deftypefnx {} {[@var{t}, @var{y}] =} ode23 (@var{fun}, @var{trange}, @var{init}, @var{ode_opt}) 1140@deftypefnx {} {[@var{t}, @var{y}, @var{te}, @var{ye}, @var{ie}] =} ode23 (@dots{}) 1141@deftypefnx {} {@var{solution} =} ode23 (@dots{}) 1142@deftypefnx {} {} ode23 (@dots{}) 1143 1144Solve a set of non-stiff Ordinary Differential Equations (non-stiff ODEs) 1145with the well known explicit @nospell{Bogacki-Shampine} method of order 3. 1146 1147@var{fun} is a function handle, inline function, or string containing the 1148name of the function that defines the ODE: @code{y' = f(t,y)}. The function 1149must accept two inputs where the first is time @var{t} and the second is a 1150column vector of unknowns @var{y}. 1151 1152@var{trange} specifies the time interval over which the ODE will be 1153evaluated. Typically, it is a two-element vector specifying the initial and 1154final times (@code{[tinit, tfinal]}). If there are more than two elements 1155then the solution will also be evaluated at these intermediate time 1156instances. 1157 1158By default, @code{ode23} uses an adaptive timestep with the 1159@code{integrate_adaptive} algorithm. The tolerance for the timestep 1160computation may be changed by using the options @qcode{"RelTol"} and 1161@qcode{"AbsTol"}. 1162 1163@var{init} contains the initial value for the unknowns. If it is a row 1164vector then the solution @var{y} will be a matrix in which each column is 1165the solution for the corresponding initial value in @var{init}. 1166 1167The optional fourth argument @var{ode_opt} specifies non-default options to 1168the ODE solver. It is a structure generated by @code{odeset}. 1169 1170The function typically returns two outputs. Variable @var{t} is a 1171column vector and contains the times where the solution was found. The 1172output @var{y} is a matrix in which each column refers to a different 1173unknown of the problem and each row corresponds to a time in @var{t}. 1174 1175The output can also be returned as a structure @var{solution} which has a 1176field @var{x} containing a row vector of times where the solution was 1177evaluated and a field @var{y} containing the solution matrix such that each 1178column corresponds to a time in @var{x}. Use 1179@w{@code{fieldnames (@var{solution})}} to see the other fields and 1180additional information returned. 1181 1182If no output arguments are requested, and no @qcode{"OutputFcn"} is 1183specified in @var{ode_opt}, then the @qcode{"OutputFcn"} is set to 1184@code{odeplot} and the results of the solver are plotted immediately. 1185 1186If using the @qcode{"Events"} option then three additional outputs may be 1187returned. @var{te} holds the time when an Event function returned a zero. 1188@var{ye} holds the value of the solution at time @var{te}. @var{ie} 1189contains an index indicating which Event function was triggered in the case 1190of multiple Event functions. 1191 1192Example: Solve the @nospell{Van der Pol} equation 1193 1194@example 1195@group 1196fvdp = @@(@var{t},@var{y}) [@var{y}(2); (1 - @var{y}(1)^2) * @var{y}(2) - @var{y}(1)]; 1197[@var{t},@var{y}] = ode23 (fvdp, [0, 20], [2, 0]); 1198@end group 1199@end example 1200 1201Reference: For the definition of this method see 1202@url{https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods}. 1203@xseealso{@ref{XREFodeset,,odeset}, @ref{XREFodeget,,odeget}, @ref{XREFode45,,ode45}, @ref{XREFode15s,,ode15s}} 1204@end deftypefn 1205 1206 1207@c ode23s scripts/ode/ode23s.m 1208@anchor{XREFode23s} 1209@deftypefn {} {[@var{t}, @var{y}] =} ode23s (@var{fun}, @var{trange}, @var{init}) 1210@deftypefnx {} {[@var{t}, @var{y}] =} ode23s (@var{fun}, @var{trange}, @var{init}, @var{ode_opt}) 1211@deftypefnx {} {[@var{t}, @var{y}] =} ode23s (@dots{}, @var{par1}, @var{par2}, @dots{}) 1212@deftypefnx {} {[@var{t}, @var{y}, @var{te}, @var{ye}, @var{ie}] =} ode23s (@dots{}) 1213@deftypefnx {} {@var{solution} =} ode23s (@dots{}) 1214 1215Solve a set of stiff Ordinary Differential Equations (stiff ODEs) with a 1216@nospell{Rosenbrock} method of order (2,3). 1217 1218@var{fun} is a function handle, inline function, or string containing the 1219name of the function that defines the ODE: @code{M y' = f(t,y)}. The 1220function must accept two inputs where the first is time @var{t} and the 1221second is a column vector of unknowns @var{y}. @var{M} is a constant mass 1222matrix, non-singular and possibly sparse. Set the field @qcode{"Mass"} in 1223@var{odeopts} using @var{odeset} to specify a mass matrix. 1224 1225@var{trange} specifies the time interval over which the ODE will be 1226evaluated. Typically, it is a two-element vector specifying the initial 1227and final times (@code{[tinit, tfinal]}). If there are more than two 1228elements then the solution will also be evaluated at these intermediate 1229time instances using an interpolation procedure of the same order as the 1230one of the solver. 1231 1232By default, @code{ode23s} uses an adaptive timestep with the 1233@code{integrate_adaptive} algorithm. The tolerance for the timestep 1234computation may be changed by using the options @qcode{"RelTol"} and 1235@qcode{"AbsTol"}. 1236 1237@var{init} contains the initial value for the unknowns. If it is a row 1238vector then the solution @var{y} will be a matrix in which each column is 1239the solution for the corresponding initial value in @var{init}. 1240 1241The optional fourth argument @var{ode_opt} specifies non-default options to 1242the ODE solver. It is a structure generated by @code{odeset}. 1243@code{ode23s} will ignore the following options: @qcode{"BDF"}, 1244@qcode{"InitialSlope"}, @qcode{"MassSingular"}, @qcode{"MStateDependence"}, 1245@qcode{"MvPattern"}, @qcode{"MaxOrder"}, @qcode{"Non-negative"}. 1246 1247The function typically returns two outputs. Variable @var{t} is a 1248column vector and contains the times where the solution was found. The 1249output @var{y} is a matrix in which each column refers to a different 1250unknown of the problem and each row corresponds to a time in @var{t}. If 1251@var{trange} specifies intermediate time steps, only those will be returned. 1252 1253The output can also be returned as a structure @var{solution} which has a 1254field @var{x} containing a row vector of times where the solution was 1255evaluated and a field @var{y} containing the solution matrix such that each 1256column corresponds to a time in @var{x}. Use 1257@w{@code{fieldnames (@var{solution})}} to see the other fields and 1258additional information returned. 1259 1260If using the @qcode{"Events"} option then three additional outputs may be 1261returned. @var{te} holds the time when an Event function returned a zero. 1262@var{ye} holds the value of the solution at time @var{te}. @var{ie} 1263contains an index indicating which Event function was triggered in the case 1264of multiple Event functions. 1265 1266Example: Solve the stiff @nospell{Van der Pol} equation 1267 1268@example 1269@group 1270f = @@(@var{t},@var{y}) [@var{y}(2); 1000*(1 - @var{y}(1)^2) * @var{y}(2) - @var{y}(1)]; 1271opt = odeset ('Mass', [1 0; 0 1], 'MaxStep', 1e-1); 1272[vt, vy] = ode23s (f, [0 2000], [2 0], opt); 1273@end group 1274@end example 1275@xseealso{@ref{XREFodeset,,odeset}, @ref{XREFdaspk,,daspk}, @ref{XREFdassl,,dassl}} 1276@end deftypefn 1277 1278 1279@c ode15s scripts/ode/ode15s.m 1280@anchor{XREFode15s} 1281@deftypefn {} {[@var{t}, @var{y}] =} ode15s (@var{fun}, @var{trange}, @var{y0}) 1282@deftypefnx {} {[@var{t}, @var{y}] =} ode15s (@var{fun}, @var{trange}, @var{y0}, @var{ode_opt}) 1283@deftypefnx {} {[@var{t}, @var{y}, @var{te}, @var{ye}, @var{ie}] =} ode15s (@dots{}) 1284@deftypefnx {} {@var{solution} =} ode15s (@dots{}) 1285@deftypefnx {} {} ode15s (@dots{}) 1286Solve a set of stiff Ordinary Differential Equations (ODEs) or stiff 1287semi-explicit index 1 Differential Algebraic Equations (DAEs). 1288 1289@code{ode15s} uses a variable step, variable order BDF (Backward 1290Differentiation Formula) method that ranges from order 1 to 5. 1291 1292@var{fun} is a function handle, inline function, or string containing the 1293name of the function that defines the ODE: @code{y' = f(t,y)}. The function 1294must accept two inputs where the first is time @var{t} and the second is a 1295column vector of unknowns @var{y}. 1296 1297@var{trange} specifies the time interval over which the ODE will be 1298evaluated. Typically, it is a two-element vector specifying the initial and 1299final times (@code{[tinit, tfinal]}). If there are more than two elements 1300then the solution will also be evaluated at these intermediate time 1301instances. 1302 1303@var{init} contains the initial value for the unknowns. If it is a row 1304vector then the solution @var{y} will be a matrix in which each column is 1305the solution for the corresponding initial value in @var{init}. 1306 1307The optional fourth argument @var{ode_opt} specifies non-default options to 1308the ODE solver. It is a structure generated by @code{odeset}. 1309 1310The function typically returns two outputs. Variable @var{t} is a 1311column vector and contains the times where the solution was found. The 1312output @var{y} is a matrix in which each column refers to a different 1313unknown of the problem and each row corresponds to a time in @var{t}. 1314 1315The output can also be returned as a structure @var{solution} which has a 1316field @var{x} containing a row vector of times where the solution was 1317evaluated and a field @var{y} containing the solution matrix such that each 1318column corresponds to a time in @var{x}. Use 1319@w{@code{fieldnames (@var{solution})}} to see the other fields and 1320additional information returned. 1321 1322If no output arguments are requested, and no @qcode{"OutputFcn"} is 1323specified in @var{ode_opt}, then the @qcode{"OutputFcn"} is set to 1324@code{odeplot} and the results of the solver are plotted immediately. 1325 1326If using the @qcode{"Events"} option then three additional outputs may be 1327returned. @var{te} holds the time when an Event function returned a zero. 1328@var{ye} holds the value of the solution at time @var{te}. @var{ie} 1329contains an index indicating which Event function was triggered in the case 1330of multiple Event functions. 1331 1332Example: Solve @nospell{Robertson's} equations: 1333 1334@smallexample 1335@group 1336function r = robertson_dae (@var{t}, @var{y}) 1337 r = [ -0.04*@var{y}(1) + 1e4*@var{y}(2)*@var{y}(3) 1338 0.04*@var{y}(1) - 1e4*@var{y}(2)*@var{y}(3) - 3e7*@var{y}(2)^2 1339@var{y}(1) + @var{y}(2) + @var{y}(3) - 1 ]; 1340endfunction 1341opt = odeset ("Mass", [1 0 0; 0 1 0; 0 0 0], "MStateDependence", "none"); 1342[@var{t},@var{y}] = ode15s (@@robertson_dae, [0, 1e3], [1; 0; 0], opt); 1343@end group 1344@end smallexample 1345@xseealso{@ref{XREFdecic,,decic}, @ref{XREFodeset,,odeset}, @ref{XREFodeget,,odeget}, @ref{XREFode23,,ode23}, @ref{XREFode45,,ode45}} 1346@end deftypefn 1347 1348 1349@c ode15i scripts/ode/ode15i.m 1350@anchor{XREFode15i} 1351@deftypefn {} {[@var{t}, @var{y}] =} ode15i (@var{fun}, @var{trange}, @var{y0}, @var{yp0}) 1352@deftypefnx {} {[@var{t}, @var{y}] =} ode15i (@var{fun}, @var{trange}, @var{y0}, @var{yp0}, @var{ode_opt}) 1353@deftypefnx {} {[@var{t}, @var{y}, @var{te}, @var{ye}, @var{ie}] =} ode15i (@dots{}) 1354@deftypefnx {} {@var{solution} =} ode15i (@dots{}) 1355@deftypefnx {} {} ode15i (@dots{}) 1356Solve a set of fully-implicit Ordinary Differential Equations (ODEs) or 1357index 1 Differential Algebraic Equations (DAEs). 1358 1359@code{ode15i} uses a variable step, variable order BDF (Backward 1360Differentiation Formula) method that ranges from order 1 to 5. 1361 1362@var{fun} is a function handle, inline function, or string containing the 1363name of the function that defines the ODE: @code{0 = f(t,y,yp)}. The 1364function must accept three inputs where the first is time @var{t}, the 1365second is the function value @var{y} (a column vector), and the third 1366is the derivative value @var{yp} (a column vector). 1367 1368@var{trange} specifies the time interval over which the ODE will be 1369evaluated. Typically, it is a two-element vector specifying the initial and 1370final times (@code{[tinit, tfinal]}). If there are more than two elements 1371then the solution will also be evaluated at these intermediate time 1372instances. 1373 1374@var{y0} and @var{yp0} contain the initial values for the unknowns @var{y} 1375and @var{yp}. If they are row vectors then the solution @var{y} will be a 1376matrix in which each column is the solution for the corresponding initial 1377value in @var{y0} and @var{yp0}. 1378 1379@var{y0} and @var{yp0} must be consistent initial conditions, meaning that 1380@code{f(t,y0,yp0) = 0} is satisfied. The function @code{decic} may be used 1381to compute consistent initial conditions given initial guesses. 1382 1383The optional fifth argument @var{ode_opt} specifies non-default options to 1384the ODE solver. It is a structure generated by @code{odeset}. 1385 1386The function typically returns two outputs. Variable @var{t} is a 1387column vector and contains the times where the solution was found. The 1388output @var{y} is a matrix in which each column refers to a different 1389unknown of the problem and each row corresponds to a time in @var{t}. 1390 1391The output can also be returned as a structure @var{solution} which has a 1392field @var{x} containing a row vector of times where the solution was 1393evaluated and a field @var{y} containing the solution matrix such that each 1394column corresponds to a time in @var{x}. Use 1395@w{@code{fieldnames (@var{solution})}} to see the other fields and 1396additional information returned. 1397 1398If no output arguments are requested, and no @qcode{"OutputFcn"} is 1399specified in @var{ode_opt}, then the @qcode{"OutputFcn"} is set to 1400@code{odeplot} and the results of the solver are plotted immediately. 1401 1402If using the @qcode{"Events"} option then three additional outputs may be 1403returned. @var{te} holds the time when an Event function returned a zero. 1404@var{ye} holds the value of the solution at time @var{te}. @var{ie} 1405contains an index indicating which Event function was triggered in the case 1406of multiple Event functions. 1407 1408Example: Solve @nospell{Robertson's} equations: 1409 1410@smallexample 1411@group 1412function r = robertson_dae (@var{t}, @var{y}, @var{yp}) 1413 r = [ -(@var{yp}(1) + 0.04*@var{y}(1) - 1e4*@var{y}(2)*@var{y}(3)) 1414 -(@var{yp}(2) - 0.04*@var{y}(1) + 1e4*@var{y}(2)*@var{y}(3) + 3e7*@var{y}(2)^2) 1415@var{y}(1) + @var{y}(2) + @var{y}(3) - 1 ]; 1416endfunction 1417[@var{t},@var{y}] = ode15i (@@robertson_dae, [0, 1e3], [1; 0; 0], [-1e-4; 1e-4; 0]); 1418@end group 1419@end smallexample 1420@xseealso{@ref{XREFdecic,,decic}, @ref{XREFodeset,,odeset}, @ref{XREFodeget,,odeget}} 1421@end deftypefn 1422 1423 1424@c decic scripts/ode/decic.m 1425@anchor{XREFdecic} 1426@deftypefn {} {[@var{y0_new}, @var{yp0_new}] =} decic (@var{fun}, @var{t0}, @var{y0}, @var{fixed_y0}, @var{yp0}, @var{fixed_yp0}) 1427@deftypefnx {} {[@var{y0_new}, @var{yp0_new}] =} decic (@var{fun}, @var{t0}, @var{y0}, @var{fixed_y0}, @var{yp0}, @var{fixed_yp0}, @var{options}) 1428@deftypefnx {} {[@var{y0_new}, @var{yp0_new}, @var{resnorm}] =} decic (@dots{}) 1429 1430Compute consistent implicit ODE initial conditions @var{y0_new} and 1431@var{yp0_new} given initial guesses @var{y0} and @var{yp0}. 1432 1433A maximum of @code{length (@var{y0})} components between @var{fixed_y0} and 1434@var{fixed_yp0} may be chosen as fixed values. 1435 1436@var{fun} is a function handle. The function must accept three inputs where 1437the first is time @var{t}, the second is a column vector of unknowns 1438@var{y}, and the third is a column vector of unknowns @var{yp}. 1439 1440@var{t0} is the initial time such that 1441@code{@var{fun}(@var{t0}, @var{y0_new}, @var{yp0_new}) = 0}, specified as a 1442scalar. 1443 1444@var{y0} is a vector used as the initial guess for @var{y}. 1445 1446@var{fixed_y0} is a vector which specifies the components of @var{y0} to 1447hold fixed. Choose a maximum of @code{length (@var{y0})} components between 1448@var{fixed_y0} and @var{fixed_yp0} as fixed values. 1449Set @var{fixed_y0}(i) component to 1 if you want to fix the value of 1450@var{y0}(i). 1451Set @var{fixed_y0}(i) component to 0 if you want to allow the value of 1452@var{y0}(i) to change. 1453 1454@var{yp0} is a vector used as the initial guess for @var{yp}. 1455 1456@var{fixed_yp0} is a vector which specifies the components of @var{yp0} to 1457hold fixed. Choose a maximum of @code{length (@var{yp0})} components 1458between @var{fixed_y0} and @var{fixed_yp0} as fixed values. 1459Set @var{fixed_yp0}(i) component to 1 if you want to fix the value of 1460@var{yp0}(i). 1461Set @var{fixed_yp0}(i) component to 0 if you want to allow the value of 1462@var{yp0}(i) to change. 1463 1464The optional seventh argument @var{options} is a structure array. Use 1465@code{odeset} to generate this structure. The relevant options are 1466@code{RelTol} and @code{AbsTol} which specify the error thresholds used to 1467compute the initial conditions. 1468 1469The function typically returns two outputs. Variable @var{y0_new} is a 1470column vector and contains the consistent initial value of @var{y}. The 1471output @var{yp0_new} is a column vector and contains the consistent initial 1472value of @var{yp}. 1473 1474The optional third output @var{resnorm} is the norm of the vector of 1475residuals. If @var{resnorm} is small, @code{decic} has successfully 1476computed the initial conditions. If the value of @var{resnorm} is large, 1477use @code{RelTol} and @code{AbsTol} to adjust it. 1478 1479Example: Compute initial conditions for @nospell{Robertson's} equations: 1480 1481@smallexample 1482@group 1483function r = robertson_dae (@var{t}, @var{y}, @var{yp}) 1484 r = [ -(@var{yp}(1) + 0.04*@var{y}(1) - 1e4*@var{y}(2)*@var{y}(3)) 1485 -(@var{yp}(2) - 0.04*@var{y}(1) + 1e4*@var{y}(2)*@var{y}(3) + 3e7*@var{y}(2)^2) 1486@var{y}(1) + @var{y}(2) + @var{y}(3) - 1 ]; 1487endfunction 1488@end group 1489[@var{y0_new},@var{yp0_new}] = decic (@@robertson_dae, 0, [1; 0; 0], [1; 1; 0], 1490[-1e-4; 1; 0], [0; 0; 0]); 1491@end smallexample 1492@xseealso{@ref{XREFode15i,,ode15i}, @ref{XREFodeset,,odeset}} 1493@end deftypefn 1494 1495 1496@c odeset scripts/ode/odeset.m 1497@anchor{XREFodeset} 1498@deftypefn {} {@var{odestruct} =} odeset () 1499@deftypefnx {} {@var{odestruct} =} odeset (@var{"field1"}, @var{value1}, @var{"field2"}, @var{value2}, @dots{}) 1500@deftypefnx {} {@var{odestruct} =} odeset (@var{oldstruct}, @var{"field1"}, @var{value1}, @var{"field2"}, @var{value2}, @dots{}) 1501@deftypefnx {} {@var{odestruct} =} odeset (@var{oldstruct}, @var{newstruct}) 1502@deftypefnx {} {} odeset () 1503 1504Create or modify an ODE options structure. 1505 1506When called with no input argument and one output argument, return a new ODE 1507options structure that contains all possible fields initialized to their 1508default values. If no output argument is requested, display a list of 1509the common ODE solver options along with their default value. 1510 1511If called with name-value input argument pairs @var{"field1"}, 1512@var{"value1"}, @var{"field2"}, @var{"value2"}, @dots{} return a new 1513ODE options structure with all the most common option fields 1514initialized, @strong{and} set the values of the fields @var{"field1"}, 1515@var{"field2"}, @dots{} to the values @var{value1}, @var{value2}, 1516@dots{}. 1517 1518If called with an input structure @var{oldstruct} then overwrite the 1519values of the options @var{"field1"}, @var{"field2"}, @dots{} with 1520new values @var{value1}, @var{value2}, @dots{} and return the 1521modified structure. 1522 1523When called with two input ODE options structures @var{oldstruct} and 1524@var{newstruct} overwrite all values from the structure 1525@var{oldstruct} with new values from the structure @var{newstruct}. 1526Empty values in @var{newstruct} will not overwrite values in 1527@var{oldstruct}. 1528 1529The most commonly used ODE options, which are always assigned a value 1530by @code{odeset}, are the following: 1531 1532@table @asis 1533@item @code{AbsTol}: positive scalar | vector, def. @code{1e-6} 1534Absolute error tolerance. 1535 1536@item @code{BDF}: @{@qcode{"off"}@} | @qcode{"on"} 1537Use BDF formulas in implicit multistep methods. 1538@emph{Note}: This option is not yet implemented. 1539 1540@item @code{Events}: function_handle 1541Event function. An event function must have the form 1542@code{[value, isterminal, direction] = my_events_f (t, y)} 1543 1544@item @code{InitialSlope}: vector 1545Consistent initial slope vector for DAE solvers. 1546 1547@item @code{InitialStep}: positive scalar 1548Initial time step size. 1549 1550@item @code{Jacobian}: matrix | function_handle 1551Jacobian matrix, specified as a constant matrix or a function of 1552time and state. 1553 1554@item @code{JConstant}: @{@qcode{"off"}@} | @qcode{"on"} 1555Specify whether the Jacobian is a constant matrix or depends on the 1556state. 1557 1558@item @code{JPattern}: sparse matrix 1559If the Jacobian matrix is sparse and non-constant but maintains a 1560constant sparsity pattern, specify the sparsity pattern. 1561 1562@item @code{Mass}: matrix | function_handle 1563Mass matrix, specified as a constant matrix or a function of 1564time and state. 1565 1566@item @code{MassSingular}: @{@qcode{"maybe"}@} | @qcode{"yes"} | @qcode{"on"} 1567Specify whether the mass matrix is singular. 1568 1569@item @code{MaxOrder}: @{@qcode{5}@} | @qcode{4} | @qcode{3} | @qcode{2} | @qcode{1} 1570Maximum order of formula. 1571 1572@item @code{MaxStep}: positive scalar 1573Maximum time step value. 1574 1575@item @code{MStateDependence}: @{@qcode{"weak"}@} | @qcode{"none"} | @qcode{"strong"} 1576Specify whether the mass matrix depends on the state or only on time. 1577 1578@item @code{MvPattern}: sparse matrix 1579If the mass matrix is sparse and non-constant but maintains a 1580constant sparsity pattern, specify the sparsity pattern. 1581@emph{Note}: This option is not yet implemented. 1582 1583@item @code{NonNegative}: scalar | vector 1584Specify elements of the state vector that are expected to remain 1585non-negative during the simulation. 1586 1587@item @code{NormControl}: @{@qcode{"off"}@} | @qcode{"on"} 1588Control error relative to the 2-norm of the solution, rather than its 1589absolute value. 1590 1591@item @code{OutputFcn}: function_handle 1592Function to monitor the state during the simulation. For the form of 1593the function to use see @code{odeplot}. 1594 1595@item @code{OutputSel}: scalar | vector 1596Indices of elements of the state vector to be passed to the output 1597monitoring function. 1598 1599@item @code{Refine}: positive scalar 1600Specify whether output should be returned only at the end of each 1601time step or also at intermediate time instances. The value should be 1602a scalar indicating the number of equally spaced time points to use 1603within each timestep at which to return output. 1604@emph{Note}: This option is not yet implemented. 1605 1606@item @code{RelTol}: positive scalar 1607Relative error tolerance. 1608 1609@item @code{Stats}: @{@qcode{"off"}@} | @qcode{"on"} 1610Print solver statistics after simulation. 1611 1612@item @code{Vectorized}: @{@qcode{"off"}@} | @qcode{"on"} 1613Specify whether @code{odefun} can be passed multiple values of the 1614state at once. 1615 1616@end table 1617 1618Field names that are not in the above list are also accepted and 1619added to the result structure. 1620 1621@xseealso{@ref{XREFodeget,,odeget}} 1622@end deftypefn 1623 1624 1625@c odeget scripts/ode/odeget.m 1626@anchor{XREFodeget} 1627@deftypefn {} {@var{val} =} odeget (@var{ode_opt}, @var{field}) 1628@deftypefnx {} {@var{val} =} odeget (@var{ode_opt}, @var{field}, @var{default}) 1629 1630Query the value of the property @var{field} in the ODE options structure 1631@var{ode_opt}. 1632 1633If called with two input arguments and the first input argument 1634@var{ode_opt} is an ODE option structure and the second input argument 1635@var{field} is a string specifying an option name, then return the option 1636value @var{val} corresponding to @var{field} from @var{ode_opt}. 1637 1638If called with an optional third input argument, and @var{field} is 1639not set in the structure @var{ode_opt}, then return the default value 1640@var{default} instead. 1641@xseealso{@ref{XREFodeset,,odeset}} 1642@end deftypefn 1643 1644 1645@c odeplot scripts/ode/odeplot.m 1646@anchor{XREFodeplot} 1647@deftypefn {} {@var{stop_solve} =} odeplot (@var{t}, @var{y}, @var{flag}) 1648 1649Open a new figure window and plot the solution of an ode problem at each 1650time step during the integration. 1651 1652The types and values of the input parameters @var{t} and @var{y} depend on 1653the input @var{flag} that is of type string. Valid values of @var{flag} 1654are: 1655 1656@table @option 1657@item @qcode{"init"} 1658The input @var{t} must be a column vector of length 2 with the first and 1659last time step (@code{[@var{tfirst} @var{tlast}]}. The input @var{y} 1660contains the initial conditions for the ode problem (@var{y0}). 1661 1662@item @qcode{""} 1663The input @var{t} must be a scalar double specifying the time for which 1664the solution in input @var{y} was calculated. 1665 1666@item @qcode{"done"} 1667The inputs should be empty, but are ignored if they are present. 1668@end table 1669 1670@code{odeplot} always returns false, i.e., don't stop the ode solver. 1671 1672Example: solve an anonymous implementation of the 1673@nospell{@qcode{"Van der Pol"}} equation and display the results while 1674solving. 1675 1676@example 1677@group 1678fvdp = @@(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; 1679 1680opt = odeset ("OutputFcn", @@odeplot, "RelTol", 1e-6); 1681sol = ode45 (fvdp, [0 20], [2 0], opt); 1682@end group 1683@end example 1684 1685Background Information: 1686This function is called by an ode solver function if it was specified in 1687the @qcode{"OutputFcn"} property of an options structure created with 1688@code{odeset}. The ode solver will initially call the function with the 1689syntax @code{odeplot ([@var{tfirst}, @var{tlast}], @var{y0}, "init")}. The 1690function initializes internal variables, creates a new figure window, and 1691sets the x limits of the plot. Subsequently, at each time step during the 1692integration the ode solver calls @code{odeplot (@var{t}, @var{y}, [])}. 1693At the end of the solution the ode solver calls 1694@code{odeplot ([], [], "done")} so that odeplot can perform any clean-up 1695actions required. 1696@xseealso{@ref{XREFodeset,,odeset}, @ref{XREFodeget,,odeget}, @ref{XREFode23,,ode23}, @ref{XREFode45,,ode45}} 1697@end deftypefn 1698 1699