share/contrib/Eulix/Eulix.html

   1 <?xml version="1.0" encoding="utf-8" ?>
   2 <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
   3 <html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
   4 <head>
   5 <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
   6 <meta name="generator" content="Docutils 0.13.1: http://docutils.sourceforge.net/" />
   7 <title></title>
   8 <style type="text/css">
   9
  10 /*
  11 :Author: Your Name
  12 :Contact: Your Email Address
  13 :Copyright: This stylesheet has been placed in the public domain.
  14
  15 Stylesheet for use with Docutils.
  16   background-color: #efefef }
  17 */
  18
  19 @import url(html4css1.css);
  20
  21 /* Your customizations go here.  For example: */
  22
  23 pre.literal-block, pre.doctest-block {
  24   margin-left: 2em ;
  25   margin-right: 2em ;
  26   background-color: #efefef ;
  27   font-weight: bold ;
  28   font-size: 12pt ;
  29   font-family: monospace }
  30 tt.docutils {
  31   background-color: #e0e0ff ;
  32   font-weight: bold ;
  33   font-size: 10pt }
  34
  35
  36 </style>
  37 </head>
  38 <body>
  39 <div class="document">
  40
  41
  42 <div class="section" id="eulix-a-versatile-numerical-ode-solver-of-maxima">
  43 <h1>Eulix - a versatile numerical ODE solver of Maxima</h1>
  44 <div class="section" id="properties">
  45 <h2>Properties</h2>
  46 <p>Eulix is an extrapolated, linearly implicit (by default) Euler method.
  47 It can approximate the solution of an initial value problem
  48 for a system of ordinary differential equations or for
  49 differential algebraic equations (DAEs) of index 1.
  50 More generally, it can solve an implicit system like  M y'(t)= f(t,y)
  51 where M is a constant mass matrix.</p>
  52 <p>It uses variable step sizes and variable orders.
  53 It offers dense output and root finding, i.e., an implicitly defined
  54 final time. Furthermore, it can deliver an accurate spline approximation
  55 of the solution.</p>
  56 <p>The linearly implicit version is  A-stable up to order 2  and
  57 A(alpha) stable with alpha &gt;= 89.81 up to order 8 and probably above,
  58 where A-stable = A(alpha)-stable with alpha = 90 degrees</p>
  59 <p><cite>Eulix</cite> is written in <cite>pure Maxima</cite> and all of its coefficients are
  60 rational numbers. Therefore, it can approximate the solution to arbitrary
  61 precision if <strong>fpprec</strong> is set sufficiently high <strong>and</strong> if all floating point
  62 constants in the right-hand-sides and in the initial values are specified as
  63 <cite>bfloat</cite> constants.</p>
  64 </div>
  65 <div class="section" id="credits">
  66 <h2>Credits</h2>
  67 <p>An essential part of <cite>Eulix</cite> consists of a simplified version of the FORTRAN code SEULEX which has been
  68 published in:</p>
  69 <pre class="literal-block">
  70 E. Hairer, G.Wanner
  71 Solving Ordinary Differential Equations II
  72 Stiff and Differential-Algebraic Problems
  73 Springer 1991
  74 </pre>
  75 </div>
  76 <div class="section" id="features">
  77 <h2>Features</h2>
  78 <p><cite>Eulix</cite> offers four APIs - one high level one (<cite>Eulix</cite>) which is an extension of the
  79 interface provided by <cite>rk</cite> and <cite>rkf45</cite>,
  80 and two intermediate level interfaces, <a class="reference internal" href="#eulix-table">Eulix_Table</a> and <a class="reference internal" href="#eulix-spline">Eulix_Spline</a>.
  81 Lastly, there is a low-level interface <a class="reference internal" href="#eulix-step">Eulix_Step</a> which computes a single time step.
  82 <br/> All four interfaces allow for a number of options, see the chapter on <a class="reference internal" href="#options">Options</a> below.</p>
  83 </div>
  84 <div class="section" id="eulix">
  85 <h2>Eulix</h2>
  86 <p><cite>Eulix</cite> is called as:</p>
  87 <pre class="literal-block">
  88 Eulix(Expressions,Vars,initials,range,[options])
  89 </pre>
  90 <table class="docutils field-list" frame="void" rules="none">
  91 <col class="field-name" />
  92 <col class="field-body" />
  93 <tbody valign="top">
  94 <tr class="field"><th class="field-name">Expressions:</th><td class="field-body"><p class="first">is a list of expressions (or a single expression)
  95 which give the r.h.s of the ODE system. <br/>
  96 As a special case there might be
  97 an equation with r.h.s. or l.h.s equal to 0 to designate an algebraic equation.</p>
  98 </td>
  99 </tr>
 100 <tr class="field"><th class="field-name">Vars:</th><td class="field-body"><p class="first">is a list of (dependent) variables which occur in these <cite>expressions</cite>.
 101 For each variable there has to be an initial value specified by the list <cite>initials</cite>.</p>
 102 </td>
 103 </tr>
 104 <tr class="field"><th class="field-name">initials:</th><td class="field-body"><p class="first">initial values for each of the independent variables.</p>
 105 </td>
 106 </tr>
 107 <tr class="field"><th class="field-name">range:</th><td class="field-body"><p class="first">is a list of 3 to 5 elements.</p>
 108 <table class="last docutils field-list" frame="void" rules="none">
 109 <col class="field-name" />
 110 <col class="field-body" />
 111 <tbody valign="top">
 112 <tr class="field"><th class="field-name"><tt class="docutils literal">range[1]</tt>:</th><td class="field-body">specifies the independent variable which will be called <cite>time</cite>
 113 in the following.</td>
 114 </tr>
 115 <tr class="field"><th class="field-name"><tt class="docutils literal">range[2]</tt>:</th><td class="field-body">and <tt class="docutils literal">range[3]</tt> specify the initial and final time, respectively. If <tt class="docutils literal">range[3]</tt> &lt; <tt class="docutils literal">range[2]</tt>
 116 the ODE is integrated <strong>backwards</strong> in time.</td>
 117 </tr>
 118 <tr class="field"><th class="field-name"><tt class="docutils literal">range[4]</tt>:</th><td class="field-body">is called <tt class="docutils literal">delta_t</tt> (if present). Depending on an option
 119 given below, this determines the exact or maximum distance of time values in the
 120 list of results which <cite>Eulix</cite> returns. A small value of delta_t isn't expensive
 121 since it has to evaluate a polynomial only.
 122 If <tt class="docutils literal">delta_t=0</tt> only the  natural time steps are returned.</td>
 123 </tr>
 124 <tr class="field"><th class="field-name"><tt class="docutils literal">range[5]</tt>:</th><td class="field-body">is a single  <cite>root condition</cite> or a list of <cite>root conditions</cite> (if present).
 125 <br/> A <cite>root condition</cite> is a real-valued expression which depends on the dependent variables
 126 (see <cite>Vars</cite>) as well as on the independent variable <tt class="docutils literal">range[1]</tt>.
 127 The integration stops if either one of the <cite>root conditions</cite> evaluates to zero
 128 or if the final time <tt class="docutils literal">range[3]</tt> is reached. A root condition implies <cite>dense output</cite>.
 129 If <tt class="docutils literal">delta_t</tt> (<tt class="docutils literal">range[4]</tt>) is zero it is reset to <tt class="docutils literal"><span class="pre">(range[3]-range[2])/100</span></tt> .
 130 <br/> Hint: A <cite>root condition</cite> might need quite a small value of <tt class="docutils literal">delta_t</tt> since <tt class="docutils literal">Eulix_Table</tt> might miss
 131 the sign change in a <cite>root condition</cite> otherwise.</td>
 132 </tr>
 133 </tbody>
 134 </table>
 135 </td>
 136 </tr>
 137 </tbody>
 138 </table>
 139 <p><cite>Eulix</cite> has four different return values.</p>
 140 <p>If there are <strong>no</strong> <cite>root conditions</cite>, <cite>Eulix</cite> returns a list of lists. Each of these
 141 lists consist of a time value and the values of the dependent variables at this
 142 point of time - except if the option <tt class="docutils literal">tabular=none</tt> (see <a class="reference internal" href="#options">Options</a> below).
 143 <br/> If the option <cite>combined_t_y_list=false</cite> is given, <cite>Eulix</cite>
 144 returns a list of two lists, the first of which is a list of time values, the second
 145 of which is a list of lists of the values of the dependent variables at the corresponding times.</p>
 146 <p>If there <strong>are</strong> <cite>root conditions</cite>, <cite>Eulix</cite> returns a list of a <cite>solution list</cite>
 147 and what is returned if there are no <cite>root conditions</cite>.
 148 The <cite>solution list</cite> contains three items: the index of the <cite>root condition</cite> (starting with 1)
 149 which has triggered (i.e. went to zero), the time at which this has happened and
 150 the state (a list of the values of the dependent variables) at this point of time.
 151 <br/> If <cite>root conditions</cite> have been specified but none of which has triggered before
 152 the final time <tt class="docutils literal">range[3]</tt> has been reached, the <cite>solution list</cite> is a list
 153 of three zeros.</p>
 154 <p>In the following we call a matrix with a single column a <cite>vector</cite>.</p>
 155 <p>For the remaining three interfaces you have to specify the right-hand-sides, as well as the root
 156 conditions as vector-valued functiona of two arguments: time and the state vector. Furthermore,
 157 you have to give a vector-valued function which returns the value of the derivative of the
 158 right-hand-side w.r.t. time. In addition, you have to give a Jacobian matrix of the
 159 right-hand-side; this might be an approximation only. The mass-matrix (even if trivial) has to
 160 specified by the option <tt class="docutils literal">mass_matrix</tt> (see <a class="reference internal" href="#options">Options</a> below).</p>
 161 </div>
 162 <div class="section" id="id1">
 163 <span id="eulix-table"></span><h2>Eulix_Table</h2>
 164 <p><cite>Eulix_Table</cite> is called as:</p>
 165 <pre class="literal-block">
 166 Eulix_Table(t0,y0,delta_t,t_end,Rhs,Rhs_t,Rhs_jac,Roots,[options])
 167 </pre>
 168 <table class="docutils field-list" frame="void" rules="none">
 169 <col class="field-name" />
 170 <col class="field-body" />
 171 <tbody valign="top">
 172 <tr class="field"><th class="field-name">t0:</th><td class="field-body">initial time</td>
 173 </tr>
 174 <tr class="field"><th class="field-name">y0:</th><td class="field-body">initial state (a <cite>vector</cite>)</td>
 175 </tr>
 176 <tr class="field"><th class="field-name">delta_t:</th><td class="field-body">time increment for the generated list of solutions values. Due to <cite>dense output</cite>
 177 this is independent of the time step used by the integrator. A small value of delta_t
 178 isn't expensive since it has to evaluate a polynomial only.
 179 If <cite>delta_t</cite> is zero, only   the natural time steps or even only the final state are output,
 180 depending on the  option <tt class="docutils literal">tabular</tt> (see  <a class="reference internal" href="#options">Options</a> below).</td>
 181 </tr>
 182 <tr class="field"><th class="field-name">t_end:</th><td class="field-body">final time. If <tt class="docutils literal">t_end</tt> &lt; <tt class="docutils literal">t0</tt> a negative time step is used. In that case <tt class="docutils literal">delta_t</tt>
 183 must be negative, as well.</td>
 184 </tr>
 185 <tr class="field"><th class="field-name">Rhs:</th><td class="field-body">a vector-valued function of (t,y), where y is the state vector. It defines the right-hand-side
 186 of the ODE system</td>
 187 </tr>
 188 <tr class="field"><th class="field-name">Rhs_t:</th><td class="field-body">a vector-valued function of (t,y), which gives the partial derivative w.r.t. t</td>
 189 </tr>
 190 <tr class="field"><th class="field-name">Rhs_jac:</th><td class="field-body">a matrix-valued function of (t,y), which gives the exact or approximate Jacobian
 191 of the right-hand-side w.r.t. the state vector y</td>
 192 </tr>
 193 <tr class="field"><th class="field-name">Roots:</th><td class="field-body">a list of scalar-valued functions R(t,y). If R(t,y(t))=0 the integration stops.
 194 A <cite>root condition</cite> implies dense output. If <tt class="docutils literal">delta_t</tt> isn't specified or is set to zero,
 195 it is reset to <tt class="docutils literal"><span class="pre">(range[3]-range[2])/100</span></tt> .
 196 <br/> Hint: A <cite>root condition</cite> might need quite a small value of <tt class="docutils literal">delta_t</tt> since Eulix_Table might miss
 197 the sign change in a <cite>root condition</cite> otherwise.</td>
 198 </tr>
 199 <tr class="field"><th class="field-name">options:</th><td class="field-body">see <a class="reference internal" href="#options">Options</a> below</td>
 200 </tr>
 201 <tr class="field"><th class="field-name">Output:</th><td class="field-body">if <cite>Roots</cite> is an empty list, <cite>Eulix_Table</cite> returns a list of two lists,
 202 <tt class="docutils literal">tlist</tt> and <tt class="docutils literal">ylist</tt> which contain the points of time and the state vectors at these times, resp.
 203 - except if the option <tt class="docutils literal">tabular=none</tt>, when only <strong>y_final</strong> is returned.
 204 <br/> if <cite>Roots</cite> is <strong>not</strong> empty, it returns a list <strong>[Solution,[tlist,ylist]]</strong> where
 205 <cite>Solution</cite> is a list of three values: the index of the <cite>root condition</cite> which triggered the event,
 206 the time and state of this event. If none of the <cite>root conditions</cite> has triggered, the first element
 207 of <cite>Solution</cite> is zero - except if the option <tt class="docutils literal">tabular=none</tt> when only <strong>[Solution,y_final]</strong>
 208 is returned</td>
 209 </tr>
 210 </tbody>
 211 </table>
 212 </div>
 213 <div class="section" id="id2">
 214 <span id="eulix-spline"></span><h2>Eulix_Spline</h2>
 215 <p><cite>Eulix_Spline</cite> is called as:</p>
 216 <pre class="literal-block">
 217 Eulix_Spline(t0,y0,t_end,Rhs,Rhs_t,Rhs_jac,[options])
 218 </pre>
 219 <table class="docutils field-list" frame="void" rules="none">
 220 <col class="field-name" />
 221 <col class="field-body" />
 222 <tbody valign="top">
 223 <tr class="field"><th class="field-name">parameters:</th><td class="field-body">of the same name are the same as those in <cite>Eulix_Table</cite> above.</td>
 224 </tr>
 225 <tr class="field"><th class="field-name">options:</th><td class="field-body">see <a class="reference internal" href="#options">Options</a> below</td>
 226 </tr>
 227 <tr class="field"><th class="field-name">Output:</th><td class="field-body"><cite>Eulix_Spline</cite> returns a list of objects of struct <tt class="docutils literal">Eulix_Spline_Type</tt>. This is defined as
 228 <br/> <tt class="docutils literal">defstruct(Eulix_Spline_Type(t_left,t_right,NewtonPoly))</tt>
 229 <br/> Each such object defines a polynomial (in <tt class="docutils literal">t</tt>) which represents the solution of the
 230 ODE / DAE in the time interval <tt class="docutils literal">[t_left,t_right]</tt>.
 231 <br/> <cite>Eulix</cite> provides a function &nbsp; <tt class="docutils literal">Eulix_SplineEval(t,Spls)</tt> &nbsp; which evaluates the spline at a point <tt class="docutils literal">t</tt>.
 232 Here, <tt class="docutils literal">Spls</tt> is the list returned by <cite>Eulix_Spline</cite>. It's possible to evaluates the
 233 derivative(s) of this spline w.r.t. time. For this, one should write a function similar to
 234 <cite>Eulix_SplineEval</cite> - see for example the function <cite>Eulix_Dense_Deriv</cite>, below.</td>
 235 </tr>
 236 </tbody>
 237 </table>
 238 </div>
 239 <div class="section" id="id3">
 240 <span id="eulix-step"></span><h2>Eulix_Step</h2>
 241 <p><cite>Eulix_Step</cite> is called as:</p>
 242 <pre class="literal-block">
 243 [DO_DQ,tn,hn,me,failed]: Eulix_Step(y,t,Rhs,Rhs_t,Rhs_jac,h,me,[options])
 244 </pre>
 245 <table class="docutils field-list" frame="void" rules="none">
 246 <col class="field-name" />
 247 <col class="field-body" />
 248 <tbody valign="top">
 249 <tr class="field"><th class="field-name">parameters:</th><td class="field-body">of the same name are the same as those in <cite>Eulix_Table</cite> above.</td>
 250 </tr>
 251 <tr class="field"><th class="field-name">h:</th><td class="field-body">step size, which must be non-zero. A negative step size is allowed.
 252 For the very first step the user must estimate a suitable step size.
 253 One might use the technique which is applied by <cite>Eulix_Table</cite> - see
 254 <cite>Automatic selection of the initial step size</cite> in that function.
 255 For any step except the first one, one should supply the value <tt class="docutils literal">hn</tt>
 256 as new step size except for the very last step (to hit a desired final time)</td>
 257 </tr>
 258 <tr class="field"><th class="field-name">me:</th><td class="field-body">width of extrapolation tableau (= order).
 259 For the very first step the user must estimate a suitable initial order,
 260 e.g. 6. For any step except the first one, one should supply the value
 261 <tt class="docutils literal">me</tt> which is returned by <cite>Eulix_Step</cite></td>
 262 </tr>
 263 <tr class="field"><th class="field-name">DO_DQ:</th><td class="field-body">this value depends on the <cite>dense_output</cite> option, see <a class="reference internal" href="#options">Options</a> below.
 264 If this is false, <tt class="docutils literal">DO_DQ[1]</tt> contains the new state vector at time <tt class="docutils literal">tn</tt>.
 265 <br/> Otherwiese <tt class="docutils literal">DO_DQ</tt> contains the dense output information.
 266 This can be passed to
 267 &nbsp; <tt class="docutils literal">Eulix_Dense_Out(dt,DO_DQ)</tt> &nbsp; or &nbsp; <tt class="docutils literal">Eulix_Dense_Deriv(dt,DO_DQ)</tt>.
 268 <br/> Both of these take a time difference <tt class="docutils literal">dt</tt> relative to the
 269 <strong>new</strong> time <tt class="docutils literal">tn</tt>. Therefore, to evaluate the solution
 270 at time <tt class="docutils literal">t1</tt> in <tt class="docutils literal">[t,tn]</tt> with <tt class="docutils literal">tn</tt> &gt; <tt class="docutils literal">t</tt>, the value
 271 <tt class="docutils literal">dt</tt> = <tt class="docutils literal">t1</tt>-<tt class="docutils literal">tn</tt> will be negative. When integrating backwards in time,
 272 i.e., <tt class="docutils literal">tn</tt> &lt; <tt class="docutils literal">t</tt>, this value will be positive.
 273 <br/> <cite>Eulix_Dense_Out</cite> evaluates <tt class="docutils literal">y(t)</tt> while <cite>Eulix_Dense_Deriv</cite>
 274 delivers <tt class="docutils literal"><span class="pre">y'(t)</span></tt>. This can be used in Newton's method as is done in
 275 <cite>Eulix_Table</cite> for determining zeros of the <cite>root conditions</cite>.</td>
 276 </tr>
 277 <tr class="field"><th class="field-name">tn:</th><td class="field-body">is the new time</td>
 278 </tr>
 279 <tr class="field"><th class="field-name">hn:</th><td class="field-body">is the new optimal size for the next step. This should be passed
 280 as parameter <tt class="docutils literal">h</tt> to <cite>Eulix_Step</cite> unless this would exceed the desired
 281 final time.</td>
 282 </tr>
 283 <tr class="field"><th class="field-name">me:</th><td class="field-body">is the new optimal tableau width for the next step. This should be
 284 passed unchanged to <cite>Eulix_Step</cite> for the next call.</td>
 285 </tr>
 286 <tr class="field"><th class="field-name">failed:</th><td class="field-body">If this is <cite>true</cite> the output parameter <cite>DO_DQ</cite> doesn't contain
 287 useful information. <cite>Eulix_Step</cite> had to reject a step more than 10 times
 288 in succession. This should only occur if the right-hand-size is
 289 discontinuous or if the values delivered by <cite>Rhs_t</cite> or <cite>Rhs_jac</cite> are wrong.</td>
 290 </tr>
 291 </tbody>
 292 </table>
 293 </div>
 294 <div class="section" id="id4">
 295 <span id="options"></span><h2>Options</h2>
 296 <p>Some options are only applicable to some of the interface functions above.
 297 The following abbreviations are used to indicate for which function an
 298 option is effective: <cite>Eul</cite> for <cite>Eulix</cite>, <cite>Tab</cite> for <cite>Eulix_Table</cite>, <cite>Spl</cite> for
 299 <cite>Eulix_Spline</cite> and <cite>Stp</cite> for <cite>Eulix_Step</cite>.
 300 <br/> The options are given in the form <cite>Option Name</cite> = <cite>default</cite>.</p>
 301 <table class="docutils field-list" frame="void" rules="none">
 302 <col class="field-name" />
 303 <col class="field-body" />
 304 <tbody valign="top">
 305 <tr class="field"><th class="field-name">absolute_tolerance:</th><td class="field-body"><p class="first">= 1E-6 &nbsp; <tt class="docutils literal">[Eul,Tab,Stp,Spl]</tt></p>
 306 </td>
 307 </tr>
 308 <tr class="field"><th class="field-name">relative_tolerance:</th><td class="field-body"><p class="first">= 1E-6 &nbsp; <tt class="docutils literal">[Eul,Tab,Stp,Spl]</tt>
 309 <br/> The (local) <strong>required tolerance</strong> is computed as (with the state vector <tt class="docutils literal">y</tt>):</p>
 310 <pre class="literal-block">
 311 absolute_tolerance + relative_tolerance *  mat_norm(y,'inf)
 312 </pre>
 313 </td>
 314 </tr>
 315 <tr class="field"><th class="field-name">root_tolerance:</th><td class="field-body"><p class="first">= 0 &nbsp; <tt class="docutils literal">[Eul,Tab]</tt> This defines the tolerance for Newton's method to
 316 determine a zero of a <cite>root condition</cite>.
 317 <br/> If this is negative or zero (equivalent to -1E-6), the tolerance is
 318 a relative tolerance according to:</p>
 319 <pre class="literal-block">
 320 abs(root_eps)*max(mat_norm(y_left,'inf),mat_norm(y_right,'inf))
 321 </pre>
 322 <p>Here, <tt class="docutils literal">y_left</tt> and <tt class="docutils literal">y_right</tt> are the state vectors at the left (right, rsp.)
 323 start interval for Newton's method.</p>
 324 </td>
 325 </tr>
 326 <tr class="field"><th class="field-name">maximum_order:</th><td class="field-body"><p class="first">= 20 &nbsp; <tt class="docutils literal">[Eul,Tab,Stp]</tt> &nbsp; This is the maximum tableau width
 327 that is used by <cite>Eulix_Step</cite></p>
 328 </td>
 329 </tr>
 330 <tr class="field"><th class="field-name">dense_output:</th><td class="field-body"><p class="first">= true &nbsp; <tt class="docutils literal">[Stp,Spl]</tt> &nbsp; This makes <cite>Eulix_Step</cite> compute the
 331 dense output information <tt class="docutils literal">DO_DQ</tt>. This can be triggered in a call to <cite>Eulix</cite>
 332 and <cite>Eulix_Table</cite> by other means.</p>
 333 </td>
 334 </tr>
 335 <tr class="field"><th class="field-name">work_estimates:</th><td class="field-body"><p class="first">= [1,1,0] &nbsp; <tt class="docutils literal">[Eul,Tab,Stp,Spl]</tt> &nbsp; For computing
 336 the optimal step size and order for the next step, <cite>Eulix_Step</cite> needs
 337 a relative amount of work to evaluate the right-hand-side <cite>Rhs</cite> and
 338 <cite>Rhs_t</cite>,   to compute the Jacobian <cite>Rhs_jac</cite> and to decompose the Jacobian
 339 into a <cite>Q-R-decomposition</cite>.</p>
 340 </td>
 341 </tr>
 342 <tr class="field"><th class="field-name">check_parameters:</th><td class="field-body"><p class="first">= true &nbsp; <tt class="docutils literal">[Eul,Tab,Stp,Spl]</tt> &nbsp; If this is true,
 343 several checks to the parameters are performed. The functions
 344 <tt class="docutils literal">Eulix_Table</tt> and <tt class="docutils literal">Eulix_Spline</tt> set this to false after the first
 345 call to <tt class="docutils literal">Eulix_Step</tt>.</p>
 346 </td>
 347 </tr>
 348 <tr class="field"><th class="field-name">mass_matrix:</th><td class="field-body"><p class="first">= <tt class="docutils literal">ident(dim)</tt> &nbsp; <tt class="docutils literal">[Stp,Spl]</tt> This parameter should take a matrix
 349 of the size of the differential system or <tt class="docutils literal">false</tt>.
 350 <br/> If it has the value <tt class="docutils literal">false</tt>, an extrapolated <strong>explicit</strong> Euler scheme
 351 is used. In this case, the function <cite>Rhs_jac</cite> will not be called.
 352 <br/> If it has the value of a (square) matrix <tt class="docutils literal">M</tt>, the (possibly)
 353 implicit ODE &nbsp; <tt class="docutils literal">M <span class="pre">y'(t)</span> = f(t,y(t))</tt> &nbsp; is solved. This matrix need not
 354 be regular. Indeed, for an algebraic equation, the corressponding row of
 355 <tt class="docutils literal">M</tt> should be zero. This can be triggered on a call to <cite>Eulix</cite>, as well.
 356 <br/> If <cite>mass_matrix</cite> has the value of a matrix, the <strong>implicit</strong> (extrapolated)
 357 Euler scheme will be used.</p>
 358 </td>
 359 </tr>
 360 <tr class="field"><th class="field-name">initial_step_size:</th><td class="field-body"><p class="first">= 0 &nbsp; <tt class="docutils literal">[Eul,Tab,Spl]</tt> defines
 361 the very first step size; for integration backwards in time, this should
 362 have a negative value. If it has the (default) value 0, an estimate
 363 is computed. This depends on the right-hand-side as well as on the option
 364 <tt class="docutils literal">absolute_tolerance</tt> - see
 365 <cite>Automatic selection of the initial step size</cite> in that function <cite>Eulix_Step</cite>.</p>
 366 </td>
 367 </tr>
 368 <tr class="field"><th class="field-name">initial_order:</th><td class="field-body"><p class="first">= 6 &nbsp; <tt class="docutils literal">[Eul,Tab,Spl]</tt> defines the initial width of the extrapolation
 369 tableau.</p>
 370 </td>
 371 </tr>
 372 <tr class="field"><th class="field-name">logging:</th><td class="field-body"><p class="first">= <tt class="docutils literal">false</tt> &nbsp; <tt class="docutils literal">[Eul,Tab,Stp,Spl]</tt> &nbsp; If this is true,
 373 some internal information is printed, e.g. rejected steps, the initial
 374 step size, options which are in effect and for each time step the values
 375 of the current time, step size and tableau width (order).</p>
 376 </td>
 377 </tr>
 378 <tr class="field"><th class="field-name">tabular:</th><td class="field-body"><p class="first">= <tt class="docutils literal">all</tt> &nbsp; <tt class="docutils literal">[Eul,Tab]</tt> &nbsp; This option can have the values
 379 <tt class="docutils literal">all</tt> (default), <tt class="docutils literal">none</tt> and <tt class="docutils literal">tabular</tt>.
 380 <br/> If it has the value <tt class="docutils literal">none</tt>, only the final state is returned
 381 or a list <strong>[Solution,y_final]</strong> when <cite>root conditions</cite> have been specified.
 382 <br/> If it has the value <tt class="docutils literal">tabular</tt> a list of solutions at <strong>equidistant</strong> times
 383 is produced. For the value <tt class="docutils literal">all</tt>, the natural time steps are intermixed, as well.</p>
 384 </td>
 385 </tr>
 386 <tr class="field"><th class="field-name">combined_t_y_list:</th><td class="field-body"><p class="first last">= <tt class="docutils literal">true</tt> &nbsp; <tt class="docutils literal">[Eul]</tt> &nbsp; If this is true,
 387 <cite>Eulix</cite> returns a <strong>single list</strong> of lists containing the time value and the
 388 corresponding solution components at this point of time.
 389 <br/> If this is false, <cite>Eulix</cite> returns &nbsp; <tt class="docutils literal">[tlist,ylist]</tt> &nbsp;, where
 390 <tt class="docutils literal">ylist</tt> is a list of lists containing the solution components at the
 391 corresponding time in <tt class="docutils literal">tlist</tt>.</p>
 392 </td>
 393 </tr>
 394 </tbody>
 395 </table>
 396 </div>
 397 <div class="section" id="examples">
 398 <h2>Examples</h2>
 399 <p>The example files are all named <tt class="docutils literal"><span class="pre">Eulix_*.mac</span></tt>. In the following list, we only
 400 give the suffix which the star stands for.</p>
 401 <table class="docutils field-list" frame="void" rules="none">
 402 <col class="field-name" />
 403 <col class="field-body" />
 404 <tbody valign="top">
 405 <tr class="field"><th class="field-name">T1:</th><td class="field-body">This is a simple 1-d ODE which requires frequent changes of the step size
 406 and order. It's partially stiff. The exact solution is known.</td>
 407 </tr>
 408 <tr class="field"><th class="field-name">T1R:</th><td class="field-body">This solves the same ODE <strong>backwards</strong> in time</td>
 409 </tr>
 410 <tr class="field"><th class="field-name">Step_T1:</th><td class="field-body">This solves the same ODE using the low level <cite>Eulix_Step</cite> interface</td>
 411 </tr>
 412 <tr class="field"><th class="field-name">T2:</th><td class="field-body">This is a simple 1-d ODE which simulates arbitrary stiffness.
 413 The exact solution is known.</td>
 414 </tr>
 415 <tr class="field"><th class="field-name">T3:</th><td class="field-body">This is the Volterra-Lotka (2-d) system which is rather demanding.</td>
 416 </tr>
 417 <tr class="field"><th class="field-name">T3N:</th><td class="field-body">This demonstrates the option <tt class="docutils literal">tabular=none</tt> for the same ODE</td>
 418 </tr>
 419 <tr class="field"><th class="field-name">Spline_T3:</th><td class="field-body">This demonstrates the <cite>Eulix_Spline</cite> interface for this ODE</td>
 420 </tr>
 421 <tr class="field"><th class="field-name">T4:</th><td class="field-body">Here we solve a differential-algebraic equation (DAE) for a pendulum</td>
 422 </tr>
 423 <tr class="field"><th class="field-name">T5:</th><td class="field-body">This is a rather demanding AIDS model which requires a good error estimator.
 424 Here we demonstrate the <cite>root finding</cite> feature.</td>
 425 </tr>
 426 <tr class="field"><th class="field-name">T5N:</th><td class="field-body">Same as above except with <tt class="docutils literal">tabular=none</tt></td>
 427 </tr>
 428 <tr class="field"><th class="field-name">Table_T5:</th><td class="field-body">This same model as <tt class="docutils literal">T5</tt> but this time using the <cite>Eulix_Table</cite>
 429 interface. Here we demonstrate the <cite>root finding</cite> feature.</td>
 430 </tr>
 431 <tr class="field"><th class="field-name">T6:</th><td class="field-body">Simple test with a highly non-linear right hand side</td>
 432 </tr>
 433 <tr class="field"><th class="field-name">T7:</th><td class="field-body">Simple test of the quality of dense output using the exact solution exp(-t*t).</td>
 434 </tr>
 435 <tr class="field"><th class="field-name">TP:</th><td class="field-body">Here we solve a system of 2 equations with a tolerance of 1E-30.
 436 Since the exact solution is known, we can check the tolerance achieved.</td>
 437 </tr>
 438 <tr class="field"><th class="field-name">TS:</th><td class="field-body">This is an extremely stiff ODE (<cite>Shampine's ball of flame</cite>)</td>
 439 </tr>
 440 <tr class="field"><th class="field-name">Step_TP:</th><td class="field-body">A 2-d ODE with a known solution, demonstrating high precision
 441 solving using the low level <cite>Eulix_Step</cite> interface.</td>
 442 </tr>
 443 <tr class="field"><th class="field-name">Table_TC:</th><td class="field-body">A very stiff system modelling a chemical reaction using the <cite>Eulix_Table</cite> interface,
 444 see Hairer/Wanner Solving ODE II .</td>
 445 </tr>
 446 </tbody>
 447 </table>
 448 <p>The examples <cite>Eulix_T1.mac</cite>, <cite>Eulix_T2.mac</cite>, <cite>Eulix_T3.mac</cite> and <cite>Eulix_TS.mac</cite>
 449 contain comparisons with <cite>rkf45</cite>.</p>
 450 </div>
 451 </div>
 452 </div>
 453 <div class="footer">
 454 <hr class="footer" />
 455 <a class="reference external" href="Eulix.rst">View document source</a>.
 456 Generated on: 2017-08-27 10:17 UTC.
 457 Generated by <a class="reference external" href="http://docutils.sourceforge.net/">Docutils</a> from <a class="reference external" href="http://docutils.sourceforge.net/rst.html">reStructuredText</a> source.
 458
 459 </div>
 460 </body>
 461 </html>