d03faf solves the Helmholtz equation in Cartesian coordinates in three dimensions using the standard seven-point finite difference approximation. This routine is designed to be particularly efficient on vector processors.

2 Specification

Fortran Interface

Subroutine d03faf (

xs, xf, l, lbdcnd, bdxs, bdxf, ys, yf, m, mbdcnd, bdys, bdyf, zs, zf, n, nbdcnd, bdzs, bdzf, lambda, ldf, ldf2, f, pertrb, w, lwrk, ifail)

Integer, Intent (In)	::	l, lbdcnd, m, mbdcnd, n, nbdcnd, ldf, ldf2, lwrk
Integer, Intent (Inout)	::	ifail
Real (Kind=nag_wp), Intent (In)	::	xs, xf, bdxs(ldf2,n+1), bdxf(ldf2,n+1), ys, yf, bdys(ldf,n+1), bdyf(ldf,n+1), zs, zf, bdzs(ldf,m+1), bdzf(ldf,m+1), lambda
Real (Kind=nag_wp), Intent (Inout)	::	f(ldf,ldf2,n+1)
Real (Kind=nag_wp), Intent (Out)	::	pertrb, w(lwrk)

C Header Interface

#include <nag.h>

void

d03faf_ (const double *xs, const double *xf, const Integer *l, const Integer *lbdcnd, const double bdxs[], const double bdxf[], const double *ys, const double *yf, const Integer *m, const Integer *mbdcnd, const double bdys[], const double bdyf[], const double *zs, const double *zf, const Integer *n, const Integer *nbdcnd, const double bdzs[], const double bdzf[], const double *lambda, const Integer *ldf, const Integer *ldf2, double f[], double *pertrb, double w[], const Integer *lwrk, Integer *ifail)

The routine may be called by the names d03faf or nagf_pde_dim3_ellip_helmholtz.

3 Description

d03faf solves the three-dimensional Helmholtz equation in Cartesian coordinates:

\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}} + \frac{\partial^{2} u}{\partial z^{2}} + λ u = f (x, y, z) .

This subroutine forms the system of linear equations resulting from the standard seven-point finite difference equations, and then solves the system using a method based on the fast Fourier transform (FFT) described by Swarztrauber (1984). This subroutine is based on the routine HW3CRT from FISHPACK (see Swarztrauber and Sweet (1979)).

More precisely, the routine replaces all the second derivatives by second-order central difference approximations, resulting in a block tridiagonal system of linear equations. The equations are modified to allow for the prescribed boundary conditions. Either the solution or the derivative of the solution may be specified on any of the boundaries, or the solution may be specified to be periodic in any of the three dimensions. By taking the discrete Fourier transform in the

x

- and

y

-directions, the equations are reduced to sets of tridiagonal systems of equations. The Fourier transforms required are computed using the multiple FFT routines found in Chapter C06.

4 References

Swarztrauber P N (1984) Fast Poisson solvers Studies in Numerical Analysis (ed G H Golub) Mathematical Association of America

Swarztrauber P N and Sweet R A (1979) Efficient Fortran subprograms for the solution of separable elliptic partial differential equations ACM Trans. Math. Software 5 352–364

5 Arguments

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

On entry: the lower bound of the range of

x

, i.e.,

xs \leq x \leq xf

Constraint:

xs < xf

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

On entry: the upper bound of the range of

x

, i.e.,

xs \leq x \leq xf

Constraint:

xs < xf

3: $l$ – Integer Input

On entry: the number of panels into which the interval

(xs, xf)

is subdivided. Hence, there will be

l + 1

grid points in the

x

-direction given by

x_{i} = xs + (i - 1) \times δ x

, for

i = 1, 2, \dots, l + 1

, where

δ x = (xf - xs) / l

is the panel width.

Constraint:

l \geq 5

4: $lbdcnd$ – Integer Input

On entry: indicates the type of boundary conditions at

x = xs

and

x = xf

$lbdcnd = 0$: If the solution is periodic in $x$ , i.e., $u (xs, y, z) = u (xf, y, z)$ .
$lbdcnd = 1$: If the solution is specified at $x = xs$ and $x = xf$ .
$lbdcnd = 2$: If the solution is specified at $x = xs$ and the derivative of the solution with respect to $x$ is specified at $x = xf$ .
$lbdcnd = 3$: If the derivative of the solution with respect to $x$ is specified at $x = xs$ and $x = xf$ .
$lbdcnd = 4$: If the derivative of the solution with respect to $x$ is specified at $x = xs$ and the solution is specified at $x = xf$ .

Constraint:

0 \leq lbdcnd \leq 4

5: $bdxs (ldf2, n + 1)$ – Real (Kind=nag_wp) array Input

On entry: the values of the derivative of the solution with respect to

x

x = xs

. When

lbdcnd = 3

4

bdxs (j, k) = u_{x} (xs, y_{j}, z_{k})

, for

j = 1, 2, \dots, m + 1

and

k = 1, 2, \dots, n + 1

When lbdcnd has any other value, bdxs is not referenced.

6: $bdxf (ldf2, n + 1)$ – Real (Kind=nag_wp) array Input

On entry: the values of the derivative of the solution with respect to

x

x = xf

. When

lbdcnd = 2

3

bdxf (j, k) = u_{x} (xf, y_{j}, z_{k})

, for

j = 1, 2, \dots, m + 1

and

k = 1, 2, \dots, n + 1

When lbdcnd has any other value, bdxf is not referenced.

7: $ys$ – Real (Kind=nag_wp) Input

On entry: the lower bound of the range of

y

, i.e.,

ys \leq y \leq yf

Constraint:

ys < yf

8: $yf$ – Real (Kind=nag_wp) Input

On entry: the upper bound of the range of

y

, i.e.,

ys \leq y \leq yf

Constraint:

ys < yf

9: $m$ – Integer Input

On entry: the number of panels into which the interval

(ys, yf)

is subdivided. Hence, there will be

m + 1

grid points in the

y

-direction given by

y_{j} = ys + (j - 1) \times δ y

, for

j = 1, 2, \dots, m + 1

, where

δ y = (yf - ys) / m

is the panel width.

Constraint:

m \geq 5

10: $mbdcnd$ – Integer Input

On entry: indicates the type of boundary conditions at

y = ys

and

y = yf

$mbdcnd = 0$: If the solution is periodic in $y$ , i.e., $u (x, yf, z) = u (x, ys, z)$ .
$mbdcnd = 1$: If the solution is specified at $y = ys$ and $y = yf$ .
$mbdcnd = 2$: If the solution is specified at $y = ys$ and the derivative of the solution with respect to $y$ is specified at $y = yf$ .
$mbdcnd = 3$: If the derivative of the solution with respect to $y$ is specified at $y = ys$ and $y = yf$ .
$mbdcnd = 4$: If the derivative of the solution with respect to $y$ is specified at $y = ys$ and the solution is specified at $y = yf$ .

Constraint:

0 \leq mbdcnd \leq 4

11: $bdys (ldf, n + 1)$ – Real (Kind=nag_wp) array Input

On entry: the values of the derivative of the solution with respect to

y

y = ys

. When

mbdcnd = 3

4

bdys (i, k) = u_{y} (x_{i}, ys, z_{k})

, for

i = 1, 2, \dots, l + 1

and

k = 1, 2, \dots, n + 1

When mbdcnd has any other value, bdys is not referenced.

12: $bdyf (ldf, n + 1)$ – Real (Kind=nag_wp) array Input

On entry: the values of the derivative of the solution with respect to

y

y = yf

. When

mbdcnd = 2

3

bdyf (i, k) = u_{y} (x_{i}, yf, z_{k})

, for

i = 1, 2, \dots, l + 1

and

k = 1, 2, \dots, n + 1

When mbdcnd has any other value, bdyf is not referenced.

13: $zs$ – Real (Kind=nag_wp) Input

On entry: the lower bound of the range of

z

, i.e.,

zs \leq z \leq zf

Constraint:

zs < zf

14: $zf$ – Real (Kind=nag_wp) Input

On entry: the upper bound of the range of

z

, i.e.,

zs \leq z \leq zf

Constraint:

zs < zf

15: $n$ – Integer Input

On entry: the number of panels into which the interval

(zs, zf)

is subdivided. Hence, there will be

n + 1

grid points in the

z

-direction given by

z_{k} = zs + (k - 1) \times δ z

, for

k = 1, 2, \dots, n + 1

, where

δ z = (zf - zs) / n

is the panel width.

Constraint:

n \geq 5

16: $nbdcnd$ – Integer Input

On entry: specifies the type of boundary conditions at

z = zs

and

z = zf

$nbdcnd = 0$: if the solution is periodic in $z$ , i.e., $u (x, y, zf) = u (x, y, zs)$ .
$nbdcnd = 1$: if the solution is specified at $z = zs$ and $z = zf$ .
$nbdcnd = 2$: if the solution is specified at $z = zs$ and the derivative of the solution with respect to $z$ is specified at $z = zf$ .
$nbdcnd = 3$: if the derivative of the solution with respect to $z$ is specified at $z = zs$ and $z = zf$ .
$nbdcnd = 4$: if the derivative of the solution with respect to $z$ is specified at $z = zs$ and the solution is specified at $z = zf$ .

Constraint:

0 \leq nbdcnd \leq 4

17: $bdzs (ldf, m + 1)$ – Real (Kind=nag_wp) array Input

On entry: the values of the derivative of the solution with respect to

z

z = zs

. When

nbdcnd = 3

4

bdzs (i, j) = u_{z} (x_{i}, y_{j}, zs)

, for

i = 1, 2, \dots, l + 1

and

j = 1, 2, \dots, m + 1

When nbdcnd has any other value, bdzs is not referenced.

18: $bdzf (ldf, m + 1)$ – Real (Kind=nag_wp) array Input

On entry: the values of the derivative of the solution with respect to

z

z = zf

. When

nbdcnd = 2

3

bdzf (i, j) = u_{z} (x_{i}, y_{j}, zf)

, for

i = 1, 2, \dots, l + 1

and

j = 1, 2, \dots, m + 1

When nbdcnd has any other value, bdzf is not referenced.

19: $lambda$ – Real (Kind=nag_wp) Input

On entry: the constant

λ

in the Helmholtz equation. For certain positive values of

λ

a solution to the differential equation may not exist, and close to these values the solution of the discretized problem will be extremely ill-conditioned. If

λ > 0

, d03faf will set

ifail = 3

, but will still attempt to find a solution. However, since in general the values of

λ

for which no solution exists cannot be predicted a priori, you are advised to treat any results computed with

λ > 0

with great caution.

20: $ldf$ – Integer Input

On entry: the first dimension of the arrays f, bdys, bdyf, bdzs and bdzf as declared in the (sub)program from which d03faf is called.

Constraint:

ldf \geq l + 1

21: $ldf2$ – Integer Input

On entry: the second dimension of the array f and the first dimension of the arrays bdxs and bdxf as declared in the (sub)program from which d03faf is called.

Constraint:

ldf2 \geq m + 1

22: $f (ldf, ldf2, n + 1)$ – Real (Kind=nag_wp) array Input/Output

On entry: the values of the right-side of the Helmholtz equation and boundary values (if any).

f (i, j, k) = f (x_{i}, y_{j}, z_{k}) i = 2, 3, \dots, l, j = 2, 3, \dots, m ​ and ​ k = 2, 3, \dots, n .

On the boundaries f is defined by

lbdcnd	$f (1, j, k)$	$f (l + 1, j, k)$
0	$f (xs, y_{j}, z_{k})$	$f (xs, y_{j}, z_{k})$
1	$u (xs, y_{j}, z_{k})$	$u (xf, y_{j}, z_{k})$
2	$u (xs, y_{j}, z_{k})$	$f (xf, y_{j}, z_{k})$	$j = 1, 2, \dots, m + 1$
3	$f (xs, y_{j}, z_{k})$	$f (xf, y_{j}, z_{k})$	$k = 1, 2, \dots, n + 1$
4	$f (xs, y_{j}, z_{k})$	$u (xf, y_{j}, z_{k})$

mbdcnd	$f (i, 1, k)$	$f (i, m + 1, k)$
0	$f (x_{i}, ys, z_{k})$	$f (x_{i}, ys, z_{k})$
1	$u (x_{i}, ys, z_{k})$	$u (x_{i}, yf, z_{k})$
2	$u (x_{i}, ys, z_{k})$	$f (x_{i}, yf, z_{k})$	$i = 1, 2, \dots, l + 1$
3	$f (x_{i}, ys, z_{k})$	$f (x_{i}, yf, z_{k})$	$k = 1, 2, \dots, n + 1$
4	$f (x_{i}, ys, z_{k})$	$u (x_{i}, yf, z_{k})$

nbdcnd	$f (i, j, 1)$	$f (i, j, n + 1)$
0	$f (x_{i}, y_{j}, zs)$	$f (x_{i}, y_{j}, zs)$
1	$u (x_{i}, y_{j}, zs)$	$u (x_{i}, y_{j}, zf)$
2	$u (x_{i}, y_{j}, zs)$	$f (x_{i}, y_{j}, zf)$	$i = 1, 2, \dots, l + 1$
3	$f (x_{i}, y_{j}, zs)$	$f (x_{i}, y_{j}, zf)$	$j = 1, 2, \dots, m + 1$
4	$f (x_{i}, y_{j}, zs)$	$u (x_{i}, y_{j}, zf)$

Note: if the table calls for both the solution

u

and the right-hand side

f

on a boundary, the solution must be specified.

On exit: contains the solution

u (i, j, k)

of the finite difference approximation for the grid point

(x_{i}, y_{j}, z_{k})

, for

i = 1, 2, \dots, l + 1

j = 1, 2, \dots, m + 1

and

k = 1, 2, \dots, n + 1

23: $pertrb$ – Real (Kind=nag_wp) Output

On exit:

pertrb = 0

, unless a solution to Poisson's equation

(λ = 0)

is required with a combination of periodic or derivative boundary conditions (lbdcnd, mbdcnd and

nbdcnd = 0

3

). In this case a solution may not exist. pertrb is a constant, calculated and subtracted from the array f, which ensures that a solution exists. d03faf then computes this solution, which is a least squares solution to the original approximation. This solution is not unique and is unnormalized. The value of pertrb should be small compared to the right-hand side f, otherwise a solution has been obtained to an essentially different problem. This comparison should always be made to ensure that a meaningful solution has been obtained.

24: $w (lwrk)$ – Real (Kind=nag_wp) array Output

25: $lwrk$ – Integer Input

On entry: the dimension of the array w as declared in the (sub)program from which d03faf is called. w is no longer used as workspace, lwrk can, therefore, be safely set to zero.

26: $ifail$ – Integer Input/Output

On entry: ifail must be set to

0

−1

1

to set behaviour on detection of an error; these values have no effect when no error is detected.

A value of

0

causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of

−1

means that an error message is printed while a value of

1

means that it is not.

If halting is not appropriate, the value

−1

1

is recommended. If message printing is undesirable, then the value

1

is recommended. Otherwise, the value

0

is recommended. 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:

$ifail = 1$: On entry, $l = ⟨ value ⟩$ .
Constraint: $l \geq 5$ .

On entry, $lbdcnd = ⟨ value ⟩$ .
Constraint: $lbdcnd \geq 0$ and $lbdcnd \leq 4$ .

On entry, $ldf2 = ⟨ value ⟩$ and $m = ⟨ value ⟩$ .
Constraint: $ldf2 \geq m + 1$ .

On entry, $ldf = ⟨ value ⟩$ and $l = ⟨ value ⟩$ .
Constraint: $ldf \geq l + 1$ .

On entry, $m = ⟨ value ⟩$ .
Constraint: $m \geq 5$ .

On entry, $mbdcnd = ⟨ value ⟩$ .
Constraint: $mbdcnd \geq 0$ and $mbdcnd \leq 4$ .

On entry, $n = ⟨ value ⟩$ .
Constraint: $n \geq 5$ .

On entry, $nbdcnd = ⟨ value ⟩$ .
Constraint: $nbdcnd \geq 0$ and $nbdcnd \leq 4$ .

On entry, $xs = ⟨ value ⟩$ and $xf = ⟨ value ⟩$ .
Constraint: $xs < xf$ .

On entry, $ys = ⟨ value ⟩$ and $yf = ⟨ value ⟩$ .
Constraint: $ys < yf$

On entry, $zs = ⟨ value ⟩$ and $zf = ⟨ value ⟩$ .
Constraint: $zs < zf$ .

$ifail = 3$: $λ > 0.0$ in Helmholtz's equation – a solution may not exist.

$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

d03faf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.

d03faf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.

Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

The execution time is roughly proportional to

l \times m \times n \times (\log_{2} l + \log_{2} m + 5)

, but also depends on input arguments lbdcnd and mbdcnd.

10 Example

This example solves the Helmholz equation

\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}} + \frac{\partial^{2} u}{\partial z^{2}} + λ u = f (x, y, z)

for

(x, y, z) \in [0, 1] \times [0, 2 π] \times [0, \frac{π}{2}]

, where

λ = −2

, and

f (x, y, z)

is derived from the exact solution

u (x, y, z) = x^{4} \sin y \cos z .

The equation is subject to the following boundary conditions, again derived from the exact solution given above.

$u (0, y, z)$ and $u (1, y, z)$ are prescribed (i.e., $lbdcnd = 1$ ).
$u (x, 0, z) = u (x, 2 π, z)$ (i.e., $mbdcnd = 0$ ).
$u (x, y, 0)$ and $u_{z} (x, y, \frac{π}{2})$ are prescribed (i.e., $nbdcnd = 2$ ).

d03fa: FL CL CPP AD

NAG FL Interfaced03faf (dim3_​ellip_​helmholtz)

▸▿ 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
d03faf (dim3_ellip_helmholtz)