| | |
-
- Dynamics
class Dynamics() |
| |
Class for solving the initial value problem for systems of ordinary
differential equations. For uniformity all methods return - besides a
list/stack or dict of solutions at each time point - (a suggestion for)
the next time point in the time stepping, including those methods that are
dependent on externally controlled time stepping and do not have their
own time stepping. The methods of this class support state vectors which
are either dicts or lists/stacks or 'd' arrays. Method "matrixexp" is a
bit special and takes (=requires) a column vector from the misclib.Matrix
class as the input state vector.
Instances are punched out with a construct like (breakpoints is optional,
cf tiptoe function):
solution = Dynamics(model, state [, breakpoints])
where "model" is the name of the function of the time and the state vector
that returns the values of the derivatives. "state" can be a dict, a stack,
a list, a 'd' array or a Matrix column vector (the only solver method
capable of handling the latter is matrixexp). The different solvers
(methods of solution) are called like
t = tstart
while t <= tfin:
tnext = t + deltat
t, state = solution.method(t, tnext)
where "method" is the name of the solver method; euler, runge_kutta4 etc.
If the state vector is the result put out by the program using this class,
the corresponding time to put out is "tnext" - if derivatives are the
result put out, then "t" is more correct, or it may even depend on whether
the solution method is explicit or implicit...
The advantage in using dicts is - besides more clarity due to the
possibility of using text strings rather than indices to identify the
elements of the derivatives vector and the state vector - is that the
order in which the different elements of the vectors are treated in the
function called by the methods in Dynamics is not crucial!
Even if some of the solver methods belonging to this class handles
reasonably stiff problems very well, you might consider the StiffDynamics
class for problems which are very stiff. The solver methods of the
StiffDynamics class are rather slow, though.
The user if referred to Dahlquist, Bjorck & Anderson Ch. 8 for the theory
behind the majority of the methods used in this class.
|
| |
Methods defined here:
- __init__(self, model, state, breakpoints=None)
- Initiate the model and the state dict/list/array, and initiate
history lists used by some of the methods.
- abm2(self, t, tnext, tolf=1.4901161193847656e-08, maxniter=4)
- Adams-Bashforth-Moulton 2:nd order, a predictor-corrector method. The
algorithm is claimed to be more accurate for a reasonably large number
of iterations but is also claimed to be more stable for a smaller number
of iterations. Different orders of Adams-Bashforth-Moulton may have
different accuracy and stability properties - this is the reason for
three abmX methods being here. abm2 uses one computation of the
model/derivative function per iteration.
'tolf' is the maximum allowed fractional difference between the sum of
the absolute values of the state vector variables for a prediction and
a correction. 'maxniter' is the maximum number of iterations made.
The computation is NOT stopped if convergence is not reached, but a
warning will be issued to stdout.
This method uses the first time step twice to build its history at
the very outset. After that the historic points will differ from
one another.
- abm3(self, t, tnext, tolf=1.4901161193847656e-08, maxniter=4)
- Adams-Bashforth-Moulton 3:rd order, a predictor-corrector method. The
algorithm is claimed to be more accurate for a reasonably large number
of iterations but is also claimed to be more stable for a smaller number
of iterations. Different orders of Adams-Bashforth-Moulton may have
different accuracy and stability properties - this is the reason for
three abmX methods being here. abm3 uses one computation of the
model/derivative function per iteration.
'tolf' is the maximum allowed fractional difference between the sum of
the absolute values of the state vector variables for a prediction and
a correction. 'maxniter' is the maximum number of iterations made.
The computation is NOT stopped if convergence is not reached, but a
warning will be issued to stdout.
This method uses abm2 to build up its history before it has created
its own history. Cf. abm2 for details.
- abm4(self, t, tnext, tolf=1.4901161193847656e-08, maxniter=4)
- Adams-Bashforth-Moulton 4:th order, a predictor-corrector method. The
algorithm is claimed to be more accurate for a reasonably large number
of iterations but is also claimed to be more stable for a smaller number
of iterations. Different orders of Adams-Bashforth-Moulton may have
different accuracy and stability properties - this is the reason for
three abmX methods being here. abm4 uses one computation of the
model/derivative function per iteration.
'tolf' is the maximum allowed fractional difference between the sum of
the absolute values of the state vector variables for a prediction and
a correction. 'maxniter' is the maximum number of iterations made.
The computation is NOT stopped if convergence is not reached, but a
warning will be issued to stdout.
This method uses abm3 to build up its history before it has created
its own history. Cf. abm3 for details.
- adams_bashforth2(self, t, tnext)
- Adams-Bashforth 2:nd order, a linear multistep method. The method
uses the derivative-function-value history and is fast since it only
uses one computation of the model/derivative function per time step.
Some of the discretization error problems are reduced. The method
should work well for smooth problems, but adams_bashforth4 would
normally be recommended (this method is used by other methods, a
fact that is the main reason for it being here in the first place).
This method uses the first time step twice to build its history at
the very outset. After that the historic points will differ from
one another.
- adams_bashforth3(self, t, tnext)
- Adams-Bashforth 3:rd order, a linear multistep method. The method
uses the derivative-function-value history and is fast since it only
uses one computation of the model/derivative function per time step.
Some of the discretization error problems are reduced. The method
should work well for smooth problems, but adams_bashforth4 would
normally be recommended (this method is used by adams_bashforth4, a
fact that is the main reason for it being here in the first place).
This method uses adams_bashforth2 to build up its history before
it has created its own history. Cf. adams_bashforth2 for details.
- adams_bashforth4(self, t, tnext)
- Adams-Bashforth 4:th order, a linear multistep method. The method
uses the derivative-function-value history and is fast since it only
uses one computation of the model/derivative function per time step.
Some of the discretization error problems are reduced. The method
should work well for smooth problems.
This method uses adams_bashforth3 to build up its history before
it has created its own history. Cf. adams_bashforth3 for details.
- euler(self, t, tnext)
- Simple forward Euler stepping with constant time step. It is fast
since it only uses uses one computation of the model/derivative
function per time step, but the numerics may be problematic due to
the large discretization error - the Euler scheme is only first-order
correct. Trying to overcome this by using very short time steps may
produce a large accumulated round-off error. But it is fast, and
sometimes useful.
- euler_ex(self, t, tnext)
- An enhanced Euler solver with one Richardson extrapolation step.
The solution is first computed using the full time step, then twice
consequtively using the halved time step. The the solutions are
combined using the Richardson extrapolation scheme. The model
/derivative function is called twice.
- euler_plus(self, t, tnext)
- The implicit trapezoidal method
ynp1 - yn - 0.5*h*(f(tn, yn) + f(tnp1, ynp1)) = 0
used in a single-step predictor-corrector scheme using the simple
Euler scheme for the predictor. It uses two computations of the
model/derivative function per time step.
- heun(self, t, tnext)
- Heun's method is really second-order Runge-Kutta. It uses two
computations of the model/derivative function per time step.
- matrixexp(self, t, tnext, tolf=4.4408920985006262e-16, maxterms=128)
- Solver using the matrix exponential method.
This solver only works for systems of odes which are linear in the
state vector. The state vector must be entered as a matrix column
vector and coeffmatrix as a matrix containing the present (possibly
time dependent) coefficients, both objects belonging to the
misclib.Matrix class. The differential equation is: dY/deltat = C + M*Y
where Y is the state column vector, C is a matrix column vector of
constants (constant between t and t+deltat), and M is a matrix of
constant coefficients (constant between t and t+deltat).
The solution is Y(t+deltat) = C*deltat + exp(M*Y(t)*deltat), where
the exponential is expanded in a McLaurin series based on t as the
zero point in time, using 'tolf' as the fractional tolerance (based
on the fractional difference between the present sum and the previous)
and 'maxterms' as the absolute maximum number of terms.
AND: matrixexp handles stiff systems well and with good accuracy !!!!
- modmidpoint(self, t, tnext, nsteps=4)
- W.B. Gragg's modified midpoint method. It uses nsteps computations of
the model/derivative function for each time step.
- modmidpoint_ex(self, t, tnext, nsteps=2)
- The modified midpoint method using one Richardson extrapolation
step, so it really uses nsteps + 2*nsteps computations of the
model/derivative function for each time step!
- rkck45(self, t, tnext, individ=False, tolf=3.637978807091713e-12, factor=0.25, maxniter=10)
- This solver uses the Runge-Kutta-Cash-&-Karp algorithm (cf.
http://en.wikipedia.org/wiki/Cash-Karp): fourth and fifth order
Runge-Kutta with Cash & Karps's coefficients is used in an adaptive
step size scheme that uses the difference between fourth and fifth
order solutions to control the stepping. Six computations of the
model/derivative function are used in each iteration. Very
efficient, accurate and also stable for most problems.
When individ==False, 'tolf' is the maximum allowed fractional difference
between the sum of the absolute values of the state vector variables
for the fourth and fifth order solutions. When individ==True, the
comparison between the fourth and fifth order solutions is made for
each of the stage variables. 'factor' is the factor by which the
present time step is multiplied when the fractional difference is
greater than tolf. 'maxniter' is the maximum number of iterations made.
The computation is NOT stopped if convergence is not reached, but a
warning will be issued to stdout.
NB. This method uses adaptive step size control which means that it may
compute the model/derivative function for times prior to what was just
used in the same single call. This might cause problems if there are
setups in the model that uses values of the derivatives or states from
prior calls. Short time steps might make this a minor problem, but
"ten cuidado"...
- rke4(self, t, tnext)
- Runge-Kutta - England 4th order version. The algorithm is the same
as the one used in GASP IV using two halved steps, with the difference
that adaptive step size control is not used. The model/derivative
function is computed eight times per full time step.
- rkf45(self, t, tnext, individ=False, tolf=2.9103830456733704e-11, factor=0.25, maxniter=10)
- The Runge-Kutta-Fehlberg algorithm (inspired by Forsythe-Malcolm-Moler):
fourth and fifth order Runge-Kutta with Fehlberg's coefficients is used
in an adaptive step size scheme that uses the difference between fourth
and fifth order solutions to control the stepping. Six computations of
the model/derivative function are used in each iteration. Efficient,
accurate and also stable for most problems.
When individ==False, 'tolf' is the maximum allowed fractional difference
between the sum of the absolute values of the state vector variables
for the fourth and fifth order solutions. When individ==True, the
comparison between the fourth and fifth order solutions is made for
each of the stage variables. 'factor' is the factor by which the
present time step is multiplied when the fractional difference is
greater than tolf. 'maxniter' is the maximum number of iterations made.
The computation is NOT stopped if convergence is not reached, but a
warning will be issued to stdout.
NB. This method uses adaptive step size control which means that it may
compute the model/derivative function for times prior to what was just
used in the same single call. This might cause problems if there are
setups in the model that uses values of the derivatives or states from
prior calls. Short time steps might make this a minor problem, but "ten
cuidado"...
- runge_kutta4(self, t, tnext)
- Standard fourth-order Runge-Kutta. It uses four computations of
the model/derivative function per time step.
- tiptoe(self, t, tstep)
- Method used to get around known points of discontinuity or points at
which there is continuity but where the first derivative is not defined
in the equations defining a set of odes, such as step functions or
Heaviside functions for the coefficients. It avoids those points by
allowing the solution to progress in time to a point exactly on the
singularity and sets off the continued progress from a point just after
the discontinuity, using he solution at the point of the discontinuity
as a new starting point just after the discontinuity. "Just after" is
expressed in terms of machine epsilon. The method is used like this:
t, tnext = solution.tiptoe(t, tstep)
t, state = solution.method(t, tnext)
For tiptoe to be used, breakpoints must be given to the instance object
punched out from Dynamics. breakpoints may be a single float, a list or
a tuple. The points will be deleted from the list when they are passed.
- trapezoid(self, t, tnext, tolf=1.4901161193847656e-08, maxniter=4)
- The implicit trapezoidal method:
ynp1 - yn - 0.5*h*(f(tn, yn) + f(tnp1, ynp1)) = 0
used in a predictor-corrector scheme using the simple Euler
scheme for the predictor.
'tolf' is the maximum allowed fractional difference between the sum of
the absolute values of the state vector variables for a prediction and
a correction. 'maxniter' is the maximum number of iterations made.
The computation is NOT stopped if convergence is not reached, but a
warning will be issued to stdout.
Data descriptors defined here:
- __dict__
- dictionary for instance variables (if defined)
- __weakref__
- list of weak references to the object (if defined)
| |