naginterfaces.library.pde.dim1_parab_coll¶

naginterfaces.library.pde.dim1_parab_coll(m, ts, tout, pdedef, bndary, u, xbkpts, npoly, uinit, acc, comm, itask, itrace, ind, data=None, io_manager=None, spiked_sorder='C')[source]¶

dim1_parab_coll integrates a system of linear or nonlinear parabolic partial differential equations (PDEs) in one space variable. The spatial discretization is performed using a Chebyshev $C^{0}$ collocation method, and the method of lines is employed to reduce the PDEs to a system of ordinary differential equations (ODEs). The resulting system is solved using a backward differentiation formula method.

For full information please refer to the NAG Library document for d03pd

https://support.nag.com/numeric/nl/nagdoc_30.3/flhtml/d03/d03pdf.html

Parameters

mint

The coordinate system used:

$m = 0$

Indicates Cartesian coordinates.

$m = 1$

Indicates cylindrical polar coordinates.

$m = 2$

Indicates spherical polar coordinates.

tsfloat

The initial value of the independent variable $t$ .

toutfloat

The final value of $t$ to which the integration is to be carried out.

pdedefcallable (p, q, r, ires) = pdedef(t, x, u, ux, ires, data=None)

$p d e d e f$ must compute the values of the functions $P_{i, j}$ , $Q_{i}$ and $R_{i}$ which define the system of PDEs.

The functions may depend on $x$ , $t$ , $U$ and $U_{x}$ and must be evaluated at a set of points.

Parameters

tfloat: The current value of the independent variable $t$ .
xfloat, ndarray, shape $(nptl)$: Contains a set of mesh points at which $P_{i, j}$ , $Q_{i}$ and $R_{i}$ are to be evaluated. $x [0]$ and $x [nptl - 1]$ contain successive user-supplied break-points and the elements of the array will satisfy $x [0] < x [1] < \dots < x [nptl - 1]$ .
ufloat, ndarray, shape $(npde, nptl)$: $u [i - 1, j - 1]$ contains the value of the component $U_{i} (x, t)$ where $x = x [j - 1]$ , for $j = 1, 2, \dots, nptl$ , for $i = 1, 2, \dots, n p d e$ .
uxfloat, ndarray, shape $(npde, nptl)$: $u x [i - 1, j - 1]$ contains the value of the component $\frac{\partial U_{i} (x, t)}{\partial x}$ where $x = x [j - 1]$ , for $j = 1, 2, \dots, nptl$ , for $i = 1, 2, \dots, n p d e$ .
iresint: Set to $- 1$ or $1$ .
dataarbitrary, optional, modifiable in place: User-communication data for callback functions.

Returns

pfloat, array-like, shape $(npde, npde, nptl)$

$p [n p d e \times n p d e \times (k - 1) + n p d e \times (j - 1) + (i - 1)]$ must be set to the value of $P_{i, j} (x, t, U, U_{x})$ where $x = x [k - 1]$ , for $k = 1, 2, \dots, nptl$ , for $j = 1, 2, \dots, n p d e$ , for $i = 1, 2, \dots, n p d e$ .

qfloat, array-like, shape $(npde, nptl)$

$q [i - 1, j - 1]$ must be set to the value of $Q_{i} (x, t, U, U_{x})$ where $x = x [j - 1]$ , for $j = 1, 2, \dots, nptl$ , for $i = 1, 2, \dots, n p d e$ .

rfloat, array-like, shape $(npde, nptl)$

$r [i - 1, j - 1]$ must be set to the value of $R_{i} (x, t, U, U_{x})$ where $x = x [j - 1]$ , for $j = 1, 2, \dots, nptl$ , for $i = 1, 2, \dots, n p d e$ .

iresint

Should usually remain unchanged. However, you may set $i r e s$ to force the integration function to take certain actions as described below:

$i r e s = 2$

Indicates to the integrator that control should be passed back immediately to the calling function with the error indicator set to $e r r n o$ = 6.

$i r e s = 3$

Indicates to the integrator that the current time step should be abandoned and a smaller time step used instead. You may wish to set $i r e s = 3$ when a physically meaningless input or output value has been generated. If you consecutively set $i r e s = 3$ , dim1_parab_coll returns to the calling function with the error indicator set to $e r r n o$ = 4.

bndarycallable (beta, gamma, ires) = bndary(t, u, ux, ibnd, ires, data=None)

$b n d a r y$ must compute the functions $β_{i}$ and $γ_{i}$ which define the boundary conditions as in equation (3).

Parameters

tfloat

The current value of the independent variable $t$ .

ufloat, ndarray, shape $(npde)$

$u [i - 1]$ contains the value of the component $U_{i} (x, t)$ at the boundary specified by $i b n d$ , for $i = 1, 2, \dots, n p d e$ .

uxfloat, ndarray, shape $(npde)$

$u x [i - 1]$ contains the value of the component $\frac{\partial U_{i} (x, t)}{\partial x}$ at the boundary specified by $i b n d$ , for $i = 1, 2, \dots, n p d e$ .

ibndint

Specifies which boundary conditions are to be evaluated.

$i b n d = 0$

$b n d a r y$ must set up the coefficients of the left-hand boundary, $x = a$ .

$i b n d \neq 0$

$b n d a r y$ must set up the coefficients of the right-hand boundary, $x = b$ .

iresint

Set to $- 1$ or $1$ .

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

betafloat, array-like, shape $(npde)$

$b e t a [i - 1]$ must be set to the value of $β_{i} (x, t)$ at the boundary specified by $i b n d$ , for $i = 1, 2, \dots, n p d e$ .

gammafloat, array-like, shape $(npde)$

$g a m m a [i - 1]$ must be set to the value of $γ_{i} (x, t, U, U_{x})$ at the boundary specified by $i b n d$ , for $i = 1, 2, \dots, n p d e$ .

iresint

Should usually remain unchanged. However, you may set $i r e s$ to force the integration function to take certain actions as described below:

$i r e s = 2$

Indicates to the integrator that control should be passed back immediately to the calling function with the error indicator set to $e r r n o$ = 6.

$i r e s = 3$

Indicates to the integrator that the current time step should be abandoned and a smaller time step used instead. You may wish to set $i r e s = 3$ when a physically meaningless input or output value has been generated. If you consecutively set $i r e s = 3$ , dim1_parab_coll returns to the calling function with the error indicator set to $e r r n o$ = 4.

ufloat, array-like, shape $(npde, (nbkpts - 1) \times n p o l y + 1)$

If $i n d = 1$ the value of $u$ must be unchanged from the previous call.

xbkptsfloat, array-like, shape $(nbkpts)$

The values of the break-points in the space direction. $x b k p t s [0]$ must specify the left-hand boundary, $a$ , and $x b k p t s [nbkpts - 1]$ must specify the right-hand boundary, $b$ .

npolyint

The degree of the Chebyshev polynomial to be used in approximating the PDE solution between each pair of break-points.

uinitcallable u = uinit(npde, x, data=None)

$u i n i t$ must compute the initial values of the PDE components $U_{i} (x_{j}, t_{0})$ , for $j = 1, 2, \dots, npts$ , for $i = 1, 2, \dots, n p d e$ .

Parameters

npdeint: The number of PDEs in the system.
xfloat, ndarray, shape $(npts)$: $x [j - 1]$ , contains the values of the $j$ th mesh point, for $j = 1, 2, \dots, npts$ .
dataarbitrary, optional, modifiable in place: User-communication data for callback functions.

Returns

ufloat, array-like, shape $(n p d e, npts)$: $u [i - 1, j - 1]$ must be set to the initial value $U_{i} (x_{j}, t_{0})$ , for $j = 1, 2, \dots, npts$ , for $i = 1, 2, \dots, n p d e$ .

accfloat

A positive quantity for controlling the local error estimate in the time integration. If $E (i, j)$ is the estimated error for $U_{i}$ at the $j$ th mesh point, the error test is:

| E (i, j) | = a c c \times (1.0 + | u [i - 1, j - 1] |) .

commdict, communication object, modified in place

Communication structure.

On initial entry: need not be set.

itaskint

Specifies the task to be performed by the ODE integrator.

$i t a s k = 1$

Normal computation of output values $u$ at $t = t o u t$ .

$i t a s k = 2$

One step and return.

$i t a s k = 3$

Stop at first internal integration point at or beyond $t = t o u t$ .

itraceint

The level of trace information required from dim1_parab_coll and the underlying ODE solver. $i t r a c e$ may take the value $- 1$ , $0$ , $1$ , $2$ or $3$ .

$i t r a c e = - 1$

No output is generated.

$i t r a c e = 0$

Only warning messages from the PDE solver are printed.

$i t r a c e > 0$

Output from the underlying ODE solver is printed. This output contains details of Jacobian entries, the nonlinear iteration and the time integration during the computation of the ODE system.

If $i t r a c e < - 1$ , $- 1$ is assumed and similarly if $i t r a c e > 3$ , $3$ is assumed.

The advisory messages are given in greater detail as $i t r a c e$ increases. You are advised to set $i t r a c e = 0$ , unless you are experienced with submodule ode.

indint

Indicates whether this is a continuation call or a new integration.

$i n d = 0$

Starts or restarts the integration in time.

$i n d = 1$

Continues the integration after an earlier exit from the function. In this case, only the argument $t o u t$ should be reset between calls to dim1_parab_coll.

dataarbitrary, optional

User-communication data for callback functions.

io_managerFileObjManager, optional

Manager for I/O in this routine.

spiked_sorderstr, optional

If $u$ is spiked (i.e., has unit extent in all but one dimension, or has size $1$ ), $s p i k e d_s o r d e r$ selects the storage order to associate with it in the NAG Engine:

spiked_sorder = $'C'$: row-major storage will be used;
spiked_sorder = $'F'$: column-major storage will be used.

Two-dimensional arrays returned from callback functions in this routine must then use the same storage order.

Returns

tsfloat: The value of $t$ corresponding to the solution values in $u$ . Normally $t s = t o u t$ .
ufloat, ndarray, shape $(npde, (nbkpts - 1) \times n p o l y + 1)$: $u [i - 1, j - 1]$ will contain the computed solution at $t = t s$ .
xfloat, ndarray, shape $((nbkpts - 1) \times n p o l y + 1)$: The mesh points chosen by dim1_parab_coll in the spatial direction. The values of $x$ will satisfy $x [0] < x [1] < \dots < x [npts - 1]$ .
indint: $i n d = 1$ .

Raises

NagValueError

(errno $1$ )

On entry, on initial entry $i n d = 1$ .

Constraint: on initial entry $i n d = 0$ .

(errno $1$ )

On entry, $i n d = ⟨ v a l u e ⟩$ .

Constraint: $i n d = 0$ or $1$ .

(errno $1$ )

On entry, $i = ⟨ v a l u e ⟩$ , $x b k p t s [i - 1] = ⟨ v a l u e ⟩$ , $j = ⟨ v a l u e ⟩$ and $x b k p t s [j - 1] = ⟨ v a l u e ⟩$ .

Constraint: $x b k p t s [0] < x b k p t s [1] < \dots < x b k p t s [nbkpts - 1]$ .

(errno $1$ )

On entry, $a c c = ⟨ v a l u e ⟩$ .

Constraint: $a c c > 0.0$ .

(errno $1$ )

On entry, $npts = ⟨ v a l u e ⟩$ , $nbkpts = ⟨ v a l u e ⟩$ and $n p o l y = ⟨ v a l u e ⟩$ .

Constraint: $npts = (nbkpts - 1) \times n p o l y + 1$ .

(errno $1$ )

On entry, $n p o l y = ⟨ v a l u e ⟩$ .

Constraint: $n p o l y \leq 49$ .

(errno $1$ )

On entry, $n p o l y = ⟨ v a l u e ⟩$ .

Constraint: $n p o l y \geq 1$ .

(errno $1$ )

On entry, $npde = ⟨ v a l u e ⟩$ .

Constraint: $npde \geq 1$ .

(errno $1$ )

On entry, $nbkpts = ⟨ v a l u e ⟩$ .

Constraint: $nbkpts \geq 2$ .

(errno $1$ )

On entry, $m = ⟨ v a l u e ⟩$ and $x b k p t s [0] = ⟨ v a l u e ⟩$ .

Constraint: $m \leq 0$ or $x b k p t s [0] \geq 0.0$

(errno $1$ )

On entry, $m = ⟨ v a l u e ⟩$ .

Constraint: $m = 0$ , $1$ or $2$ .

(errno $1$ )

On entry, $i t a s k = ⟨ v a l u e ⟩$ .

Constraint: $i t a s k = 1$ , $2$ or $3$ .

(errno $1$ )

On entry, $t o u t - t s$ is too small: $t o u t = ⟨ v a l u e ⟩$ and $t s = ⟨ v a l u e ⟩$ .

(errno $1$ )

On entry, $t o u t = ⟨ v a l u e ⟩$ and $t s = ⟨ v a l u e ⟩$ .

Constraint: $t o u t > t s$ .

(errno $4$ )

In setting up the ODE system an internal auxiliary was unable to initialize the derivative. This could be due to your setting $i r e s = 3$ in $p d e d e f$ or $b n d a r y$ .

(errno $5$ )

Singular Jacobian of ODE system. Check problem formulation.

(errno $7$ )

$a c c$ was too small to start integration: $a c c = ⟨ v a l u e ⟩$ .

(errno $8$ )

$i r e s$ set to an invalid value in call to $p d e d e f$ or $b n d a r y$ .

(errno $9$ )

Serious error in internal call to an auxiliary. Increase $i t r a c e$ for further details.

(errno $11$ )

Error during Jacobian formulation for ODE system. Increase $i t r a c e$ for further details.

(errno $14$ )

Flux function appears to depend on time derivatives.

Warns

NagAlgorithmicWarning

(errno $2$ ): Underlying ODE solver cannot make further progress from the point $t s$ with the supplied value of $a c c$ . $t s = ⟨ v a l u e ⟩$ , $a c c = ⟨ v a l u e ⟩$ .
(errno $3$ ): Repeated errors in an attempted step of underlying ODE solver. Integration was successful as far as $t s$ : $t s = ⟨ v a l u e ⟩$ .
(errno $6$ ): In evaluating residual of ODE system, $i r e s = 2$ has been set in $p d e d e f$ or $b n d a r y$ . Integration is successful as far as $t s$ : $t s = ⟨ v a l u e ⟩$ .
(errno $10$ ): Integration completed, but a small change in $a c c$ is unlikely to result in a changed solution. $a c c = ⟨ v a l u e ⟩$ .

Notes

dim1_parab_coll integrates the system of parabolic equations:

npde \sum j = 1 P_{i, j} \frac{\partial U_{j}}{\partial t} + Q_{i} = x^{- m} \frac{\partial}{\partial x} (x^{m} R_{i}), i = 1, 2, \dots, npde, a \leq x \leq b, t \geq t_{0},

where $P_{i, j}$ , $Q_{i}$ and $R_{i}$ depend on $x$ , $t$ , $U$ , $U_{x}$ and the vector $U$ is the set of solution values

U (x, t) = {[U_{1} (x, t), \dots, U_{npde} (x, t)]}^{T},

and the vector $U_{x}$ is its partial derivative with respect to $x$ . Note that $P_{i, j}$ , $Q_{i}$ and $R_{i}$ must not depend on $\frac{\partial U}{\partial t}$ .

The integration in time is from $t_{0}$ to $t_{o u t}$ , over the space interval $a \leq x \leq b$ , where $a = x_{1}$ and $b = x_{nbkpts}$ are the leftmost and rightmost of a user-defined set of break-points $x_{1}, x_{2}, \dots, x_{nbkpts}$ . The coordinate system in space is defined by the value of $m$ ; $m = 0$ for Cartesian coordinates, $m = 1$ for cylindrical polar coordinates and $m = 2$ for spherical polar coordinates.

The system is defined by the functions $P_{i, j}$ , $Q_{i}$ and $R_{i}$ which must be specified in $p d e d e f$ .

The initial values of the functions $U (x, t)$ must be given at $t = t_{0}$ , and must be specified in $u i n i t$ .

The functions $R_{i}$ , for $i = 1, 2, \dots, npde$ , which may be thought of as fluxes, are also used in the definition of the boundary conditions for each equation. The boundary conditions must have the form

β_{i} (x, t) R_{i} (x, t, U, U_{x}) = γ_{i} (x, t, U, U_{x}), i = 1, 2, \dots, npde,

where $x = a$ or $x = b$ .

The boundary conditions must be specified in $b n d a r y$ . Thus, the problem is subject to the following restrictions:

$t_{0} < t_{o u t}$ , so that integration is in the forward direction;
$P_{i, j}$ , $Q_{i}$ and the flux $R_{i}$ must not depend on any time derivatives;
the evaluation of the functions $P_{i, j}$ , $Q_{i}$ and $R_{i}$ is done at both the break-points and internally selected points for each element in turn, that is $P_{i, j}$ , $Q_{i}$ and $R_{i}$ are evaluated twice at each break-point. Any discontinuities in these functions must, therefore, be at one or more of the break-points $x_{1}, x_{2}, \dots, x_{nbkpts}$ ;
at least one of the functions $P_{i, j}$ must be nonzero so that there is a time derivative present in the problem;
if $m > 0$ and $x_{1} = 0.0$ , which is the left boundary point, then it must be ensured that the PDE solution is bounded at this point. This can be done by either specifying the solution at $x = 0.0$ or by specifying a zero flux there, that is $β_{i} = 1.0$ and $γ_{i} = 0.0$ . See also Further Comments.

The parabolic equations are approximated by a system of ODEs in time for the values of $U_{i}$ at the mesh points. This ODE system is obtained by approximating the PDE solution between each pair of break-points by a Chebyshev polynomial of degree $n p o l y$ . The interval between each pair of break-points is treated by dim1_parab_coll as an element, and on this element, a polynomial and its space and time derivatives are made to satisfy the system of PDEs at $n p o l y - 1$ spatial points, which are chosen internally by the code and the break-points. In the case of just one element, the break-points are the boundaries. The user-defined break-points and the internally selected points together define the mesh. The smallest value that $n p o l y$ can take is one, in which case, the solution is approximated by piecewise linear polynomials between consecutive break-points and the method is similar to an ordinary finite element method.

In total there are $(nbkpts - 1) \times n p o l y + 1$ mesh points in the spatial direction, and $npde \times ((nbkpts - 1) \times n p o l y + 1)$ ODEs in the time direction; one ODE at each break-point for each PDE component and ( $n p o l y - 1$ ) ODEs for each PDE component between each pair of break-points. The system is then integrated forwards in time using a backward differentiation formula method.

References

Berzins, M, 1990, Developments in the NAG Library software for parabolic equations, Scientific Software Systems, (eds J C Mason and M G Cox), 59–72, Chapman and Hall

Berzins, M and Dew, P M, 1991, Algorithm 690: Chebyshev polynomial software for elliptic-parabolic systems of PDEs, ACM Trans. Math. Software (17), 178–206

Zaturska, N B, Drazin, P G and Banks, W H H, 1988, On the flow of a viscous fluid driven along a channel by a suction at porous walls, Fluid Dynamics Research (4)

NAG and Python

Return to Front

naginterfaces.library.pde.dim1_parab_coll¶

naginterfaces.library.pde.dim1_​parab_​coll¶

naginterfaces.library.pde.dim1_parab_coll¶