IVPAG
Solves an initial-value problem for ordinary differential equations using either Adams-Moulton’s or Gear’s BDF method.
Required Arguments
IDO — Flag indicating the state of the computation. (Input/Output)
IDO | State |
|---|---|
1 | Initial entry |
2 | Normal re-entry |
3 | Final call to release workspace |
4 | Return because of interrupt 1 |
5 | Return because of interrupt 2 with step accepted |
6 | Return because of interrupt 2 with step rejected |
7 | Return for new value of matrix A. |
Normally, the initial call is made with IDO = 1. The routine then sets IDO = 2, and this value is then used for all but the last call that is made with IDO = 3. This final call is only used to release workspace, which was automatically allocated by the initial call with IDO = 1. See Comment 5 for a description of the interrupts.
When IDO = 7, the matrix A at t must be recomputed and IVPAG/DIVPAG called again. No other argument (including IDO) should be changed. This value of IDO is returned only if PARAM(19) = 2.
FCN — User-supplied subroutine to evaluate functions. The usage is CALL FCN (N, T, Y, YPRIME), where
N – Number of equations. (Input)
T – Independent variable, t. (Input)
Y – Array of size N containing the dependent variable values, y. (Input)
YPRIME – Array of size N containing the values of the vector yʹ evaluated at (t, y). (Output)
See Comment 3.
See Comment 3.
FCN must be declared EXTERNAL in the calling program.
FCNJ — User-supplied subroutine to compute the Jacobian. The usage is CALL FCNJ (N, T, Y, DYPDY) where
N – Number of equations. (Input)
T – Independent variable, t. (Input)
Y – Array of size N containing the dependent variable values, y(t). (Input)
DYPDY – An array, with data structure and type determined by PARAM(14) = MTYPE, containing the required partial derivatives ∂fi∕∂yj. (Output)
These derivatives are to be evaluated at the current values of (t, y). When the Jacobian is dense, MTYPE = 0 or MTYPE = 2, the leading dimension of DYPDY has the value N. When the Jacobian matrix is banded, MTYPE = 1, and the leading dimension of DYPDY has the value 2 * NLC + NUC + 1. If the matrix is banded positive definite symmetric, MTYPE = 3, and the leading dimension of DYPDY has the value NUC + 1.
FCNJ must be declared EXTERNAL in the calling program. If PARAM(19) = IATYPE is nonzero, then FCNJ should compute the Jacobian of the righthand side of the equation Ayʹ = f(t, y). The subroutine FCNJ is used only if PARAM(13) = MITER = 1.
T — Independent variable, t. (Input/Output)
On input, T contains the initial independent variable value. On output, T is replaced by TEND unless error or other normal conditions arise. See IDO for details.
On input, T contains the initial independent variable value. On output, T is replaced by TEND unless error or other normal conditions arise. See IDO for details.
TEND — Value of t = tend where the solution is required. (Input)
The value tend may be less than the initial value of t.
The value tend may be less than the initial value of t.
Y — Array of size NEQ of dependent variables, y(t). (Input/Output)
On input, Y contains the initial values, y(t0). On output, Y contains the approximate solution, y(t).
On input, Y contains the initial values, y(t0). On output, Y contains the approximate solution, y(t).
Optional Arguments
NEQ— Number of differential equations. (Input)
Default: NEQ = size (Y,1)
Default: NEQ = size (Y,1)
A — Matrix structure used when the system is implicit. (Input)
The matrix A is referenced only if PARAM(19) = IATYPE is nonzero. Its data structure is determined by PARAM(14) = MTYPE. The matrix A must be nonsingular and MITER must be 1 or 2. See Comment 3.
The matrix A is referenced only if PARAM(19) = IATYPE is nonzero. Its data structure is determined by PARAM(14) = MTYPE. The matrix A must be nonsingular and MITER must be 1 or 2. See Comment 3.
TOL — Tolerance for error control. (Input)
An attempt is made to control the norm of the local error such that the global error is proportional to TOL.
Default: TOL = .001
An attempt is made to control the norm of the local error such that the global error is proportional to TOL.
Default: TOL = .001
PARAM — A floating-point array of size 50 containing optional parameters. (Input/Output)
If a parameter is zero, then the default value is used. These default values are given below. Parameters that concern values of the step size are applied in the direction of integration. The following parameters may be set by the user:
If a parameter is zero, then the default value is used. These default values are given below. Parameters that concern values of the step size are applied in the direction of integration. The following parameters may be set by the user:
PARAM | Meaning | |
|---|---|---|
1 | HINIT | Initial value of the step size H. Always nonnegative. Default: 0.001|tend − t0|. |
2 | HMIN | Minimum value of the step size H. Default: 0.0. |
3 | HMAX | Maximum value of the step size H. Default: No limit, beyond the machine scale, is imposed on the step size. |
4 | MXSTEP | Maximum number of steps allowed. Default: 500. |
5 | MXFCN | Maximum number of function evaluations allowed. Default: No enforced limit. |
6 | MAXORD | Maximum order of the method. Default: If Adams-Moulton method is used, then 12. If Gear’s or BDF method is used, then 5. The defaults are the maximum values allowed. |
7 | INTRP1 | If this value is set nonzero, the subroutine will return before every step with IDO = 4. See Comment 5. Default: 0. |
8 | INTRP2 | If this value is nonzero, the subroutine will return after every successful step with IDO = 5 and return with IDO = 6 after every unsuccessful step. See Comment 5. Default: 0 |
9 | SCALE | A measure of the scale of the problem, such as an approximation to the average value of a norm of the Jacobian along the solution. Default: 1.0 |
10 | INORM | Switch determining error norm. In the following, ei is the absolute value of an estimate of the error in yi(t). Default: 0. 0 — min(absolute error, relative error) = max(eiwi); i = 1, …, N, where wi = max(∣yi(t)∣, 1.0). 1 — absolute error = max(ei), i = 1, …, NEQ. 2 — max(ei / wi), i = 1, …, N where wi = max(∣yi(t)∣, FLOOR), and FLOOR is the value PARAM(11). 3 — Scaled Euclidean norm defined as where wi = max(∣yi(t)∣, 1.0). Other definitions of YMAX can be specified by the user, as explained in Comment 1. |
11 | FLOOR | Used in the norm computation associated the parameter INORM. Default: 1.0. |
12 | METH | Integration method indicator. 1 = METH selects the Adams-Moulton method. 2 = METH selects Gear’s BDF method. Default: 1. |
13 | MITER | Nonlinear solver method indicator. Note: If the problem is stiff and a chord or modified Newton method is most efficient, use MITER = 1 or = 2. 0 = MITER selects functional iteration. The value IATYPE must be set to zero with this option. 1 = MITER selects a chord method with a user-provided Jacobian. 2 = MITER selects a chord method with a divided-difference Jacobian. 3 = MITER selects a chord method with the Jacobian replaced by a diagonal matrix based on a directional derivative. The value IATYPE must be set to zero with this option. Default: 0. |
14 | MTYPE | Matrix type for A (if used) and the Jacobian (if MITER = 1 or MITER = 2). When both are used, A and the Jacobian must be of the same type. 0 = MTYPE selects full matrices. 1 = MTYPE selects banded matrices. 2 = MTYPE selects symmetric positive definite matrices. 3 = MTYPE selects banded symmetric positive definite matrices. Default: 0. |
15 | NLC | Number of lower codiagonals, used if MTYPE = 1. Default: 0. |
16 | NUC | Number of upper codiagonals, used if MTYPE = 1 or MTYPE = 3. Default: 0. |
17 | Not used. | |
18 | EPSJ | Relative tolerance used in computing divided difference Jacobians. Default: SQRT(AMACH(4)) . |
19 | IATYPE | Type of the matrix A. 0 = IATYPE implies A is not used (the system is explicit). 1 = IATYPE if A is a constant matrix. 2 = IATYPE if A depends on t. Default: 0. |
20 | LDA | Leading dimension of array A exactly as specified in the dimension statement in the calling program. Used if IATYPE is not zero. Default: N if MTYPE = 0 or = 2 NUC + NLC + 1 if MTYPE = 1 NUC + 1 if MTYPE = 3 |
21−30 | Not used. |
The following entries in the array PARAM are set by the program:
PARAM | Meaning | |
|---|---|---|
31 | HTRIAL | Current trial step size. |
32 | HMINC | Computed minimum step size. |
33 | HMAXC | Computed maximum step size. |
34 | NSTEP | Number of steps taken. |
35 | NFCN | Number of function evaluations used. |
36 | NJE | Number of Jacobian evaluations. |
37−50 | Not used. |
FORTRAN 90 Interface
Generic: CALL IVPAG (IDO, FCN, FCNJ, T, TEND, Y [, …])
Specific: The specific interface names are S_IVPAG and D_IVPAG.
FORTRAN 77 Interface
Single: CALL IVPAG (IDO, NEQ, FCN, FCNJ, A, T, TEND, TOL, PARAM, Y)
Double: The double precision name is DIVPAG.
Description
The routine IVPAG solves a system of first-order ordinary differential equations of the form yʹ = f (t, y) or Ayʹ = f (t, y) with initial conditions where A is a square nonsingular matrix of order N. Two classes of implicit linear multistep methods are available. The first is the implicit Adams-Moulton method (up to order twelve); the second uses the backward differentiation formulas BDF (up to order five). The BDF method is often called Gear’s stiff method. In both cases, because basic formulas are implicit, a system of nonlinear equations must be solved at each step. The deriviative matrix in this system has the form L = A + ηJ where η is a small number computed by IVPAG and J is the Jacobian. When it is used, this matrix is computed in the user-supplied routine FCNJ or else it is approximated by divided differences as a default. Using defaults, A is the identity matrix. The data structure for the matrix L may be identified to be real general, real banded, symmetric positive definite, or banded symmetric positive definite. The default structure for L is real general.
Comments
1. Workspace and a user-supplied error norm subroutine may be explicitly provided, if desired, by use of I2PAG/DI2PAG. The reference is:
CALL I2PAG (IDO, NEQ, FCN, FCNJ, A, T, TEND, TOL, PARAM, Y, YTEMP, YMAX, ERROR, SAVE1, SAVE2, PW, IPVT, VNORM)
None of the additional array arguments should be changed from the first call with IDO = 1 until after the final call with IDO = 3. The additional arguments are as follows:
YTEMP — Array of size NMETH. (Workspace)
YMAX — Array of size NEQ containing the maximum Y-values computed so far. (Output)
ERROR — Array of size NEQ containing error estimates for each component of Y. (Output)
SAVE1 — Array of size NEQ. (Workspace)
SAVE2 — Array of size NEQ. (Workspace)
PW — Array of size NPW. (Workspace)
IPVT — Array of size NEQ. (Workspace)
VNORM — A Fortran subroutine to compute the norm of the error. (Input)
The routine may be provided by the user, or the IMSL routine I3PRK/DI3PRK may be used. In either case, the name must be declared in a Fortran EXTERNAL statement. If usage of the IMSL routine is intended, then the name I3PRK/DI3PRK should be specified. The usage of the error norm routine is CALL VNORM (NEQ, V, Y, YMAX, ENORM) where
The routine may be provided by the user, or the IMSL routine I3PRK/DI3PRK may be used. In either case, the name must be declared in a Fortran EXTERNAL statement. If usage of the IMSL routine is intended, then the name I3PRK/DI3PRK should be specified. The usage of the error norm routine is CALL VNORM (NEQ, V, Y, YMAX, ENORM) where
Arg | Definition |
|---|---|
NEQ | Number of equations. (Input). |
V | Array of size N containing the vector whose norm is to be computed. (Input) |
Y | Array of size N containing the values of the dependent variable. (Input) |
YMAX | Array of size N containing the maximum values of |y(t)|. (Input). |
ENORM | Norm of the vector V. (Output). |
VNORM must be declared EXTERNAL in the calling program.
2. Informational errors
Type | Code | Description |
|---|---|---|
4 | 1 | After some initial success, the integration was halted by repeated error-test failures. |
4 | 2 | The maximum number of function evaluations have been used. |
4 | 3 | The maximum number of steps allowed have been used. The problem may be stiff. |
4 | 4 | On the next step T + H will equal T. Either TOL is too small, or the problem is stiff. Note: If the Adams-Moulton method is the one used in the integration, then users can switch to the BDF methods. If the BDF methods are being used, then these comments are gratuitous and indicate that the problem is too stiff for this combination of method and value of TOL. |
4 | 5 | After some initial success, the integration was halted by a test on TOL. |
4 | 6 | Integration was halted after failing to pass the error test even after dividing the initial step size by a factor of 1.0E + 10. The value TOL may be too small. |
4 | 7 | Integration was halted after failing to achieve corrector convergence even after dividing the initial step size by a factor of 1.0E + 10. The value TOL may be too small. |
4 | 8 | IATYPE is nonzero and the input matrix A multiplying y′ is singular. |
3. Both explicit systems, of the form yʹ = f (t, y), and implicit systems, Ayʹ = f (t, y), can be solved. If the system is explicit, then PARAM(19) = 0; and the matrix A is not referenced. If the system is implicit, then PARAM(14) determines the data structure of the array A. If PARAM(19) = 1, then A is assumed to be a constant matrix. The value of A used on the first call (with IDO = 1) is saved until after a call with IDO = 3. The value of A must not be changed between these calls. If PARAM(19) = 2, then the matrix is assumed to be a function of t.
4. If MTYPE is greater than zero, then MITER must equal 1 or 2.
5. If PARAM(7) is nonzero, the subroutine returns with IDO= 4 and will resume calculation at the point of interruption if re-entered with IDO = 4. If PARAM(8) is nonzero, the subroutine will interrupt immediately after decides to accept the result of the most recent trial step. The value IDO = 5 is returned if the routine plans to accept, or IDO = 6 if it plans to reject. The value IDO may be changed by the user (by changing IDO from 6 to 5) to force acceptance of a step that would otherwise be rejected. Relevant parameters to observe after return from an interrupt are IDO, HTRIAL, NSTEP, NFCN, NJE, T and Y. The array Y contains the newly computed trial value y(t).
Examples
Example 1
Euler’s equation for the motion of a rigid body not subject to external forces is
Its solution is, in terms of Jacobi elliptic functions, y(t) = sn(t; k), y2 (t) = cn(t; k), y3 (t) = dn(t; k) where k2 = 0.51. The Adams-Moulton method of IVPAG is used to solve this system, since this is the default. All parameters are set to defaults.
The last call to IVPAG with IDO = 3 releases IMSL workspace that was reserved on the first call to IVPAG. It is not necessary to release the workspace in this example because the program ends after solving a single problem. The call to release workspace is made as a model of what would be needed if the program included further calls to IMSL routines.
Because PARAM(13) = MITER = 0, functional iteration is used and so subroutine FCNJ is never called. It is included only because the calling sequence for IVPAG requires it.
USE IVPAG_INT
USE UMACH_INT
IMPLICIT NONE
INTEGER N, NPARAM
PARAMETER (N=3, NPARAM=50)
! SPECIFICATIONS FOR LOCAL VARIABLES
INTEGER IDO, IEND, NOUT
REAL A(1,1), T, TEND, TOL, Y(N)
! SPECIFICATIONS FOR SUBROUTINES
! SPECIFICATIONS FOR FUNCTIONS
EXTERNAL FCN, FCNJ
! Initialize
!
IDO = 1
T = 0.0
Y(1) = 0.0
Y(2) = 1.0
Y(3) = 1.0
TOL = 1.0E-6
! Write title
CALL UMACH (2, NOUT)
WRITE (NOUT,99998)
! Integrate ODE
IEND = 0
10 CONTINUE
IEND = IEND + 1
TEND = IEND
! The array a(*,*) is not used.
CALL IVPAG (IDO, FCN, FCNJ, T, TEND, Y, TOL=TOL)
IF (IEND .LE. 10) THEN
WRITE (NOUT,99999) T, Y
! Finish up
IF (IEND .EQ. 10) IDO = 3
GO TO 10
END IF
99998 FORMAT (11X, 'T', 14X, 'Y(1)', 11X, 'Y(2)', 11X, 'Y(3)')
99999 FORMAT (4F15.5)
END
!
SUBROUTINE FCN (N, X, Y, YPRIME)
! SPECIFICATIONS FOR ARGUMENTS
INTEGER N
REAL X, Y(N), YPRIME(N)
!
YPRIME(1) = Y(2)*Y(3)
YPRIME(2) = -Y(1)*Y(3)
YPRIME(3) = -0.51*Y(1)*Y(2)
RETURN
END
!
SUBROUTINE FCNJ (N, X, Y, DYPDY)
! SPECIFICATIONS FOR ARGUMENTS
INTEGER N
REAL X, Y(N), DYPDY(N,*)
! This subroutine is never called
RETURN
END
Output
T Y(1) Y(2) Y(3)
1.00000 0.80220 0.59705 0.81963
2.00000 0.99537 -0.09615 0.70336
3.00000 0.64141 -0.76720 0.88892
4.00000 -0.26961 -0.96296 0.98129
5.00000 -0.91173 -0.41079 0.75899
6.00000 -0.95751 0.28841 0.72967
7.00000 -0.42877 0.90342 0.95197
8.00000 0.51092 0.85963 0.93106
9.00000 0.97567 0.21926 0.71730
10.00000 0.87790 -0.47884 0.77906
Example 2
The BDF method of IVPAG is used to solve Example 2 of IVPRK. We set PARAM(12) = 2 to designate the BDF method. A chord or modified Newton method, with the Jacobian computed by divided differences, is used to solve the nonlinear equations. Thus, we set PARAM(13) = 2. The number of evaluations of yʹ is printed after the last output point, showing the efficiency gained when using a stiff solver compared to using IVPRK on this problem. The number of evaluations may vary, depending on the accuracy and other arithmetic characteristics of the computer.
USE IVPAG_INT
USE UMACH_INT
IMPLICIT NONE
INTEGER MXPARM, N
PARAMETER (MXPARM=50, N=2)
! SPECIFICATIONS FOR PARAMETERS
INTEGER MABSE, MBDF, MSOLVE
PARAMETER (MABSE=1, MBDF=2, MSOLVE=2)
! SPECIFICATIONS FOR LOCAL VARIABLES
INTEGER IDO, ISTEP, NOUT
REAL A(1,1), PARAM(MXPARM), T, TEND, TOL, Y(N)
! SPECIFICATIONS FOR SUBROUTINES
! SPECIFICATIONS FOR FUNCTIONS
EXTERNAL FCN, FCNJ
!
CALL UMACH (2, NOUT)
! Set initial conditions
T = 0.0
Y(1) = 1.0
Y(2) = 0.0
! Set error tolerance
TOL = 0.001
! Set PARAM to defaults
PARAM = 0.0E0
!
PARAM(10) = MABSE
! Select BDF method
PARAM(12) = MBDF
! Select chord method and
! a divided difference Jacobian.
PARAM(13) = MSOLVE
! Print header
WRITE (NOUT,99998)
IDO = 1
ISTEP = 0
10 CONTINUE
ISTEP = ISTEP + 24
TEND = ISTEP
! The array a(*,*) is not used.
CALL IVPAG (IDO, FCN, FCNJ, T, TEND, Y, TOL=TOL, &
PARAM=PARAM)
IF (ISTEP .LE. 240) THEN
WRITE (NOUT,'(I6,3F12.3)') ISTEP/24, T, Y
! Final call to release workspace
IF (ISTEP .EQ. 240) IDO = 3
GO TO 10
END IF
! Show number of function calls.
WRITE (NOUT,99999) PARAM(35)
99998 FORMAT (4X, 'ISTEP', 5X, 'Time', 9X, 'Y1', 11X, 'Y2')
99999 FORMAT (4X, 'Number of fcn calls with IVPAG =', F6.0)
END
SUBROUTINE FCN (N, T, Y, YPRIME)
! SPECIFICATIONS FOR ARGUMENTS
INTEGER N
REAL T, Y(N), YPRIME(N)
! SPECIFICATIONS FOR SAVE VARIABLES
REAL AK1, AK2, AK3
SAVE AK1, AK2, AK3
!
DATA AK1, AK2, AK3/294.0E0, 3.0E0, 0.01020408E0/
!
YPRIME(1) = -Y(1) - Y(1)*Y(2) + AK1*Y(2)
YPRIME(2) = -AK2*Y(2) + AK3*(1.0E0-Y(2))*Y(1)
RETURN
END
SUBROUTINE FCNJ (N, T, Y, DYPDY)
! SPECIFICATIONS FOR ARGUMENTS
INTEGER N
REAL T, Y(N), DYPDY(N,*)
!
RETURN
END
Output
ISTEP Time Y1 Y2
1 24.000 0.689 0.002
2 48.000 0.636 0.002
3 72.000 0.590 0.002
4 96.000 0.550 0.002
5 120.000 0.515 0.002
6 144.000 0.485 0.002
7 168.000 0.458 0.002
8 192.000 0.434 0.001
9 216.000 0.412 0.001
10 240.000 0.392 0.001
Number of fcn calls with IVPAG = 73.
Example 3
The BDF method of IVPAG is used to solve the so-called Robertson problem:
Output is obtained after each unit of the independent variable. A user-provided subroutine for the Jacobian matrix is used. An absolute error tolerance of 10-5 is required.
USE IVPAG_INT
USE UMACH_INT
IMPLICIT NONE
INTEGER MXPARM, N
PARAMETER (MXPARM=50, N=3)
! SPECIFICATIONS FOR PARAMETERS
INTEGER MABSE, MBDF, MSOLVE
PARAMETER (MABSE=1, MBDF=2, MSOLVE=1)
! SPECIFICATIONS FOR LOCAL VARIABLES
INTEGER IDO, ISTEP, NOUT
REAL A(1,1), PARAM(MXPARM), T, TEND, TOL, Y(N)
! SPECIFICATIONS FOR SUBROUTINES
! SPECIFICATIONS FOR FUNCTIONS
EXTERNAL FCN, FCNJ
!
CALL UMACH (2, NOUT)
! Set initial conditions
T = 0.0
Y(1) = 1.0
Y(2) = 0.0
Y(3) = 0.0
! Set error tolerance
TOL = 1.0E-5
! Set PARAM to defaults
PARAM = 0.0E0
! Select absolute error control
PARAM(10) = MABSE
! Select BDF method
PARAM(12) = MBDF
! Select chord method and
! a user-provided Jacobian.
PARAM(13) = MSOLVE
! Print header
WRITE (NOUT,99998)
IDO = 1
ISTEP = 0
10 CONTINUE
ISTEP = ISTEP + 1
TEND = ISTEP
! The array a(*,*) is not used.
CALL IVPAG (IDO, FCN, FCNJ, T, TEND, Y, TOL=TOL, PARAM=PARAM)
IF (ISTEP .LE. 10) THEN
WRITE (NOUT,'(I6,F12.2,3F13.5)') ISTEP, T, Y
! Final call to release workspace
IF (ISTEP .EQ. 10) IDO = 3
GO TO 10
END IF
99998 FORMAT (4X, 'ISTEP', 5X, 'Time', 9X, 'Y1', 11X, 'Y2', 11X, &
'Y3')
END
SUBROUTINE FCN (N, T, Y, YPRIME)
! SPECIFICATIONS FOR ARGUMENTS
INTEGER N
REAL T, Y(N), YPRIME(N)
! SPECIFICATIONS FOR SAVE VARIABLES
REAL C1, C2, C3
SAVE C1, C2, C3
!
DATA C1, C2, C3/0.04E0, 1.0E4, 3.0E7/
!
YPRIME(1) = -C1*Y(1) + C2*Y(2)*Y(3)
YPRIME(3) = C3*Y(2)**2
YPRIME(2) = -YPRIME(1) - YPRIME(3)
RETURN
END
SUBROUTINE FCNJ (N, T, Y, DYPDY)
! SPECIFICATIONS FOR ARGUMENTS
INTEGER N
REAL T, Y(N), DYPDY(N,*)
! SPECIFICATIONS FOR SAVE VARIABLES
REAL C1, C2, C3
SAVE C1, C2, C3
! SPECIFICATIONS FOR SUBROUTINES
EXTERNAL SSET
!
DATA C1, C2, C3/0.04E0, 1.0E4, 3.0E7/
! Clear array to zero
CALL SSET (N**2, 0.0, DYPDY, 1)
! Compute partials
DYPDY(1,1) = -C1
DYPDY(1,2) = C2*Y(3)
DYPDY(1,3) = C2*Y(2)
DYPDY(3,2) = 2.0*C3*Y(2)
DYPDY(2,1) = -DYPDY(1,1)
DYPDY(2,2) = -DYPDY(1,2) - DYPDY(3,2)
DYPDY(2,3) = -DYPDY(1,3)
RETURN
END
Output
ISTEP Time Y1 Y2 Y3
1 1.00 0.96647 0.00003 0.03350
2 2.00 0.94164 0.00003 0.05834
3 3.00 0.92191 0.00002 0.07806
4 4.00 0.90555 0.00002 0.09443
5 5.00 0.89153 0.00002 0.10845
6 6.00 0.87928 0.00002 0.12070
7 7.00 0.86838 0.00002 0.13160
8 8.00 0.85855 0.00002 0.14143
9 9.00 0.84959 0.00002 0.15039
10 10.00 0.84136 0.00002 0.15862
Example 4
Solve the partial differential equation
with the initial condition
u(t = 0, x) = sin x
and the boundary conditions
u(t, x = 0) = u(t, x = π) = 0
on the square [0, 1] × [0, π], using the method of lines with a piecewise-linear Galerkin discretization. The exact solution is u(t, x) = exp(1 − et) sin x. The interval [0, π] is divided into equal intervals by choosing breakpoints xk = kπ/(N + 1) for k = 0, …, N + 1. The unknown function u(t, x) is approximated by
where ɸk (x) is the piecewise linear function that equals 1 at xk and is zero at all of the other breakpoints. We approximate the partial differential equation by a system of N ordinary differential equations, A dc/dt = Rc where A and R are matrices of order N. The matrix A is given by
where h = 1/(N + 1) is the mesh spacing. The matrix R is given by
The integrals involving
are assigned the values of the integrals on the right-hand side, by using the boundary values and integration by parts. Because this system may be stiff, Gear’s BDF method is used.
In the following program, the array Y(1:N) corresponds to the vector of coefficients, c. Note that Y contains N + 2 elements; Y(0) and Y(N + 1) are used to store the boundary values. The matrix A depends on t so we set PARAM(19) = 2 and evaluate A when IVPAG returns with IDO = 7. The subroutine FCN computes the vector Rc, and the subroutine FCNJ computes R. The matrices A and R are stored as band-symmetric positive-definite structures having one upper co-diagonal.
USE IVPAG_INT
USE CONST_INT
USE WRRRN_INT
USE SSET_INT
IMPLICIT NONE
INTEGER LDA, N, NPARAM, NUC
PARAMETER (N=9, NPARAM=50, NUC=1, LDA=NUC+1)
! SPECIFICATIONS FOR PARAMETERS
INTEGER NSTEP
PARAMETER (NSTEP=4)
! SPECIFICATIONS FOR LOCAL VARIABLES
INTEGER I, IATYPE, IDO, IMETH, INORM, ISTEP, MITER, MTYPE
REAL A(LDA,N), C, HINIT, PARAM(NPARAM), PI, T, TEND, TMAX, &
TOL, XPOINT(0:N+1), Y(0:N+1)
CHARACTER TITLE*10
! SPECIFICATIONS FOR COMMON /COMHX/
COMMON /COMHX/ HX
REAL HX
! SPECIFICATIONS FOR INTRINSICS
INTRINSIC EXP, REAL, SIN
REAL EXP, REAL, SIN
! SPECIFICATIONS FOR SUBROUTINES
! SPECIFICATIONS FOR FUNCTIONS
EXTERNAL FCN, FCNJ
! Initialize PARAM
HINIT = 1.0E-3
INORM = 1
IMETH = 2
MITER = 1
MTYPE = 3
IATYPE = 2
PARAM = 0.0E0
PARAM(1) = HINIT
PARAM(10) = INORM
PARAM(12) = IMETH
PARAM(13) = MITER
PARAM(14) = MTYPE
PARAM(16) = NUC
PARAM(19) = IATYPE
! Initialize other arguments
PI = CONST('PI')
HX = PI/REAL(N+1)
CALL SSET (N-1, HX/6., A(1:,2), LDA)
CALL SSET (N, 2.*HX/3., A(2:,1), LDA)
DO 10 I=0, N + 1
XPOINT(I) = I*HX
Y(I) = SIN(XPOINT(I))
10 CONTINUE
TOL = 1.0E-6
T = 0.0
TMAX = 1.0
! Integrate ODE
IDO = 1
ISTEP = 0
20 CONTINUE
ISTEP = ISTEP + 1
TEND = TMAX*REAL(ISTEP)/REAL(NSTEP)
30 CALL IVPAG (IDO, FCN, FCNJ, T, TEND, Y(1:), NEQ=N, A=A, &
TOL=TOL, PARAM=PARAM)
! Set matrix A
IF (IDO .EQ. 7) THEN
C = EXP(-T)
CALL SSET (N-1, C*HX/6., A(1:,2), LDA)
CALL SSET (N, 2.*C*HX/3., A(2:,1), LDA)
GO TO 30
END IF
IF (ISTEP .LE. NSTEP) THEN
! Print solution
WRITE (TITLE,'(A,F5.3,A)') 'U(T=', T, ')'
CALL WRRRN (TITLE, Y, 1, N+2, 1)
! Final call to release workspace
IF (ISTEP .EQ. NSTEP) IDO = 3
GO TO 20
END IF
END
!
SUBROUTINE FCN (N, T, Y, YPRIME)
! SPECIFICATIONS FOR ARGUMENTS
INTEGER N
REAL T, Y(*), YPRIME(N)
! SPECIFICATIONS FOR LOCAL VARIABLES
INTEGER I
! SPECIFICATIONS FOR COMMON /COMHX/
COMMON /COMHX/ HX
REAL HX
! SPECIFICATIONS FOR SUBROUTINES
EXTERNAL SSCAL
!
YPRIME(1) = -2.0*Y(1) + Y(2)
DO 10 I=2, N - 1
YPRIME(I) = -2.0*Y(I) + Y(I-1) + Y(I+1)
10 CONTINUE
YPRIME(N) = -2.0*Y(N) + Y(N-1)
CALL SSCAL (N, 1.0/HX, YPRIME, 1)
RETURN
END
!
SUBROUTINE FCNJ (N, T, Y, DYPDY)
! SPECIFICATIONS FOR ARGUMENTS
INTEGER N
REAL T, Y(*), DYPDY(2,*)
! SPECIFICATIONS FOR COMMON /COMHX/
COMMON /COMHX/ HX
REAL HX
! SPECIFICATIONS FOR SUBROUTINES
EXTERNAL SSET
!
CALL SSET (N-1, 1.0/HX, DYPDY(1,2), 2)
CALL SSET (N, -2.0/HX, DYPDY(2,1), 2)
RETURN
END
Output
U(T=0.250)
1 2 3 4 5 6 7 8
0.0000 0.2321 0.4414 0.6076 0.7142 0.7510 0.7142 0.6076
9 10 11
0.4414 0.2321 0.0000
U(T=0.500)
1 2 3 4 5 6 7 8
0.0000 0.1607 0.3056 0.4206 0.4945 0.5199 0.4945 0.4206
9 10 11
0.3056 0.1607 0.0000
U(T=0.750)
1 2 3 4 5 6 7 8
0.0000 0.1002 0.1906 0.2623 0.3084 0.3243 0.3084 0.2623
9 10 11
0.1906 0.1002 0.0000
U(T=1.000)
1 2 3 4 5 6 7 8
0.0000 0.0546 0.1039 0.1431 0.1682 0.1768 0.1682 0.1431
9 10 11
0.1039 0.0546 0.0000
