NAG Toolbox: nag_pde_3d_ellip_fd (d03ec)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

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

The number of nodes in the first coordinate direction, $n_1$.

Constraint: $\mathbf{n1}>1$.

2:     $\mathrm{n2}$ – int64int32nag_int scalar

The number of nodes in the second coordinate direction, $n_2$.

Constraint: $\mathbf{n2}>1$.

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

The number of nodes in the third coordinate direction, $n_3$.

Constraint: $\mathbf{n3}>1$.

4:     $\mathrm{a}\left(\mathit{lda},\mathbf{sda},\mathbf{n3}\right)$ – double array

lda, the first dimension of the array, must satisfy the constraint $\mathit{lda}\ge \mathbf{n1}$.

$\mathbf{a}\left(\mathit{i},\mathit{j},\mathit{k}\right)$ must contain the coefficient of $t_{\mathit{i}\mathit{j},\mathit{k}-1}$ in the $\left(\mathit{i},\mathit{j},\mathit{k}\right)$th equation of the system (1), for $\mathit{i}=1,2,\dots ,\mathbf{n1}$, $\mathit{j}=1,2,\dots ,\mathbf{n2}$ and $\mathit{k}=1,2,\dots ,\mathbf{n3}$. The elements of a, for $k=1$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the box.

5:     $\mathrm{b}\left(\mathit{lda},\mathbf{sda},\mathbf{n3}\right)$ – double array

lda, the first dimension of the array, must satisfy the constraint $\mathit{lda}\ge \mathbf{n1}$.

$\mathbf{b}\left(\mathit{i},\mathit{j},\mathit{k}\right)$ must contain the coefficient of $t_{\mathit{i},\mathit{j}-1,\mathit{k}}$ in the $\left(\mathit{i},\mathit{j},\mathit{k}\right)$th equation of the system (1), for $\mathit{i}=1,2,\dots ,\mathbf{n1}$, $\mathit{j}=1,2,\dots ,\mathbf{n2}$ and $\mathit{k}=1,2,\dots ,\mathbf{n3}$. The elements of b, for $j=1$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the box.

6:     $\mathrm{c}\left(\mathit{lda},\mathbf{sda},\mathbf{n3}\right)$ – double array

lda, the first dimension of the array, must satisfy the constraint $\mathit{lda}\ge \mathbf{n1}$.

$\mathbf{c}\left(\mathit{i},\mathit{j},\mathit{k}\right)$ must contain the coefficient of $t_{\mathit{i}-1,\mathit{j}\mathit{k}}$ in the $\left(\mathit{i},\mathit{j},\mathit{k}\right)$th equation of the system (1), for $\mathit{i}=1,2,\dots ,\mathbf{n1}$, $\mathit{j}=1,2,\dots ,\mathbf{n2}$ and $\mathit{k}=1,2,\dots ,\mathbf{n3}$. The elements of c, for $i=1$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the box.

7:     $\mathrm{d}\left(\mathit{lda},\mathbf{sda},\mathbf{n3}\right)$ – double array

lda, the first dimension of the array, must satisfy the constraint $\mathit{lda}\ge \mathbf{n1}$.

$\mathbf{d}\left(\mathit{i},\mathit{j},\mathit{k}\right)$ must contain the coefficient of $t_{\mathit{i}\mathit{j}\mathit{k}}$ (the ‘central’ term) in the $\left(\mathit{i},\mathit{j},\mathit{k}\right)$th equation of the system (1), for $\mathit{i}=1,2,\dots ,\mathbf{n1}$, $\mathit{j}=1,2,\dots ,\mathbf{n2}$ and $\mathit{k}=1,2,\dots ,\mathbf{n3}$. The elements of d are checked to ensure that they are nonzero. If any element is found to be zero, the corresponding algebraic equation is assumed to be $t_{ijk}=q_{ijk}$. This feature can be used to define the equations for nodes at which, for example, Dirichlet boundary conditions are applied, or for nodes external to the problem of interest. Setting $\mathbf{d}\left(i,j,k\right)=0.0$ at appropriate points, and the corresponding value of $\mathbf{q}\left(i,j,k\right)$ to the appropriate value, namely the prescribed value of $\mathbf{t}\left(i,j,k\right)$ in the Dirichlet case, or to zero at an external point.

8:     $\mathrm{e}\left(\mathit{lda},\mathbf{sda},\mathbf{n3}\right)$ – double array

lda, the first dimension of the array, must satisfy the constraint $\mathit{lda}\ge \mathbf{n1}$.

$\mathbf{e}\left(\mathit{i},\mathit{j},\mathit{k}\right)$ must contain the coefficient of $t_{\mathit{i}+1,\mathit{j}\mathit{k}}$ in the $\left(\mathit{i},\mathit{j},\mathit{k}\right)$th equation of the system (1), for $\mathit{i}=1,2,\dots ,\mathbf{n1}$, $\mathit{j}=1,2,\dots ,\mathbf{n2}$ and $\mathit{k}=1,2,\dots ,\mathbf{n3}$. The elements of e, for $i=\mathbf{n1}$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the box.

9:     $\mathrm{f}\left(\mathit{lda},\mathbf{sda},\mathbf{n3}\right)$ – double array

lda, the first dimension of the array, must satisfy the constraint $\mathit{lda}\ge \mathbf{n1}$.

$\mathbf{f}\left(\mathit{i},\mathit{j},\mathit{k}\right)$ must contain the coefficient of $t_{\mathit{i},\mathit{j}+1,\mathit{k}}$ in the $\left(\mathit{i},\mathit{j},\mathit{k}\right)$th equation of the system (1), , for $\mathit{i}=1,2,\dots ,\mathbf{n1}$, $\mathit{j}=1,2,\dots ,\mathbf{n2}$ and $\mathit{k}=1,2,\dots ,\mathbf{n3}$. The elements of f, for $j=\mathbf{n2}$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the box.

10:   $\mathrm{g}\left(\mathit{lda},\mathbf{sda},\mathbf{n3}\right)$ – double array

lda, the first dimension of the array, must satisfy the constraint $\mathit{lda}\ge \mathbf{n1}$.

$\mathbf{g}\left(\mathit{i},\mathit{j},\mathit{k}\right)$ must contain the coefficient of $t_{\mathit{i}\mathit{j},\mathit{k}+1}$ in the $\left(\mathit{i},\mathit{j},\mathit{k}\right)$th equation of the system (1), for $\mathit{i}=1,2,\dots ,\mathbf{n1}$, $\mathit{j}=1,2,\dots ,\mathbf{n2}$ and $\mathit{k}=1,2,\dots ,\mathbf{n3}$. The elements of g, for $k=\mathbf{n3}$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the box.

11:   $\mathrm{q}\left(\mathit{lda},\mathbf{sda},\mathbf{n3}\right)$ – double array

lda, the first dimension of the array, must satisfy the constraint $\mathit{lda}\ge \mathbf{n1}$.

$\mathbf{q}\left(\mathit{i},\mathit{j},\mathit{k}\right)$ must contain $q_{\mathit{i}\mathit{j}\mathit{k}}$, for $\mathit{i}=1,2,\dots ,\mathbf{n1}$, $\mathit{j}=1,2,\dots ,\mathbf{n2}$ and $\mathit{k}=1,2,\dots ,\mathbf{n3}$, i.e., the source-term values at the nodal points of the system (1).

12:   $\mathrm{t}\left(\mathit{lda},\mathbf{sda},\mathbf{n3}\right)$ – double array

lda, the first dimension of the array, must satisfy the constraint $\mathit{lda}\ge \mathbf{n1}$.

$\mathbf{t}\left(\mathit{i},\mathit{j},\mathit{k}\right)$ must contain the element $t_{\mathit{i}\mathit{j}\mathit{k}}$ of an approximate solution to the equations, for $\mathit{i}=1,2,\dots ,\mathbf{n1}$, $\mathit{j}=1,2,\dots ,\mathbf{n2}$ and $\mathit{k}=1,2,\dots ,\mathbf{n3}$.

If no better approximation is known, an array of zeros can be used.

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

The iteration acceleration factor. A value of $1.0$ is adequate for most typical problems. However, if convergence is slow, the value can be reduced, typically to $0.2$ or $0.1$. If divergence is obtained, the value can be increased, typically to $2.0$, $5.0$ or $10.0$.

Constraint: $0.0<\mathbf{aparam}\le \left({\left(\mathbf{n1}-1\right)}^2+{\left(\mathbf{n2}-1\right)}^2+{\left(\mathbf{n3}-1\right)}^2\right)/3.0$.

14:   $\mathrm{itmax}$ – int64int32nag_int scalar

The maximum number of iterations to be used by the function in seeking the solution. A reasonable value might be $20$ for a problem with $3000$ nodes and convergence criteria of about ${10}^{-3}$ of the original residual and change.

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

On the first call of nag_pde_3d_ellip_fd (d03ec), itcoun must be set to $0$. On subsequent entries, its value must be unchanged from the previous call.

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

Indicates whether or not the system of equations has a unique solution. For systems which have a unique solution, ndir must be set to any nonzero value. For systems derived from problems to which an arbitrary constant can be added to the solution, for example Poisson's equation with all Neumann boundary conditions, ndir should be set to $0$ and the values of the next three arguments must be specified. For such problems the function subtracts the value of the function derived at the node (ixn, iyn, izn) from the whole solution after each iteration to reduce the possibility of large rounding errors. You must also ensure for such problems that the appropriate compatibility condition on the source terms q is satisfied. See the comments at the end of Description.

17:   $\mathrm{ixn}$ – int64int32nag_int scalar

Is ignored unless ndir is equal to zero, in which case it must specify the first index of the nodal point at which the solution is to be set to zero. The node should not correspond to a corner node, or to a node external to the region of interest.

18:   $\mathrm{iyn}$ – int64int32nag_int scalar

Is ignored unless ndir is equal to zero, in which case it must specify the second index of the nodal point at which the solution is to be set to zero. The node should not correspond to a corner node, or to a node external to the region of interest.

19:   $\mathrm{izn}$ – int64int32nag_int scalar

Is ignored unless ndir is equal to zero, in which case it must specify the third index of the nodal point at which the solution is to be set to zero. The node should not correspond to a corner node, or to a node external to the region of interest.

20:   $\mathrm{conres}$ – double scalar

The convergence criterion to be used on the maximum absolute value of the normalized residual vector components. The latter is defined as the residual of the algebraic equation divided by the central coefficient when the latter is not equal to $0.0$, and defined as the residual when the central coefficient is zero.

conres should not be less than a reasonable multiple of the machine precision.

21:   $\mathrm{conchn}$ – double scalar

The convergence criterion to be used on the maximum absolute value of the change made at each iteration to the elements of the array t, namely the dependent variable. conchn should not be less than a reasonable multiple of the machine accuracy multiplied by the maximum value of t attained.

Convergence is achieved when both the convergence criteria are satisfied. You can therefore set convergence on either the residual or on the change, or (as is recommended) on a requirement that both are below prescribed limits.

Optional Input Parameters

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

Default: the second dimension of the arrays a, b, c, d, e, f, g, q, t. (An error is raised if these dimensions are not equal.)

The second dimension of the arrays a, b, c, d, e, f, g, q, t, wrksp1, wrksp2, wrksp3 and wrksp4.

Constraint: $\mathbf{sda}\ge \mathbf{n2}$.

Output Parameters

1:     $\mathrm{t}\left(\mathit{lda},\mathbf{sda},\mathbf{n3}\right)$ – double array

The solution derived by the function.

2:     $\mathrm{itcoun}$ – int64int32nag_int scalar

Its value is increased by the number of iterations used on this call (namely itused). It therefore stores the accumulated number of iterations actually used.

For subsequent calls for the same problem, i.e., with the same n1, n2 and n3 but possibly different coefficients and/or source terms, as occur with nonlinear systems or with time-dependent systems, itcoun should not be reset, i.e., it must contain the accumulated number of iterations. In this way a suitable cycling of the sequence of iteration arguments is obtained in the calls to nag_pde_3d_ellip_fd_iter (d03ub).

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

The number of iterations actually used on that call.

4:     $\mathrm{resids}\left(\mathbf{itmax}\right)$ – double array

The maximum absolute value of the residuals calculated at the $\mathit{i}$th iteration, for $\mathit{i}=1,2,\dots ,\mathbf{itused}$. If the residual of the solution is sought you must calculate this in the function from which nag_pde_3d_ellip_fd (d03ec) is called. The sequence of values resids indicates the rate of convergence.

5:     $\mathrm{chngs}\left(\mathbf{itmax}\right)$ – double array

The maximum absolute value of the changes made to the components of the dependent variable t at the $\mathit{i}$th iteration, for $\mathit{i}=1,2,\dots ,\mathbf{itused}$. The sequence of values chngs indicates the rate of convergence.

6:     $\mathrm{ifail}$ – int64int32nag_int scalar

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

Error Indicators and Warnings

   $\mathbf{ifail}=1$

On entry,	$\mathbf{n1}<2$,
or	$\mathbf{n2}<2$,
or	$\mathbf{n3}<2$.

   $\mathbf{ifail}=2$

On entry,	$\mathit{lda}<\mathbf{n1}$,
or	$\mathbf{sda}<\mathbf{n2}$.

   $\mathbf{ifail}=3$

On entry,

$\mathbf{aparam}\le 0.0$.

   $\mathbf{ifail}=4$

On entry,

$\mathbf{aparam}>\left({\left(\mathbf{n1}-1\right)}^2+{\left(\mathbf{n2}-1\right)}^2+{\left(\mathbf{n3}-1\right)}^2\right)/3.0$.

   $\mathbf{ifail}=5$

Convergence was not achieved after itmax iterations.

   $\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: d03ec_example

function d03ec_example


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

% Solve Laplace equation in a rectangular box using d03ec. 
% Diricchlet boundary conditions are given as eact soltion on boundary.

% Mesh size
n1 = int64(16);
n2 = int64(20);
n3 = int64(24);

% Allocate arrays and set up initial values.
a = zeros(n1, n2, n3);
b = a; c = a; c = a; e = a; f = a; g = a; q = a; u0 = a;
aparam = 1;
itmax = int64(100); itcoun = int64(0); ndir = int64(1);
ixn = int64(0); iyn = ixn; izn = ixn;
conres = 1e-06; conchn = 1e-06;

% Set up mesh coordinates.
delta = 12.0/double(n1*(n1-1));
for i = 1:n1    
  x(i) = delta*double(i*(i-1))/2;        
end
delta = 20.0/double(n2*(n2-1));
for i = 1:n2
  y(i) = delta*double(i*(i-1))/2;
end
delta = 30.0/double(n3*(n3-1));
for i = 1:n3
  z(i) = delta*double(i*(i-1))/2;
end

% Set up difference equation coefficients, source terms and initial
% approximation.
for k = 2:n3-1
  a(2:n1-1,2:n2-1,k) = 2/((z(k)-z(k-1))*(z(k+1)-z(k-1)));
  g(2:n1-1,2:n2-1,k) = 2/((z(k+1)-z(k))*(z(k+1)-z(k-1)));
end
for j = 2:n2-1
  b(2:n1-1,j,2:n3-1) = 2/((y(j)-y(j-1))*(y(j+1)-y(j-1)));
  f(2:n1-1,j,2:n3-1) = 2/((y(j+1)-y(j))*(y(j+1)-y(j-1)));
end
for i = 2:n1-1
  c(i,2:n2-1,2:n3-1) = 2/((x(i)-x(i-1))*(x(i+1)-x(i-1)));
  e(i,2:n2-1,2:n3-1) = 2/((x(i+1)-x(i))*(x(i+1)-x(i-1)));
end
d = - a - b - c - e - f - g;    

scale = 1/y(n2);
exact = @(x,y,z)(exp((x+1)*scale)* cos(sqrt(2)*y*scale)* exp((-z-1)*scale));
q(1,   1:n2, 1:n3) = exact(x(1), y',   z);
q(n1,  1:n2, 1:n3) = exact(x(n1),y',   z);
q(1:n1,1   , 1:n3) = exact(x',   y(1), z);
q(1:n1,n2  , 1:n3) = exact(x',   y(n2),z);
q(1:n1,1:n2, 1)    = exact(x',   y,    z(1));
q(1:n1,1:n2, n3)   = exact(x',   y,    z(n3));

[u, itcoun, itused, resids, chngs, ifail] = ...
d03ec( ...
       n1, n2, n3, a, b, c, d, e, f, g, q, u0, aparam, itmax, itcoun, ...
       ndir, ixn, iyn, izn, conres, conchn);

% Output solution statistics.
fprintf('Number of iterations  %4d\n',itused);
fprintf('Maximum residual       %8.3e\n',resids(itused));

% Plot results using slices on the (x,y)-plane.
fig1 = figure;
[xs,ys,zs] = meshgrid(y,x,z);
slice(xs,ys,zs,u,[],[],double([1 3 5 7]));
axis([y(1) y(end) x(1) x(end) 0 7]);
xlabel('x'); ylabel('y'); zlabel('z');
title('Laplace''s Equation in a rectangular box');
view(16, 14);
h = colorbar;
ylabel(h, 'U(x,y,z)')

d03ec example results

Number of iterations    12
Maximum residual       1.193e-07