NAG Toolbox: nag_pde_3d_ellip_helmholtz (d03fa)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{xs}$ – double scalar

The lower bound of the range of $x$, i.e., $\mathbf{xs}\le x\le \mathbf{xf}$.

Constraint: $\mathbf{xs}<\mathbf{xf}$.

2:     $\mathrm{xf}$ – double scalar

The upper bound of the range of $x$, i.e., $\mathbf{xs}\le x\le \mathbf{xf}$.

Constraint: $\mathbf{xs}<\mathbf{xf}$.

3:     $\mathrm{l}$ – int64int32nag_int scalar

The number of panels into which the interval (xs,xf) is subdivided. Hence, there will be $\mathbf{l}+1$ grid points in the $x$-direction given by $x_{\mathit{i}}=\mathbf{xs}+\left(\mathit{i}-1\right)\times \delta x$, for $\mathit{i}=1,2,\dots ,\mathbf{l}+1$, where $\delta x=\left(\mathbf{xf}-\mathbf{xs}\right)/\mathbf{l}$ is the panel width.

Constraint: $\mathbf{l}\ge 5$.

4:     $\mathrm{lbdcnd}$ – int64int32nag_int scalar

Indicates the type of boundary conditions at $x=\mathbf{xs}$ and $x=\mathbf{xf}$.

$\mathbf{lbdcnd}=0$

If the solution is periodic in $x$, i.e., $u\left(\mathbf{xs},y,z\right)=u\left(\mathbf{xf},y,z\right)$.

$\mathbf{lbdcnd}=1$

If the solution is specified at $x=\mathbf{xs}$ and $x=\mathbf{xf}$.

$\mathbf{lbdcnd}=2$

If the solution is specified at $x=\mathbf{xs}$ and the derivative of the solution with respect to $x$ is specified at $x=\mathbf{xf}$.

$\mathbf{lbdcnd}=3$

If the derivative of the solution with respect to $x$ is specified at $x=\mathbf{xs}$ and $x=\mathbf{xf}$.

$\mathbf{lbdcnd}=4$

If the derivative of the solution with respect to $x$ is specified at $x=\mathbf{xs}$ and the solution is specified at $x=\mathbf{xf}$.

Constraint: $0\le \mathbf{lbdcnd}\le 4$.

5:     $\mathrm{bdxs}\left(\mathit{ldf2},\mathbf{n}+1\right)$ – double array

ldf2, the first dimension of the array, must satisfy the constraint $\mathit{ldf2}\ge \mathbf{m}+1$.

The values of the derivative of the solution with respect to $x$ at $x=\mathbf{xs}$. When $\mathbf{lbdcnd}=3$ or $4$, $\mathbf{bdxs}\left(\mathit{j},\mathit{k}\right)=u_x\left(\mathbf{xs},y_j,z_k\right)$, for $\mathit{j}=1,2,\dots ,\mathbf{m}+1$ and $\mathit{k}=1,2,\dots ,\mathbf{n}+1$.

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

6:     $\mathrm{bdxf}\left(\mathit{ldf2},\mathbf{n}+1\right)$ – double array

ldf2, the first dimension of the array, must satisfy the constraint $\mathit{ldf2}\ge \mathbf{m}+1$.

The values of the derivative of the solution with respect to $x$ at $x=\mathbf{xf}$. When $\mathbf{lbdcnd}=2$ or $3$, $\mathbf{bdxf}\left(\mathit{j},\mathit{k}\right)=u_x\left(\mathbf{xf},y_{\mathit{j}},z_{\mathit{k}}\right)$, for $\mathit{j}=1,2,\dots ,\mathbf{m}+1$ and $\mathit{k}=1,2,\dots ,\mathbf{n}+1$.

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

7:     $\mathrm{ys}$ – double scalar

The lower bound of the range of $y$, i.e., $\mathbf{ys}\le y\le \mathbf{yf}$.

Constraint: $\mathbf{ys}<\mathbf{yf}$.

8:     $\mathrm{yf}$ – double scalar

The upper bound of the range of $y$, i.e., $\mathbf{ys}\le y\le \mathbf{yf}$.

Constraint: $\mathbf{ys}<\mathbf{yf}$.

9:     $\mathrm{m}$ – int64int32nag_int scalar

The number of panels into which the interval (ys,yf) is subdivided. Hence, there will be $\mathbf{m}+1$ grid points in the $y$-direction given by $y_{\mathit{j}}=\mathbf{ys}+\left(\mathit{j}-1\right)\times \delta y$, for $\mathit{j}=1,2,\dots ,\mathbf{m}+1$, where $\delta y=\left(\mathbf{yf}-\mathbf{ys}\right)/\mathbf{m}$ is the panel width.

Constraint: $\mathbf{m}\ge 5$.

10:   $\mathrm{mbdcnd}$ – int64int32nag_int scalar

Indicates the type of boundary conditions at $y=\mathbf{ys}$ and $y=\mathbf{yf}$.

$\mathbf{mbdcnd}=0$

If the solution is periodic in $y$, i.e., $u\left(x,\mathbf{yf},z\right)=u\left(x,\mathbf{ys},z\right)$.

$\mathbf{mbdcnd}=1$

If the solution is specified at $y=\mathbf{ys}$ and $y=\mathbf{yf}$.

$\mathbf{mbdcnd}=2$

If the solution is specified at $y=\mathbf{ys}$ and the derivative of the solution with respect to $y$ is specified at $y=\mathbf{yf}$.

$\mathbf{mbdcnd}=3$

If the derivative of the solution with respect to $y$ is specified at $y=\mathbf{ys}$ and $y=\mathbf{yf}$.

$\mathbf{mbdcnd}=4$

If the derivative of the solution with respect to $y$ is specified at $y=\mathbf{ys}$ and the solution is specified at $y=\mathbf{yf}$.

Constraint: $0\le \mathbf{mbdcnd}\le 4$.

11:   $\mathrm{bdys}\left(\mathit{ldf},\mathbf{n}+1\right)$ – double array

ldf, the first dimension of the array, must satisfy the constraint $\mathit{ldf}\ge \mathbf{l}+1$.

The values of the derivative of the solution with respect to $y$ at $y=\mathbf{ys}$. When $\mathbf{mbdcnd}=3$ or $4$, $\mathbf{bdys}\left(\mathit{i},\mathit{k}\right)=u_y\left(x_{\mathit{i}},\mathbf{ys},z_{\mathit{k}}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{l}+1$ and $\mathit{k}=1,2,\dots ,\mathbf{n}+1$.

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

12:   $\mathrm{bdyf}\left(\mathit{ldf},\mathbf{n}+1\right)$ – double array

ldf, the first dimension of the array, must satisfy the constraint $\mathit{ldf}\ge \mathbf{l}+1$.

The values of the derivative of the solution with respect to $y$ at $y=\mathbf{yf}$. When $\mathbf{mbdcnd}=2$ or $3$, $\mathbf{bdyf}\left(\mathit{i},\mathit{k}\right)=u_y\left(x_{\mathit{i}},\mathbf{yf},z_{\mathit{k}}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{l}+1$ and $\mathit{k}=1,2,\dots ,\mathbf{n}+1$.

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

13:   $\mathrm{zs}$ – double scalar

The lower bound of the range of $z$, i.e., $\mathbf{zs}\le z\le \mathbf{zf}$.

Constraint: $\mathbf{zs}<\mathbf{zf}$.

14:   $\mathrm{zf}$ – double scalar

The upper bound of the range of $z$, i.e., $\mathbf{zs}\le z\le \mathbf{zf}$.

Constraint: $\mathbf{zs}<\mathbf{zf}$.

15:   $\mathrm{n}$ – int64int32nag_int scalar

The number of panels into which the interval (zs,zf) is subdivided. Hence, there will be $\mathbf{n}+1$ grid points in the $z$-direction given by $z_{\mathit{k}}=\mathbf{zs}+\left(\mathit{k}-1\right)\times \delta z$, for $\mathit{k}=1,2,\dots ,\mathbf{n}+1$, where $\delta z=\left(\mathbf{zf}-\mathbf{zs}\right)/\mathbf{n}$ is the panel width.

Constraint: $\mathbf{n}\ge 5$.

16:   $\mathrm{nbdcnd}$ – int64int32nag_int scalar

Specifies the type of boundary conditions at $z=\mathbf{zs}$ and $z=\mathbf{zf}$.

$\mathbf{nbdcnd}=0$

if the solution is periodic in $z$, i.e., $u\left(x,y,\mathbf{zf}\right)=u\left(x,y,\mathbf{zs}\right)$.

$\mathbf{nbdcnd}=1$

if the solution is specified at $z=\mathbf{zs}$ and $z=\mathbf{zf}$.

$\mathbf{nbdcnd}=2$

if the solution is specified at $z=\mathbf{zs}$ and the derivative of the solution with respect to $z$ is specified at $z=\mathbf{zf}$.

$\mathbf{nbdcnd}=3$

if the derivative of the solution with respect to $z$ is specified at $z=\mathbf{zs}$ and $z=\mathbf{zf}$.

$\mathbf{nbdcnd}=4$

if the derivative of the solution with respect to $z$ is specified at $z=\mathbf{zs}$ and the solution is specified at $z=\mathbf{zf}$.

Constraint: $0\le \mathbf{nbdcnd}\le 4$.

17:   $\mathrm{bdzs}\left(\mathit{ldf},\mathbf{m}+1\right)$ – double array

ldf, the first dimension of the array, must satisfy the constraint $\mathit{ldf}\ge \mathbf{l}+1$.

The values of the derivative of the solution with respect to $z$ at $z=\mathbf{zs}$. When $\mathbf{nbdcnd}=3$ or $4$, $\mathbf{bdzs}\left(\mathit{i},\mathit{j}\right)=u_z\left(x_i,y_j,\mathbf{zs}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{l}+1$ and $\mathit{j}=1,2,\dots ,\mathbf{m}+1$.

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

18:   $\mathrm{bdzf}\left(\mathit{ldf},\mathbf{m}+1\right)$ – double array

ldf, the first dimension of the array, must satisfy the constraint $\mathit{ldf}\ge \mathbf{l}+1$.

The values of the derivative of the solution with respect to $z$ at $z=\mathbf{zf}$. When $\mathbf{nbdcnd}=2$ or $3$, $\mathbf{bdzf}\left(\mathit{i},\mathit{j}\right)=u_z\left(x_i,y_j,\mathbf{zf}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{l}+1$ and $\mathit{j}=1,2,\dots ,\mathbf{m}+1$.

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

19:   $\mathrm{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 $\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 $\lambda >0$ with great caution.

20:   $\mathrm{f}\left(\mathit{ldf},\mathit{ldf2},\mathbf{n}+1\right)$ – double array

ldf, the first dimension of the array, must satisfy the constraint $\mathit{ldf}\ge \mathbf{l}+1$.

The values of the right-side of the Helmholtz equation and boundary values (if any).

$$\mathbf{f}\left(i,j,k\right)=f\left(x_i,y_j,z_k\right)i=2,3,\dots ,\mathbf{l},j=2,3,\dots ,\mathbf{m}\text{ and }k=2,3,\dots ,\mathbf{n}\text{.}$$

On the boundaries f is defined by

lbdcnd	$\mathbf{f}\left(1,j,k\right)$	$\mathbf{f}\left(\mathbf{l}+1,j,k\right)$
0	$f\left(\mathbf{xs},y_j,z_k\right)$	$f\left(\mathbf{xs},y_j,z_k\right)$
1	$u\left(\mathbf{xs},y_j,z_k\right)$	$u\left(\mathbf{xf},y_j,z_k\right)$
2	$u\left(\mathbf{xs},y_j,z_k\right)$	$f\left(\mathbf{xf},y_j,z_k\right)$	$j=1,2,\dots ,\mathbf{m}+1$
3	$f\left(\mathbf{xs},y_j,z_k\right)$	$f\left(\mathbf{xf},y_j,z_k\right)$	$k=1,2,\dots ,\mathbf{n}+1$
4	$f\left(\mathbf{xs},y_j,z_k\right)$	$u\left(\mathbf{xf},y_j,z_k\right)$
mbdcnd	$\mathbf{f}\left(i,1,k\right)$	$\mathbf{f}\left(i,\mathbf{m}+1,k\right)$
0	$f\left(x_i,\mathbf{ys},z_k\right)$	$f\left(x_i,\mathbf{ys},z_k\right)$
1	$u\left(\mathbf{ys},x_i,z_k\right)$	$u\left(\mathbf{yf},x_i,z_k\right)$
2	$u\left(x_i,\mathbf{ys},z_k\right)$	$f\left(x_i,\mathbf{yf},z_k\right)$	$i=1,2,\dots ,\mathbf{l}+1$
3	$f\left(x_i,\mathbf{ys},z_k\right)$	$f\left(x_i,\mathbf{yf},z_k\right)$	$k=1,2,\dots ,\mathbf{n}+1$
4	$f\left(x_i,\mathbf{ys},z_k\right)$	$u\left(x_i,\mathbf{yf},z_k\right)$
nbdcnd	$\mathbf{f}\left(i,j,1\right)$	$\mathbf{f}\left(i,j,\mathbf{n}+1\right)$
0	$f\left(x_i,y_j,\mathbf{zs}\right)$	$f\left(x_i,y_j,\mathbf{zs}\right)$
1	$u\left(x_i,y_j,\mathbf{zs}\right)$	$u\left(x_i,y_j,\mathbf{zf}\right)$
2	$u\left(x_i,y_j,\mathbf{zs}\right)$	$f\left(x_i,y_j,\mathbf{zf}\right)$	$i=1,2,\dots ,\mathbf{l}+1$
3	$f\left(x_i,y_j,\mathbf{zs}\right)$	$f\left(x_i,y_j,\mathbf{zf}\right)$	$j=1,2,\dots ,\mathbf{m}+1$
4	$f\left(x_i,y_j,\mathbf{zs}\right)$	$u\left(x_i,y_j,\mathbf{zf}\right)$

Note: if the table calls for both the solution $u$ and the right-hand side $f$ on a boundary, then the solution must be specified.

Optional Input Parameters

1:     $\mathrm{lwrk}$ – int64int32nag_int scalar

Default: $2\times \left(\mathbf{n}+1\right)\times \max \left(\mathbf{l},\mathbf{m}\right)+3\times \mathbf{l}+3\times \mathbf{m}+4\times \mathbf{n}+6$

The dimension of the array w. $2\times \left(\mathbf{n}+1\right)\times \max \left(\mathbf{l},\mathbf{m}\right)+3\times \mathbf{l}+3\times \mathbf{m}+4\times \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 $\mathbf{ifail}=0$ or $-1$, a message is output giving the exact value of lwrk required to solve the current problem.

Output Parameters

1:     $\mathrm{f}\left(\mathit{ldf},\mathit{ldf2},\mathbf{n}+1\right)$ – double array

$\mathit{ldf2}=\mathbf{m}+1$.

Contains the solution $u\left(\mathit{i},\mathit{j},\mathit{k}\right)$ of the finite difference approximation for the grid point $\left(x_{\mathit{i}},y_{\mathit{j}},z_{\mathit{k}}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{l}+1$, $\mathit{j}=1,2,\dots ,\mathbf{m}+1$ and $\mathit{k}=1,2,\dots ,\mathbf{n}+1$.

2:     $\mathrm{pertrb}$ – double scalar

$\mathbf{pertrb}=0$, unless a solution to Poisson's equation $\left(\lambda =0\right)$ is required with a combination of periodic or derivative boundary conditions (lbdcnd, mbdcnd and $\mathbf{nbdcnd}=0$ or $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:     $\mathrm{ifail}$ – int64int32nag_int scalar

$\mathbf{ifail}=\mathbf{0}$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

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

   $\mathbf{ifail}=1$

On entry,	$\mathbf{xs}\ge \mathbf{xf}$,
or	$\mathbf{l}<5$,
or	$\mathbf{lbdcnd}<0$,
or	$\mathbf{lbdcnd}>4$,
or	$\mathbf{ys}\ge \mathbf{yf}$,
or	$\mathbf{m}<5$,
or	$\mathbf{mbdcnd}<0$,
or	$\mathbf{mbdcnd}>4$,
or	$\mathbf{zs}\ge \mathbf{zf}$,
or	$\mathbf{n}<5$,
or	$\mathbf{nbdcnd}<0$,
or	$\mathbf{nbdcnd}>4$,
or	$\mathit{ldf}<\mathbf{l}+1$,
or	$\mathit{ldf2}<\mathbf{m}+1$.

   $\mathbf{ifail}=2$

On entry,

lwrk is too small.

W  $\mathbf{ifail}=3$

On entry,

$\lambda >0$.

   $\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

   $\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

   $\mathbf{ifail}=-999$

Dynamic memory allocation failed.

Accuracy

Further Comments

Example

Open in the MATLAB editor: d03fa_example

function d03fa_example


fprintf('d03fa example results\n\n');

xs = 0; xf = 1;    l = 16;
ys = 0; yf = 2*pi; m = 32;
zs = 0; zf = pi/2; n = 20;

bdxs = zeros(m+1, n+1); bdxf = zeros(m+1, n+1);
bdys = zeros(l+1, n+1); bdyf = zeros(l+1, n+1);
bdzs = zeros(l+1, m+1); bdzf = zeros(l+1, m+1);

lbdcnd = int64(1); 
mbdcnd = int64(0); 
nbdcnd = int64(2); 
lambda = -2;

% Define the grid points for later use.
dx = (xf-xs)/l;  x  = [xs:dx:xf];
dy = (yf-ys)/m;  y  = [ys:dy:yf];
dz = (zf-zs)/n;  z  = [zs:dz:zf];

% Define the array of derivative boundary values (z only here).
XY = [x.^4]'*sin(y);
bdzf = -XY;

% Define the f array including boundary conditions.
YZ = [sin(y)]'*cos(z);
for i = 1:l
  f(i,:,:) = (4*x(i)^2*(3-x(i)^2))*YZ;
  g(i,:,:) = x(i)^4*YZ;
end
g(l+1,:,:) = YZ;
f(1,:,:)   = 0;
f(l+1,:,:) = YZ;
f(:,:,1)   = XY;

[u, 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);
% calculate error
maxerr   = max(max(max(abs(u - g))));
fprintf('Maximum error in computed solution = %10.3e\n',maxerr);

fig1 = figure;
[xs,ys,zs] = meshgrid(y,x,z);
slice(xs,ys,zs,u,[1.8 5],[0.95],[0.78]);
axis([y(1) y(end) x(1) x(end) 0 2]);
xlabel('y'); ylabel('x'); zlabel('z');
title('Helholtz Equation solution in a box');
view(111, 26);
h = colorbar;
ylabel(h, 'u(x,y,z)')

d03fa example results

Maximum error in computed solution =  5.177e-04