|
|
IDL Analyst Reference Guide: Differential Equations |
|
The IMSL_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.
| Note This routine requires an IDL Analyst license. For more information, contact your ITT Visual Information Solutions sales or technical support representative. |
Result = IMSL_PDE_MOL(t, y, xbreak, f_ut, f_bc [, /DOUBLE] [, DERIV_INIT=array] [, HINIT=value] [, TOLERANCE=value] )
Three-dimensional array of size npde by nx by N_ELEMENTS(t) containing the approximate solutions for each specified value of the independent variable.
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.
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.
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.
Scalar string specifying an user-supplied function to evaluate ut. Function f_ut accepts the following arguments:
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 IMSL_PDE_MOL are:
| Note Users must supply the values ak and bk, which determine the values gk. Since gk can depend on t values, gk' also are required. |
Named variable into which an one-dimensional array of length npde containing the derivatives is stored. (Output):
If present and nonzero, double precision is used.
Two-dimensional array that supplies the derivative values ux(x, t(0)). This derivative information is input as:
Default: Derivatives are computed using cubic spline interpolation
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
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*e, where e is machine epsilon.
Let M = npde, N = nx and xi = xbreak(i). The routine IMSL_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 fi(x) and yi(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 2M N 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 gk(t) must be continuous and have a smooth derivative, or the boundary conditions will not be properly imposed for t > t0.
If ak = bk = 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.
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 x M diagonal matrix. The element a1 is replaced by diag(a1,1, ..., a1,M ). The elements aN, b1 and bN are handled in the same manner. The fi(pj) and yi(pj) elements are replaced by fi(pj)IM and yi(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 IDL Analyst cubic spline routine IMSL_CSINTERP to construct functions:
such that:
The IDL Analyst routine IMSL_SPVALUE is used to approximate the values:
There is an optional use of IMSL_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 IMSL_ODE. 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 IMSL_PDE_MOL. The examples are small and not indicative of the complexities that most practitioners will face in their applications.
This equation solves the normalized linear diffusion PDE, ut = uxx, 0 £ x £ 1, t > t0. 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 = td = 0.09.
Due to restrictions in the type of boundary conditions successfully processed by IMSL_PDE_MOL, it is necessary to provide the derivative boundary value function g¢ at x = 0 and at x = 1. The function g at x = 0 makes a smooth transition from the value u0 at t = t0 to the value u1 at t = td. The transition phase for g¢ is computed by evaluating a cubic interpolating polynomial. For this purpose, the function subprogram IMSL_SPVALUE. The interpolation is performed as a first step in the user-supplied procedure f_bc. The function and derivative values g(t0) = u0, g¢(t0) = 0, g(td) = u1, and g¢(td) = 0, are used as input to routine IMSL_CSINTERP to obtain the coefficients evaluated by IMSL_SPVALUE. Notice that g¢(t) = 0, t > td. The evaluation routine IMSL_SPVALUE will not yield this value so logic in the procedure f_bc assigns g¢(t) = 0, t > td.
Save the following code as pde_mol_example1, then compile and run:
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 = IMSL_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) = $ IMSL_SPVALUE(t, ppoly, xderiv = 1) END ELSE BEGIN ; These are for x = 1 alpha(0) = 0.0 beta(0) = 1.0 gammap(0) = 0.0 END RETURN END PRO pde_mol_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 = IMSL_PDE_MOL(t, y, xbreak, 'f_ut', 'f_bc') num = INDGEN(8) + 1 FOR i = 1, 10 DO BEGIN PRINT, 'solution at t = ', t(i) PRINT, num, FORMAT = '(8I7)' PM, res(0, *, i), FORMAT = '(8F7.4)' ENDFOR END
IDL Prints:
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
This example solves Problem C 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 8-6. The problem is defined as:
Save the following code as pde_mol_example2, then compile and run:
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
END ELSE BEGIN
d = 1.0
v = 1.0
END
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
RETURN
END
PRO pde_mol_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 = IMSL_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 = IMSL_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
|
|
In this example, using IMSL_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)px], 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.
Save the following code as pde_mol_example3, then compile and run:
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_example3
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 = IMSL_MACHINE(/FLOAT)
tol = SQRT(mach.MAX_REL_SPACE)
t = [FINDGEN(nstep + 1)*(nstep*0.001)/(nstep)]
; Solve the problem
res = IMSL_PDE_MOL(t, y, xbreak, 'f_ut', 'f_bc', $
Tolerance = tol, Deriv_Init = di)
; Print results at every other ti=ti+1
FOR i = 2, 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)
PRINT
ENDFOR
END
IDL Prints:
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
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(px) 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(px) cos(pt). 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 IMSL_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.
Save the following code as pde_mol_example4, then compile and run:
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 RETURN END PRO pde_mol_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 = IMSL_MACHINE(/FLOAT) tol = SQRT(mach.MAX_REL_SPACE) t = [0.0, 0.01 + FINDGEN(nstep)*2.0/(nstep)] ; Solve the problem u = IMSL_PDE_MOL(t, y, xbreak, 'f_ut', 'f_bc', $ Tolerance = tol, Deriv_Init = di) err = 0.0 pde_error = FLTARR(nstep) FOR j = 1, N_ELEMENTS(t) - 1 DO BEGIN FOR i = 0, 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
IDL Prints: