hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
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 xx- and yy-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 xx, i.e., xsxxfxsxxf.
Constraint: xs < xfxs<xf.
2:     xf – double scalar
The upper bound of the range of xx, i.e., xsxxfxsxxf.
Constraint: xs < xfxs<xf.
3:     l – int64int32nag_int scalar
The number of panels into which the interval (xs,xf) is subdivided. Hence, there will be l + 1l+1 grid points in the xx-direction given by xi = xs + (i1) × δxxi=xs+(i-1)×δx, for i = 1,2,,l + 1i=1,2,,l+1, where δx = (xfxs) / lδx=(xf-xs)/l is the panel width.
Constraint: l5l5.
4:     lbdcnd – int64int32nag_int scalar
Indicates the type of boundary conditions at x = xsx=xs and x = xfx=xf.
lbdcnd = 0lbdcnd=0
If the solution is periodic in xx, i.e., u(xs,y,z) = u(xf,y,z)u(xs,y,z)=u(xf,y,z).
lbdcnd = 1lbdcnd=1
If the solution is specified at x = xsx=xs and x = xfx=xf.
lbdcnd = 2lbdcnd=2
If the solution is specified at x = xsx=xs and the derivative of the solution with respect to xx is specified at x = xfx=xf.
lbdcnd = 3lbdcnd=3
If the derivative of the solution with respect to xx is specified at x = xsx=xs and x = xfx=xf.
lbdcnd = 4lbdcnd=4
If the derivative of the solution with respect to xx is specified at x = xsx=xs and the solution is specified at x = xfx=xf.
Constraint: 0lbdcnd40lbdcnd4.
5:     bdxs(ldf2,n + 1n+1) – double array
ldf2, the first dimension of the array, must satisfy the constraint ldf2m + 1ldf2m+1.
The values of the derivative of the solution with respect to xx at x = xsx=xs. When lbdcnd = 3lbdcnd=3 or 44, bdxs(j,k) = ux(xs,yj,zk)bdxsjk=ux(xs,yj,zk), for j = 1,2,,m + 1j=1,2,,m+1 and k = 1,2,,n + 1k=1,2,,n+1.
When lbdcnd has any other value, bdxs is not referenced.
6:     bdxf(ldf2,n + 1n+1) – double array
ldf2, the first dimension of the array, must satisfy the constraint ldf2m + 1ldf2m+1.
The values of the derivative of the solution with respect to xx at x = xfx=xf. When lbdcnd = 2lbdcnd=2 or 33, bdxf(j,k) = ux(xf,yj,zk)bdxfjk=ux(xf,yj,zk), for j = 1,2,,m + 1j=1,2,,m+1 and k = 1,2,,n + 1k=1,2,,n+1.
When lbdcnd has any other value, bdxf is not referenced.
7:     ys – double scalar
The lower bound of the range of yy, i.e., ysyyfysyyf.
Constraint: ys < yfys<yf.
8:     yf – double scalar
The upper bound of the range of yy, i.e., ysyyfysyyf.
Constraint: ys < yfys<yf.
9:     m – int64int32nag_int scalar
The number of panels into which the interval (ys,yf) is subdivided. Hence, there will be m + 1m+1 grid points in the yy-direction given by yj = ys + (j1) × δyyj=ys+(j-1)×δy, for j = 1,2,,m + 1j=1,2,,m+1, where δy = (yfys) / mδy=(yf-ys)/m is the panel width.
Constraint: m5m5.
10:   mbdcnd – int64int32nag_int scalar
Indicates the type of boundary conditions at y = ysy=ys and y = yfy=yf.
mbdcnd = 0mbdcnd=0
If the solution is periodic in yy, i.e., u(x,yf,z) = u(x,ys,z)u(x,yf,z)=u(x,ys,z).
mbdcnd = 1mbdcnd=1
If the solution is specified at y = ysy=ys and y = yfy=yf.
mbdcnd = 2mbdcnd=2
If the solution is specified at y = ysy=ys and the derivative of the solution with respect to yy is specified at y = yfy=yf.
mbdcnd = 3mbdcnd=3
If the derivative of the solution with respect to yy is specified at y = ysy=ys and y = yfy=yf.
mbdcnd = 4mbdcnd=4
If the derivative of the solution with respect to yy is specified at y = ysy=ys and the solution is specified at y = yfy=yf.
Constraint: 0mbdcnd40mbdcnd4.
11:   bdys(ldf,n + 1n+1) – double array
ldf, the first dimension of the array, must satisfy the constraint ldfl + 1ldfl+1.
The values of the derivative of the solution with respect to yy at y = ysy=ys. When mbdcnd = 3mbdcnd=3 or 44, bdys(i,k) = uy(xi,ys,zk)bdysik=uy(xi,ys,zk), for i = 1,2,,l + 1i=1,2,,l+1 and k = 1,2,,n + 1k=1,2,,n+1.
When mbdcnd has any other value, bdys is not referenced.
12:   bdyf(ldf,n + 1n+1) – double array
ldf, the first dimension of the array, must satisfy the constraint ldfl + 1ldfl+1.
The values of the derivative of the solution with respect to yy at y = yfy=yf. When mbdcnd = 2mbdcnd=2 or 33, bdyf(i,k) = uy(xi,yf,zk)bdyfik=uy(xi,yf,zk), for i = 1,2,,l + 1i=1,2,,l+1 and k = 1,2,,n + 1k=1,2,,n+1.
When mbdcnd has any other value, bdyf is not referenced.
13:   zs – double scalar
The lower bound of the range of zz, i.e., zszzfzszzf.
Constraint: zs < zfzs<zf.
14:   zf – double scalar
The upper bound of the range of zz, i.e., zszzfzszzf.
Constraint: zs < zfzs<zf.
15:   n – int64int32nag_int scalar
The number of panels into which the interval (zs,zf) is subdivided. Hence, there will be n + 1n+1 grid points in the zz-direction given by zk = zs + (k1) × δzzk=zs+(k-1)×δz, for k = 1,2,,n + 1k=1,2,,n+1, where δz = (zfzs) / nδz=(zf-zs)/n is the panel width.
Constraint: n5n5.
16:   nbdcnd – int64int32nag_int scalar
Specifies the type of boundary conditions at z = zsz=zs and z = zfz=zf.
nbdcnd = 0nbdcnd=0
if the solution is periodic in zz, i.e., u(x,y,zf) = u(x,y,zs)u(x,y,zf)=u(x,y,zs).
nbdcnd = 1nbdcnd=1
if the solution is specified at z = zsz=zs and z = zfz=zf.
nbdcnd = 2nbdcnd=2
if the solution is specified at z = zsz=zs and the derivative of the solution with respect to zz is specified at z = zfz=zf.
nbdcnd = 3nbdcnd=3
if the derivative of the solution with respect to zz is specified at z = zsz=zs and z = zfz=zf.
nbdcnd = 4nbdcnd=4
if the derivative of the solution with respect to zz is specified at z = zsz=zs and the solution is specified at z = zfz=zf.
Constraint: 0nbdcnd40nbdcnd4.
17:   bdzs(ldf,m + 1m+1) – double array
ldf, the first dimension of the array, must satisfy the constraint ldfl + 1ldfl+1.
The values of the derivative of the solution with respect to zz at z = zsz=zs. When nbdcnd = 3nbdcnd=3 or 44, bdzs(i,j) = uz(xi,yj,zs)bdzsij=uz(xi,yj,zs), for i = 1,2,,l + 1i=1,2,,l+1 and j = 1,2,,m + 1j=1,2,,m+1.
When nbdcnd has any other value, bdzs is not referenced.
18:   bdzf(ldf,m + 1m+1) – double array
ldf, the first dimension of the array, must satisfy the constraint ldfl + 1ldfl+1.
The values of the derivative of the solution with respect to zz at z = zfz=zf. When nbdcnd = 2nbdcnd=2 or 33, bdzf(i,j) = uz(xi,yj,zf)bdzfij=uz(xi,yj,zf), for i = 1,2,,l + 1i=1,2,,l+1 and j = 1,2,,m + 1j=1,2,,m+1.
When nbdcnd has any other value, bdzf is not referenced.
19:   lambda – double scalar
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λ>0, then nag_pde_3d_ellip_helmholtz (d03fa) will set ifail = 3ifail=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λ>0 with great caution.
20:   f(ldf,ldf2,n + 1n+1) – double array
ldf, the first dimension of the array, must satisfy the constraint ldfl + 1ldfl+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)f(1,j,k) f(l + 1,j,k)f(l+1,j,k)
0 f(xs,yj,zk)f(xs,yj,zk) f(xs,yj,zk)f(xs,yj,zk)  
1 u(xs,yj,zk)u(xs,yj,zk) u(xf,yj,zk)u(xf,yj,zk)
2 u(xs,yj,zk)u(xs,yj,zk) f(xf,yj,zk)f(xf,yj,zk) j = 1,2,,m + 1j=1,2,,m+1
3 f(xs,yj,zk)f(xs,yj,zk) f(xf,yj,zk)f(xf,yj,zk) k = 1,2,,n + 1k=1,2,,n+1
4 f(xs,yj,zk)f(xs,yj,zk) u(xf,yj,zk)u(xf,yj,zk)
       
mbdcnd f(i,1,k)f(i,1,k) f(i,m + 1,k)f(i,m+1,k)
0 f(xi,ys,zk)f(xi,ys,zk) f(xi,ys,zk)f(xi,ys,zk)
1 u(ys,xi,zk)u(ys,xi,zk) u(yf,xi,zk)u(yf,xi,zk)
2 u(xi,ys,zk)u(xi,ys,zk) f(xi,yf,zk)f(xi,yf,zk) i = 1,2,,l + 1i=1,2,,l+1
3 f(xi,ys,zk)f(xi,ys,zk) f(xi,yf,zk)f(xi,yf,zk) k = 1,2,,n + 1k=1,2,,n+1
4 f(xi,ys,zk)f(xi,ys,zk) u(xi,yf,zk)u(xi,yf,zk)
       
nbdcnd f(i,j,1)f(i,j,1) f(i,j,n + 1)f(i,j,n+1)
0 f(xi,yj,zs)f(xi,yj,zs) f(xi,yj,zs)f(xi,yj,zs)
1 u(xi,yj,zs)u(xi,yj,zs) u(xi,yj,zf)u(xi,yj,zf)
2 u(xi,yj,zs)u(xi,yj,zs) f(xi,yj,zf)f(xi,yj,zf) i = 1,2,,l + 1i=1,2,,l+1
3 f(xi,yj,zs)f(xi,yj,zs) f(xi,yj,zf)f(xi,yj,zf) j = 1,2,,m + 1j=1,2,,m+1
4 f(xi,yj,zs)f(xi,yj,zs) u(xi,yj,zf)u(xi,yj,zf)
Note: if the table calls for both the solution uu and the right-hand side ff 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 + 62×(n+1)×max(l,m)+3×l+3×m+4×n+6 is an upper bound on the required size of w. If lwrk is too small, the function exits with ifail = 2ifail=2, and if on entry ifail = 0ifail=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 + 62×(n+1)×max(l,m)+3×l+3×m+4×n+6

Input Parameters Omitted from the MATLAB Interface

ldf ldf2 w

Output Parameters

1:     f(ldf,ldf2,n + 1n+1) – double array
ldfl + 1ldfl+1.
ldf2m + 1ldf2m+1.
Contains the solution u(i,j,k)u(i,j,k) of the finite difference approximation for the grid point (xi,yj,zk)(xi,yj,zk), for i = 1,2,,l + 1i=1,2,,l+1, j = 1,2,,m + 1j=1,2,,m+1 and k = 1,2,,n + 1k=1,2,,n+1.
2:     pertrb – double scalar
pertrb = 0pertrb=0, unless a solution to Poisson's equation (λ = 0)(λ=0) is required with a combination of periodic or derivative boundary conditions (lbdcnd, mbdcnd and nbdcnd = 0nbdcnd=0 or 33). 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
ifail = 0ifail=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 = 1ifail=1
On entry,xsxfxsxf,
orl < 5l<5,
orlbdcnd < 0lbdcnd<0,
orlbdcnd > 4lbdcnd>4,
orysyfysyf,
orm < 5m<5,
ormbdcnd < 0mbdcnd<0,
ormbdcnd > 4mbdcnd>4,
orzszfzszf,
orn < 5n<5,
ornbdcnd < 0nbdcnd<0,
ornbdcnd > 4nbdcnd>4,
orldf < l + 1ldf<l+1,
orldf2 < m + 1ldf2<m+1.
  ifail = 2ifail=2
On entry,lwrk is too small.
W ifail = 3ifail=3
On entry,λ > 0λ>0.

Accuracy

Not applicable.

Further Comments

The execution time is roughly proportional to l × m × n × (log2l + log2m + 5) l×m×n× ( log2l + log2m +5 ) , 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



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013