MMOLCH

   more...
Solves a system of partial differential equations of the form ut = f(xtuuxuxx) using the method of lines. The solution is represented with cubic Hermite polynomials.
Note: MMOLCH replaces deprecated function MOLCH.
Required Arguments
IDO — Flag indicating the state of the computation. (Input/Output)
IDO
State
1
Initial entry
2
Normal reentry
3
Final call, release workspace
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.
FCNUT — User-supplied subroutine to evaluate the function ut. The usage is CALL FCNUT (NPDESXT, UUXUXXUT[,…]) where
Required Arguments
NPDES — Number of equations. (Input)
X — Space variable, x. (Input)
T — Time variable, t. (Input)
U — Array of length NPDES containing the dependent variable values, u. (Input)
UX — Array of length NPDES containing the first derivatives ux. (Input)
UXX — Array of length NPDES containing the second derivative uxx. (Input)
UT — Array of length NPDES containing the computed derivatives, ut. (Output)
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 detailed description of this argument see FCN_DATA below. (Input/Output)
FCNUT must be declared EXTERNAL in the calling program.
FCNBC — User-supplied subroutine to evaluate the boundary conditions. The boundary conditions accepted by MMOLCH are αkuk + βkux = k(t). Users must supply the values αk and βk, and functions k(t). The usage is CALL FCNBC (NPDES, X, T, ALPHA, BETA, GAMMA []), where
Required Arguments
NPDES – Number of equations. (Input)
X — Space variable, x. This value directs which boundary condition to compute. (Input)
T — Time variable, t. (Input)
ALPHA — Array of length NPDES containing the αk values. (Output)
BETA — Array of length NPDES containing the βk values. (Output)
GAMMA — Array of length NPDES containing the values of k(t). (Output)
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 detailed description of this argument see FCN_DATA below. (Input/Output)
FCNBC must be declared EXTERNAL in the calling program.
T — Independent variable, t. (Input/Output)
On input, T supplies the initial time, t0. On output, T is set to the value to which the integration has been updated. Normally, this new value is TEND.
TEND — Value of t = tend at which the solution is desired. (Input)
XBREAK — Array of length NX containing the break points for the cubic Hermite splines used in the x discretization. (Input)
The points in the array XBREAK must be strictly increasing. The values XBREAK(1) and XBREAK(NX) are the endpoints of the interval.
Y — Array of size NPDES by NX containing the solution. (Input/Output)
The array Y contains the solution as Y(ki) = uk(xt) at x = XBREAK(i). On input, Y contains the initial values. On output, Y contains the computed solution. The user can optionally supply the derivative values, ux(xt0). The user allocates twice the space for Y to pass this information. The optional derivative information is input as
at x = XBREAK(i). The array Y contains the optional derivative values as output:
at x = XBREAK(i). To signal that this information is provided, set INPDER = 1.
Optional Arguments
NPDES — Number of differential equations. (Input)
Default: NPDES = size (Y,1).
NX — Number of mesh points or lines. (Input)
Default: NX = size (XBREAK,1).
TOL — Differential equation error tolerance. (Input)
An attempt is made to control the local error in such a way that the global error is proportional to TOL.
Default: TOL = 100 * machine precision.
HINIT — Initial step size in the t integration. (Input)
This value must be nonnegative. If HINIT is zero, an initial step size of 0.001tend-t0 will be arbitrarily used. The step will be applied in the direction of integration.
Default: HINIT = 0.0.
INPDER — Set INPDER = 1 if the user is supplying the derivative values, ux(xt0), in the array Y. (Input)
Default: INPDER = 0.
FCN_DATA — A derived type, s_fcn_data, which may be used to pass additional information to/from the user-supplied function. (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.
FORTRAN 90 Interface
Generic: CALL MMOLCH (IDO, FCNUT, FCNBC, T, TEND, XBREAK, Y [])
Specific: The specific interface names are S_MMOLCH and D_MMOLCH.
Description
Let M = NPDES, N = NX and xi = XBREAK(I). The routine MMOLCH uses the method of lines to solve the partial differential equation system
with the initial conditions
and the boundary conditions
for k = 1, M.
Cubic Hermite polynomials are used in the x variable approximation so that the trial solution is expanded in the series
where ɸi(x) and Ψi(x) are the standard basis functions for the cubic Hermite polynomials with the knots x1 < x2 <  < xN. These are piecewise cubic polynomials with continuous first derivatives. At the breakpoints, they satisfy
According to the collocation method, the coefficients of the approximation are obtained so that the trial solution satisfies the differential equations at the two Gaussian points in each subinterval,
for j = 1, N. The collocation approximation to the differential equation is
for k = 1, M and j = 1, , 2(N  1).
This is a system of 2M(N  1) ordinary differential equations in 2MN unknown coefficient functions, aik and bik. This system can be written in the matrix-vector form as with c(t0) = c0 where c is a vector of coefficients of length 2MN and c0 holds the initial values of the coefficients. The last 2M equations are obtained from the boundary conditions.
If αk = βk = 0, it is assumed that no boundary condition is desired for the k-th unknown at the left endpoint. A similar comment holds for the right endpoint. Thus, collocation is done at the endpoint. This is generally a useful feature for systems of first-order partial differential equations.
The input/output array Y contains the values of the ai,k. The initial values of the bi,k are obtained by using the IMSL cubic spline routine CSINT (see Chapter 3, “Interpolation and Approximation”) to construct functions
such that
The IMSL routine CSDER, (see Chapter 3, “Interpolation and Approximation”), is used to approximate the values
If INPDER = 1, the user should provide the initial values of bi,k.
The order of matrix A is 2MN and its maximum bandwidth is 6M-1. The band structure of the Jacobian of F with respect to c is the same as the band structure of A. This system is solved using a modified version of IVPAG. Numerical Jacobians are used exclusively. The algorithm is unchanged. Gear’s BDF method is used as the default because the system is typically stiff. For more details, see Sewell (1982).
We now present three examples of PDEs that illustrate how users can interface their problems with IMSL PDE solving software. The examples are small and not indicative of the complexities that most practitioners will face in their applications. A set of seven sample application problems, some of them with more than one equation, is given in Sincovec and Madsen (1975). Two further examples are given in Madsen and Sincovec (1979).
Comments
Informational errors
Type
Code
Description
4
1
After some initial success, the integration was halted by repeated error test failures.
4
2
On the next step, X + H will equal X. Either TOL is too small or the problem is stiff.
4
3
After some initial success, the integration was halted by a test on TOL.
4
4
Integration was halted after failing to pass the error test even after reducing the step size by a factor of 1.0E + 10. TOL may be too small.
4
5
Integration was halted after failing to achieve corrector convergence even after reducing the step size by a factor of 1.0E + 10. TOL may be too small.
Examples
Example 1
The normalized linear diffusion PDE, ut = uxx, 0 ≤ x  1, t > 0, is solved. The initial values are u(x0) = u0 = 1. There is a “zero-flux” boundary condition at x = 1, namely ux(1, t) = 0, (t > 0). The boundary value of u(0, t) is abruptly changed from u0 to the value 0, for t > 0.
When the boundary conditions are discontinuous, or incompatible with the initial conditions such as in this example, it may be important to use double precision.
 
USE MMOLCH_INT
USE WRRRN_INT
IMPLICIT NONE
 
INTEGER, PARAMETER :: NPDES=1, NX=8
 
INTEGER :: I, IDO, J, NSTEP
REAL :: HINIT, T, TEND, TOL
REAL :: XBREAK(NX), Y(NPDES,NX), U0
CHARACTER :: TITLE*19
 
EXTERNAL FCNBC, FCNUT
! SET BREAKPOINTS AND INITIAL CONDITIONS
U0 = 1.0
DO I=1,NX
XBREAK(I) = FLOAT(I-1)/FLOAT(NX-1)
Y(1,I) = U0
END DO
! SET PARAMETERS FOR MMOLCH
TOL = 10.e-4
HINIT = 0.01*TOL
T = 0.0
IDO = 1
NSTEP = 10
DO J=1,NSTEP
TEND = FLOAT(J)/FLOAT(NSTEP)
! SOLVE THE PROBLEM
CALL MMOLCH (IDO, FCNUT, FCNBC, T, TEND, XBREAK, Y, TOL=TOL, &
HINIT=HINIT)
! PRINT RESULTS
WRITE (TITLE,'(A,F4.2)') 'Solution at T =', TEND
CALL WRRRN (TITLE, Y)
END DO
! LAST CALL, TO RELEASE WORKSPACE
IDO = 3
CALL MMOLCH (IDO, FCNUT, FCNBC, T, TEND, XBREAK, Y, TOL=TOL, &
HINIT=HINIT)
STOP
END
 
SUBROUTINE FCNUT (NPDES, X, T, U, UX, UXX, UT)
 
INTEGER NPDES
REAL X, T, U(*), UX(*), UXX(*), UT(*)
! DEFINE THE PDE
UT(1) = UXX(1)
RETURN
END
 
SUBROUTINE FCNBC (NPDES, X, T, ALPHA, BTA, GAM)
 
INTEGER NPDES
REAL X, T, ALPHA(*), BTA(*), GAM(*)
! DEFINE THE BOUNDARY CONDITIONS
IF (X .EQ. 0.0) THEN
! THESE ARE FOR X=0
ALPHA(1) = 1.0
BTA(1) = 0.0
GAM(1) = 0.0
ELSE
! THESE ARE FOR X=1
ALPHA(1) = 0.0
BTA(1) = 1.0
GAM(1) = 0.0
END IF
RETURN
END
Output
 
Solution at T =0.10
1 2 3 4 5 6 7 8
0.0000 0.2507 0.4771 0.6617 0.7972 0.8857 0.9341 0.9493
Solution at T =0.20
1 2 3 4 5 6 7 8
0.0000 0.1762 0.3424 0.4893 0.6100 0.6992 0.7538 0.7721
Solution at T =0.30
1 2 3 4 5 6 7 8
0.0000 0.1356 0.2642 0.3793 0.4751 0.5471 0.5916 0.6067
Solution at T =0.40
1 2 3 4 5 6 7 8
0.0000 0.1057 0.2060 0.2960 0.3711 0.4276 0.4626 0.4745
Solution at T =0.50
1 2 3 4 5 6 7 8
0.0000 0.0825 0.1610 0.2313 0.2900 0.3341 0.3616 0.3708
Solution at T =0.60
1 2 3 4 5 6 7 8
0.0000 0.0645 0.1258 0.1808 0.2267 0.2612 0.2826 0.2899
Solution at T =0.70
1 2 3 4 5 6 7 8
0.0000 0.0504 0.0983 0.1413 0.1772 0.2041 0.2209 0.2266
Solution at T =0.80
1 2 3 4 5 6 7 8
0.0000 0.0394 0.0769 0.1105 0.1385 0.1597 0.1728 0.1772
Solution at T =0.90
1 2 3 4 5 6 7 8
0.0000 0.0309 0.0602 0.0865 0.1084 0.1249 0.1352 0.1387
Solution at T =1.00
1 2 3 4 5 6 7 8
0.0000 0.0242 0.0471 0.0677 0.0849 0.0979 0.1059 0.1086
Example 2
In this example, using MMOLCH, we solve the linear normalized diffusion PDE ut = uxx but with an optional usage that provides values of the derivatives, ux, of the initial data. Due to errors in the numerical derivatives computed by spline interpolation, more precise derivative values are required when the initial data is u(x, 0) = 1 + cos[(2n-1)πx], n > 1. The boundary conditions are “zero flux” conditions ux(0, t) = ux(1, t) = 0 for t > 0.
 
USE MMOLCH_INT
USE CONST_INT
USE WRRRN_INT
USE PGOPT_INT
IMPLICIT NONE
 
INTEGER, PARAMETER :: NPDES=1, NX=10
INTEGER :: I, IDO, J, NSTEP, N, IPAGE
REAL :: HINIT, T, TEND, TOL, XBREAK(NX)
REAL :: Y(NPDES,2*NX), PI, ARG1
CHARACTER :: TITLE*36
 
EXTERNAL FCNBC, FCNUT
REAL FLOAT
 
N = 5
PI = CONST('pi')
DO I=1,NX
XBREAK(I) = FLOAT(I-1)/FLOAT(NX-1)
ARG1 = (2.*N-1)*PI
! SET FUNCTION VALUES
Y(1,I) = 1. + COS(ARG1*XBREAK(I))
! SET FIRST DERIVATIVE VALUES
Y(1,I+NX) = -ARG1*SIN(ARG1*XBREAK(I))
END DO
! SET PARAMETERS FOR MMOLCH
TOL = 10.0e-4
HINIT = 0.01*TOL
! OUTPUT AT STEPS OF 0.001
TEND = 0.
T = 0.0
IDO = 1
NSTEP = 10
DO J=1,NSTEP
TEND = TEND + 0.001
! SOLVE THE PROBLEM
CALL MMOLCH (IDO, FCNUT, FCNBC, T, TEND, XBREAK, Y, NPDES=NPDES, &
NX=NX, HINIT=HINIT, TOL=TOL, INPDER=1)
! PRINT RESULTS
IPAGE = 70
CALL PGOPT(-1, IPAGE)
WRITE (TITLE,'(A,F5.3)') 'Solution and derivatives at T =', T
CALL WRRRN (TITLE, Y)
END DO
! LAST CALL, TO RELEASE WORKSPACE
IDO = 3
CALL MMOLCH (IDO, FCNUT, FCNBC, T, TEND, XBREAK, Y, NPDES=NPDES, &
NX=NX, HINIT=HINIT, TOL=TOL, INPDER=1)
END
 
SUBROUTINE FCNUT (NPDES, X, T, U, UX, UXX, UT)
 
INTEGER NPDES
REAL X, T, U(*), UX(*), UXX(*), UT(*)
! DEFINE THE PDE
UT(1) = UXX(1)
RETURN
END
 
SUBROUTINE FCNBC (NPDES, X, T, ALPHA, BTA, GAM)
 
INTEGER NPDES
REAL X, T, ALPHA(*), BTA(*), GAM(*)
! DEFINE THE BOUNDARY CONDITIONS
ALPHA(1) = 0.0
BTA(1) = 1.0
GAM(1) = 0.0
RETURN
END
Output
 
Solution and derivatives at T =0.001
1 2 3 4 5 6 7 8 9
1.482 0.518 1.482 0.518 1.482 0.518 1.482 0.518 1.482
10 11 12 13 14 15 16 17 18
0.518 0.000 0.000 0.000 0.000 -0.000 0.000 -0.000 0.000
19 20
-0.000 -0.000
Solution and derivatives at T =0.002
1 2 3 4 5 6 7 8 9
1.235 0.765 1.235 0.765 1.235 0.765 1.235 0.765 1.235
10 11 12 13 14 15 16 17 18
0.765 0.000 0.000 0.000 0.000 -0.000 0.000 -0.000 0.000
19 20
-0.000 0.000
Solution and derivatives at T =0.003
1 2 3 4 5 6 7 8 9
1.114 0.886 1.114 0.886 1.114 0.886 1.114 0.886 1.114
10 11 12 13 14 15 16 17 18
0.886 0.000 0.000 0.000 0.000 -0.000 0.000 -0.000 0.000
19 20
-0.000 -0.000
Solution and derivatives at T =0.004
1 2 3 4 5 6 7 8 9
1.055 0.945 1.055 0.945 1.055 0.945 1.055 0.945 1.055
10 11 12 13 14 15 16 17 18
0.945 0.000 0.000 0.000 0.000 -0.000 -0.000 0.000 0.000
19 20
-0.000 -0.000
Solution and derivatives at T =0.005
1 2 3 4 5 6 7 8 9
1.027 0.973 1.027 0.973 1.027 0.973 1.027 0.973 1.027
10 11 12 13 14 15 16 17 18
0.973 0.000 -0.000 0.000 0.000 0.000 -0.000 -0.000 0.000
19 20
-0.000 -0.000
Solution and derivatives at T =0.006
1 2 3 4 5 6 7 8 9
1.013 0.987 1.013 0.987 1.013 0.987 1.013 0.987 1.013
10 11 12 13 14 15 16 17 18
0.987 0.000 0.000 0.000 -0.000 0.000 0.000 -0.000 0.000
19 20
-0.000 -0.000
Solution and derivatives at T =0.007
1 2 3 4 5 6 7 8 9
1.006 0.994 1.006 0.994 1.006 0.994 1.006 0.994 1.006
10 11 12 13 14 15 16 17 18
0.994 0.000 0.000 0.000 -0.000 0.000 0.000 -0.000 -0.000
19 20
-0.000 -0.000
Solution and derivatives at T =0.008
1 2 3 4 5 6 7 8 9
1.003 0.997 1.003 0.997 1.003 0.997 1.003 0.997 1.003
10 11 12 13 14 15 16 17 18
0.997 0.000 0.000 0.000 -0.000 -0.000 0.000 0.000 -0.000
19 20
-0.000 -0.000
Solution and derivatives at T =0.009
1 2 3 4 5 6 7 8 9
1.002 0.998 1.002 0.998 1.002 0.998 1.002 0.998 1.002
10 11 12 13 14 15 16 17 18
0.998 0.000 0.000 0.000 -0.000 -0.000 -0.000 0.000 -0.000
19 20
-0.000 0.000
Solution and derivatives at T =0.010
1 2 3 4 5 6 7 8 9
1.001 0.999 1.001 0.999 1.001 0.999 1.001 0.999 1.001
10 11 12 13 14 15 16 17 18
0.999 0.000 0.000 0.000 -0.000 -0.000 -0.000 0.000 -0.000
19 20
-0.000 -0.000
 
Example 3
In this example, we consider the linear normalized hyperbolic PDE, utt = uxx, the “vibrating string” equation. This naturally leads to a system of first order PDEs. Define a new dependent variable ut = v. Then, vt = uxx is the second equation in the system. We take as initial data u(x,0) = sin(πx) and ut(x, 0) = v(x, 0) = 0. The ends of the string are fixed so u(0, t) = u(1, t) = v(0, t) = v(1, t) = 0. The exact solution to this problem is u(xt) = sin(πx) cos(πt). Residuals are computed at the output values of t for 0 < t2. Output is obtained at 200 steps in increments of 0.01.
Even though the sample code MMOLCH gives satisfactory results for this PDE, users should be aware that for nonlinear problems, “shocks” can develop in the solution. The appearance of shocks may cause the code to fail in unpredictable ways. See Courant and Hilbert (1962), pages 488-490, for an introductory discussion of shocks in hyperbolic systems.
 
USE MMOLCH_INT
USE UMACH_INT
USE CONST_INT
 
IMPLICIT NONE
 
INTEGER, PARAMETER :: NPDES=2, NX=10
INTEGER :: I, IDO, J, NOUT, NSTEP
REAL :: HINIT, T, TEND, TOL, XBREAK(NX)
REAL :: Y(NPDES,2*NX), PI, ERROR, ERRU
CHARACTER :: TITLE*36
 
EXTERNAL FCNBC, FCNUT
REAL FLOAT
CALL UMACH (2,NOUT)
! SET BREAKPOINTS AND INITIAL CONDITIONS
PI = CONST('pi')
DO I=1,NX
XBREAK(I) = FLOAT(I-1)/FLOAT(NX-1)
! SET FUNCTION VALUES
Y(1,I) = SIN(PI*XBREAK(I))
Y(2,I) = 0.
! SET FIRST DERIVATIVE VALUES
Y(1,I+NX) = PI*COS(PI*XBREAK(I))
Y(2,I+NX) = 0.
END DO
! SET PARAMETERS FOR MMOLCH
TOL = 10.0e-4
HINIT = 0.01*TOL
! OUTPUT AT STEPS OF 0.01
TEND = 0.
T = 0.0
IDO = 1
NSTEP = 200
DO J=1,NSTEP
TEND = TEND + 0.01
! SOLVE THE PROBLEM
CALL MMOLCH (IDO, FCNUT, FCNBC, T, TEND, XBREAK, Y, &
HINIT=HINIT, TOL=TOL, INPDER=1)
! COMPUTE MAXIMUM ERROR
ERRU = 0.0
DO I=1,NX
ERROR = Y(1,I) - SIN(PI*XBREAK(I))*COS(PI*TEND)
ERRU = AMAX1(ERRU,ABS(ERROR))
END DO
END DO
! PRINT ERROR
WRITE (NOUT, *) ' Maximum error in u(x,t): ', ERRU
! LAST CALL, TO RELEASE WORKSPACE
IDO = 3
CALL MMOLCH (IDO, FCNUT, FCNBC, T, TEND, XBREAK, Y, &
HINIT=HINIT, TOL=TOL, INPDER=1)
END
 
SUBROUTINE FCNUT (NPDES, X, T, U, UX, UXX, UT)
 
INTEGER NPDES
REAL X, T, U(*), UX(*), UXX(*), UT(*)
! DEFINE THE PDEs
UT(1) = U(2)
UT(2) = UXX(1)
RETURN
END
 
SUBROUTINE FCNBC (NPDES, X, T, ALPHA, BTA, GAM)
 
INTEGER NPDES
REAL X, T, ALPHA(*), BTA(*), GAM(*)
! DEFINE THE BOUNDARY CONDITIONS
ALPHA(1) = 1.0
BTA(1) = 0.0
GAM(1) = 0.0
ALPHA(2) = 1.0
BTA(2) = 0.0
GAM(2) = 0.0
RETURN
END
 
Output
 
Maximum error in u(x,t): 5.49525E-3