PDE_MOL Function
Solves a system of partial differential equations of the form ut = f(x, t, u, ux, uxx) using the method of lines. The solution is represented with cubic Hermite polynomials.
Usage
result = PDE_MOL(t, y, xbreak, f_ut, f_bc)
Input Parameters
t—One-dimensional array containing values of independent variable. Element t(0) should contain the initial independent variable value (the initial time, t0) and the remaining elements of t should be values of the independent variable at which a solution is desired.
y—Two-dimensional array of size npde by nx containing the initial values, where npde is the number of differential equations and nx is the number of mesh points or lines. It must satisfy the boundary conditions.
xbreak—One-dimensional array of length nx containing the breakpoints for the cubic Hermite splines used in the x discretization. The points in xbreak must be strictly increasing. The values xbreak(0) and xbreak(nx − 1) are the endpoints of the interval.
f_ut—Scalar string specifying an user-supplied function to evaluate ut. Function f_ut accepts the following input parameters:
npde—Number of equations.x—Space variable, x.t—Time variable, t.u—One-dimensional array of length npde containing the dependent values, u.ux—One-dimensional array of length npde containing the first derivatives, ux.uxx—One-dimensional array of length npde containing the second derivative, uxx.The return value of this function is an one-dimensional array of length npde containing the computed derivatives ut.
f_bc—Scalar string specifying user-supplied procedure to evaluate boundary conditions. The boundary conditions accepted by PDE_MOL are:
note | Users must supply the values αk and βk, which determine the values γk. Since γk can depend on t values, γk' also are required. |
npde—Number of equations. (Input)x—Space variable, x. (Input)t—Time variable, t. (Input)alpha—Named variable into which an one-dimensional array of length npde containing the αk values is stored. (Output)beta—Named variable into which an one-dimensional array of length npde containing the βk values is stored. (Output)gammap—Named variable into which an one-dimensional array of length npde containing the derivatives:
is stored. (Output)
Returned Value
result—Three-dimensional array of size npde by nx by N_ELEMENTS(t) containing the approximate solutions for each specified value of the independent variable.
Input Keywords
Double—If present and nonzero, double precision is used.
Tolerance—Differential equation error tolerance. An attempt to control the local error in such a way that the global relative error is proportional to Tolerance. Default: Tolerance = 100.0*ε, where ε is machine epsilon.
Hinit—Initial step size in the t integration. This value must be nonnegative. If Hinit is zero, an initial step size of 0.001|ti+1 - ti| will be arbitrarily used. The step will be applied in the direction of integration. Default: Hinit = 0.0
Deriv_Init—Two-dimensional array that supplies the derivative values
ux(x, t(0)). This derivative information is input as:
ux(x, t(0)). This derivative information is input as:
Default: Derivatives are computed using cubic spline interpolation
Discussion
Let M = npde, N = nx and xi = xbreak(i). The routine PDE_MOL uses the method of lines to solve the partial differential equation system:
with the initial conditions:
uk = uk(x, t) at t = t0, where t0 = t(0)
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 standard basis functions for 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 equation 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, ai,k and bi,k. This system can be written in the matrix−vector form as A dc/dt = F (t, y) with c(t0) = c0 where c is a vector of coefficients of length 2M N and c0 holds the initial values of the coefficients. The last 2M equations are obtained by differentiating the boundary conditions:
for k = 1, ..., M.
The initial conditions uk(x, t0) must satisfy the boundary conditions. Also, the γk(t) must be continuous and have a smooth derivative, or the boundary conditions will not be properly imposed for t > t0.
If αk = βk = 0, it is assumed that no boundary condition is desired for the kth 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.
If the number of partial differential equations is M = 1 and the number of breakpoints is N = 4, then:
The vector c is:
c = [a1, b1, a2, b2, a3, b3, a4, b4]T
and the right-side F is:
If M > 1, then each entry in the above matrix is replaced by an M × M diagonal matrix. The element α1 is replaced by diag(α1,1, ..., α1,M ). The elements αN, β1 and βN are handled in the same manner. The φi(pj) and ψi(pj) elements are replaced by φi(pj)IM and ψi(pj)IM where IM is the identity matrix of order M. See Madsen and Sincovec (1979) for further details about discretization errors and Jacobian matrix structure.
The input array y contains the values of the ak,i. The initial values of the bk,i are obtained by using the PV-WAVE cubic spline routine CSINTERP (Chapter 4: Interpolation and Approximation) to construct functions:
such that:
The PV-WAVE routine SPVALUE, Chapter 4: Interpolation and Approximation is used to approximate the values:
There is an optional use of PDE_MOL that allows the user to provide the initial values of bk,i.
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 the ODE Function. Some of the linear solvers were removed. Numerical Jacobians are used exclusively. The algorithm is unchanged. Gear’s BDF method is used as the default because the system is typically stiff.
Four examples of PDEs are now presented that illustrate how users can interface their problems with PDE_MOL. The examples are small and not indicative of the complexities that most practitioners will face in their applications.
Example 1
The normalized linear diffusion PDE, ut = uxx, 0 ≤ x ≤ 1, t > t0, is solved. The initial values are t0 = 0, u(x, t0) = u0 = 1. There is a “zero-flux” boundary condition at x = 1, namely ux(1, t) = 0, (t > t0). The boundary value of u(0, t) is abruptly changed from u0 to the value u1 = 0.1. This transition is completed by t = tδ = 0.09.
Due to restrictions in the type of boundary conditions successfully processed by PDE_MOL, it is necessary to provide the derivative boundary value function γ' at x = 0 and x = 1. The function γ at x = 0 makes a smooth transition from the value u0 at t = t0 to the value u1 at t = tδ. The transition phase for γ' is computed by evaluating a cubic interpolating polynomial. For this purpose, the function subprogram SPVALUE, Chapter 4: Interpolation and Approximation is used. The interpolation is performed as a first step in the user-supplied procedure f_bc. The function and derivative values γ(t0) = u0, γ'(t0) = 0, γ(tδ) = u1, and γ'(tδ) = 0, are used as input to routine CSINTERP to obtain the coefficients evaluated by SPVALUE. Notice that γ'(t) = 0, t > tδ. The evaluation routine SPVALUE will not yield this value so logic in the procedure f_bc assigns γ'(t) = 0, t > tδ.
; PDE_MOL Function (PV-WAVE Advantage)
FUNCTION f_ut, npde, x, t, u, ux, uxx
; Define the PDE
ut = uxx
RETURN, ut
END
PRO f_bc, npde, x, t, alpha, beta, gammap
COMMON ex1_pde, first, ppoly
first = 1
alpha = FLTARR(npde)
beta = FLTARR(npde)
gammap = FLTARR(npde)
delta = 0.09
; Compute interpolant first time only
IF (first EQ 1) THEN BEGIN
first = 0
ppoly = CSINTERP([0.0, delta], [1.0, 0.1], $
ileft = 1, left = 0.0, iright = 1, right = 0.0)
ENDIF
; Define the boundary conditions.
IF (x EQ 0.0) THEN BEGIN
alpha(0) = 1.0
beta(0) = 0.0
gammap(0) = 0.0
; If in the boundary layer, compute nonzero gamma prime
IF (t LE delta) THEN gammap(0) = $
SPVALUE(t, ppoly, xderiv = 1)
ENDIF ELSE BEGIN
; These are for x = 1
alpha(0) = 0.0
beta(0) = 1.0
gammap(0) = 0.0
ENDELSE
END
PRO EXAMPLE1
COMMON ex1_pde, first, ppoly
npde = 1
nx = 8
nstep = 10
; Set breakpoints and initial conditions
xbreak = FINDGEN(nx)/(nx - 1)
y = FLTARR(npde, nx)
y(*) = 1.0
; Initialize the solver
t = FINDGEN(nstep)/(nstep) + 0.1
t = [0.0, t*t]
; Solve the problem
res = PDE_MOL(t, y, xbreak, 'f_ut', 'f_bc')
num = INDGEN(8) + 1
FOR i=1L, 10 DO BEGIN
PRINT, 'solution at t = ', t(i)
PRINT, num, Format = '(8I7)'
PM, res(0, *, i), Format = '(8F7.4)'
ENDFOR
END
This results in the following output:
solution at t = 0.0100000
1 2 3 4 5 6 7 8
0.9691 0.9972 0.9999 1.0000 1.0000 1.0000 1.0000 1.0000
solution at t = 0.0400000
1 2 3 4 5 6 7 8
0.6247 0.8708 0.9624 0.9908 0.9981 0.9997 1.0000 1.0000
solution at t = 0.0900000
1 2 3 4 5 6 7 8
0.1000 0.4602 0.7169 0.8671 0.9436 0.9781 0.9917 0.9951
solution at t = 0.160000
1 2 3 4 5 6 7 8
0.1000 0.3130 0.5071 0.6681 0.7893 0.8708 0.9168 0.9315
solution at t = 0.250000
1 2 3 4 5 6 7 8
0.1000 0.2567 0.4045 0.5354 0.6428 0.7224 0.7710 0.7874
solution at t = 0.360000
1 2 3 4 5 6 7 8
0.1000 0.2176 0.3292 0.4292 0.5125 0.5751 0.6139 0.6270
solution at t = 0.490000
1 2 3 4 5 6 7 8
0.1000 0.1852 0.2661 0.3386 0.3992 0.4448 0.4731 0.4827
solution at t = 0.640000
1 2 3 4 5 6 7 8
0.1000 0.1588 0.2147 0.2648 0.3066 0.3381 0.3577 0.3643
solution at t = 0.810000
1 2 3 4 5 6 7 8
0.1000 0.1387 0.1754 0.2083 0.2358 0.2565 0.2694 0.2738
solution at t = 1.00000
1 2 3 4 5 6 7 8
0.1000 0.1242 0.1472 0.1678 0.1850 0.1980 0.2060 0.2087
Example 2
Here, Problem C is solved from Sincovec and Madsen (1975). The equation is of diffusion-convection type with discontinuous coefficients. This problem illustrates a simple method for programming the evaluation routine for the derivative, ut. Note that the weak discontinuities at x = 0.5 are not evaluated in the expression for ut. The results are shown in Figure 6-7: Diffusion-Convection Type with Discontinuous Coefficients. The problem is defined as:
FUNCTION f_ut, npde, x, t, u, ux, uxx
; Define the PDE
ut = FLTARR(npde)
IF (x LE 0.5) THEN BEGIN
d = 5.0
v = 1000.0
ENDIF ELSE BEGIN
d = 1.0
v = 1.0
ENDELSE
ut(0) = d*uxx(0) - v*ux(0)
RETURN, ut
END
PRO f_bc, npde, x, t, alpha, beta, gammap
; Define the Boundary Conditions
alpha = FLTARR(npde)
beta = FLTARR(npde)
gammap = FLTARR(npde)
alpha(0) = 1.0
beta(0) = 0.0
gammap(0) = 0.0
END
PRO EXAMPLE2
npde = 1
nx = 100
nstep = 10
; Set breakpoints and initial conditions
xbreak = FINDGEN(nx)/(nx - 1)
y = FLTARR(npde, 100)
y(*) = 0.0
y(0) = 1.0
; Initialize the solver
mach = MACHINE(/FLOAT)
tol = SQRT(mach.MAX_REL_SPACE)
hinit = 0.01*tol
PRINT, 'tol = ', tol, ' and hinit = ', hinit
t = [0.0, FINDGEN(nstep)/(nstep) + 0.1]
; Solve the problem
res = PDE_MOL(t, y, xbreak, 'f_ut', 'f_bc', $
tolerance = tol, hinit = hinit)
; Plot results at current ti=ti+1
PLOT, xbreak, res(0,*,10), psym = 3, yrange=[0 , 1.25], $
title = 'Solution at t = 1.0'
END
 |
Example 3
In this example, using PDE_MOL, the linear normalized diffusion PDE ut = uxx is solved but with an optional use 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. Note that the initial data is compatible with these end conditions since the derivative function:
vanishes at x = 0 and x = 1.
This optional usage signals that the derivative of the initial data is passed by the user.
FUNCTION f_ut, npde, x, t, u, ux, uxx
; Define the PDE
ut = fltARR(npde)
ut(0) = uxx(0)
RETURN, ut
END
PRO f_bc, npde, x, t, alpha, beta, gammap
; Define the boundary conditions
alpha = FLTARR(npde)
beta = FLTARR(npde)
gammap = FLTARR(npde)
alpha(0) = 0.0
beta(0) = 1.0
gammap(0) = 0.0
RETURN
END
PRO pde_mol_ex3
npde = 1
nx = 10
nstep = 10
arg = 9.0*!Pi
; Set breakpoints and initial conditions
xbreak = FINDGEN(nx)/(nx - 1)
y = FLTARR(npde, nx)
y(0, *) = 1.0 + COS(arg*xbreak)
di = y
di(0, *) = -arg*SIN(arg*xbreak)
; Initialize the solver
mach = MACHINE(/FLOAT)
tol = SQRT(mach.MAX_REL_SPACE)
t = [FINDGEN(nstep + 1)*(nstep*0.001)/(nstep)]
; Solve the problem
res = PDE_MOL(t, y, xbreak, 'f_ut', 'f_bc', $
Tolerance=tol, Deriv_Init=di)
; Print results at every other ti=ti+1
FOR i=2L, 10, 2 DO BEGIN
PRINT, 'solution at t = ', t(i)
PM, res(0, *, i), Format = '(10F10.4)'
PRINT, 'derivative at t = ', t(i)
PM, di(0, *, i)
ENDFOR
END
This results in the following output:
solution at t = 0.00200000 1.2329 0.7671 1.2329 0.7671 1.2329
0.7671 1.2329 0.7671 1.2329 0.7671
derivative at t = 0.00200000
0.00000 9.58505e-07 7.96148e-09 1.25302e-06
-1.61002e-07 1.91968e-06 -1.60244e-06 3.85856e-06
-4.83314e-06 2.02301e-06
solution at t = 0.00400000
1.0537 0.9463 1.0537 0.9463 1.0537
0.9463 1.0537 0.9463 1.0537 0.9463
derivative at t = 0.00400000
0.00000 6.64098e-07 -5.12883e-07 8.55131e-07
-6.11177e-07 -2.76893e-06 7.84288e-08 2.97113e-06
-2.32777e-07 2.02301e-06
solution at t = 0.00600000
1.0121 0.9879 1.0121 0.9879 1.0121
0.9879 1.0121 0.9879 1.0121 0.9879
derivative at t = 0.00600000
0.00000 7.42109e-07 -5.29244e-08 -1.98559e-07
-1.19702e-06 -8.66795e-07 1.17180e-07 7.09625e-07
4.31432e-07 2.02301e-06
solution at t = 0.00800000
1.0027 0.9973 1.0027 0.9973 1.0027
0.9973 1.0027 0.9973 1.0027 0.9973
derivative at t = 0.00800000
0.00000 3.56892e-07 -3.80790e-07 -9.99308e-07
-1.96765e-07 7.72356e-07 8.50576e-08 1.11979e-07
4.74838e-07 2.02301e-06
solution at t = 0.0100000
1.0008 0.9992 1.0008 0.9992 1.0008
0.9992 1.0008 0.9992 1.0008 0.9992
derivative at t = 0.0100000
0.00000 2.40533e-07 -4.27171e-07 -1.25933e-06
3.60702e-08 6.42627e-07 -1.00818e-07 2.08207e-07
1.12973e-06 2.02301e-06
Example 4
In this example, 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. 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(x, t) = sin(πx) cos(πt). Residuals are computed at the output values of t for 0 < t ≤ 2. Output is obtained at 200 steps in increments of 0.01.
Even though the sample code PDE_MOL 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), pp 488-490, for an introductory discussion of shocks in hyperbolic systems.
FUNCTION f_ut, npde, x, t, u, ux, uxx
; Define the PDE
ut = FLTARR(npde)
ut(0) = u(1)
ut(1) = uxx(0)
RETURN, ut
END
PRO f_bc, npde, x, t, alpha, beta, gammap
; Define the boundary conditions
alpha = FLTARR(npde)
beta = FLTARR(npde)
gammap = FLTARR(npde)
alpha(0) = 1
alpha(1) = 1
beta(0) = 0
beta(1) = 0
gammap(0) = 0
gammap(1) = 0
END
PRO EXAMPLE4
npde = 2
nx = 10
nstep = 200
; Set breakpoints and initial conditions
xbreak = FINDGEN(nx)/(nx - 1)
y = FLTARR(npde, nx)
y(0, *) = SIN(!Pi*xbreak)
y(1, *) = 0
di = y
di(0, *) = !Pi*COS(!Pi*xbreak)
di(1, *) = 0.0
; Initialize the solver
mach = MACHINE(/FLOAT)
tol = SQRT(mach.MAX_REL_SPACE)
t = [0.0, 0.01 + FINDGEN(nstep)*2.0/(nstep)]
; Solve the problem
u = PDE_MOL(t, y, xbreak, 'f_ut', 'f_bc', $
Tolerance = tol, Deriv_Init = di)
err = 0.0
pde_error = FLTARR(nstep)
FOR j=1L, N_ELEMENTS(t) - 1 DO BEGIN
FOR i=0L, nx - 1 DO BEGIN
err = (err) > (u(0, i, j) - $
SIN(!Pi*xbreak(i))*COS(!Pi*t(j)))
ENDFOR
ENDFOR
PRINT, 'Maximum error in u(x, t) = ', err
END
