odeAdamsKrogh¶
Solves an initial-value problem for a system of ordinary differential equations of order one or two using a variable order Adams method.
Synopsis¶
odeAdamsKrogh (t, tend, ido, y, hidrvs, fcn)
Required Arguments¶
- float
t
(Input/Output) - On input,
t
contains the initial independent variable value. On output,t
is replaced bytend
unless error conditions arise. Seeido
for details. - float
tend
(Input) - Value of t = tend where the solution is required.
- int
ido
(Input/Output) Flag indicating the state of the computation.
ido
State 1 Initial entry input value. 2 Normal re-entry input value. On output, if ido
= 2 then the integration is finished. If the integrator is called with a new value fortend
, the integration continues. If the integrator is called withtend
unchanged, an error message is issued.3 Input value to use on final call to release workspace. >3 Output value that indicates that a fatal error has occurred. The initial call is made with
ido
= 1. The function then setsido
= 2, and this value is used for all but the last call that is made withido
= 3. This final call is only used to release workspace which was automatically allocated by the initial call withido
= 1.- float
y[]
(Input/Output) - An array of length k containing the dependent variables, y(t), and
first derivatives, if any. k will be the sum of the orders of the
equations in the system of equations to solve, that is, the sum of the
elements of
eqOrder
. On input,y
contains the initial values, \(y(t_0)\) and \(y'(t_0)\) (if needed). On output,y
contains the approximate solution, \(y(t)\). For example, for a system of first order equations,y
[i‑1] is the i‑th dependent variable. For a system of second order equations,y
[2i‑2] is the i‑th dependent variable andy
[2i‑1] is the derivative of the i‑th dependent variable. For systems of equations in which one or more equations is of order 2, optional argumenteqOrder
must be used to denote the order of each equation so that the derivatives iny
can be identified. By default it is assumed that all equations are of order 1 andy
contains only dependent variables. - float
hidrvs[]
(Output) - An array of length
neq
containing the highest order derivatives at the pointy
. - void
fcn
(neq
,ido
,t
,y[]
,hidrvs
) (Input) - User-supplied function to evaluate derivatives.
Arguments
- int
neq
(Input) - Number of differential equations in the system of equations to solve.
- int
ido
(Input) - Flag indicating the state of the computation. This flag corresponds to
the
ido
argument described above. Iffcn
has complicated subexpressions, which depend only weakly or not at all ony
then these subexpressions need only be computed whenido
= 1 and their values then reused whenido
=
2. - float
t
(Input) - Independent variable, t.
- float
y[]
(Input) - An array of length k containing the dependent variable values, y, and first derivatives, if any. k will be the sum of the orders of the equations in the system of equations to solve.
- float
hidrvs[]
(Output) - An array of length
neq
containing the values of the highest order derivatives evaluated at (t, y).
Optional Arguments¶
eqOrder
, int (Input)An array of length
neq
specifying the orders of the equations in the system of equations to solve. The elements ofeqOrder
can be 1 or 2.eqOrder
must be used with argumenty
to define systems of mixed or higher order.Default:
eqOrder
= [1,1,1,…,1].eqErr
, float (Input)An array of length
neq
specifying the error tolerance for each equation. Let e(i) be the error tolerance for equation i for i = 0,…,neq
‑1. ThenValue Explanation \(e(i) > 0\) Implies an absolute error tolerance of \(e(i)\) is to be used for equation i. \(e(i) = 0\) Implies no error checking is to be performed for equation i. \(e(i) < 0\) Implies a relative error test is to be performed for equation i. In this case, the base error tolerance used will be \(|e(i)|\) and the relative error factor used will be \((15/16 * |e(i)|)\). Thus the actual absolute error tolerance used will be \(|e(i)| \times (15/16 \times |e(i)|)\). Default: An absolute error tolerance of 1.e‑5 is used for single precision and 1.e‑10 for double precision for all equations.
stepsizeInc
, float (Input)Factor used for increasing the stepsize. One should set
stepsizeInc
such that 9/8 <=stepsizeInc
<= 4.Default:
stepsizeInc
= 2.0.stepsizeDec
, float (Input)Factor used for decreasing the stepsize. One should set
stepsizeDec
such that 1/4 <=stepsizeDec
<= 7/8.Default:
stepsizeDec
= 0.5.minStepsize
, float (Input)Absolute value of the minimum stepsize permitted.
Default:
minStepsize
= 10.0/machine
(2).maxStepsize
, float (Input)Absolute value of the maximum stepsize permitted.
Default:
maxStepsize
=machine
(2).
Description¶
odeAdamsKrogh
is based on the JPL Library routine SIVA
.
odeAdamsKrogh
uses a variable order Adams method to solve the initial
value problem
or more generally
where y is the vector
\(z_i^{(k)}\) is the kth derivative of \(z_i\) with respect to t, \(d_i\) is the order of the ith differential equation, and η is a vector with the same dimension as y.
Note that the systems of equations solved by odeAdamsKrogh
can be of
order one, order two, or mixed order one and two.
See “Changing Stepsize in the Integration of Differential Equations Using Modified Divided Differences,” Krogh (1974).
Examples¶
Example 1¶
In this example a system of two equations of order two is solved.
The initial conditions are
Since the system is of order two, optional argument eqOrder
must be used
to specify the orders of the equations. Also, because the system is of order
two, y
[0] contains the first dependent variable, y[1]
contains
the derivative of the first dependent variable, y[
2]
contains
the second dependent variable, and y[
3]
contains the derivative
of the second dependent variable.
from __future__ import print_function
from numpy import *
from pyimsl.math.odeAdamsKrogh import odeAdamsKrogh
from pyimsl.math.constant import constant
def fcn(neq, ido, t, y, hidrvs):
tp = y[0] * y[0] + y[2] * y[2]
tp = 1.0e0 / (tp * sqrt(tp))
hidrvs[0] = -y[0] * tp
hidrvs[1] = -y[2] * tp
neq = 2
ido = [1]
t = [0.0]
y = [1.0, 0, 0, 1.0]
hidrvs = []
korder = [2, 2]
# Write Title
print(" T Y1/Y2 Y1P/Y2P Y1PP/Y2PP")
# Integrate ODE
iend = 0
delta = 2.0 * constant("PI")
for k in range(5):
iend += 1
tend = t[0] + delta
if (tend > 20.0):
tend = 20.0
odeAdamsKrogh(neq, t, tend, ido, y, hidrvs, fcn,
eqOrder=korder)
if (iend < 5):
print("%8.4f %10.4f %15.4f %15.4f" %
(t[0], y[0], y[1], hidrvs[0]))
print(" %12.4f %15.4f %15.4f" %
(y[2], y[3], hidrvs[1]))
if (iend == 4):
ido[0] = 3
Output¶
T Y1/Y2 Y1P/Y2P Y1PP/Y2PP
6.2832 1.0000 -0.0000 -1.0000
0.0000 1.0000 -0.0000
12.5664 1.0000 -0.0000 -1.0000
0.0000 1.0000 -0.0000
18.8496 1.0000 -0.0000 -1.0000
0.0000 1.0000 -0.0000
20.0000 0.4081 -0.9129 -0.4081
0.9129 0.4081 -0.9129
Example 2¶
This contrived example illustrates how to use odeAdamsKrogh
to solve a
system of equations of mixed order.
The height, y(t), of an object of mass m above the surface of the Earth can be modeled using Newton’s second law as:
or
where ‑mg is the downward force of gravity and ‑ky’ is the force due to air resistance, in a direction opposing the velocity. If the object is a meteor, the mass, m, and air resistance, k, will decrease as the meteor burns up in the atmosphere. The mass is proportional to \(r^3\) (r = radius) and the air resistance, presumably dependent on the surface area, may be assumed to be proportional to \(r^2\), so that \(k/m = k_0/r\). The rate at which the meteor’s radius decreases as it burns up may depend on r, on the velocity y’, and, since the density of the atmosphere depends on y, on y itself. However, we will construct a very simple model where the rate is just proportional to the square of the velocity,
We solve (1) and (2), with \(k_0 = 0.005\), \(c_0 = 10^{-8}\), \(g = 9.8\) and initial conditions \(y(0) = 100,000 \text{ meters}\), \(y'(0) = ‑1000 \text{ meters/second}\), \(r(0) = 1 \text{ meter}\).
from __future__ import print_function
from numpy import *
from pyimsl.math.odeAdamsKrogh import odeAdamsKrogh
def fcn(neq, ido, t, y, hidrvs):
hidrvs[0] = -9.8 - .005 / y[2] * y[1]
hidrvs[1] = -1.e-8 * y[1] * y[1]
neq = 2
ido = [1]
t = [0.0]
y = [100000.0, -1000.0, 1]
korder = [2, 1]
eqnerr = [.003, .003]
hidrvs = []
# Write Title
print(" T Y1/Y2 Y1P/Y2P Y1PP/Y2PP")
# Integrate ODE
iend = 0
delta = 10.0
for k in range(6):
iend += 1
tend = t[0] + delta
if (tend > 50.0):
tend = 50.0
odeAdamsKrogh(neq, t, tend, ido, y, hidrvs, fcn,
eqOrder=korder,
eqErr=eqnerr)
if (iend < 6):
print("%8.4f %10.4f %12.4f %15.4f" %
(t[0], y[0], y[1], hidrvs[0]))
print(" %15.4f %15.4f" %
(y[2], hidrvs[1]))
if (iend == 5):
ido[0] = 3
Output¶
T Y1/Y2 Y1P/Y2P Y1PP/Y2PP
10.0000 89773.0459 -1044.0096 -3.9701
0.8954 -0.0109
20.0000 79150.9980 -1078.6333 -2.9083
0.7826 -0.0116
30.0000 68240.9574 -1101.0377 -1.5031
0.6635 -0.0121
40.0000 57184.9223 -1106.9633 0.4253
0.5413 -0.0123
50.0000 46178.1509 -1089.8294 3.1700
0.4201 -0.0119
Warning Errors¶
IMSL_
TOLERANCE_TOO_SMALL |
The requested error tolerance, # is too small. Using # instead. |
IMSL_ RESTART |
The stepsize has been reduced too rapidly The integrator is going to do a restart. |
Fatal Errors¶
IMSL_ADJUST_STEPSIZE1 |
The current step length = #, is
less than the minimum steplength,
“minStepsize ” = #, at the
conclusion of the starting phase of
the integration. Decreasing
“minStepsize ” to a value
less than or equal to # may help. |
IMSL_ADJUST_STEPSIZE2 |
The integrator needs to take a step
smaller than # in order to maintain
the requested local error.
Decreasing “minStepsize ” to
a value less than or equal to # may
help. |
IMSL_INCORRECT_TEND |
Either the new output point
precedes the last one or it has the
same value. “tend ” = #. |
IMSL_ADJUST_ERROR_TOLERANCE |
The step length, H = #, is so small
that when Tn + H is formed, the
result will be the same as Tn ,
where Tn is the base value of
the independent variable. If this
problem is not due to a
nonintegrable singularity, it can
probably be corrected by
translating “t ” so that it
is closer to 0. Reducing the error
tolerance for the equations through
argument “eqErr ” may also
help with this problem. |
IMSL_ERROR_TOLERANCE |
A local error tolerance of zero has been requested. |
IMSL_ERROR_PREVIOUS |
A fatal error has occurred because of the error reported in the previous error message. |
IMSL_STOP_USER_FCN |
Request from user supplied function to stop algorithm. User flag = “#”. |