doc/info/odepack.texi

   1 @c -*- Mode: texinfo -*-
   2
   3 @menu
   4 * Introduction to ODEPACK::
   5 * Getting Started with ODEPACK::
   6 * Functions and Variables for odepack::
   7 @end menu
   8
   9 @node Introduction to ODEPACK, Getting Started with ODEPACK, Package numericalio, Package odepack
  10 @section Introduction to ODEPACK
  11
  12 @quotation
  13 ODEPACK is a collection of Fortran solvers for the initial value
  14 problem for ordinary differential equation systems.  It consists of nine
  15 solvers, namely a basic solver called LSODE and eight variants of it --
  16 LSODES, LSODA, LSODAR, LSODPK, LSODKR, LSODI, LSOIBT, and LSODIS.
  17 The collection is suitable for both stiff and nonstiff systems.  It
  18 includes solvers for systems given in explicit form, dy/dt = f(t,y),
  19 and also solvers for systems given in linearly implicit form,
  20 A(t,y) dy/dt = g(t,y).  Two of the solvers use general sparse matrix
  21 solvers for the linear systems that arise.  Two others use iterative
  22 (preconditioned Krylov) methods instead of direct methods for these
  23 linear systems.  The most recent addition is LSODIS, which solves
  24 implicit problems with general sparse treatment of all matrices involved.
  25 @end quotation
  26 @footnote{From @url{https://www.netlib.org/odepack/opkd-sum}}
  27
  28 References:
  29 [1] Fortran Code is from @url{https://www.netlib.org/odepack/}
  30
  31 @opencatbox{Categories:}
  32 @category{Numerical methods}
  33 @category{Share packages}
  34 @category{Package odepack}
  35 @closecatbox
  36
  37 @node Getting Started with ODEPACK, Functions and Variables for odepack, Introduction to ODEPACK, Package odepack
  38 @subsection Getting Started with ODEPACK
  39
  40 Of the eight variants of the solver, Maxima currently only has an
  41 interface to @code{dlsode}.
  42
  43 Let's say we have this system of equations to solve:
  44 @example
  45   f1 = -.04d0*y1 + 1d4*y2*y3
  46   f3 = 3d7*y2*y2
  47   dy1/dt = f1
  48   dy2/dt = -f1 - f3
  49   dy3/dt = f3
  50 @end example
  51 The independent variable is @math{t}; the  dependent variables are @math{y1}, @math{y2},
  52 and @math{y3},
  53
  54 To start the solution, set up the differential equations to solved:
  55 @example
  56 load("dlsode");
  57 f1: -.04d0*y1 + 1d4*y2*y3$
  58 f3: 3d7*y2*y2$
  59 f2: -f1 - f3$
  60 fex: [f1, f2, f3];
  61 @end example
  62
  63 Initialize the solver, where we have selected method 21:
  64 @example
  65 (%i6) state : dlsode_init(fex, ['t,y1,y2,y3], 21);
  66 (%o6) [[f, #<Function "LAMBDA ($T $Y1 $Y2 $Y3)" @{49DAC061@}>],
  67 [vars, [t, y1, y2, y3]], [mf, 21], [neq, 3], [lrw, 58], [liw, 23], [rwork, @{Li\
  68 sp Array: #(0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
  69                0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
  70                0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
  71                0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0)@}],
  72 [iwork, @{Lisp Array: #(0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0)@}],
  73 [fjac, #<Function "LAMBDA ($T $Y1 $Y2 $Y3)" @{49D52AC9@}>]]
  74 @end example
  75 The arrays rwork and iwork carry state between calls to
  76 @mref{dlsode_step}, so they should not be modified by the user.  In
  77 fact, this state should not be modified by the user at all.
  78
  79 Now that the algorithm has been initialized we can compute solutions
  80 to the differential equation, using the @var{state} returned above.
  81
  82 For this example, we want to compute the solution at times
  83 @code{0.4*10^k} for @math{k} from 0 to 11, with the initial values of 1, 0, 0
  84 for the dependent variables and with a relative tolerance of 1d-4 and
  85 absolute tolerances of 1e-6, 1e-10, and 1d-6 for the dependent
  86 variables.
  87
  88 Then
  89 @example
  90 y: [1d0, 0d0, 0d0];
  91 t: 0d0;
  92 rtol : 1d-4;
  93 atol: [1d-6, 1d-10, 1d-6];
  94 istate: 1;
  95 t:0d0;
  96 tout:.4d0;
  97
  98 for k : 1 thru 12 do
  99   block([],
 100     result: dlsode_step(y, t, tout, rtol, atol, istate, state),
 101     printf(true, "At t = ~12,4,2e   y = ~@{~14,6,2e~@}~%", result[1], result[2]),
 102     istate : result[3],
 103     tout : tout * 10);
 104 @end example
 105
 106 This produces the output:
 107 @example
 108 At t =   4.0000e-01   y =   9.851726e-01  3.386406e-05  1.479357e-02
 109 At t =   4.0000e+00   y =   9.055142e-01  2.240418e-05  9.446344e-02
 110 At t =   4.0000e+01   y =   7.158050e-01  9.184616e-06  2.841858e-01
 111 At t =   4.0000e+02   y =   4.504846e-01  3.222434e-06  5.495122e-01
 112 At t =   4.0000e+03   y =   1.831701e-01  8.940379e-07  8.168290e-01
 113 At t =   4.0000e+04   y =   3.897016e-02  1.621193e-07  9.610297e-01
 114 At t =   4.0000e+05   y =   4.935213e-03  1.983756e-08  9.950648e-01
 115 At t =   4.0000e+06   y =   5.159269e-04  2.064759e-09  9.994841e-01
 116 At t =   4.0000e+07   y =   5.306413e-05  2.122677e-10  9.999469e-01
 117 At t =   4.0000e+08   y =   5.494530e-06  2.197824e-11  9.999945e-01
 118 At t =   4.0000e+09   y =   5.129458e-07  2.051784e-12  9.999995e-01
 119 At t =   4.0000e+10   y =  -7.170563e-08 -2.868225e-13  1.000000e+00
 120 @end example
 121
 122
 123 @node Functions and Variables for odepack, , Getting Started with ODEPACK, Package odepack
 124 @section Functions and Variables for odepack
 125
 126 @anchor{dlsode_init}
 127 @deffn {Function} dlsode_init (@var{fex}, @var{vars}, @var{method})
 128
 129 This must be called before running the solver.  This function returns
 130 a state object for use in the solver.  The user must not modify the
 131 state.
 132
 133 The ODE to be solved is given in @var{fex}, which is a list of the
 134 equations.  @var{vars} is a list of independent variable and the
 135 dependent variables.  The list of dependent variables must be in the
 136 same order as the equations if @var{fex}.  Finally, @var{method}
 137 indicates the method to be used by the solver:
 138
 139 @table @code
 140 @item 10
 141 Nonstiff (Adams) method, no Jacobian used.
 142 @item 21
 143 Stiff (BDF) method, user-supplied full Jacobian.
 144 @item 22
 145 Stiff method, internally generated full Jacobian.
 146 @end table
 147
 148 The returned state object is a list of lists.  The sublist is a list
 149 of two elements:
 150 @table @code
 151 @item f
 152 The compiled function for the ODE.
 153 @item vars
 154 The list independent and dependent variables (@var{vars}).
 155 @item mf
 156 The method to be used (@var{method}).
 157 @item neq
 158 The number of equations.
 159 @item lrw
 160 Length of the work vector for real values.
 161 @item liw
 162 Length of the work vector for integer values.
 163 @item rwork
 164 Lisp array holding the real-valued work vector.
 165 @item iwork
 166 Lisp array holding the integer-valued work vector.
 167 @item fjac
 168 Compiled analytical Jacobian of the equations
 169 @end table
 170
 171 See also @mrefdot{dlsode_step}  @xref{Getting Started with ODEPACK} for
 172 an example of usage.
 173
 174 @opencatbox{Categories:}
 175 @category{Package odepack}
 176 @closecatbox
 177
 178 @end deffn
 179
 180 @anchor{dlsode_step}
 181 @deffn {Function} dlsode_step (@var{inity}, @var{t}, @var{tout}, @var{rtol}, @var{atol}, @var{istate}, @var{state})
 182
 183 Performs one step of the solver, returning the values of the
 184 independent and dependent variables, a success or error code.
 185
 186 @table @code
 187 @item inity
 188 For the first call (when istate = 1), the initial values
 189 @item t
 190 Current value of the independent value
 191 @item tout
 192 Next point where output is desired which must not be equal to @var{t}.
 193 @item rtol
 194 relative tolerance parameter
 195 @item atol
 196 Absolute tolerance parameter, scalar of vector.  If
 197 scalar, it applies to all dependent variables.
 198 Otherwise it must be the tolerance for each dependent
 199 variable.
 200
 201 Use @var{rtol} = 0 for pure absolute error and use @var{atol} = 0
 202 for pure relative error.
 203
 204 @item istate
 205 1 for the first call to dlsode, 2 for subsequent calls.
 206 @item state
 207 state returned by dlsode_init.
 208 @end table
 209
 210 The output is a list of the following items:
 211 @table @code
 212 @item t
 213 independent variable value
 214 @item y
 215 list of values of the dependent variables at time t.
 216 @item istate
 217 Integration status:
 218 @table @code
 219 @item 1
 220 no work because tout = tt
 221 @item 2
 222 successful result
 223 @item -1
 224 Excess work done on this call
 225 @item -2
 226 Excess accuracy requested
 227 @item -3
 228 Illegal input detected
 229 @item -4
 230  Repeated error test failures
 231 @item -5
 232 Repeated convergence failures (perhaps bad Jacobian or wrong choice of
 233 mf or tolerances)
 234 @item -6
 235 Error weight because zero during problem (solution component is
 236 vanished and atol(i) = 0.
 237 @end table
 238 @item info
 239 association list of various bits of information:
 240 @table @code
 241 @item n_steps
 242 total steps taken thus far
 243 @item n_f_eval
 244 total number of function evals
 245 @item n_j_eval
 246 total number of Jacobian evals
 247 @item method_order
 248 method order
 249 @item len_rwork
 250 Actual length used for real work array
 251 @item len_iwork
 252 Actual length used for integer work array
 253 @end table
 254 @end table
 255
 256 See also @mrefdot{dlsode_init}  @xref{Getting Started with ODEPACK} for
 257 an example of usage.
 258
 259 @opencatbox{Categories:}
 260 @category{Package odepack}
 261 @closecatbox
 262
 263 @end deffn