Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_pde_3d_ellip_helmholtz (d03fa)

## Purpose

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

## Syntax

[f, pertrb, ifail] = d03fa(xs, xf, l, lbdcnd, bdxs, bdxf, ys, yf, m, mbdcnd, bdys, bdyf, zs, zf, n, nbdcnd, bdzs, bdzf, lambda, f, 'lwrk', lwrk)
[f, pertrb, ifail] = nag_pde_3d_ellip_helmholtz(xs, xf, l, lbdcnd, bdxs, bdxf, ys, yf, m, mbdcnd, bdys, bdyf, zs, zf, n, nbdcnd, bdzs, bdzf, lambda, f, 'lwrk', lwrk)
Note: the interface to this routine has changed since earlier releases of the toolbox:
Mark 23: lwrk added as optional parameter
.

## Description

nag_pde_3d_ellip_helmholtz (d03fa) solves the three-dimensional Helmholtz equation in Cartesian coordinates:
 (∂2u)/( ∂ x2) + (∂2u)/( ∂ y2) + (∂2u)/( ∂ z2) + λu = f(x,y,z). $∂2u ∂x2 + ∂2u ∂y2 + ∂2u ∂z2 +λu=f(x,y,z).$
This function 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 function is based on the function HW3CRT from FISHPACK (see Swarztrauber and Sweet (1979)).
More precisely, the function 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$x$- and y$y$-directions, the equations are reduced to sets of tridiagonal systems of equations. The Fourier transforms required are computed using the multiple FFT functions found in Chapter C06.

## 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

## Parameters

### Compulsory Input Parameters

1:     xs – double scalar
The lower bound of the range of x$x$, i.e., xsxxf${\mathbf{xs}}\le x\le {\mathbf{xf}}$.
Constraint: xs < xf${\mathbf{xs}}<{\mathbf{xf}}$.
2:     xf – double scalar
The upper bound of the range of x$x$, i.e., xsxxf${\mathbf{xs}}\le x\le {\mathbf{xf}}$.
Constraint: xs < xf${\mathbf{xs}}<{\mathbf{xf}}$.
3:     l – int64int32nag_int scalar
The number of panels into which the interval (xs,xf) is subdivided. Hence, there will be l + 1${\mathbf{l}}+1$ grid points in the x$x$-direction given by xi = xs + (i1) × δx${x}_{\mathit{i}}={\mathbf{xs}}+\left(\mathit{i}-1\right)×\delta x$, for i = 1,2,,l + 1$\mathit{i}=1,2,\dots ,{\mathbf{l}}+1$, where δx = (xfxs) / l$\delta x=\left({\mathbf{xf}}-{\mathbf{xs}}\right)/{\mathbf{l}}$ is the panel width.
Constraint: l5${\mathbf{l}}\ge 5$.
4:     lbdcnd – int64int32nag_int scalar
Indicates the type of boundary conditions at x = xs$x={\mathbf{xs}}$ and x = xf$x={\mathbf{xf}}$.
lbdcnd = 0${\mathbf{lbdcnd}}=0$
If the solution is periodic in x$x$, i.e., u(xs,y,z) = u(xf,y,z)$u\left({\mathbf{xs}},y,z\right)=u\left({\mathbf{xf}},y,z\right)$.
lbdcnd = 1${\mathbf{lbdcnd}}=1$
If the solution is specified at x = xs$x={\mathbf{xs}}$ and x = xf$x={\mathbf{xf}}$.
lbdcnd = 2${\mathbf{lbdcnd}}=2$
If the solution is specified at x = xs$x={\mathbf{xs}}$ and the derivative of the solution with respect to x$x$ is specified at x = xf$x={\mathbf{xf}}$.
lbdcnd = 3${\mathbf{lbdcnd}}=3$
If the derivative of the solution with respect to x$x$ is specified at x = xs$x={\mathbf{xs}}$ and x = xf$x={\mathbf{xf}}$.
lbdcnd = 4${\mathbf{lbdcnd}}=4$
If the derivative of the solution with respect to x$x$ is specified at x = xs$x={\mathbf{xs}}$ and the solution is specified at x = xf$x={\mathbf{xf}}$.
Constraint: 0lbdcnd4$0\le {\mathbf{lbdcnd}}\le 4$.
5:     bdxs(ldf2,n + 1${\mathbf{n}}+1$) – double array
ldf2, the first dimension of the array, must satisfy the constraint ldf2m + 1$\mathit{ldf2}\ge {\mathbf{m}}+1$.
The values of the derivative of the solution with respect to x$x$ at x = xs$x={\mathbf{xs}}$. When lbdcnd = 3${\mathbf{lbdcnd}}=3$ or 4$4$, bdxs(j,k) = ux(xs,yj,zk)${\mathbf{bdxs}}\left(\mathit{j},\mathit{k}\right)={u}_{x}\left({\mathbf{xs}},{y}_{j},{z}_{k}\right)$, for j = 1,2,,m + 1$\mathit{j}=1,2,\dots ,{\mathbf{m}}+1$ and k = 1,2,,n + 1$\mathit{k}=1,2,\dots ,{\mathbf{n}}+1$.
When lbdcnd has any other value, bdxs is not referenced.
6:     bdxf(ldf2,n + 1${\mathbf{n}}+1$) – double array
ldf2, the first dimension of the array, must satisfy the constraint ldf2m + 1$\mathit{ldf2}\ge {\mathbf{m}}+1$.
The values of the derivative of the solution with respect to x$x$ at x = xf$x={\mathbf{xf}}$. When lbdcnd = 2${\mathbf{lbdcnd}}=2$ or 3$3$, bdxf(j,k) = ux(xf,yj,zk)${\mathbf{bdxf}}\left(\mathit{j},\mathit{k}\right)={u}_{x}\left({\mathbf{xf}},{y}_{\mathit{j}},{z}_{\mathit{k}}\right)$, for j = 1,2,,m + 1$\mathit{j}=1,2,\dots ,{\mathbf{m}}+1$ and k = 1,2,,n + 1$\mathit{k}=1,2,\dots ,{\mathbf{n}}+1$.
When lbdcnd has any other value, bdxf is not referenced.
7:     ys – double scalar
The lower bound of the range of y$y$, i.e., ysyyf${\mathbf{ys}}\le y\le {\mathbf{yf}}$.
Constraint: ys < yf${\mathbf{ys}}<{\mathbf{yf}}$.
8:     yf – double scalar
The upper bound of the range of y$y$, i.e., ysyyf${\mathbf{ys}}\le y\le {\mathbf{yf}}$.
Constraint: ys < yf${\mathbf{ys}}<{\mathbf{yf}}$.
9:     m – int64int32nag_int scalar
The number of panels into which the interval (ys,yf) is subdivided. Hence, there will be m + 1${\mathbf{m}}+1$ grid points in the y$y$-direction given by yj = ys + (j1) × δy${y}_{\mathit{j}}={\mathbf{ys}}+\left(\mathit{j}-1\right)×\delta y$, for j = 1,2,,m + 1$\mathit{j}=1,2,\dots ,{\mathbf{m}}+1$, where δy = (yfys) / m$\delta y=\left({\mathbf{yf}}-{\mathbf{ys}}\right)/{\mathbf{m}}$ is the panel width.
Constraint: m5${\mathbf{m}}\ge 5$.
10:   mbdcnd – int64int32nag_int scalar
Indicates the type of boundary conditions at y = ys$y={\mathbf{ys}}$ and y = yf$y={\mathbf{yf}}$.
mbdcnd = 0${\mathbf{mbdcnd}}=0$
If the solution is periodic in y$y$, i.e., u(x,yf,z) = u(x,ys,z)$u\left(x,{\mathbf{yf}},z\right)=u\left(x,{\mathbf{ys}},z\right)$.
mbdcnd = 1${\mathbf{mbdcnd}}=1$
If the solution is specified at y = ys$y={\mathbf{ys}}$ and y = yf$y={\mathbf{yf}}$.
mbdcnd = 2${\mathbf{mbdcnd}}=2$
If the solution is specified at y = ys$y={\mathbf{ys}}$ and the derivative of the solution with respect to y$y$ is specified at y = yf$y={\mathbf{yf}}$.
mbdcnd = 3${\mathbf{mbdcnd}}=3$
If the derivative of the solution with respect to y$y$ is specified at y = ys$y={\mathbf{ys}}$ and y = yf$y={\mathbf{yf}}$.
mbdcnd = 4${\mathbf{mbdcnd}}=4$
If the derivative of the solution with respect to y$y$ is specified at y = ys$y={\mathbf{ys}}$ and the solution is specified at y = yf$y={\mathbf{yf}}$.
Constraint: 0mbdcnd4$0\le {\mathbf{mbdcnd}}\le 4$.
11:   bdys(ldf,n + 1${\mathbf{n}}+1$) – double array
ldf, the first dimension of the array, must satisfy the constraint ldfl + 1$\mathit{ldf}\ge {\mathbf{l}}+1$.
The values of the derivative of the solution with respect to y$y$ at y = ys$y={\mathbf{ys}}$. When mbdcnd = 3${\mathbf{mbdcnd}}=3$ or 4$4$, bdys(i,k) = uy(xi,ys,zk)${\mathbf{bdys}}\left(\mathit{i},\mathit{k}\right)={u}_{y}\left({x}_{\mathit{i}},{\mathbf{ys}},{z}_{\mathit{k}}\right)$, for i = 1,2,,l + 1$\mathit{i}=1,2,\dots ,{\mathbf{l}}+1$ and k = 1,2,,n + 1$\mathit{k}=1,2,\dots ,{\mathbf{n}}+1$.
When mbdcnd has any other value, bdys is not referenced.
12:   bdyf(ldf,n + 1${\mathbf{n}}+1$) – double array
ldf, the first dimension of the array, must satisfy the constraint ldfl + 1$\mathit{ldf}\ge {\mathbf{l}}+1$.
The values of the derivative of the solution with respect to y$y$ at y = yf$y={\mathbf{yf}}$. When mbdcnd = 2${\mathbf{mbdcnd}}=2$ or 3$3$, bdyf(i,k) = uy(xi,yf,zk)${\mathbf{bdyf}}\left(\mathit{i},\mathit{k}\right)={u}_{y}\left({x}_{\mathit{i}},{\mathbf{yf}},{z}_{\mathit{k}}\right)$, for i = 1,2,,l + 1$\mathit{i}=1,2,\dots ,{\mathbf{l}}+1$ and k = 1,2,,n + 1$\mathit{k}=1,2,\dots ,{\mathbf{n}}+1$.
When mbdcnd has any other value, bdyf is not referenced.
13:   zs – double scalar
The lower bound of the range of z$z$, i.e., zszzf${\mathbf{zs}}\le z\le {\mathbf{zf}}$.
Constraint: zs < zf${\mathbf{zs}}<{\mathbf{zf}}$.
14:   zf – double scalar
The upper bound of the range of z$z$, i.e., zszzf${\mathbf{zs}}\le z\le {\mathbf{zf}}$.
Constraint: zs < zf${\mathbf{zs}}<{\mathbf{zf}}$.
15:   n – int64int32nag_int scalar
The number of panels into which the interval (zs,zf) is subdivided. Hence, there will be n + 1${\mathbf{n}}+1$ grid points in the z$z$-direction given by zk = zs + (k1) × δz${z}_{\mathit{k}}={\mathbf{zs}}+\left(\mathit{k}-1\right)×\delta z$, for k = 1,2,,n + 1$\mathit{k}=1,2,\dots ,{\mathbf{n}}+1$, where δz = (zfzs) / n$\delta z=\left({\mathbf{zf}}-{\mathbf{zs}}\right)/{\mathbf{n}}$ is the panel width.
Constraint: n5${\mathbf{n}}\ge 5$.
16:   nbdcnd – int64int32nag_int scalar
Specifies the type of boundary conditions at z = zs$z={\mathbf{zs}}$ and z = zf$z={\mathbf{zf}}$.
nbdcnd = 0${\mathbf{nbdcnd}}=0$
if the solution is periodic in z$z$, i.e., u(x,y,zf) = u(x,y,zs)$u\left(x,y,{\mathbf{zf}}\right)=u\left(x,y,{\mathbf{zs}}\right)$.
nbdcnd = 1${\mathbf{nbdcnd}}=1$
if the solution is specified at z = zs$z={\mathbf{zs}}$ and z = zf$z={\mathbf{zf}}$.
nbdcnd = 2${\mathbf{nbdcnd}}=2$
if the solution is specified at z = zs$z={\mathbf{zs}}$ and the derivative of the solution with respect to z$z$ is specified at z = zf$z={\mathbf{zf}}$.
nbdcnd = 3${\mathbf{nbdcnd}}=3$
if the derivative of the solution with respect to z$z$ is specified at z = zs$z={\mathbf{zs}}$ and z = zf$z={\mathbf{zf}}$.
nbdcnd = 4${\mathbf{nbdcnd}}=4$
if the derivative of the solution with respect to z$z$ is specified at z = zs$z={\mathbf{zs}}$ and the solution is specified at z = zf$z={\mathbf{zf}}$.
Constraint: 0nbdcnd4$0\le {\mathbf{nbdcnd}}\le 4$.
17:   bdzs(ldf,m + 1${\mathbf{m}}+1$) – double array
ldf, the first dimension of the array, must satisfy the constraint ldfl + 1$\mathit{ldf}\ge {\mathbf{l}}+1$.
The values of the derivative of the solution with respect to z$z$ at z = zs$z={\mathbf{zs}}$. When nbdcnd = 3${\mathbf{nbdcnd}}=3$ or 4$4$, bdzs(i,j) = uz(xi,yj,zs)${\mathbf{bdzs}}\left(\mathit{i},\mathit{j}\right)={u}_{z}\left({x}_{i},{y}_{j},{\mathbf{zs}}\right)$, for i = 1,2,,l + 1$\mathit{i}=1,2,\dots ,{\mathbf{l}}+1$ and j = 1,2,,m + 1$\mathit{j}=1,2,\dots ,{\mathbf{m}}+1$.
When nbdcnd has any other value, bdzs is not referenced.
18:   bdzf(ldf,m + 1${\mathbf{m}}+1$) – double array
ldf, the first dimension of the array, must satisfy the constraint ldfl + 1$\mathit{ldf}\ge {\mathbf{l}}+1$.
The values of the derivative of the solution with respect to z$z$ at z = zf$z={\mathbf{zf}}$. When nbdcnd = 2${\mathbf{nbdcnd}}=2$ or 3$3$, bdzf(i,j) = uz(xi,yj,zf)${\mathbf{bdzf}}\left(\mathit{i},\mathit{j}\right)={u}_{z}\left({x}_{i},{y}_{j},{\mathbf{zf}}\right)$, for i = 1,2,,l + 1$\mathit{i}=1,2,\dots ,{\mathbf{l}}+1$ and j = 1,2,,m + 1$\mathit{j}=1,2,\dots ,{\mathbf{m}}+1$.
When nbdcnd has any other value, bdzf is not referenced.
19:   lambda – double scalar
The constant λ$\lambda$ in the Helmholtz equation. For certain positive values of λ$\lambda$ 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$\lambda >0$, then nag_pde_3d_ellip_helmholtz (d03fa) will set ${\mathbf{ifail}}={\mathbf{3}}$, but will still attempt to find a solution. However, since in general the values of λ$\lambda$ for which no solution exists cannot be predicted a priori, you are advised to treat any results computed with λ > 0$\lambda >0$ with great caution.
20:   f(ldf,ldf2,n + 1${\mathbf{n}}+1$) – double array
ldf, the first dimension of the array, must satisfy the constraint ldfl + 1$\mathit{ldf}\ge {\mathbf{l}}+1$.
The values of the right-side of the Helmholtz equation and boundary values (if any).
 f(i,j,k) = f(xi,yj,zk) i = 2,3, … ,l, j = 2,3, … ,m ​ and ​ k = 2,3, … ,n. $f(i,j,k) = f(xi,yj,zk) i=2,3,…,l, j=2,3,…,m ​ and ​ k=2,3,…,n.$
On the boundaries f is defined by
 lbdcnd f(1,j,k)${\mathbf{f}}\left(1,j,k\right)$ f(l + 1,j,k)${\mathbf{f}}\left({\mathbf{l}}+1,j,k\right)$ 0 f(xs,yj,zk)$f\left({\mathbf{xs}},{y}_{j},{z}_{k}\right)$ f(xs,yj,zk)$f\left({\mathbf{xs}},{y}_{j},{z}_{k}\right)$ 1 u(xs,yj,zk)$u\left({\mathbf{xs}},{y}_{j},{z}_{k}\right)$ u(xf,yj,zk)$u\left({\mathbf{xf}},{y}_{j},{z}_{k}\right)$ 2 u(xs,yj,zk)$u\left({\mathbf{xs}},{y}_{j},{z}_{k}\right)$ f(xf,yj,zk)$f\left({\mathbf{xf}},{y}_{j},{z}_{k}\right)$ j = 1,2, … ,m + 1$j=1,2,\dots ,{\mathbf{m}}+1$ 3 f(xs,yj,zk)$f\left({\mathbf{xs}},{y}_{j},{z}_{k}\right)$ f(xf,yj,zk)$f\left({\mathbf{xf}},{y}_{j},{z}_{k}\right)$ k = 1,2, … ,n + 1$k=1,2,\dots ,{\mathbf{n}}+1$ 4 f(xs,yj,zk)$f\left({\mathbf{xs}},{y}_{j},{z}_{k}\right)$ u(xf,yj,zk)$u\left({\mathbf{xf}},{y}_{j},{z}_{k}\right)$ mbdcnd f(i,1,k)${\mathbf{f}}\left(i,1,k\right)$ f(i,m + 1,k)${\mathbf{f}}\left(i,{\mathbf{m}}+1,k\right)$ 0 f(xi,ys,zk)$f\left({x}_{i},{\mathbf{ys}},{z}_{k}\right)$ f(xi,ys,zk)$f\left({x}_{i},{\mathbf{ys}},{z}_{k}\right)$ 1 u(ys,xi,zk)$u\left({\mathbf{ys}},{x}_{i},{z}_{k}\right)$ u(yf,xi,zk)$u\left({\mathbf{yf}},{x}_{i},{z}_{k}\right)$ 2 u(xi,ys,zk)$u\left({x}_{i},{\mathbf{ys}},{z}_{k}\right)$ f(xi,yf,zk)$f\left({x}_{i},{\mathbf{yf}},{z}_{k}\right)$ i = 1,2, … ,l + 1$i=1,2,\dots ,{\mathbf{l}}+1$ 3 f(xi,ys,zk)$f\left({x}_{i},{\mathbf{ys}},{z}_{k}\right)$ f(xi,yf,zk)$f\left({x}_{i},{\mathbf{yf}},{z}_{k}\right)$ k = 1,2, … ,n + 1$k=1,2,\dots ,{\mathbf{n}}+1$ 4 f(xi,ys,zk)$f\left({x}_{i},{\mathbf{ys}},{z}_{k}\right)$ u(xi,yf,zk)$u\left({x}_{i},{\mathbf{yf}},{z}_{k}\right)$ nbdcnd f(i,j,1)${\mathbf{f}}\left(i,j,1\right)$ f(i,j,n + 1)${\mathbf{f}}\left(i,j,{\mathbf{n}}+1\right)$ 0 f(xi,yj,zs)$f\left({x}_{i},{y}_{j},{\mathbf{zs}}\right)$ f(xi,yj,zs)$f\left({x}_{i},{y}_{j},{\mathbf{zs}}\right)$ 1 u(xi,yj,zs)$u\left({x}_{i},{y}_{j},{\mathbf{zs}}\right)$ u(xi,yj,zf)$u\left({x}_{i},{y}_{j},{\mathbf{zf}}\right)$ 2 u(xi,yj,zs)$u\left({x}_{i},{y}_{j},{\mathbf{zs}}\right)$ f(xi,yj,zf)$f\left({x}_{i},{y}_{j},{\mathbf{zf}}\right)$ i = 1,2, … ,l + 1$i=1,2,\dots ,{\mathbf{l}}+1$ 3 f(xi,yj,zs)$f\left({x}_{i},{y}_{j},{\mathbf{zs}}\right)$ f(xi,yj,zf)$f\left({x}_{i},{y}_{j},{\mathbf{zf}}\right)$ j = 1,2, … ,m + 1$j=1,2,\dots ,{\mathbf{m}}+1$ 4 f(xi,yj,zs)$f\left({x}_{i},{y}_{j},{\mathbf{zs}}\right)$ u(xi,yj,zf)$u\left({x}_{i},{y}_{j},{\mathbf{zf}}\right)$
Note: if the table calls for both the solution u$u$ and the right-hand side f$f$ on a boundary, then the solution must be specified.

### Optional Input Parameters

1:     lwrk – int64int32nag_int scalar
The dimension of the array w as declared in the (sub)program from which nag_pde_3d_ellip_helmholtz (d03fa) is called. 2 × (n + 1) × max (l,m) + 3 × l + 3 × m + 4 × n + 6$2×\left({\mathbf{n}}+1\right)×\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{l}},{\mathbf{m}}\right)+3×{\mathbf{l}}+3×{\mathbf{m}}+4×{\mathbf{n}}+6$ is an upper bound on the required size of w. If lwrk is too small, the function exits with ${\mathbf{ifail}}={\mathbf{2}}$, and if on entry ifail = 0${\mathbf{ifail}}=0$ or 1$-1$, a message is output giving the exact value of lwrk required to solve the current problem.
Default: 2 × (n + 1) × max (l,m) + 3 × l + 3 × m + 4 × n + 6$2×\left({\mathbf{n}}+1\right)×\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{l}},{\mathbf{m}}\right)+3×{\mathbf{l}}+3×{\mathbf{m}}+4×{\mathbf{n}}+6$

ldf ldf2 w

### Output Parameters

1:     f(ldf,ldf2,n + 1${\mathbf{n}}+1$) – double array
ldfl + 1$\mathit{ldf}\ge {\mathbf{l}}+1$.
ldf2m + 1$\mathit{ldf2}\ge {\mathbf{m}}+1$.
Contains the solution u(i,j,k)$u\left(\mathit{i},\mathit{j},\mathit{k}\right)$ of the finite difference approximation for the grid point (xi,yj,zk)$\left({x}_{\mathit{i}},{y}_{\mathit{j}},{z}_{\mathit{k}}\right)$, for i = 1,2,,l + 1$\mathit{i}=1,2,\dots ,{\mathbf{l}}+1$, j = 1,2,,m + 1$\mathit{j}=1,2,\dots ,{\mathbf{m}}+1$ and k = 1,2,,n + 1$\mathit{k}=1,2,\dots ,{\mathbf{n}}+1$.
2:     pertrb – double scalar
pertrb = 0${\mathbf{pertrb}}=0$, unless a solution to Poisson's equation (λ = 0)$\left(\lambda =0\right)$ is required with a combination of periodic or derivative boundary conditions (lbdcnd, mbdcnd and nbdcnd = 0${\mathbf{nbdcnd}}=0$ or 3$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. nag_pde_3d_ellip_helmholtz (d03fa) 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.
3:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

ifail = 1${\mathbf{ifail}}=1$
 On entry, xs ≥ xf${\mathbf{xs}}\ge {\mathbf{xf}}$, or l < 5${\mathbf{l}}<5$, or lbdcnd < 0${\mathbf{lbdcnd}}<0$, or lbdcnd > 4${\mathbf{lbdcnd}}>4$, or ys ≥ yf${\mathbf{ys}}\ge {\mathbf{yf}}$, or m < 5${\mathbf{m}}<5$, or mbdcnd < 0${\mathbf{mbdcnd}}<0$, or mbdcnd > 4${\mathbf{mbdcnd}}>4$, or zs ≥ zf${\mathbf{zs}}\ge {\mathbf{zf}}$, or n < 5${\mathbf{n}}<5$, or nbdcnd < 0${\mathbf{nbdcnd}}<0$, or nbdcnd > 4${\mathbf{nbdcnd}}>4$, or ldf < l + 1$\mathit{ldf}<{\mathbf{l}}+1$, or ldf2 < m + 1$\mathit{ldf2}<{\mathbf{m}}+1$.
ifail = 2${\mathbf{ifail}}=2$
 On entry, lwrk is too small.
W ifail = 3${\mathbf{ifail}}=3$
 On entry, λ > 0$\lambda >0$.

## Accuracy

Not applicable.

The execution time is roughly proportional to l × m × n × (log2l + log2m + 5) ${\mathbf{l}}×{\mathbf{m}}×{\mathbf{n}}×\left({\mathrm{log}}_{2}{\mathbf{l}}+{\mathrm{log}}_{2}{\mathbf{m}}+5\right)$, but also depends on input parameters lbdcnd and mbdcnd.

## Example

```function nag_pde_3d_ellip_helmholtz_example
xs = 0;
xf = 1;
l = 16;
m = 32;
n = 20;
lbdcnd = int64(1);
bdxs = zeros(33, 21);
bdxf = zeros(33, 21);
ys = 0;
yf = 2*pi;
mbdcnd = int64(0);
bdys = zeros(17, 21);
bdyf = zeros(17, 21);
zs = 0;
zf = pi/2;
nbdcnd = int64(2);
bdzs = zeros(17, 33);
bdzf = zeros(17, 33);
lambda = -2;
f = zeros(17, 33, 21);

x = zeros(l+1, 1);
y = zeros(m+1, 1);
z = zeros(n+1, 1);

% Define the grid points for later use.
dx = (xf-xs)/double(l);
for i = 1:l+1
x(i) = xs + (i-1)*dx;
end
dy = (yf-ys)/m;
for j = 1:m+1
y(j) = ys + (j-1)*dy;
end
dz = (zf-zs)/n;
for k = 1:n+1
z(k) = zs + (k-1)*dz;
end

% Define the array of derivative boundary values.
for j = 1:m+1
for i = 1:l+1
bdzf(i,j) = -x(i)^4*sin(y(j));
end
end

% Note that for this example all other boundary arrays are
% dummy variables.
%
% We define the function boundary values in the f array.
for k = 1:n+1
for j = 1:m+1
f(l+1,j,k) = sin(y(j))*cos(z(k));
end
end
for j = 1:m+1
for i = 1:l+1
f(i,j,1) = x(i)^4*sin(y(j));
end
end

% Define the values of the right hand side of the Helmholtz
% equation.
for k = 2:n+1
for j = 1:m+1
for i = 2:l
f(i,j,k) = 4*x(i)^2*(3-x(i)^2)*sin(y(j))*cos(z(k));
end
end
end

for k=1:n+1
for j=1:m+1
for i=1:l+1
end
end
end
[fOut, pertrb, ifail] = ...
nag_pde_3d_ellip_helmholtz(xs, xf, int64(l), lbdcnd, bdxs, bdxf, ys, yf, int64(m), mbdcnd, bdys, bdyf, ...
zs, zf, int64(n), nbdcnd, bdzs, bdzf, lambda, f);
pertrb
ifail
```
```

pertrb =

0

ifail =

0

```
```function d03fa_example
xs = 0;
xf = 1;
l = 16;
m = 32;
n = 20;
lbdcnd = int64(1);
bdxs = zeros(33, 21);
bdxf = zeros(33, 21);
ys = 0;
yf = 2*pi;
mbdcnd = int64(0);
bdys = zeros(17, 21);
bdyf = zeros(17, 21);
zs = 0;
zf = pi/2;
nbdcnd = int64(2);
bdzs = zeros(17, 33);
bdzf = zeros(17, 33);
lambda = -2;
f = zeros(17, 33, 21);

x = zeros(l+1, 1);
y = zeros(m+1, 1);
z = zeros(n+1, 1);

% Define the grid points for later use.
dx = (xf-xs)/double(l);
for i = 1:l+1
x(i) = xs + (i-1)*dx;
end
dy = (yf-ys)/m;
for j = 1:m+1
y(j) = ys + (j-1)*dy;
end
dz = (zf-zs)/n;
for k = 1:n+1
z(k) = zs + (k-1)*dz;
end

% Define the array of derivative boundary values.
for j = 1:m+1
for i = 1:l+1
bdzf(i,j) = -x(i)^4*sin(y(j));
end
end

% Note that for this example all other boundary arrays are
% dummy variables.
%
% We define the function boundary values in the f array.
for k = 1:n+1
for j = 1:m+1
f(l+1,j,k) = sin(y(j))*cos(z(k));
end
end
for j = 1:m+1
for i = 1:l+1
f(i,j,1) = x(i)^4*sin(y(j));
end
end

% Define the values of the right hand side of the Helmholtz
% equation.
for k = 2:n+1
for j = 1:m+1
for i = 2:l
f(i,j,k) = 4*x(i)^2*(3-x(i)^2)*sin(y(j))*cos(z(k));
end
end
end

for k=1:n+1
for j=1:m+1
for i=1:l+1
end
end
end
[fOut, pertrb, ifail] = ...
d03fa(xs, xf, int64(l), lbdcnd, bdxs, bdxf, ys, yf, int64(m), mbdcnd, bdys, bdyf, ...
zs, zf, int64(n), nbdcnd, bdzs, bdzf, lambda, f);
pertrb
ifail
```
```

pertrb =

0

ifail =

0

```