Solves a first order differential-algebraic system of equations, g(t, y, yʹ) = 0, with optional additional constraints and user-defined linear system solver.
Note: DAESL replaces deprecated routine DASPG.
T Independent
variable, t. (Input/Output)
Set T to the starting
value t0 at the first step. On output, T is set to the value
to which the integration has advanced. Normally, this new value is TEND.
TEND Final
value of the independent variable. (Input)
Update this value when
re-entering after output with IDO = 2.
IDO Flag indicating the state of the computation. (Input/Output)
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 = t0. The routine then sets IDO =2, and this value is used for all but the last entry, which is made with IDO = 3.
Y Array of size
NEQ containing
the dependent variable values, y. (Input/Output)
On
input, Y must
contain initial values. On output, Y contains the
computed solution at TEND.
YPRIME Array of
size NEQ
containing derivative values, yʹ.
(Input/Output)
This array must contain initial values, but they
need not be such that g(t, y, yʹ) = 0
at
t= t0. See the description of parameter
IYPR for more information.
GCN
User-supplied subroutine 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.
The routine GCN is also used to evaluate the NCON additional algebraic constraints
The usage is CALL GCN (T, Y, YPRIME, DELTA, D, LDD, IRES [, ]) where
Required Arguments
T Integration variable t. (Input)
Y Array of NEQ dependent variables, y. (Input)
YPRIME Array of NEQ derivative values, yʹ. (Input)
DELTA Output array of length MAX(NEQ, NCON) containing residuals. See parameter IRES for definition. (Input/Output)
D Output array dimensioned D(LDD,NEQ), containing partial derivatives. See parameter IRES for definition. (Input/Output)
LDD Leading dimension of D. (Input)
IRES Flag
indicating what is to be calculated in the user routine, GCN.
(Input/Output)
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 |
Explanation |
0 |
Do initializations, if any are required. |
1 |
Compute DELTA(i) = |
2 |
(Required only if IUJAC=1 and MATSTR = 0 or 1). |
3 |
(Required only if IUJAC =1 and MATSTR = 0 or 1). |
4 |
(Required only if IYPR=2). |
5 |
(Required only if NCON > 0). |
(Required only if ISOLVE = 1.) where cj = DELTA (1), and save this matrix in any user-defined
format. This is for later use when If MATSTR = 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 Note: For MATSTR = 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 ISOLVE = 1.) The user must solve |
Optional Arguments
FCN_DATA A derived type, s_fcn_data, which may be used to pass additional integer or floating point information to or from the user-supplied subroutine. For a description of this argument see FCN_DATA below. (Input/Output)
GCN must be declared EXTERNAL in the calling program.
NEQ Number of
dependent variables, and number of differential/algebraic equations, not
counting any additional constraints. (Input)
Default: NEQ = size (Y).
NCON Number of
additional constraints. (Input)
Default: NCON = 0.
IUJAC Jacobian
calculation option. (Input)
Value |
Description |
0 |
Calculates using finite difference approximations. |
1 |
User supplies the Jacobian matrices of partial
derivatives of |
Default: IUJAC = 0 for MATSTR = 0 or
1.
IUJAC = 1 for MATSTR = 2.
IYPR Initial
yʹ
calculation method. (Input)
Value |
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 = t0. Any constraints must be satisfied at t = t0. |
1 |
Consistent values of YPRIME are calculated by Petzolds 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 IUJAC= 1 and ISOLVE = 0, and additional derivatives corresponding to IRES= 4 are to be calculated in GCN. |
Default: IYPR = 1.
MATSTR
Parameter specifying the Jacobian matrix structure
(Input)
Set to:
Value |
Description |
0 |
The Jacobian matrices (whether IUJAC = 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 IUJAC= 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 IRES = 6 or 7 for more details). If MATSTR = 2, ISOLVE and IUJAC are set to 1 internally. |
Default: MATSTR = 0.
ISOLVE Solve method. (Input)
Value |
Description |
0 |
DAESL solves the linear systems. |
1 |
The user wishes to solve the linear system in routine GCN. See parameter GCN for details. |
Default: ISOLVE = 0 for MATSTR = 0 or
1.
ISOLVE = 1 for MATSTR = 2.
ML Number of
non-zero diagonals below the main diagonal in the Jacobian matrices when band
storage mode is used. (Input)
ML is ignored if MATSTR ≠
1.
Default: ML = NEQ-1.
MU Number of
non-zero diagonals above the main diagonal in the Jacobian matrices when band
storage mode is used. (Input)
MU is ignored if MATSTR ≠
1.
Default: MU = NEQ-1.
RTOL Relative
error tolerance for solver. (Input)
The program attempts to
maintain a local error in Y(i) less than
RTOL*│Y(i)│ +
ATOL(i).
Default:
RTOL = , where
is machine
precision.
ATOL Array of
size NEQ
containing absolute error tolerances. (Input)
See description of
RTOL.
Default:
ATOL(i) =
0.
H0 Initial
stepsize used by the solver. (Input)
If H0
= 0, the routine
defines the initial stepsize.
Default: H0 = 0.
HMAX Maximum
stepsize used by the solver. (Input)
If HMAX=0, the routine
defines the maximum stepsize.
Default: HMAX =
0.
MAXORD Maximum
order of the backward difference formulas used. (Input).
1 ≤
MAXORD ≤
5.
Default: MAXORD = 5.
MAXSTEPS
Maximum number of steps taken from T to TEND.
(Input).
Default: MAXSTEPS = 500.
TSTOP
Integration limit point. (Input)
For efficiency reasons, the code
sometimes integrates past TEND and interpolates
a solution at TEND. If a value for
TSTOP is
specified, the code will never integrate past T=TSTOP.
Default: No
TSTOP value is
specified.
FMAG
Order-of-magnitude estimate. (Input)
FMAG is used as an
order-of-magnitude estimate of the magnitude of the functions Fi (see
description of GCN), for convergence
testing, if IYPR=2. FMAG is ignored if
IYPR=0 or 1.
Default: FMAG = 1.
FCN_DATA A derived type, s_fcn_data, which may
be used to pass additional information to/from the user-supplied
subroutine. (Input/Output)
The derived type, s_fcn_data, is defined
as:
type s_fcn_data
real(kind(1e0)), pointer, dimension(:) :: rdata
integer, pointer, dimension(:) :: idata
end type
in module mp_types. The double precision counterpart to s_fcn_data is named d_fcn_data. The user must include a use mp_types statement in the calling program to define this derived type.
Note that if this optional argument is present then FCN_DATA must also be defined as an optional argument in the user-supplied subroutine.
Generic: CALL DAESL (T, TEND, IDO, Y, YPRIME, GCN [, ])
Specific: The specific interface names are S_DAESL and D_DAESL.
Routine DAESL
finds an approximation to the solution of a system of differential-algebraic
equations with given initial data
for
and
. The routine uses BDF formulas, which are appropriate for
stiff systems. DAESL
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
values consistent with
the given initial
values. This
occurs when the initial
are not known.
A differential-algebraic system of equations is said to
have index 0 if the Jacobian matrix of partial derivatives of the with respect to the
is nonsingular. Thus it is possible to solve for all
the initial values of
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. DAESL
can generally only solve systems of index 0 or 1; other systems will usually
have to be reduced to such a form through differentiation.
This example solves the partial differential equation, with initial condition
, and boundary conditions
,
which has exact
solution
. If we approximate the
term using finite
differences, where
, and
, we get:
If Yi(t) =
U(xi,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 ML
= MU
= 1. We use this and specify the option for dealing with banded matrices
in DAESL.
The parameterand the number of equations
is passed to the evaluation routine, GCN,
with the optional argument USER_DATA.
Link to example source (daesl_ex1.f90)
USE DAESL_INT
USE MP_TYPES
IMPLICIT NONE
! NEQ = Number of equations
INTEGER, PARAMETER :: NEQ=101
REAL T, Y(NEQ), YPRIME(NEQ), TEND, X, TRUE, HX, ERRMAX
INTEGER NOUT, IDO, I, NSTEPS
REAL, TARGET :: RPARAM(1)
INTEGER, TARGET :: IPARAM(1)
TYPE (S_FCN_DATA) USER_DATA
EXTERNAL GCN
! Pass NEQ, HX to GCN
HX = 1.0 / (NEQ-1)
IPARAM(1) = NEQ
RPARAM(1) = HX
USER_DATA%RDATA=>RPARAM
USER_DATA%IDATA=>IPARAM
! Initial values for y, initial guesses for y'
DO I = 1, NEQ
X = (I-1) * HX
Y(I) = 1 + X
END DO
YPRIME = 0.0
NSTEPS = 10
! Always set IDO=1 on first call
IDO = 1
DO I = 1, NSTEPS
! Output solution at T=0.1,0.2,...,1.0
T = 0.1 * (I-1)
TEND = 0.1 * I
! Set IDO = 3 on last call
IF (I == NSTEPS) IDO = 3
! User-supplied Jacobian matrix (IUJAC=1)
! Banded Jacobian (MATSTR=1)
CALL DAESL (T, TEND, IDO, Y, YPRIME, GCN, IYPR=1, IUJAC=1, &
MATSTR=1, ML=1, MU=1, RTOL=1.0E-4, FCN_DATA=USER_DATA)
END DO
ERRMAX = 0.0
DO I = 1, NEQ
X = (I-1) * HX
TRUE = (1+X) * EXP(T)
ERRMAX = MAX(ERRMAX, ABS(Y(I) - TRUE))
END DO
CALL UMACH(2, NOUT)
WRITE (NOUT, *) ' Max Error at T=1 is ', ERRMAX
END
SUBROUTINE GCN (T, Y, YPRIME, DELTA, D, LDD, IRES, FCN_DATA)
USE MP_TYPES
IMPLICIT NONE
REAL T, Y(*), YPRIME(*), DELTA(*), D(LDD,*), HX
INTEGER IRES, LDD, I, J, NEQ, MU
TYPE (S_FCN_DATA), OPTIONAL, INTENT(INOUT) :: FCN_DATA
NEQ = FCN_DATA%IDATA(1)
HX = FCN_DATA%RDATA(1)
MU = 1
SELECT CASE (IRES)
! F_I defined here
CASE(1)
DELTA(1) = Y(1) - EXP(T)
DO I = 2, NEQ-1
DELTA(I) = -YPRIME(I) + (Y(I+1) - 2.0 * Y(I) + Y(I-1)) &
/ HX**2 + Y(I)
END DO
DELTA(NEQ) = Y(NEQ) - 2.0 * EXP(T)
! D(I-J+MU+1,J) = D(F_I)/D(Y_J)
! in band storage mode
CASE(2)
D(MU+1,1) = 1.0
DO I = 2, NEQ-1
J = I-1
D(I-J+MU+1, J) = 1.0 / HX**2
J = I
D(I-J+MU+1, J) = -2.0 / HX**2 + 1.0
J = I+1
D(I-J+MU+1, J) = 1.0 / HX**2
END DO
D(MU+1, NEQ) = 1.0
! D(I-J+MU+1,J) = D(F_I)/D(YPRIME_J)
CASE(3)
DO I = 2, NEQ-1
D(MU+1, I) = -1.0
END DO
END SELECT
END
Max Error at T=1 is 5.6743621E-5
The first-order equations of motion of a point-mass m
suspended on a massless wire of length under the influence of gravity, mg, and wire
tension, λ , in Cartesian coordinates (p,q) are
The problem above has an index number equal to 3, thus it cannot be solved with DAESL 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 IYPR=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:
With initial conditions and appropriate definitions of the dependent variables, the system becomes:
The initial conditions correspond to the pendulum starting in a horizontal position.
Since we have replaced the original constraint, , 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 NCON=0),
the pendulum length drifts substantially from L at larger times. DAESL
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
. Using these two
supplementary constraints, (NCON=2),
the pendulum length is constant.
Link to example source (daesl_ex2.f90)
USE DAESL_INT
USE MP_TYPES
IMPLICIT NONE
! NEQ = Number of equations
! NCON = Number of extra constraints
INTEGER, PARAMETER :: NEQ=5, NCON = 2
REAL, PARAMETER :: MASS=1.0, LENGTH=1.1, GRAVITY=9.806650
REAL T, Y(NEQ), YPRIME(NEQ), TEND, ATOL(NEQ), TOL, LEN
INTEGER NOUT, IDO, I, NSTEPS
REAL, TARGET :: RPARAM(3)
TYPE (S_FCN_DATA) USER_DATA
EXTERNAL GCN
! Pass Mass, Pendulum length and G as parameters
RPARAM(1) = MASS
RPARAM(2) = LENGTH
RPARAM(3) = GRAVITY
USER_DATA%RDATA=>RPARAM
! Initial values for y, guesses for initial y'
Y = 0.0
Y(1) = LENGTH
YPRIME = 0.0
TOL = 1.0E-5
ATOL = TOL
CALL UMACH(2, NOUT)
WRITE (NOUT, 5)
NSTEPS = 5
! Always set IDO=1 on first call
IDO = 1
DO I = 1, NSTEPS
! Output solution at T=10,20,30,40,50
T = 10.0 * (I-1)
TEND = 10.0 * I
! Set IDO = 3 on last call
IF (I.EQ.NSTEPS) IDO = 3
! User-supplied Jacobian matrix (IUJAC=1)
! Use new algorithm to get compatible y'
CALL DAESL (T, TEND, IDO, Y, YPRIME, GCN, NCON=NCON, RTOL=TOL, &
ATOL=ATOL, IYPR=2, IUJAC=1, MAXSTEPS=50000, &
FCN_DATA=USER_DATA)
! LEN = pendulum length (should be constant)
LEN = SQRT(Y(1)**2 + Y(2)**2)
WRITE (NOUT, 10) T, Y(1), Y(2), LEN
END DO
5 FORMAT (8X,'T',14X,'Y(1)',11X,'Y(2)',11X,'Length',/)
10 FORMAT (4F15.7)
END
SUBROUTINE GCN (T, Y, YPRIME, DELTA, D, LDD, IRES, FCN_DATA)
USE MP_TYPES
IMPLICIT NONE
! Simple swinging pendulum problem
REAL T, Y(*), YPRIME(*), DELTA(*), D(LDD,*), MASS, &
LENGTH, GRAVITY, MG, LSQ
INTEGER IRES, LDD
TYPE (S_FCN_DATA), OPTIONAL, INTENT(INOUT) :: FCN_DATA
MASS = FCN_DATA%RDATA(1)
LENGTH = FCN_DATA%RDATA(2)
GRAVITY = FCN_DATA%RDATA(3)
MG = MASS * GRAVITY
LSQ = LENGTH**2
SELECT CASE (IRES)
! F_I defined here
CASE(1)
DELTA(1) = Y(3) - YPRIME(1)
DELTA(2) = Y(4) - YPRIME(2)
DELTA(3) = -Y(1) * Y(5) - MASS * YPRIME(3)
DELTA(4) = -Y(2) * Y(5) - MASS * YPRIME(4) - MG
DELTA(5) = MASS * (Y(3)**2 + Y(4)**2) - MG * Y(2) - LSQ * Y(5)
! D(I,J) = D(F_I)/D(Y_J)
CASE(2)
D(1, 3) = 1.0
D(2, 4) = 1.0
D(3, 1) = -Y(5)
D(3, 5) = -Y(1)
D(4, 2) = -Y(5)
D(4, 5) = -Y(2)
D(5, 2) = -MG
D(5, 3) = MASS * 2.0 * Y(3)
D(5, 4) = MASS * 2.0 * Y(4)
D(5, 5) = -LSQ
! D(I,J) = D(F_I)/D(YPRIME_J)
CASE(3)
D(1, 1)= -1.0
D(2, 2)= -1.0
D(3, 3)= -MASS
D(4, 4)= -MASS
! DELTA(I) = D(F_I)/DT
CASE(4)
DELTA(1:5) = 0.0
! DELTA(I) = G_I
! D(I,J) = D(G_I)/D(Y_J)
CASE(5)
DELTA(1) = Y(1)**2 + Y(2)**2 - LSQ
DELTA(2) = Y(1) * Y(3) + Y(2) * Y(4)
D(1, 1) = 2.0 * Y(1)
D(1, 2) = 2.0 * Y(2)
D(1, 3) = 0.0
D(1, 4) = 0.0
D(1, 5) = 0.0
D(2, 1) = Y(3)
D(2, 2) = Y(4)
D(2, 3) = Y(1)
D(2, 4) = Y(2)
D(2, 5) = 0.0
END SELECT
END
T Y(1) Y(2) Length
10.0000000 1.0998126 -0.0203017 1.0999999
20.0000000 1.0970103 -0.0810476 1.1000000
30.0000000 1.0850314 -0.1808525 1.1000004
40.0000000 1.0535675 -0.3162208 1.1000000
50.0000000 0.9896186 -0.4802662 1.1000003
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 1s 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 T v = 0, where vi = 1/(i-1)!, v is an eigenvector of B T corresponding to the eigenvalue 0. Thus
so v T y(t) is
constant. Since it has the value v T y(0) =
v1 = 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 (MATSTR=2). Normally, when IRES= 6, the matrix
is computed, saved and possibly factored, using a sparse matrix factorization routine of the users 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.
Link to example source (daesl_ex3.f90)
USE DAESL_INT
USE MP_TYPES
IMPLICIT NONE
! NEQ = Number of equations
INTEGER, PARAMETER :: NEQ=100
REAL T, Y(NEQ), YPRIME(NEQ), TEND, ATOL(NEQ), CON
INTEGER NOUT, IDO, I, NSTEPS
REAL, TARGET :: RPARAM(NEQ)
INTEGER, TARGET :: IPARAM(1)
TYPE (S_FCN_DATA) USER_DATA
EXTERNAL GCN
! Pass NEQ and A^T eigenvector V to GCN
IPARAM(1) = NEQ
RPARAM(1) = 1.0
DO I = 2, NEQ
RPARAM(I) = RPARAM(I-1) / FLOAT(I-1)
END DO
USER_DATA%RDATA=>RPARAM
USER_DATA%IDATA=>IPARAM
! Initial values for y, y'
Y = 0.0
Y(1) = 1.0
YPRIME = 0.0
YPRIME(1) = -1.0
YPRIME(2) = 1.0
ATOL = 1.0E-4
NSTEPS = 10
! Always set IDO=1 on first call
IDO = 1
DO I = 1, NSTEPS
! Output solution at T=1,2,...,10
T = I-1
TEND = I
! Set IDO = 3 on last call
IF (I == NSTEPS) IDO = 3
! User-defined Jacobian matrix structure (MATSTR=2)
CALL DAESL (T, TEND, IDO, Y, YPRIME, GCN, IYPR=0, MATSTR=2, &
NCON=1, ATOL=ATOL, FCN_DATA=USER_DATA)
END DO
! Check if solution satisfies constraint
CON = 0.0
DO I = 1, NEQ
CON = CON + RPARAM(I) * Y(I)
END DO
CALL UMACH(2, NOUT)
WRITE (NOUT, *) ' V dot Y =', CON
END
SUBROUTINE GCN (T, Y, YPRIME, DELTA, D, LDD, IRES, FCN_DATA)
USE MP_TYPES
IMPLICIT NONE
REAL T, Y(*), YPRIME(*), DELTA(*), D(LDD,*), CON, CJ
INTEGER IRES, LDD, I, NEQ
SAVE CJ
TYPE (S_FCN_DATA), OPTIONAL, INTENT(INOUT) :: FCN_DATA
NEQ = FCN_DATA%IDATA(1)
SELECT CASE (IRES)
! F_I defined here
CASE(1)
DELTA(1) = YPRIME(1) + Y(1)
DO I = 2, NEQ-1
DELTA(I) = YPRIME(I) - Y(I-1) + Y(I) / FLOAT(I)
END DO
DELTA(NEQ) = YPRIME(NEQ) - Y(NEQ-1)
! Constraint is V dot Y = 1
CASE(5)
CON = -1.0
DO I = 1, NEQ
CON = CON + FCN_DATA%RDATA(I) * Y(I)
D(1,I) = FCN_DATA%RDATA(I)
END DO
DELTA(1) = CON
! 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.
CASE(6)
CJ = DELTA(1)
! If CJ > 0 not close to zero, A is nonsingular,
! so set IRES = 0.
IF (CJ >= 1.0E-4) IRES = 0
! Solve A*x=DELTA and return x in DELTA.
CASE(7)
DELTA(1) = DELTA(1) / (1.0 + CJ)
DO I = 2, NEQ-1
DELTA(I) = (DELTA(I) + DELTA(I-1)) / (1.0 / FLOAT(I) + CJ)
END DO
DELTA(NEQ) = (DELTA(NEQ) + DELTA(NEQ-1)) / CJ
END SELECT
END
V dot Y = 1.0
PHONE: 713.784.3131 FAX:713.781.9260 |