| | |
- dynamics.Dynamics()
-
- StochDynamics
class StochDynamics(dynamics.Dynamics) |
| |
May be used when noise is added to a problem, i. e. for stochastic
differential equations. All responsibility for handling the stochastics
is placed on a "noise model" function that is used together with the
derivative model normally used by the deterministic (non-stochastic)
solver methods.
Instances are punched out with a construct like (the sequence of breakpoints
is optional, cf. the tiptoe method/function in the Dynamics class):
solution = StochDynamics(model, noisemodel, state [, breakpoints])
where "model" is the name of the function of time and the state vector
that returns the values of the derivatives. "noisemodel" is the noise model
function and "state" can be a dict, a stack, a list or a 'd' array. The
solvers (there is only one at the moment) 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_maruyama is the
only method available at the moment.
StochDynamics inherits from the Dynamics class, and all the deterministic
methods in the latter are available here, a feature which makes it simple to
switch between deterministic and stochastic mode in one single simulation.
|
| |
- Method resolution order:
- StochDynamics
- dynamics.Dynamics
Methods defined here:
- __init__(self, model, noisemodel, state, breakpoints=None)
- May be used when noise is added to a problem, i. e. for stochastic
differential equations. All responsibility for handling the stochastics
is placed on a "noisemodel" function that is used together with the
derivative model normally used by the solver methods of the solvers
of the Dynamics and StiffDynamics classes.
- euler_maruyama(self, t, tnext)
- A fairly general system of SDEs - where "state" is multi-dimensional
(=vector-valued) - may be formulated using the Ito formalism as
dstate = a(t, state)dt + B(t, state).dW + C(t, state).dJ
where W are Wiener processes and J generalized Poisson jump processes.
In principle B and C have matrix character and also contain information
on the correlation between the different individual Wiener processes
and the different generalized Poisson jump processes. B and C must be
part of and handled by "noisemodel" as must the handling of the
vectors of the random increments of W and J. "a" is vector-valued as
in the case of systems of non-stochastic ODEs. The deterministic part
is solved as usual and the stochastic part is added in a second step
in this method, which uses the Euler-Maruyama solver formalism.
For the jump process it must be remembered that all solvers use
discrete time stepping, implying that the number of Poisson jumps
may be > 1 during the finite deltat. Otherwise the jump process
may be prescribed in a number of ways.
NB The present method can also be used for the Milstein scheme, but
all handling related to B and C and their derivatives as well as the
stochastic processes must be placed in "noisemodel".
Methods inherited from dynamics.Dynamics:
- 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 inherited from dynamics.Dynamics:
- __dict__
- dictionary for instance variables (if defined)
- __weakref__
- list of weak references to the object (if defined)
| |