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, numericalio-pkg, odepack-pkg
  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, odepack-pkg
  38 @subsection Getting Started with ODEPACK
  39
  40 Of the eight variants of the solver, Maxima currently only has an
  41 interface the @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 t; the  dependent variables are y1, y2,
  52 and 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 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, operatingsystem-pkg, Getting Started with ODEPACK, odepack-pkg
 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}
 172
 173 @end deffn
 174
 175 @anchor{dlsode_step}
 176 @deffn {Function} dlsode_step (@var{inity}, @var{t}, @var{tout}, @var{rtol}, @var{atol}, @var{istate}, @var{state})
 177
 178 Performs one step of the solver, returning the values of the
 179 independent and dependent variables, a success or error code.
 180
 181 @table @code
 182 @item inity
 183 For the first call (when istate = 1), the initial values
 184 @item t
 185 Current value of the independent value
 186 @item tout
 187 Next point where output is desired which must not be equal to @var{t}.
 188 @item rtol
 189 relative tolerance parameter
 190 @item atol
 191 Absolute tolerance parameter, scalar of vector.  If
 192 scalar, it applies to all dependent variables.
 193 Otherwise it must be the tolerance for each dependent
 194 variable.
 195
 196 Use @var{rtol} = 0 for pure absolute error and use @var{atol} = 0
 197 for pure relative error.
 198
 199 @item istate
 200 1 for the first call to dlsode, 2 for subsequent calls.
 201 @item state
 202 state returned by dlsode_init.
 203 @end table
 204
 205 The output is a list of the following items:
 206 @table @code
 207 @item t
 208 independent variable value
 209 @item y
 210 list of values of the dependent variables at time t.
 211 @item istate
 212 Integration status:
 213 @table @code
 214 @item 1
 215 no work because tout = tt
 216 @item 2
 217 successful result
 218 @item -1
 219 Excess work done on this call
 220 @item -2
 221 Excess accuracy requested
 222 @item -3
 223 Illegal input detected
 224 @item -4
 225  Repeated error test failures
 226 @item -5
 227 Repeated convergence failures (perhaps bad Jacobian or wrong choice of
 228 mf or tolerances)
 229 @item -6
 230 Error weight because zero during problem (solution component is
 231 vanished and atol(i) = 0.
 232 @end table
 233 @item info
 234 association list of various bits of information:
 235 @table @code
 236 @item n_steps
 237 total steps taken thus far
 238 @item n_f_eval
 239 total number of function evals
 240 @item n_j_eval
 241 total number of Jacobian evals
 242 @item method_order
 243 method order
 244 @item len_rwork
 245 Actual length used for real work array
 246 @item len_iwork
 247 Actual length used for integer work array
 248 @end table
 249 @end table
 250
 251 See also @mrefdot{dlsode_init}
 252
 253 @end deffn