differentialAlgebraicEqs¶

Solves a first order differential-algebraic system of equations, \(g(t, y, y') = 0\), with optional additional constraints and user-defined linear system solver.

Synopsis¶

differentialAlgebraicEqs (t, tend, ido, y, yprime, gcn)

Required Arguments¶

float t (Input/Output)

Set t to the starting value \(t_0\) at the first step. On output, t is set to the value to which the integration has advanced. Normally, this new value is tend.

float tend (Input)

Final value of the independent variable. Update this value when re-entering after output with ido = 2.

int ido (Input/Output)

Flag indicating the state of the computation.

`ido`	State
1	Initial entry
2	Normal re-entry after obtaining output
3	Release workspace, last call

The user sets ido = 1 on the first call at t = \(t_0\). The function then sets ido =2, and this value is used for all but the last entry, which is made with ido = 3.

float y[] (Input/Output)

Array of length neq containing the dependent variable values, y. On input, y must contain initial values. On output, y contains the computed solution at tend.

float yprime[] (Input/Output)

Array of length neq containing derivative values, yʹ. This array must contain initial values, but they need not be such that \(g(t, y, y') = 0\) at t = \(t_0\). See the description of parameter yprimeMethod for more information.

void gcn (neq, t, y[], yprime[], delta[], d[[]], ldd, ires) (Input)

User-supplied function to evaluate \(g(t, y, y')\), and any constraints. Also partial derivative evaluations and optionally linear solving steps occur here. The equations \(g(t, y, y') = 0\) consist of neq differential-algebraic equations of the form.

\[F_i\left(t,y_1, \ldots, y_{\mathit{neq}}, y'_1, \ldots, y'_{\mathit{neq}}\right) \equiv F_i\left(t,y,y'\right) = 0, \phantom{...} i = 1, \ldots, \mathit{neq}\]

The function gcn is also used to evaluate the nConstraints additional algebraic constraints.

\[C_i\left(t,y_1, \ldots, y_{\mathit{neq}}\right) \equiv C_i(t,y) = 0, \phantom{...} i=1, \ldots, \mathit{ncon}, \phantom{...} \mathit{ncon} \geq 0\]

Arguments

int neq (Input)

Number of dependent variables, and number of differential‑algebraic equations, not counting any additional constraints.

float t (Input)

Integration variable t.

float y[] (Input)

Array of neq dependent variables, y.

float yprime[] (Input)

Array of neq derivative values, yʹ.

float delta[] (Input/Output)

Array of length max(neq, nConstraints) containing residuals, \(\delta\) . See parameter ires for a definition.

float d[[]] (Input/Output

Array of length ldd × neq containing partial derivatives, d See parameter ires for a definition.

int ldd (Input)

Number of rows in d.

int ires (Input/Output)

Flag indicating what is to be calculated in the user function, gcn.

Note: ires is input only, except when ires = 6. It is input/output when ires = 6. For a detailed description see the table below.

The code calls gcn with ires = 0, 1, 2, 3, 4, 5, 6, or 7, defined as follows:

`ires`	Description
0	Do initializations which may be required in later calls to `gcn`. This is a setup flag that is input to gcn just once per problem. Initializations might be computing parameters to be used internally by `gcn` or taking any other necessary steps for what may follow in terms of evaluating derivatives or linear solves. Return and do nothing if no initializations are needed.
1	Compute \(\delta_i = F_i \left( t, y, y' \right)\) , the i‑th residual, for i =1,…, `neq`.
2	(Required only if `jacobianOption` =1 and `jacobianMatrixType` = 0 or 1). Compute \(d_{i,j} = \frac{\partial F_i \left( t, y, y' \right)}{\partial y_j}\) , the partial derivative matrix. These are derivatives of \(F_i\) with respect to \(y_j\) , for i =1,…, `neq` and j = 1,…, `neq`.
3	(Required only if `jacobianOption` =1 and `jacobianMatrixType` = 0 or 1). Compute \(d_{i,j} = \frac{\partial F_i \left( t, y, y' \right)}{\partial y'_j}\) , the partial derivative of \(F_i\) with respect to \(y'_j\) , for i = 1,…, `neq` and j = 1,…, `neq`.
4	(Required only if `yprimeMethod` = 2). Compute \(\delta_i = \frac{\partial F_i \left( t, y, y' \right)}{\partial t}\) , the partial derivative of \(F_i\) with respect to t, for i =1,…, `neq`.
5	(Required only if `nConstraints` > 0). Compute \(\delta_i = C_i(t, y)\) , the i-th residual in the additional constraints, for i =1,…, `nConstraints`, and \(d_{i,j} = \frac{\partial C_i (t,y)}{\partial y_j}\) , the partial derivative of \(C_i\) with respect to \(y_j\) for i =1,…, `nConstraints` and j =1,…, `neq`.
6	(Required only if `method` = 1.) If `jacobianMatrixType` = 2, the user must compute the matrix \(A = \frac{\partial F}{\partial y} + cj \frac{\partial F}{\partial y'}\) , where \(cj = \delta_1\) , and save this matrix in any user-defined format. This is for later use when `ires` = 7. The matrix may also be factored in this step, if desired. The array `d` is not referenced if `jacobianMatrixType` = 2. If `jacobianMatrixType` = 0 or 1, the A matrix will already be defined and passed to `gcn` in the array `d`, which will be in full matrix format if `jacobianMatrixType` = 0, and band matrix format, if `jacobianMatrixType` = 1. The user may factor `d` in this step, if desired. Note: For `jacobianMatrixType` = 0, 1, or 2, the user must set `ires` = 0 to signal that A is nonsingular. If A is nearly singular, leave `ires` = 6. This results in using a smaller step-size internally.
7	(Required only if `method` = 1.) The user must solve \(Ax=b\) , where \(b\) is passed to `gcn` in the vector `delta`, and x is returned in `delta`. If `jacobianMatrixType` = 2, A is the matrix which was computed and saved at the call with `ires` = 6; if `jacobianMatrixType` = 0 or 1, A is passed to `gcn` in the array `d`. In either case, the A matrix will remain factored if the user factored it when `ires` = 6.

Optional Arguments¶

nConstraints (Input)

Number of additional constraints.

Default: nConstraints = 0.

jacobianOption, int (Input)

Jacobian calculation option.

`jacobianOption`	Description
0	Calculates using finite difference approximations.
1	User supplies the Jacobian matrices of partial derivatives of \(F_i, i = 1, \ldots, \mathit{neq}\), in the function `gcn`, when `ires` = 2 and 3.

Default: jacobianOption = 0 for jacobianMatrixType = 0 or 1. jacobianOption = 1 for jacobianMatrixType = 2.

yprimeMethod, int (Input)

Initial yʹ calculation method.

`yprimeMethod`	Description
0	The initial input values of `yprime` are already consistent with the input values of `Y`. That is \(g(t, y, y') = 0\) at \(t = t_0\). Any constraints must be satisfied at \(t = t_0\).
1	Consistent values of `yprime` are calculated by Petzold’s original `DASSL` algorithm.
2	Consistent values of `yprime` are calculated using a new algorithm [Hanson and Krogh, 2008], which is generally more robust but requires that `jacobianOption` = 1 and `method` = 0, and additional derivatives corresponding to `ires` = 4 are to be calculated in `gcn`.

Default: yprimeMethod = 1.

jacobianMatrixType, int (Input)

Parameter specifying the Jacobian matrix structure.

jacobianMatrixType Description

0 The Jacobian matrices (whether jacobianOption = 0 or 1) are to be stored in full storage mode.

1

The Jacobian matrices are to be stored in band storage mode. In this case, if jacobianOption = 1, the partial derivative matrices have their entries for row i and column j, stored as array elements

\[d_{(i-j+mu+1,j)}\]

.

This occurs when ires = 2 or 3 in gcn.

2 A user-defined matrix structure is used (see the documentation for 6 for more details). If jacobianMatrixType = 2, method and jacobianOption are set to 1 internally.

Default: jacobianMatrixType = 0.

method, int (Input)

Solve method.

`method`	Description
0	`differentialAlgebraicEqs` solves the linear systems.
1	The user wishes to solve the linear system in function `gcn`. See parameter `gcn` for details.

Default: method = 0 for jacobianMatrixType = 0 or 1. method = 1 for jacobianMatrixType = 2.

nLowerDiag, int (Input)

Number of non-zero diagonals below the main diagonal in the Jacobian matrices when band storage mode is used. nLowerDiag is ignored if jacobianMatrixType ≠ 1.

Default: nLowerDiag = neq ‑ 1.

nUpperDiag, int (Input)

Number of non-zero diagonals above the main diagonal in the Jacobian matrices when band storage mode is used. nUpperDiag is ignored if jacobianMatrixType ≠ 1.

Default: nUpperDiag = neq ‑ 1.

relativeTolerance float (Input)

Relative error tolerance for solver. The program attempts to maintain a local error in Y(i) less than relativeTolerance*∣y[i]∣ + absoluteTolerance[i].

Default: relativeTolerance = \(\sqrt{\varepsilon}\) , where ɛ is machine precision.

absoluteTolerance float (Input)

Array of size neq containing absolute error tolerances. See description of relativeTolerance.

Default: absoluteTolerance[i] = 0.0.

initialStepsize, float (Input)s

Initial stepsize used by the solver. If initialStepsize = 0.0, the function defines the initial stepsize.

Default: initialStepsize = 0.0.

maxStepsize, float (Input)

Maximum stepsize used by the solver. If maxStepsize = 0.0, the function defines the maximum stepsize.

Default: maxStepsize = 0.0.

maxOrder, int (Input)

Maximum order of the backward difference formulas used. 1 ≤ maxOrder ≤ 5.

Default: maxOrder = 5.

maxNumberSteps, int (Input)

Maximum number of steps allowed from t to tend.

Default: maxNumberSteps = 500.

integrationLimit, float (Input)

Integration limit point. For efficiency reasons, the code sometimes integrates past tend and interpolates a solution at tend. If a value for integrationLimit is specified, the code will never integrate past t=integrationLimit.

Default: No integrationLimit value is specified.

orderMagnitudeEst, float (Input)

Order-of-magnitude estimate. orderMagnitudeEst is used as an order-of-magnitude estimate of the magnitude of the functions \(F_i\) (see description of gcn), for convergence testing, if yprimeMethod=2. orderMagnitudeEst is ignored if yprimeMethod=0 or 1.

Default: orderMagnitudeEst = 1.0.

Description¶

Function differentialAlgebraicEqs finds an approximation to the solution of a system of differential-algebraic equations \(g \left( t, y, y' \right) = 0\) with given initial data for \(y\) and \(y'\) . The function uses BDF formulas, which are appropriate for stiff systems. differentialAlgebraicEqs is based on the code DASSL designed by Linda Petzold [1982], and has been modified by Hanson and Krogh [2008]Solving Constrained Differential-Algebraic Systems Using Projections to allow the inclusion of additional constraints, including conservation principles, after each time step. The modified code also provides a more robust algorithm to calculate initial \(y'\) values consistent with the given initial \(y\) values. This occurs when the initial \(y'\) are not known.

A differential-algebraic system of equations is said to have “index 0” if the Jacobian matrix of partial derivatives of the \(F_i\) with respect to the \(y'_j\) is nonsingular. Thus it is possible to solve for all the initial values of \(y'_j\) and put the system in the form of a standard ODE system. If it is possible to reduce the system to a system of index 0 by taking first derivatives of some of the equations, the system has index 1, otherwise the index is greater than 1. See Brenan [1989] for a definition of index. differentialAlgebraicEqs can generally only solve systems of index 0 or 1; other systems will usually have to be reduced to such a form through differentiation.

Examples¶

Example 1 – Method of Lines PDE Problem¶

This example solves the partial differential equation \(U_t = U_{xx} + U\) , with initial condition \(U(x, 0) = 1 + x\) , and boundary conditions \(U(0, t) = e^t\), \(U(1, t) = 2e^t\) which has exact solution \(U(x, t) = (1+x)e^t\). If we approximate the \(U_{xx}\) term using finite differences, where \(x_i = (i-1)h\) , and \(h = 1 / (n - 1)\), we get:

\[U\left(x_1,t\right) = e^t\]

\[U'(x_i,t) = \left[ U\left(x_{i+1}, t\right) - 2U\left(x_i, t\right) + U\left(x_{i-1},t\right) \right] / h^2 + U(x_i,t), i = 2, \ldots, n - 1\]

\[U\left(x_n,t\right) = 2e^t\]

If \(Y_i(t) = U(x_i,t)\), the first and last equations are algebraic and the others are differential equations, so this is a system of differential-algebraic equations. The system has \(index = 1\), since it could be transformed into an ODE system by differentiating the first and last equations. Note that the Jacobian matrices are banded (tridiagonal), with nLowerDiag = nUpperDiag = 1. We use this and specify the option for dealing with banded matrices in differentialAlgebraicEqs.

from __future__ import print_function
from numpy import *
from pyimsl.math.differentialAlgebraicEqs import differentialAlgebraicEqs
from pyimsl.math.writeMatrix import writeMatrix


def gcn(neq, t, y, yprime, delta, d, ldd, ires):
    hx = 1.0 / (neq - 1)
    mu = 1
    if ires[0] == 1:
        # f_i defined here
        delta[0] = y[0] - exp(t)
        for i in range(1, neq - 1):
            delta[i] = -yprime[i] + \
                (y[i + 1] - 2.0 * y[i] + y[i - 1]) / pow(hx, 2) + y[i]
        delta[neq - 1] = y[neq - 1] - 2.0 * exp(t)
    elif ires[0] == 2:
        # d(i-j+mu+1,j) = d(f_i)/d(y_j) in band storage mode
        d[(mu) * neq + 0] = 1.0
        for i in range(1, neq - 1):
            j = i
            d[(i - j + mu + 1) * neq + j - 1] = 1.0 / pow(hx, 2)
            j = i + 1
            d[(i - j + mu + 1) * neq + j - 1] = -2.0 / pow(hx, 2) + 1.0
            j = i + 2
            d[(i - j + mu + 1) * neq + j - 1] = 1.0 / pow(hx, 2)
        d[(mu) * neq + neq - 1] = 1.0
    elif ires[0] == 3:
        #  d(i-j+mu+1,j) = d(f_i)/d(yprime_j )
        for i in range(1, neq - 1):
            d[(mu) * neq + i] = -1.0


neq = 101
iujac = 1
iypr = 1
matstr = 1
ml = 1
mu = 1
nsteps = 10
errmax = 0.0
rtol = 1.0e-4

hx = 1.0 / (neq - 1)
y = zeros(neq, dtype=double)
yprime = zeros(neq, dtype=double)
# Initial values
for i in range(neq):
    yprime[i] = 0.0
    x = i * hx
    y[i] = 1.0 + x

# Always set ido=1 on first call
ido = [1]
for i in range(nsteps):
    # Output solution at t=0.1,0.2,...,1.0
    t = [0.1 * i]
    tend = 0.1 * (i + 1)
    # Set ido = 3 on last call
    if i == (nsteps - 1):
        ido = [3]
    # User-supplied jacobian matrix (iujac=1)
    # Banded jacobian (matstr=1) */
    differentialAlgebraicEqs(t, tend, ido, y,
                             yprime, gcn,
                             jacobianOption=iujac,
                             yprimeMethod=iypr,
                             jacobianMatrixType=matstr,
                             nLowerDiag=ml,
                             nUpperDiag=mu,
                             relativeTolerance=rtol)

for i in range(neq):
    x = i * hx
    tr = (1.0 + x) * exp(t)
    errmax = max([errmax, abs(y[i] - tr)])

print("Max Error at T=1 is %g." % errmax)

Output¶

Max Error at T=1 is 5.6854e-05.

Example 2 – Pendulum Problem¶

The first-order equations of motion of a point-mass m suspended on a massless wire of length \(L\) under the influence of gravity, mg, and wire tension, λ , in Cartesian coordinates \((p,q)\) are

\[\begin{split}\begin{array}{l} p' = u \\ q' = v \\ mu' = -p\lambda \\ mv' = -q\lambda - mg \\ p^2 + q^2 - L^2 = 0 \\ \end{array}\end{split}\]

The problem above has an index number equal to 3, thus it cannot be solved with differentialAlgebraicEqs directly. Unfortunately, the fact that the index is greater than 1 is not obvious, but an attempt to solve it will generally produce an error message stating the corrector equation did not converge, or if yprimeMethod = 2 an error message stating that the index appears to be greater than 1 should be issued. The user then differentiates the last equation, which after replacing pʹ by u and qʹ by v, gives \(pu + qv = 0\). This system still has \(index=2\) (again not obvious, the user discovers this by unsuccessfully trying to solve the new system) and the last equation must be differentiated again, to finally (after appropriate substitutions) give the equation of total energy balance:

\[m\left(u^2+v^2\right) - mgq - L^2\lambda = 0\]

With initial conditions and appropriate definitions of the dependent variables, the system becomes:

\[\begin{split}\begin{array}{c} p(0) = L,q(0) = u(0) = v(0) = \lambda(0) = 0 \hfill \\ y_1 = p \hfill \\ y_2 = q \hfill \\ y_3 = u \hfill \\ y_4 = v \hfill \\ y_5 = \lambda \hfill \\ F_1 = y_3 - y'_1 = 0 \hfill \\ F_2 = y_4 - y'_2 = 0 \hfill \\ F_3 = -y_1y_5 - my'_3 = 0 \hfill \\ F_4 = -y_2y_5 - mg - my'_4 = 0 \hfill \\ F_5 = m\left(y_3^2 + y_4^2\right)-mgy_2 - L^2y_5 = 0 \hfill \\ \end{array}\end{split}\]

The initial conditions correspond to the pendulum starting in a horizontal position.

Since we have replaced the original constraint, \(C_1 = p^2 + q^2 - L^2 = 0\) , which requires that the pendulum length be L, by differentiating it twice, this constraint is no longer explicitly enforced, and if we try to solve the above system alone (ie, with nConstraints=0), the pendulum length drifts substantially from L at larger times. differentialAlgebraicEqs therefore allows the user to add additional constraints, to be re-enforced after each time step, so we add this original constraint, as well as the intermediate constraint \(C_2 = pu + qv = 0\) . Using these two supplementary constraints, (nConstraints = 2), the pendulum length is constant.

from __future__ import print_function
from numpy import *
from pyimsl.math.differentialAlgebraicEqs import differentialAlgebraicEqs
from pyimsl.math.writeMatrix import writeMatrix
from pyimsl.util.imslConvert import toNdarray, toCtypes


def gcn(neq, t, y, yprime, delta, d, ldd, ires):

    mass = 1.0
    length = 1.1
    gravity = 9.806650

    # Convert d to an ndarray to allow easier subscripting.
    D = toNdarray(d, (ldd, neq))

    # Simple swinging pendulum problem
    mg = mass * gravity
    lsq = pow(length, 2)
    if ires[0] == 1:
        # f_i defined here
        delta[0] = y[2] - yprime[0]
        delta[1] = y[3] - yprime[1]
        delta[2] = -y[0] * y[4] - mass * yprime[2]
        delta[3] = -y[1] * y[4] - mass * yprime[3] - mg
        delta[4] = mass * (pow(y[2], 2) + pow(y[3], 2)) - \
            mg * y[1] - lsq * y[4]
    elif ires[0] == 2:
        # d(i,j) = d(f_i)/d(y_j)
        D[0][2] = 1.0
        D[1][3] = 1.0
        D[2][0] = -y[4]
        D[2][4] = -y[0]
        D[3][1] = -y[4]
        D[3][4] = -y[1]
        D[4][1] = -mg
        D[4][2] = mass * 2.0 * y[2]
        D[4][3] = mass * 2.0 * y[3]
        D[4][4] = -lsq
    elif ires[0] == 3:
        # d(i,j) = d(f_i)/d(yprime_j)
        D[0][0] = -1.0
        D[1][1] = -1.0
        D[2][2] = -mass
        D[3][3] = -mass
    elif ires[0] == 4:
        # delta(i) = d(f_i)/dt
        delta[0] = 0.0
        delta[1] = 0.0
        delta[2] = 0.0
        delta[3] = 0.0
        delta[4] = 0.0
    elif ires[0] == 5:
        # delta(i) = g_i
        # d(i,j) = d(g_i)/d(y_j)
        delta[0] = pow(y[0], 2) + pow(y[1], 2) - lsq
        delta[1] = y[0] * y[2] + y[1] * y[3]
        D[0][0] = 2.0 * y[0]
        D[0][1] = 2.0 * y[1]
        D[0][2] = 0.0
        D[0][3] = 0.0
        D[0][4] = 0.0
        D[1][0] = y[2]
        D[1][1] = y[3]
        D[1][2] = y[0]
        D[1][3] = y[1]
        D[1][4] = 0.0

    # Copy results from ndarray back into ctype.
    toCtypes(D, d)


neq = 5
ncon = 2
nsteps = 5
iypr = 2
iujac = 1
maxsteps = 50000
mass = 1.0
length = 1.1
gravity = 9.806650

# INitial values
tol = 1.0e-5
y = zeros(neq, dtype=double)
yprime = zeros(neq, dtype=double)
atol = zeros(neq, dtype=double)
atol[:] = tol
y[0] = length
print("        T               Y(0)            Y(1)            Length")

# Always set ido=1 on first call
ido = [1]
for i in range(nsteps):
    # Output solution at t=10,20,30,40,50
    t = [10.0 * i]
    tend = 10.0 * (i + 1)
    # Set ido = 3 on last call
    if i == (nsteps - 1):
        ido = [3]
    # User-supplied jacobian matrix (iujac=1)
    # Use new algorithm to get compatible y'
    differentialAlgebraicEqs(t, tend, ido, y,
                             yprime, gcn,
                             nConstraints=ncon,
                             jacobianOption=iujac,
                             yprimeMethod=iypr,
                             relativeTolerance=tol,
                             absoluteTolerance=atol,
                             maxNumberSteps=maxsteps)

    # len = pendulum length (should be constant)
    len = sqrt(pow(y[0], 2) + pow(y[1], 2))
    print("%15.7f %15.7f %15.7f %15.7f" % (t[0], y[0], y[1], len))

Output¶

        T               Y(0)            Y(1)            Length
0000000       1.0998160      -0.0201200       1.1000000
0000000       1.0970757      -0.0801564       1.1000000
0000000       1.0854514      -0.1783126       1.1000001
0000000       1.0551496      -0.3109016       1.1000002
0000000       0.9942180      -0.4706711       1.1000003

Example 3 – User Solves Linear System¶

Consider the system of ordinary differential equations, \(y' = By\), where B is the bi-diagonal matrix with \((‑1, ‑1/2, ‑1/3, ..., ‑1/(n‑1), 0)\) on the main diagonal and with 1’s along the first sub-diagonal. The initial condition is \(y(0) = (1,0,0,...,0)^T\), and since \(y'(0) = By(0) = (-1,1,0,...,0)^T\), \(y'(0)\) is also known for this problem.

Since \(B^Tv = 0\), where \(v_i = 1/(i-1)!\), v is an eigenvector of \(B^T\)corresponding to the eigenvalue 0. Thus

\[0 = v^T\left(y'-By\right) = v^Ty' - \left(B^Tv\right)^Ty = v^Ty' = \left(v^Ty\right)'\]

so \(v^T y(t)\) is constant. Since it has the value \(v^T y(0) = v_1 = 1\) at \(t = 0\), the constraint \(v^T y(t) = 1\) is satisfied for all t. This constraint is imposed in this example.

This example also illustrates how the user can solve his/her own linear systems (jacobianMatrixType=2). Normally, when ires = 6, the matrix

\[A = \frac{\partial g}{\partial y} + cj\frac{\partial g}{\partial y'}\]

is computed, saved and possibly factored, using a sparse matrix factorization function of the user’s choice. Then when ires=7, the system Ax = delta is solved, using the matrix B saved and factored earlier, and the solution is returned in delta. In this case, B is just a bidiagonal matrix, so there is no need to save or factor A when ires = 6, since a bi-diagonal system can be solved directly using forward substitution, when ires = 7.

from __future__ import print_function
from numpy import *
from pyimsl.math.differentialAlgebraicEqs import differentialAlgebraicEqs
from pyimsl.math.writeMatrix import writeMatrix
from pyimsl.util.imslConvert import toNdarray, toCtypes

cj = 0


def gcn(neq, t, y, yprime, delta, d, ldd, ires):
    global cj

    v = zeros(neq, dtype=double)
    v[0] = 1.0
    for i in range(1, neq):
        v[i] = v[i - 1] / i

    # Convert d to an ndarray to allow easier subscripting.
    D = toNdarray(d, (ldd, neq))

    if ires[0] == 1:
        # f_i defined here
        delta[0] = yprime[0] + y[0]
        for i in range(1, neq - 1):
            delta[i] = yprime[i] - y[i - 1] + y[i] / (i + 1)
        delta[neq - 1] = yprime[neq - 1] - y[neq - 2]
    elif ires[0] == 5:
        # constraint is v dot y = 1
        con = -1.0
        for i in range(neq):
            con += v[i] * y[i]
            D[0, i] = v[i]
        delta[0] = con
    elif ires[0] == 6:
        # normally, compute matrix
        # a = df/dy + cj*df/dy' = -b + cj*i
        # here. only cj needs to be saved in this case, however,
        # since b is bidiagonal, so a*x=delta can be solved (ires=7)
        # without saving or factoring b.
        cj = delta[0]
        # if cj > 0 not close to zero, a is nonsingular,
        # so set ires = 0.
        if cj >= 1.0e-4:
            ires[0] = 0
    elif ires[0] == 7:
        # solve a*x=delta and return x in delta.
        delta[0] /= 1.0 + cj
        for i in range(1, neq - 1):
            delta[i] = (delta[i] + delta[i - 1]) / (1.0 / (i + 1) + cj)
        delta[neq - 1] = (delta[neq - 1] + delta[neq - 2]) / cj

    # Copy results from ndarray back into ctype.
    toCtypes(D, d)


neq = 100
nsteps = 10
ncon = 1
iypr = 0
matstr = 2
con = 0.0

v = zeros(neq, dtype=double)
# a^t eigenvector v
v[0] = 1.0
for i in range(1, neq):
    v[i] = v[i - 1] / i
# initial values
y = zeros(neq, dtype=double)
yprime = zeros(neq, dtype=double)
y[0] = 1.0
yprime[0] = -1.0
yprime[1] = 1.0

atol = zeros(neq, dtype=double)
atol[:] = 1.0e-4

# Always set ido=1 on first call
ido = [1]

for i in range(nsteps):
    t = [i]
    tend = (i + 1)
    # Set ido = 3 on last call
    if i == (nsteps - 1):
        ido = [3]
    differentialAlgebraicEqs(t, tend, ido, y,
                             yprime, gcn,
                             nConstraints=ncon,
                             yprimeMethod=iypr,
                             jacobianMatrixType=matstr,
                             absoluteTolerance=atol)
# check if solution satisfies constraint
for i in range(neq):
    con += v[i] * y[i]
print(" V dot Y = %f" % con)

Output¶

 V dot Y = 1.000000

Fatal Errors¶

`IMSL_SYSTEM_CONVERGENCE`	The system has index # but convergence of “`yprime`” values was not achieved.
`IMSL_SYSTEM_CONVERGENCE`	The system appears to have `index` > 1.
`IMSL_SYS_INDEX_GT_ONE`	For the Pareto distribution, the Hessian cannot be calculated because the parameter estimate is 0.
`IMSL_SYS_NOT_DIFF_ALG`	This is not a differential-algebraic system.
`IMSL_ACCURACY_EXCEEDED_1`	Accuracy requested exceeds machine precision.
`IMSL_STEPS_EXCEEDED`	More than “`maxNumberSteps`”=# steps taken between output points.
`IMSL_ERROR_TEST_FAILURE_2`	The error test has failed repeatedly.
`IMSL_CORRECTOR_FAILED_3`	The corrector iteration failed repeatedly to converge.
`IMSL_SINGULAR_MATRIX_1`	The iteration matrix is singular.
`IMSL_UNABLE_TO_SOLVE_YPR`	Unable to solve for initial “`yprime`”.
`IMSL_TEND_GT_TSTOP`	“`tend`” is greater than “`integrationLimit`”.
`IMSL_TEND_CLOSE_TO_T`	“`tend`” is too close to “`t`”.
`IMSL_TSTOP_INCONSIST_T`	“`integrationLimit`” is not consistent with“`t`”.
`IMSL_CONSTRAINTS_INCONSIST`	Constraints appear inconsistent
`IMSL_STOP_USER_FCN`	Request from user supplied function to stop algorithm. User flag = “#”.