>  Docs Center  >  IDL Reference  >  Advanced Math and Stats  >  IMSL_PDE_MOL

IMSL_PDE_MOL

IMSL_PDE_MOL

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.

This routine requires an IDL Advanced Math and Stats license. For more information, contact your sales or technical support representative.

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 φ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 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 cubic spline routine IMSL_CSINTERP to construct functions:

such that:

The IDL Advanced Math and Stats 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.

Examples


Example 1

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:

Example 2

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 the figure that follows. 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

Example 3

n 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:

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(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, 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:

Maximum error in u(x, t) =	0.000626385

Syntax


Result = IMSL_PDE_MOL(T, Y, Xbreak, F_Ut, F_Bc, Gammap [, /DOUBLE] [, DERIV_INIT=array] [, HINIT=value] [, TOLERANCE=value])

Return 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.

Arguments


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 arguments:

  • 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 IMSL_PDE_MOL are:

Note: You 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):

Keywords


DOUBLE (optional)

If present and nonzero, then double precision is used.

DERIV_INIT (optional)

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

HINIT (optional)

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: 0.0

TOLERANCE (optional)

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: 100.0*ε, where ε is machine epsilon.

Version History


6.4

Introduced

See Also


IMSL_CSINTERP, IMSL_SPVALUE



© 2019 Harris Geospatial Solutions, Inc. |  Legal
My Account    |    Store    |    Contact Us