naginterfaces.library.pde.dim2_​ellip_​discret

naginterfaces.library.pde.dim2_ellip_discret(xmin, xmax, ymin, ymax, pdef, bndy, ngx, ngy, scheme, data=None)

dim2_ellip_discret discretizes a second-order elliptic partial differential equation (PDE) on a rectangular region.

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

https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d03/d03eef.html

Parameters

xminfloat

The lower and upper $x$ coordinates of the rectangular region respectively, $x_A$ and $x_B$.

xmaxfloat

The lower and upper $x$ coordinates of the rectangular region respectively, $x_A$ and $x_B$.

yminfloat

The lower and upper $y$ coordinates of the rectangular region respectively, $y_A$ and $y_B$.

ymaxfloat

The lower and upper $y$ coordinates of the rectangular region respectively, $y_A$ and $y_B$.

pdefcallable (alpha, beta, gamma, delta, epslon, phi, psi) = pdef(x, y, data=None)

$pdef$ must evaluate the functions $\alpha (x,y)$, $\beta (x,y)$, $\gamma (x,y)$, $\delta (x,y)$, $ϵ(x,y)$, $\varphi (x,y)$ and $\psi (x,y)$ which define the equation at a general point $(x,y)$.

Parameters

xfloat

The $x$ and $y$ coordinates of the point at which the coefficients of the partial differential equation are to be evaluated.

yfloat

The $x$ and $y$ coordinates of the point at which the coefficients of the partial differential equation are to be evaluated.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

alphafloat

$alpha$ must be set to the value of $\text{alpha}(x,y)$ at the point specified by $x$ and $y$.

betafloat

$beta$ must be set to the value of $\text{beta}(x,y)$ at the point specified by $x$ and $y$.

gammafloat

$gamma$ must be set to the value of $\text{gamma}(x,y)$ at the point specified by $x$ and $y$.

deltafloat

$delta$ must be set to the value of $\text{delta}(x,y)$ at the point specified by $x$ and $y$.

epslonfloat

$epslon$ must be set to the value of $\text{epislon}(x,y)$ at the point specified by $x$ and $y$.

phifloat

$phi$ must be set to the value of $\text{phi}(x,y)$ at the point specified by $x$ and $y$.

psifloat

$psi$ must be set to the value of $\text{psi}(x,y)$ at the point specified by $x$ and $y$.

bndycallable (a, b, c) = bndy(x, y, ibnd, data=None)

$bndy$ must evaluate the functions $a(x,y)$, $b(x,y)$, and $c(x,y)$ involved in the boundary conditions.

Parameters

xfloat

The $x$ and $y$ coordinates of the point at which the boundary conditions are to be evaluated.

yfloat

The $x$ and $y$ coordinates of the point at which the boundary conditions are to be evaluated.

ibndint

Specifies on which boundary the point $(x,y)$ lies. $ibnd=0$, $1$, $2$ or $3$ according as the point lies on the bottom, right, top or left boundary.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

afloat

$a$, $b$ and $c$ must be set to the values of the functions appearing in the boundary conditions.

bfloat

$a$, $b$ and $c$ must be set to the values of the functions appearing in the boundary conditions.

cfloat

$a$, $b$ and $c$ must be set to the values of the functions appearing in the boundary conditions.

ngxint

The number of interior grid points in the $x$- and $y$-directions respectively, $n_x$ and $n_y$. If the seven-diagonal equations are to be solved by dim2_ellip_mgrid(), $ngx-1$ and $ngy-1$ should preferably be divisible by as high a power of $2$ as possible.

ngyint

The number of interior grid points in the $x$- and $y$-directions respectively, $n_x$ and $n_y$. If the seven-diagonal equations are to be solved by dim2_ellip_mgrid(), $ngx-1$ and $ngy-1$ should preferably be divisible by as high a power of $2$ as possible.

schemestr, length 1

The type of approximation to be used for the first derivatives which occur in the partial differential equation.

$scheme=\text{'C'}$

Central differences are used.

$scheme=\text{'U'}$

Upwind differences are used.

dataarbitrary, optional

User-communication data for callback functions.

Returns

afloat, ndarray, shape $((4\times (ngx+1)\times (ngy+1))/3,7)$

$a[\text{i}-1,\text{j}-1]$, for $\text{j}=1,2,\dots ,7$, for $\text{i}=1,2,\dots ,ngx\times ngy$, contains the seven-diagonal linear equations produced by the discretization described above. If $\text{lda}>ngx\times ngy$, the remaining elements are not referenced by the function, but if $\text{lda}\ge (4\times (ngx+1)\times (ngy+1))/3$ then the array $a$ can be passed directly to dim2_ellip_mgrid(), where these elements are used as workspace.

rhsfloat, ndarray, shape $\left((4\times (ngx+1)\times (ngy+1))/3\right)$

The first $ngx\times ngy$ elements contain the right-hand sides of the seven-diagonal linear equations produced by the discretization described above. If $\text{lda}>ngx\times ngy$, the remaining elements are not referenced by the function, but if $\text{lda}\ge (4\times (ngy+1)\times (ngy+1))/3$ then the array $rhs$ can be passed directly to dim2_ellip_mgrid(), where these elements are used as workspace.

Raises

NagValueError

(errno $1$)

On entry, $scheme=⟨value⟩$.

Constraint: $scheme=\text{'C'}$ or $\text{'U'}$.

(errno $1$)

On entry, $ngy=⟨value⟩$.

Constraint: $ngy\ge 3$.

(errno $1$)

On entry, $ngx=⟨value⟩$.

Constraint: $ngx\ge 3$.

(errno $1$)

On entry, $ngx\times ngy=⟨value⟩$ and $\text{lda}=⟨value⟩$.

Constraint: $\text{lda}$ must be at least $ngx\times ngy$.

(errno $1$)

On entry, $ymin=⟨value⟩$ and $ymax=⟨value⟩$.

Constraint: $ymin<ymax$.

(errno $1$)

On entry, $xmin=⟨value⟩$ and $xmax=⟨value⟩$.

Constraint: $xmin<xmax$.

(errno $2$)

Mixed derivative in equation and derivative in boundary condition at top left boundary. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $2$)

Mixed derivative in equation and derivative in boundary condition at top right boundary. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $2$)

Mixed derivative in equation and derivative in boundary condition at bottom left boundary. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $2$)

Mixed derivative in equation and derivative in boundary condition at bottom right boundary. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $2$)

Mixed derivative in equation and derivative in boundary condition at right boundary. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $2$)

Mixed derivative in equation and derivative in boundary condition at left boundary. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $2$)

Mixed derivative in equation and derivative in boundary condition at top boundary. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $2$)

Mixed derivative in equation and derivative in boundary condition at bottom boundary. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $3$)

Null boundary condition at left boundary, top end. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $3$)

Null boundary condition at top boundary, left end. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $3$)

Null boundary condition at right boundary, top end. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $3$)

Null boundary condition at top boundary, right end. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $3$)

Null boundary condition at left boundary, bottom end. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $3$)

Null boundary condition at bottom boundary, left end. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $3$)

Null boundary condition at right boundary, bottom end. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $3$)

Null boundary condition at bottom boundary, right end. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $3$)

Null boundary condition at right boundary. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $3$)

Null boundary condition at left boundary. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $3$)

Null boundary condition at top boundary. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

(errno $3$)

Null boundary condition at bottom boundary. $(x,y)=(⟨value⟩\text{,}⟨value⟩)$.

Warns

NagAlgorithmicWarning

(errno $4$)

Equation not elliptic at some point.

(errno $5$)

There is no unique solution with Neumann Boundary conditions.

(errno $6$)

The linear equations were not diagonally dominant.

Notes

No equivalent traditional C interface for this routine exists in the NAG Library.

dim2_ellip_discret discretizes a second-order linear elliptic partial differential equation of the form

$$\alpha (x,y)\frac{{\partial }^{2}U}{\partial x^2}+\beta (x,y)\frac{{\partial }^{2}U}{\partial x\partial y}+\gamma (x,y)\frac{{\partial }^{2}U}{\partial y^2}+\delta (x,y)\frac{\partial U}{\partial x}+ϵ(x,y)\frac{\partial U}{\partial y}+\varphi (x,y)U=\psi (x,y)$$

on a rectangular region

$$\begin{array}{c}\begin{array}{c}{x}_A\le x\le x_B\\ y_A\le y\le y_B\end{array}\end{array}$$

subject to boundary conditions of the form

$$a(x,y)U+b(x,y)\frac{\partial U}{\partial n}=c(x,y)$$

where $\frac{\partial U}{\partial n}$ denotes the outward pointing normal derivative on the boundary. Equation (1) is said to be elliptic if

$$4\alpha (x,y)\gamma (x,y)\ge {\left(\beta (x,y)\right)}^2$$

for all points in the rectangular region. The linear equations produced are in a form suitable for passing directly to the multigrid function dim2_ellip_mgrid().

The equation is discretized on a rectangular grid, with $n_x$ grid points in the $x$-direction and $n_y$ grid points in the $y$-direction. The grid spacing used is, therefore,

$$\begin{array}{c}\begin{array}{c}{h}_x=(x_B-x_A)/(n_x-1)\\ h_y=(y_B-y_A)/(n_y-1)\end{array}\end{array}$$

and the coordinates of the grid points $(x_i,y_j)$ are

$$\begin{array}{c}\begin{array}{c}{x}_i=x_A+(i-1)h_x\text{,}i=1,2,\dots ,n_x\text{,}\\ y_j=y_A+(j-1)h_y\text{,}j=1,2,\dots ,n_y\text{.}\end{array}\end{array}$$

At each grid point $(x_i,y_j)$ six neighbouring grid points are used to approximate the partial differential equation, so that the equation is discretized on the seven-point stencil shown in Figure [label omitted].

[figure omitted]

For convenience the approximation $u_{ij}$ to the exact solution $U(x_i,y_j)$ is denoted by $u_O$, and the neighbouring approximations are labelled according to points of the compass as shown. Where numerical labels for the seven points are required, these are also shown.

The following approximations are used for the second derivatives:

$$\begin{array}{c}\begin{array}{ccc}\frac{{\partial }^{2}U}{\partial x^2}& \simeq & \frac{1}{h_x^2}(u_E-2u_O+u_W)\\ \frac{{\partial }^{2}U}{\partial y^2}& \simeq & \frac{1}{h_y^2}(u_N-2u_O+u_S)\\ \frac{{\partial }^{2}U}{\partial x\partial y}& \simeq & \frac{1}{2h_x{h}_y}(u_N-u_{NW}+u_E-2u_O+u_W-u_{SE}+u_S)\text{.}\end{array}\end{array}$$

Two possible schemes may be used to approximate the first derivatives:

Central Differences

$$\begin{array}{c}\begin{array}{ccc}\frac{\partial U}{\partial x}& \simeq & \frac{1}{2h_x}(u_E-u_W)\\ \frac{\partial U}{\partial y}& \simeq & \frac{1}{2h_y}(u_N-u_S)\end{array}\end{array}$$

Upwind Differences

$$\begin{array}{c}\begin{array}{cccc}\frac{\partial U}{\partial x}& \simeq & \frac{1}{h_x}(u_O-u_W)& \text{if}\delta (x,y)>0\\ \frac{\partial U}{\partial x}& \simeq & \frac{1}{h_x}(u_E-u_O)& \text{if}\delta (x,y)<0\\ \frac{\partial U}{\partial y}& \simeq & \frac{1}{h_y}(u_N-u_O)& \text{if}ϵ(x,y)>0\\ \frac{\partial U}{\partial y}& \simeq & \frac{1}{h_y}(u_O-u_S)& \text{if}ϵ(x,y)<0\text{.}\end{array}\end{array}$$

Central differences are more accurate than upwind differences, but upwind differences may lead to a more diagonally dominant matrix for those problems where the coefficients of the first derivatives are significantly larger than the coefficients of the second derivatives.

The approximations used for the first derivatives may be written in a more compact form as follows:

$$\begin{array}{c}\begin{array}{ccc}\frac{\partial U}{\partial x}& \simeq & \frac{1}{2h_x}((k_x-1)u_W-2k_x{u}_O+(k_x+1)u_E)\\ \frac{\partial U}{\partial y}& \simeq & \frac{1}{2h_y}((k_y-1)u_S-2k_y{u}_O+(k_y+1)u_N)\end{array}\end{array}$$

where $k_x=sign\left(\delta \right)$ and $k_y=sign\left(ϵ\right)$ for upwind differences, and $k_x=k_y=0$ for central differences.

At all points in the rectangular domain, including the boundary, the coefficients in the partial differential equation are evaluated by calling $pdef$, and applying the approximations. This leads to a seven-diagonal system of linear equations of the form:

$$\begin{array}{c}\begin{array}{ccccccc}& & A_{ij}^6{u}_{i-1,j+1}& +& A_{ij}^7{u}_{i,j+1}& & \\ & +& A_{ij}^3{u}_{i-1,j}& +& A_{ij}^4{u}_{ij}& +& A_{ij}^5{u}_{i+1,j}\\ & & & +& A_{ij}^1{u}_{i,j-1}& +& A_{ij}^2{u}_{i+1,j-1}=f_{ij}\text{,}i=1,2,\dots ,n_x\text{and}j=1,2,\dots ,n_y\text{,}\end{array}\end{array}$$

where the coefficients are given by

$$\begin{array}{c}\begin{array}{ccc}{A}_{ij}^1& =& \beta (x_i,y_j)\frac{1}{2h_x{h}_y}+\gamma (x_i,y_j)\frac{1}{h_y^2}+ϵ(x_i,y_j)\frac{1}{2h_y}(k_y-1)\\ A_{ij}^2& =& -\beta (x_i,y_j)\frac{1}{2h_x{h}_y}\\ A_{ij}^3& =& \alpha (x_i,y_j)\frac{1}{h_x^2}+\beta (x_i,y_j)\frac{1}{2h_x{h}_y}+\delta (x_i,y_j)\frac{1}{2h_x}(k_x-1)\\ A_{ij}^4& =& -\alpha (x_i,y_j)\frac{2}{h_x^2}-\beta (x_i,y_j)\frac{1}{h_x{h}_y}-\gamma (x_i,y_j)\frac{2}{h_y^2}-\delta (x_i,y_j)\frac{k_y}{h_x}-ϵ(x_i,y_j)\frac{k_y}{h_y}-\varphi (x_i,y_j)\\ A_{ij}^5& =& \alpha (x_i,y_j)\frac{1}{h_x^2}+\beta (x_i,y_j)\frac{1}{2h_x{h}_y}+\delta (x_i,y_j)\frac{1}{2h_x}(k_x+1)\\ A_{ij}^6& =& -\beta (x_i,y_j)\frac{1}{2h_x{h}_y}\\ A_{ij}^7& =& \beta (x_i,y_j)\frac{1}{2h_x{h}_y}+\gamma (x_i,y_j)\frac{1}{h_y^2}+ϵ(x_i,y_j)\frac{1}{2h_y}(k_y+1)\\ f_{ij}& =& \psi (x_i,y_j)\end{array}\end{array}$$

These equations then have to be modified to take account of the boundary conditions. These may be Dirichlet (where the solution is given), Neumann (where the derivative of the solution is given), or mixed (where a linear combination of solution and derivative is given).

If the boundary conditions are Dirichlet, there are an infinity of possible equations which may be applied:

$$\mu u_{ij}=\mu f_{ij},\mu \ne 0\text{.}$$

If dim2_ellip_mgrid() is used to solve the discretized equations, it turns out that the choice of $\mu$ can have a dramatic effect on the rate of convergence, and the obvious choice $\mu =1$ is not the best. Some choices may even cause the multigrid method to fail altogether. In practice it has been found that a value of the same order as the other diagonal elements of the matrix is best, and the following value has been found to work well in practice:

$$\mu ={min}_{ij}(-\{\frac{2}{h_x^2}+\frac{2}{h_y^2}\},A_{ij}^4)\text{.}$$

If the boundary conditions are either mixed or Neumann (i.e., $B\ne 0$ on return from $bndy$), then one of the points in the seven-point stencil lies outside the domain. In this case the normal derivative in the boundary conditions is used to eliminate the ‘fictitious’ point, $u_{\text{outside}}$:

$$\frac{\partial U}{\partial n}\simeq \frac{1}{2h}(u_{outside}-u_{inside})\text{.}$$

It should be noted that if the boundary conditions are Neumann and $\varphi (x,y)\equiv 0$, then there is no unique solution. The function returns with $errno$ = 5 in this case, and the seven-diagonal matrix is singular.

The four corners are treated separately. $bndy$ is called twice, once along each of the edges meeting at the corner. If both boundary conditions at this point are Dirichlet and the prescribed solution values agree, then this value is used in an equation of the form (2). If the prescribed solution is discontinuous at the corner, then the average of the two values is used. If one boundary condition is Dirichlet and the other is mixed, then the value prescribed by the Dirichlet condition is used in an equation of the form given above. Finally, if both conditions are mixed or Neumann, then two ‘fictitious’ points are eliminated using two equations of the form (3).

It is possible that equations for which the solution is known at all points on the boundary, have coefficients which are not defined on the boundary. Since this function calls $pdef$ at all points in the domain, including boundary points, arithmetic errors may occur in $pdef$ which this function cannot trap. If you have an equation with Dirichlet boundary conditions (i.e., $B=0$ at all points on the boundary), but with PDE coefficients which are singular on the boundary, then dim2_ellip_mgrid() could be called directly only using interior grid points at your discretization.

After the equations have been set up as described above, they are checked for diagonal dominance. That is to say,

$$\left|A_{ij}^4\right|>\sum _{k\ne 4}\left|A_{ij}^k\right|\text{,}i=1,2,\dots ,n_x\text{and}j=1,2,\dots ,n_y\text{.}$$

If this condition is not satisfied then the function returns with $errno$ = 6. The multigrid function:meth:dim2_ellip_mgrid may still converge in this case, but if the coefficients of the first derivatives in the partial differential equation are large compared with the coefficients of the second derivative, you should consider using upwind differences ($scheme=\text{'U'}$).

Since this function is designed primarily for use with dim2_ellip_mgrid(), this document should be read in conjunction with the document for that function.

References

Wesseling, P, 1982, MGD1 – a robust and efficient multigrid method, Multigrid Methods. Lecture Notes in Mathematics (960), 614–630, Springer–Verlag