Integer, Intent (In)	::	ngx, ngy, lda
Integer, Intent (Inout)	::	ifail
Real (Kind=nag_wp), Intent (In)	::	xmin, xmax, ymin, ymax
Real (Kind=nag_wp), Intent (Inout)	::	a(lda,7)
Real (Kind=nag_wp), Intent (Out)	::	rhs(lda)
Character (1), Intent (In)	::	scheme
External	::	pdef, bndy

C Header Interface

#include <nag.h>

void

d03eef_ (const double *xmin, const double *xmax, const double *ymin, const double *ymax,
void (NAG_CALL *pdef)(const double *x, const double *y, double *alpha, double *beta, double *gamma, double *delta, double *epslon, double *phi, double *psi),
void (NAG_CALL *bndy)(const double *x, const double *y, double *a, double *b, double *c, const Integer *ibnd),
const Integer *ngx, const Integer *ngy, const Integer *lda, double a[], double rhs[], const char *scheme, Integer *ifail, const Charlen length_scheme)

The routine may be called by the names d03eef or nagf_pde_dim2_ellip_discret.

3 Description

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

α (x, y) \frac{\partial^{2} U}{\partial x^{2}} + β (x, y) \frac{\partial^{2} U}{\partial x \partial y} + γ (x, y) \frac{\partial^{2} U}{\partial y^{2}} + δ (x, y) \frac{\partial U}{\partial x} + ε (x, y) \frac{\partial U}{\partial y} + ϕ (x, y) U = ψ (x, y)

(1)

on a rectangular region

\begin{matrix} x_{A} \leq x \leq x_{B} \\ y_{A} \leq y \leq y_{B} \end{matrix}

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 α (x, y) γ (x, y) \geq {(β (x, y))}^{2}

for all points in the rectangular region. The linear equations produced are in a form suitable for passing directly to the multigrid routine d03edf.

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{matrix} h_{x} = (x_{B} - x_{A}) / (n_{x} - 1) \\ h_{y} = (y_{B} - y_{A}) / (n_{y} - 1) \end{matrix}

and the coordinates of the grid points

(x_{i}, y_{j})

are

\begin{matrix} x_{i} = x_{A} + (i - 1) h_{x}, i = 1, 2, \dots, n_{x}, \\ y_{j} = y_{A} + (j - 1) h_{y}, j = 1, 2, \dots, n_{y} . \end{matrix}

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

Figure 1

For convenience the approximation

u_{i j}

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}{l} \frac{\partial^{2} U}{\partial x^{2}} & ≃ & \frac{1}{h_{x}^{2}} (u_{E} - 2 u_{O} + u_{W}) \\ \frac{\partial^{2} U}{\partial y^{2}} & ≃ & \frac{1}{h_{y}^{2}} (u_{N} - 2 u_{O} + u_{S}) \\ \frac{\partial^{2} U}{\partial x \partial y} & ≃ & \frac{1}{2 h_{x} h_{y}} (u_{N} - u_{NW} + u_{E} - 2 u_{O} + u_{W} - u_{SE} + u_{S}) . \end{array}

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

Central Differences

\begin{array}{l} \frac{\partial U}{\partial x} & ≃ & \frac{1}{2 h_{x}} (u_{E} - u_{W}) \\ \frac{\partial U}{\partial y} & ≃ & \frac{1}{2 h_{y}} (u_{N} - u_{S}) \end{array}

Upwind Differences

\begin{array}{l} \frac{\partial U}{\partial x} & ≃ & \frac{1}{h_{x}} (u_{O} - u_{W}) & if δ (x, y) > 0 \\ \frac{\partial U}{\partial x} & ≃ & \frac{1}{h_{x}} (u_{E} - u_{O}) & if δ (x, y) < 0 \\ \frac{\partial U}{\partial y} & ≃ & \frac{1}{h_{y}} (u_{N} - u_{O}) & if ε (x, y) > 0 \\ \frac{\partial U}{\partial y} & ≃ & \frac{1}{h_{y}} (u_{O} - u_{S}) & if ε (x, y) < 0 . \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}{l} \frac{\partial U}{\partial x} & ≃ & \frac{1}{2 h_{x}} ((k_{x} - 1) u_{W} - 2 k_{x} u_{O} + (k_{x} + 1) u_{E}) \\ \frac{\partial U}{\partial y} & ≃ & \frac{1}{2 h_{y}} ((k_{y} - 1) u_{S} - 2 k_{y} u_{O} + (k_{y} + 1) u_{N}) \end{array}

where

k_{x} = sign δ

and

k_{y} = sign ε

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}{l} A_{i j}^{6} u_{i - 1, j + 1} & + & A_{i j}^{7} u_{i, j + 1} \\ + & A_{i j}^{3} u_{i - 1, j} & + & A_{i j}^{4} u_{i j} & + & A_{i j}^{5} u_{i + 1, j} \\ + & A_{i j}^{1} u_{i, j - 1} & + & A_{i j}^{2} u_{i + 1, j - 1} = f_{i j}, i = 1, 2, \dots, n_{x} ​ and ​ j = 1, 2, \dots, n_{y}, \end{array}

where the coefficients are given by

\begin{array}{l} A_{i j}^{1} & = & β (x_{i}, y_{j}) \frac{1}{2 h_{x} h_{y}} + γ (x_{i}, y_{j}) \frac{1}{h_{y}^{2}} + ε (x_{i}, y_{j}) \frac{1}{2 h_{y}} (k_{y} - 1) \\ A_{i j}^{2} & = & - β (x_{i}, y_{j}) \frac{1}{2 h_{x} h_{y}} \\ A_{i j}^{3} & = & α (x_{i}, y_{j}) \frac{1}{h_{x}^{2}} + β (x_{i}, y_{j}) \frac{1}{2 h_{x} h_{y}} + δ (x_{i}, y_{j}) \frac{1}{2 h_{x}} (k_{x} - 1) \\ A_{i j}^{4} & = & - α (x_{i}, y_{j}) \frac{2}{h_{x}^{2}} - β (x_{i}, y_{j}) \frac{1}{h_{x} h_{y}} - γ (x_{i}, y_{j}) \frac{2}{h_{y}^{2}} - δ (x_{i}, y_{j}) \frac{k_{y}}{h_{x}} - ε (x_{i}, y_{j}) \frac{k_{y}}{h_{y}} - ϕ (x_{i}, y_{j}) \\ A_{i j}^{5} & = & α (x_{i}, y_{j}) \frac{1}{h_{x}^{2}} + β (x_{i}, y_{j}) \frac{1}{2 h_{x} h_{y}} + δ (x_{i}, y_{j}) \frac{1}{2 h_{x}} (k_{x} + 1) \\ A_{i j}^{6} & = & - β (x_{i}, y_{j}) \frac{1}{2 h_{x} h_{y}} \\ A_{i j}^{7} & = & β (x_{i}, y_{j}) \frac{1}{2 h_{x} h_{y}} + γ (x_{i}, y_{j}) \frac{1}{h_{y}^{2}} + ε (x_{i}, y_{j}) \frac{1}{2 h_{y}} (k_{y} + 1) \\ f_{i j} & = & ψ (x_{i}, y_{j}) \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:

μ u_{i j} = μ f_{i j}, μ \neq 0 .

(2)

If d03edf is used to solve the discretized equations, it turns out that the choice of

μ

can have a dramatic effect on the rate of convergence, and the obvious choice

μ = 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:

μ = \min_{i j} (- \{\frac{2}{h_{x}^{2}} + \frac{2}{h_{y}^{2}}\}, A_{i j}^{4}) .

If the boundary conditions are either mixed or Neumann (i.e.,

B \neq 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_{outside}

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

(3)

It should be noted that if the boundary conditions are Neumann and

ϕ (x, y) \equiv 0

, then there is no unique solution. The routine returns with

ifail = 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 routine calls pdef at all points in the domain, including boundary points, arithmetic errors may occur in pdef which this routine 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 d03edf 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,

|A_{i j}^{4}| > \sum_{k \neq 4} |A_{i j}^{k}|, i = 1, 2, \dots, n_{x} ​ and ​ j = 1, 2, \dots, n_{y} .

If this condition is not satisfied then the routine returns with

ifail = 6

. The multigrid routined03edf 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 ='U'

Since this routine is designed primarily for use with d03edf, this document should be read in conjunction with the document for that routine.

4 References

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

5 Arguments

1: $xmin$ – Real (Kind=nag_wp) Input

2: $xmax$ – Real (Kind=nag_wp) Input

On entry: the lower and upper

x

coordinates of the rectangular region respectively,

x_{A}

and

x_{B}

Constraint:

xmin < xmax

3: $ymin$ – Real (Kind=nag_wp) Input

4: $ymax$ – Real (Kind=nag_wp) Input

On entry: the lower and upper

y

coordinates of the rectangular region respectively,

y_{A}

and

y_{B}

Constraint:

ymin < ymax

5: $pdef$ – Subroutine, supplied by the user. External Procedure

pdef must evaluate the functions

α (x, y)

β (x, y)

γ (x, y)

δ (x, y)

ε (x, y)

ϕ (x, y)

and

ψ (x, y)

which define the equation at a general point

(x, y)

The specification of pdef is:

Fortran Interface

Subroutine pdef (

x, y, alpha, beta, gamma, delta, epslon, phi, psi)

Real (Kind=nag_wp), Intent (In)	::	x, y
Real (Kind=nag_wp), Intent (Out)	::	alpha, beta, gamma, delta, epslon, phi, psi

C Header Interface

void	pdef_ (const double x, const double y, double alpha, double beta, double gamma, double delta, double epslon, double phi, double *psi)

1: $x$ – Real (Kind=nag_wp) Input
2: $y$ – Real (Kind=nag_wp) Input: On entry: the $x$ and $y$ coordinates of the point at which the coefficients of the partial differential equation are to be evaluated.
3: $alpha$ – Real (Kind=nag_wp) Output
4: $beta$ – Real (Kind=nag_wp) Output
5: $gamma$ – Real (Kind=nag_wp) Output
6: $delta$ – Real (Kind=nag_wp) Output
7: $epslon$ – Real (Kind=nag_wp) Output
8: $phi$ – Real (Kind=nag_wp) Output
9: $psi$ – Real (Kind=nag_wp) Output: On exit: alpha, beta, gamma, delta, epslon, phi and psi must be set to the values of $α (x, y)$ , $β (x, y)$ , $γ (x, y)$ , $δ (x, y)$ , $ε (x, y)$ , $ϕ (x, y)$ and $ψ (x, y)$ respectively at the point specified by x and y.

pdef must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d03eef is called. Arguments denoted as Input must not be changed by this procedure.

Note: pdef should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d03eef. If your code inadvertently does return any NaNs or infinities, d03eef is likely to produce unexpected results.

6: $bndy$ – Subroutine, supplied by the user. External Procedure

bndy must evaluate the functions

a (x, y)

b (x, y)

, and

c (x, y)

involved in the boundary conditions.

The specification of bndy is:

Fortran Interface

Subroutine bndy (

x, y, a, b, c, ibnd)

Integer, Intent (In)	::	ibnd
Real (Kind=nag_wp), Intent (In)	::	x, y
Real (Kind=nag_wp), Intent (Out)	::	a, b, c

C Header Interface

void	bndy_ (const double x, const double y, double a, double b, double c, const Integer ibnd)

1: $x$ – Real (Kind=nag_wp) Input
2: $y$ – Real (Kind=nag_wp) Input: On entry: the $x$ and $y$ coordinates of the point at which the boundary conditions are to be evaluated.
3: $a$ – Real (Kind=nag_wp) Output
4: $b$ – Real (Kind=nag_wp) Output
5: $c$ – Real (Kind=nag_wp) Output: On exit: a, b and c must be set to the values of the functions appearing in the boundary conditions.
6: $ibnd$ – Integer Input: On entry: 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.

bndy must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d03eef is called. Arguments denoted as Input must not be changed by this procedure.

Note: bndy should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d03eef. If your code inadvertently does return any NaNs or infinities, d03eef is likely to produce unexpected results.

7: $ngx$ – Integer Input

8: $ngy$ – Integer Input

On entry: 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 d03edf,

ngx - 1

and

ngy - 1

should preferably be divisible by as high a power of

2

as possible.

Constraints:

$ngx \geq 3$ ;
$ngy \geq 3$ .

9: $lda$ – Integer Input

On entry: the first dimension of the array a and the dimension of the array rhs as declared in the (sub)program from which d03eef is called.

Constraint: if only the seven-diagonal equations are required,

lda \geq ngx \times ngy

. If a call to this routine is to be followed by a call to d03edf to solve the seven-diagonal linear equations,

lda \geq (4 \times (ngx + 1) \times (ngy + 1)) / 3

Note: this routine only checks the former condition. d03edf, if called, will check the latter condition.

10: $a (lda, 7)$ – Real (Kind=nag_wp) array Output

On exit:

a (i, j)

, for

i = 1, 2, \dots, ngx \times ngy

and

j = 1, 2, \dots, 7

, contains the seven-diagonal linear equations produced by the discretization described above. If

lda > ngx \times ngy

, the remaining elements are not referenced by the routine, but if

lda \geq (4 \times (ngx + 1) \times (ngy + 1)) / 3

then the array a can be passed directly to d03edf, where these elements are used as workspace.

11: $rhs (lda)$ – Real (Kind=nag_wp) array Output

On exit: the first

ngx \times ngy

elements contain the right-hand sides of the seven-diagonal linear equations produced by the discretization described above. If

lda > ngx \times ngy

, the remaining elements are not referenced by the routine, but if

lda \geq (4 \times (ngy + 1) \times (ngy + 1)) / 3

then the array rhs can be passed directly to d03edf, where these elements are used as workspace.

12: $scheme$ – Character(1) Input

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

$scheme ='C'$: Central differences are used.
$scheme ='U'$: Upwind differences are used.

Constraint:

scheme ='C'

'U'

Note: generally speaking, if at least one of the coefficients multiplying the first derivatives (delta or epslon as returned by pdef) are large compared with the coefficients multiplying the second derivatives, then upwind differences may be more appropriate. Upwind differences are less accurate than central differences, but may result in more rapid convergence for strongly convective equations. The easiest test is to try both schemes.

13: $ifail$ – Integer Input/Output

On entry: ifail must be set to

0

- 1 or 1

. If you are unfamiliar with this argument you should refer to Section 4 in the Introduction to the NAG Library FL Interface for details.

For environments where it might be inappropriate to halt program execution when an error is detected, the value

- 1 or 1

is recommended. If the output of error messages is undesirable, then the value

1

is recommended. Otherwise, because for this routine the values of the output arguments may be useful even if

ifail \neq 0

on exit, the recommended value is

- 1

. When the value $- 1 or 1$ is used it is essential to test the value of ifail on exit.

On exit:

ifail = 0

unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry

ifail = 0

- 1

, explanatory error messages are output on the current error message unit (as defined by x04aaf).

Errors or warnings detected by the routine:

Note: in some cases d03eef may return useful information.

$ifail = 1$: On entry, $ngx = 〈value〉$ .
Constraint: $ngx \geq 3$ .

On entry, $ngx \times ngy = 〈value〉$ and $lda = 〈value〉$ .
Constraint: lda must be at least $ngx \times ngy$ .

On entry, $ngy = 〈value〉$ .
Constraint: $ngy \geq 3$ .

On entry, $scheme = 〈value〉$ .
Constraint: $scheme ='C'$ or $'U'$ .

On entry, $xmin = 〈value〉$ and $xmax = 〈value〉$ .
Constraint: $xmin < xmax$ .

On entry, $ymin = 〈value〉$ and $ymax = 〈value〉$ .
Constraint: $ymin < ymax$ .

$ifail = 2$: Mixed derivative in equation and derivative in boundary condition at bottom boundary. $(x, y) = (〈value〉, 〈value〉)$ .

Mixed derivative in equation and derivative in boundary condition at bottom left boundary. $(x, y) = (〈value〉, 〈value〉)$ .

Mixed derivative in equation and derivative in boundary condition at bottom right boundary. $(x, y) = (〈value〉, 〈value〉)$ .

Mixed derivative in equation and derivative in boundary condition at left boundary. $(x, y) = (〈value〉, 〈value〉)$ .

Mixed derivative in equation and derivative in boundary condition at right boundary. $(x, y) = (〈value〉, 〈value〉)$ .

Mixed derivative in equation and derivative in boundary condition at top boundary. $(x, y) = (〈value〉, 〈value〉)$ .

Mixed derivative in equation and derivative in boundary condition at top left boundary. $(x, y) = (〈value〉, 〈value〉)$ .

Mixed derivative in equation and derivative in boundary condition at top right boundary. $(x, y) = (〈value〉, 〈value〉)$ .

$ifail = 3$: Null boundary condition at bottom boundary, left end. $(x, y) = (〈value〉, 〈value〉)$ .

Null boundary condition at bottom boundary, right end. $(x, y) = (〈value〉, 〈value〉)$ .

Null boundary condition at bottom boundary. $(x, y) = (〈value〉, 〈value〉)$ .

Null boundary condition at left boundary, bottom end. $(x, y) = (〈value〉, 〈value〉)$ .

Null boundary condition at left boundary, top end. $(x, y) = (〈value〉, 〈value〉)$ .

Null boundary condition at left boundary. $(x, y) = (〈value〉, 〈value〉)$ .

Null boundary condition at right boundary, bottom end. $(x, y) = (〈value〉, 〈value〉)$ .

Null boundary condition at right boundary, top end. $(x, y) = (〈value〉, 〈value〉)$ .

Null boundary condition at right boundary. $(x, y) = (〈value〉, 〈value〉)$ .

Null boundary condition at top boundary, left end. $(x, y) = (〈value〉, 〈value〉)$ .

Null boundary condition at top boundary, right end. $(x, y) = (〈value〉, 〈value〉)$ .

Null boundary condition at top boundary. $(x, y) = (〈value〉, 〈value〉)$ .

$ifail = 4$: Equation not elliptic at some point.

$ifail = 5$: There is no unique solution with Neumann Boundary conditions.

$ifail = 6$: The linear equations were not diagonally dominant.

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

Not applicable.

8 Parallelism and Performance

d03eef is not threaded in any implementation.

9 Further Comments

If this routine is used as a preprocessor to the multigrid routine d03edf it should be noted that the rate of convergence of that routine is strongly dependent upon the number of levels in the multigrid scheme, and thus the choice of ngx and ngy is very important.

10 Example

The program solves the elliptic partial differential equation

\frac{\partial^{2} U}{\partial x^{2}} + \frac{\partial^{2} U}{\partial y^{2}} + 50 \{\frac{\partial U}{\partial x} + \frac{\partial U}{\partial y}\} = f (x, y)

on the unit square

0 \leq x

y \leq 1

, with boundary conditions

$\frac{\partial U}{\partial n}$ given on $x = 0$ and $y = 0$ ,
$U$ given on $x = 1$ and $y = 1$ .

The function

f (x, y)

and the exact form of the boundary conditions are derived from the exact solution

U (x, y) = \sin x \sin y

The equation is first solved using central differences. Since the coefficients of the first derivatives are large, the linear equations are not diagonally dominated, and convergence is slow. The equation is solved a second time with upwind differences, showing that convergence is more rapid, but the solution is less accurate.

Interfaces: FL CL AD

NAG FL Interface Introduction

D03 (Pde) Chapter Contents

D03 (Pde) Chapter Introduction

d03ee: FL CL AD

NAG FL Interfaced03eef (dim2_​ellip_​discret)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG FL Interface
d03eef (dim2_ellip_discret)