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_fd (d03ec)

Purpose

nag_pde_3d_ellip_fd (d03ec) uses the Strongly Implicit Procedure to calculate the solution to a system of simultaneous algebraic equations of seven-point molecule form on a three-dimensional topologically-rectangular mesh. (‘Topological’ means that a polar grid, for example, can be used if it is equivalent to a rectangular box.)

Syntax

[t, itcoun, itused, resids, chngs, ifail] = d03ec(n1, n2, n3, a, b, c, d, e, f, g, q, t, aparam, itmax, itcoun, ndir, ixn, iyn, izn, conres, conchn, 'sda', sda)
[t, itcoun, itused, resids, chngs, ifail] = nag_pde_3d_ellip_fd(n1, n2, n3, a, b, c, d, e, f, g, q, t, aparam, itmax, itcoun, ndir, ixn, iyn, izn, conres, conchn, 'sda', sda)

Description

Given a set of simultaneous equations
 Mt = q $Mt=q$ (1)
(which could be nonlinear) derived, for example, from a finite difference representation of a three-dimensional elliptic partial differential equation and its boundary conditions, the function determines the values of the dependent variable t$t$. M$M$ is a square (n1 × n2 × n3)$\left({n}_{1}×{n}_{2}×{n}_{3}\right)$ by (n1 × n2 × n3)$\left({n}_{1}×{n}_{2}×{n}_{3}\right)$ matrix and q$q$ is a known vector of length (n1 × n2 × n3)$\left({n}_{1}×{n}_{2}×{n}_{3}\right)$.
The equations must be of seven-diagonal form:
 aijktij,k − 1 + bijkti,j − 1,k + cijkti − 1,jk + dijktijk + eijkti + 1,jk + fijkti,j + 1,k + gijktij,k + 1 = qijk $aijktij,k-1+bijkti,j-1,k+cijkti-1,jk+dijktijk+eijkti+1,jk+fijkti,j+1,k+gijktij,k+1=qijk$
for i = 1,2,,n1$i=1,2,\dots ,{n}_{1}$, j = 1,2,,n2$j=1,2,\dots ,{n}_{2}$ and k = 1,2,,n3$k=1,2,\dots ,{n}_{3}$, provided that dijk0.0${d}_{ijk}\ne 0.0$.
Indeed, if dijk = 0.0${d}_{ijk}=0.0$, then the equation is assumed to be:
 tijk = qijk. $tijk=qijk.$
The system is solved iteratively from a starting approximation t(1)${t}^{\left(1\right)}$ by the formulae:
 r(n) = q − Mt(n) Ms(n) = r(n) t(n + 1) = t(n) + s(n).
$r (n) = q-Mt(n) Ms(n) = r (n) t (n+1) = t (n) +s (n) .$
Thus r(n)${r}^{\left(n\right)}$ is the residual of the n$n$th approximate solution t(n)${t}^{\left(n\right)}$, and s(n)${s}^{\left(n\right)}$ is the update change vector.
The calling program supplies an initial approximation for the values of the dependent variable in the array t, the coefficients of the seven-point molecule system of equations in the arrays a, b, c, d, e, f and g, and the source terms in the array q. The function derives the residual of the latest approximate solution, and then uses the approximate LU$LU$ factorization of the Strongly Implicit Procedure with the necessary acceleration parameter adjustment by calling nag_pde_3d_ellip_fd_iter (d03ub) at each iteration. nag_pde_3d_ellip_fd (d03ec) combines the newly derived change with the old approximation to obtain the new approximate solution for t$t$. The new solution is checked for convergence against the user-supplied convergence criteria, and if these have not been satisfied, the iterative cycle is repeated. Convergence is based on both the maximum absolute normalized residuals (calculated with reference to the previous approximate solution as these are calculated at the commencement of each iteration) and on the maximum absolute change made to the values of t$t$.
Problems in topologically non-rectangular-box-shaped regions can be solved using the function by surrounding the region by a circumscribing topologically rectangular box. The equations for the nodal values external to the region of interest are set to zero (i.e., dijk = tijk = 0${d}_{ijk}={t}_{ijk}=0$) and the boundary conditions are incorporated into the equations for the appropriate nodes.
If there is no better initial approximation when starting the iterative cycle, one can use an array of zeros as the initial approximation.
The function can be used to solve linear elliptic equations in which case the arrays a, b, c, d, e, f, g and q remain constant and for which a single call provides the required solution. It can also be used to solve nonlinear elliptic equations, in which case some or all of these arrays may require updating during the progress of the iterations as more accurate solutions are derived. The function will then have to be called repeatedly in an outer iterative cycle. Dependent on the nonlinearity, some under-relaxation of the coefficients and/or source terms may be needed during their recalculation using the new estimates of the solution.
The function can also be used to solve each step of a time-dependent parabolic equation in three space dimensions. The solution at each time step can be expressed in terms of an elliptic equation if the Crank–Nicolson or other form of implicit time integration is used.
Neither diagonal dominance, nor positive-definiteness, of the matrix M$M$ formed from the arrays a, b, c, d, e, f and g is necessary to ensure convergence.
For problems in which the solution is not unique in the sense that an arbitrary constant can be added to the solution (for example Poisson's equation with all Neumann boundary conditions), a parameter is incorporated so that the solution can be rescaled. A specified nodal value is subtracted from the whole solution t$t$ after the completion of every iteration. This keeps rounding errors to a minimum for those cases when convergence is slow. For such problems there is generally an associated compatibility condition. For the example mentioned this compatibility condition equates the total net source within the region (i.e., the source integrated over the region) with the total net outflow across the boundaries defined by the Neumann conditions (i.e., the normal derivative integrated along the whole boundary). It is very important that the algebraic equations derived to model such a problem implement accurately the compatibility condition. If they do not, a net source or sink is very likely to be represented by the set of algebraic equations and no steady-state solution of the equations exists.

References

Jacobs D A H (1972) The strongly implicit procedure for the numerical solution of parabolic and elliptic partial differential equations Note RD/L/N66/72 Central Electricity Research Laboratory
Stone H L (1968) Iterative solution of implicit approximations of multi-dimensional partial differential equations SIAM J. Numer. Anal. 5 530–558
Weinstein H G, Stone H L and Kwan T V (1969) Iterative procedure for solution of systems of parabolic and elliptic equations in three dimensions Industrial and Engineering Chemistry Fundamentals 8 281–287

Parameters

Compulsory Input Parameters

1:     n1 – int64int32nag_int scalar
The number of nodes in the first coordinate direction, n1${n}_{1}$.
Constraint: n1 > 1${\mathbf{n1}}>1$.
2:     n2 – int64int32nag_int scalar
The number of nodes in the second coordinate direction, n2${n}_{2}$.
Constraint: n2 > 1${\mathbf{n2}}>1$.
3:     n3 – int64int32nag_int scalar
The number of nodes in the third coordinate direction, n3${n}_{3}$.
Constraint: n3 > 1${\mathbf{n3}}>1$.
4:     a(lda,sda,n3) – double array
lda, the first dimension of the array, must satisfy the constraint ldan1$\mathit{lda}\ge {\mathbf{n1}}$.
a(i,j,k)${\mathbf{a}}\left(\mathit{i},\mathit{j},\mathit{k}\right)$ must contain the coefficient of tij,k1${t}_{\mathit{i}\mathit{j},\mathit{k}-1}$ in the (i,j,k)$\left(\mathit{i},\mathit{j},\mathit{k}\right)$th equation of the system (1), for i = 1,2,,n1$\mathit{i}=1,2,\dots ,{\mathbf{n1}}$, j = 1,2,,n2$\mathit{j}=1,2,\dots ,{\mathbf{n2}}$ and k = 1,2,,n3$\mathit{k}=1,2,\dots ,{\mathbf{n3}}$. The elements of a, for k = 1$k=1$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the box.
5:     b(lda,sda,n3) – double array
lda, the first dimension of the array, must satisfy the constraint ldan1$\mathit{lda}\ge {\mathbf{n1}}$.
b(i,j,k)${\mathbf{b}}\left(i,j,k\right)$ must contain the coefficient of ti,j1,k${t}_{i,j-1,k}$ in the (i,j,k)$\left(i,j,k\right)$th equation of the system (1), for i = 1,2,,n1$\mathit{i}=1,2,\dots ,{\mathbf{n1}}$, j = 1,2,,n2$\mathit{j}=1,2,\dots ,{\mathbf{n2}}$ and k = 1,2,,n3$\mathit{k}=1,2,\dots ,{\mathbf{n3}}$. The elements of b, for j = 1$j=1$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the box.
6:     c(lda,sda,n3) – double array
lda, the first dimension of the array, must satisfy the constraint ldan1$\mathit{lda}\ge {\mathbf{n1}}$.
c(i,j,k)${\mathbf{c}}\left(\mathit{i},\mathit{j},\mathit{k}\right)$ must contain the coefficient of ti1,jk${t}_{\mathit{i}-1,\mathit{j}\mathit{k}}$ in the (i,j,k)$\left(\mathit{i},\mathit{j},\mathit{k}\right)$th equation of the system (1), for i = 1,2,,n1$\mathit{i}=1,2,\dots ,{\mathbf{n1}}$, j = 1,2,,n2$\mathit{j}=1,2,\dots ,{\mathbf{n2}}$ and k = 1,2,,n3$\mathit{k}=1,2,\dots ,{\mathbf{n3}}$. The elements of c, for i = 1$i=1$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the box.
7:     d(lda,sda,n3) – double array
lda, the first dimension of the array, must satisfy the constraint ldan1$\mathit{lda}\ge {\mathbf{n1}}$.
d(i,j,k)${\mathbf{d}}\left(i,j,k\right)$ must contain the coefficient of tijk${t}_{ijk}$ (the ‘central’ term) in the (i,j,k)$\left(i,j,k\right)$th equation of the system (1), for i = 1,2,,n1$\mathit{i}=1,2,\dots ,{\mathbf{n1}}$, j = 1,2,,n2$\mathit{j}=1,2,\dots ,{\mathbf{n2}}$ and k = 1,2,,n3$\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 tijk = qijk${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 d(i,j,k) = 0.0${\mathbf{d}}\left(i,j,k\right)=0.0$ at appropriate points, and the corresponding value of q(i,j,k)${\mathbf{q}}\left(i,j,k\right)$ to the appropriate value, namely the prescribed value of t(i,j,k)${\mathbf{t}}\left(i,j,k\right)$ in the Dirichlet case, or to zero at an external point.
8:     e(lda,sda,n3) – double array
lda, the first dimension of the array, must satisfy the constraint ldan1$\mathit{lda}\ge {\mathbf{n1}}$.
e(i,j,k)${\mathbf{e}}\left(\mathit{i},\mathit{j},\mathit{k}\right)$ must contain the coefficient of ti + 1,jk${t}_{\mathit{i}+1,\mathit{j}\mathit{k}}$ in the (i,j,k)$\left(\mathit{i},\mathit{j},\mathit{k}\right)$th equation of the system (1), for i = 1,2,,n1$\mathit{i}=1,2,\dots ,{\mathbf{n1}}$, j = 1,2,,n2$\mathit{j}=1,2,\dots ,{\mathbf{n2}}$ and k = 1,2,,n3$\mathit{k}=1,2,\dots ,{\mathbf{n3}}$. The elements of e, for i = n1$i={\mathbf{n1}}$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the box.
9:     f(lda,sda,n3) – double array
lda, the first dimension of the array, must satisfy the constraint ldan1$\mathit{lda}\ge {\mathbf{n1}}$.
f(i,j,k)${\mathbf{f}}\left(\mathit{i},\mathit{j},\mathit{k}\right)$ must contain the coefficient of ti,j + 1,k${t}_{\mathit{i},\mathit{j}+1,\mathit{k}}$ in the (i,j,k)$\left(\mathit{i},\mathit{j},\mathit{k}\right)$th equation of the system (1), , for i = 1,2,,n1$\mathit{i}=1,2,\dots ,{\mathbf{n1}}$, j = 1,2,,n2$\mathit{j}=1,2,\dots ,{\mathbf{n2}}$ and k = 1,2,,n3$\mathit{k}=1,2,\dots ,{\mathbf{n3}}$. The elements of f, for j = n2$j={\mathbf{n2}}$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the box.
10:   g(lda,sda,n3) – double array
lda, the first dimension of the array, must satisfy the constraint ldan1$\mathit{lda}\ge {\mathbf{n1}}$.
g(i,j,k)${\mathbf{g}}\left(\mathit{i},\mathit{j},\mathit{k}\right)$ must contain the coefficient of tij,k + 1${t}_{\mathit{i}\mathit{j},\mathit{k}+1}$ in the (i,j,k)$\left(\mathit{i},\mathit{j},\mathit{k}\right)$th equation of the system (1), for i = 1,2,,n1$\mathit{i}=1,2,\dots ,{\mathbf{n1}}$, j = 1,2,,n2$\mathit{j}=1,2,\dots ,{\mathbf{n2}}$ and k = 1,2,,n3$\mathit{k}=1,2,\dots ,{\mathbf{n3}}$. The elements of g, for k = n3$k={\mathbf{n3}}$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the box.
11:   q(lda,sda,n3) – double array
lda, the first dimension of the array, must satisfy the constraint ldan1$\mathit{lda}\ge {\mathbf{n1}}$.
q(i,j,k)${\mathbf{q}}\left(\mathit{i},\mathit{j},\mathit{k}\right)$ must contain qijk${q}_{\mathit{i}\mathit{j}\mathit{k}}$, for i = 1,2,,n1$\mathit{i}=1,2,\dots ,{\mathbf{n1}}$, j = 1,2,,n2$\mathit{j}=1,2,\dots ,{\mathbf{n2}}$ and k = 1,2,,n3$\mathit{k}=1,2,\dots ,{\mathbf{n3}}$, i.e., the source-term values at the nodal points of the system (1).
12:   t(lda,sda,n3) – double array
lda, the first dimension of the array, must satisfy the constraint ldan1$\mathit{lda}\ge {\mathbf{n1}}$.
t(i,j,k)${\mathbf{t}}\left(\mathit{i},\mathit{j},\mathit{k}\right)$ must contain the element tijk${t}_{\mathit{i}\mathit{j}\mathit{k}}$ of an approximate solution to the equations, for i = 1,2,,n1$\mathit{i}=1,2,\dots ,{\mathbf{n1}}$, j = 1,2,,n2$\mathit{j}=1,2,\dots ,{\mathbf{n2}}$ and k = 1,2,,n3$\mathit{k}=1,2,\dots ,{\mathbf{n3}}$.
If no better approximation is known, an array of zeros can be used.
13:   aparam – double scalar
The iteration acceleration factor. A value of 1.0$1.0$ is adequate for most typical problems. However, if convergence is slow, the value can be reduced, typically to 0.2$0.2$ or 0.1$0.1$. If divergence is obtained, the value can be increased, typically to 2.0$2.0$, 5.0$5.0$ or 10.0$10.0$.
Constraint: 0.0 < aparam((n11)2 + (n21)2 + (n31)2) / 3.0$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:   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$20$ for a problem with 3000$3000$ nodes and convergence criteria of about 103${10}^{-3}$ of the original residual and change.
15:   itcoun – int64int32nag_int scalar
On the first call of nag_pde_3d_ellip_fd (d03ec), itcoun must be set to 0$0$. On subsequent entries, its value must be unchanged from the previous call.
16:   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$0$ and the values of the next three parameters 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 Section [Description].
17:   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:   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:   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:   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$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:   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:     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 as declared in the (sub)program from which nag_pde_3d_ellip_fd (d03ec) is called.
Constraint: ${\mathbf{sda}}\ge {\mathbf{n2}}$.

Input Parameters Omitted from the MATLAB Interface

lda wrksp1 wrksp2 wrksp3 wrksp4

Output Parameters

1:     t(lda,sda,n3) – double array
ldan1$\mathit{lda}\ge {\mathbf{n1}}$.
The solution derived by the function.
2:     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 parameters is obtained in the calls to nag_pde_3d_ellip_fd_iter (d03ub).
3:     itused – int64int32nag_int scalar
The number of iterations actually used on that call.
4:     resids(itmax) – double array
The maximum absolute value of the residuals calculated at the i$\mathit{i}$th iteration, for i = 1,2,,itused$\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:     chngs(itmax) – double array
The maximum absolute value of the changes made to the components of the dependent variable t at the i$\mathit{i}$th iteration, for i = 1,2,,itused$\mathit{i}=1,2,\dots ,{\mathbf{itused}}$. The sequence of values chngs indicates the rate of convergence.
6:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Note: nag_pde_3d_ellip_fd (d03ec) may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
 On entry, n1 < 2${\mathbf{n1}}<2$, or n2 < 2${\mathbf{n2}}<2$, or n3 < 2${\mathbf{n3}}<2$.
ifail = 2${\mathbf{ifail}}=2$
 On entry, lda < n1$\mathit{lda}<{\mathbf{n1}}$, or ${\mathbf{sda}}<{\mathbf{n2}}$.
ifail = 3${\mathbf{ifail}}=3$
 On entry, aparam ≤ 0.0${\mathbf{aparam}}\le 0.0$.
ifail = 4${\mathbf{ifail}}=4$
 On entry, aparam > ((n1 − 1)2 + (n2 − 1)2 + (n3 − 1)2) / 3.0${\mathbf{aparam}}>\left({\left({\mathbf{n1}}-1\right)}^{2}+{\left({\mathbf{n2}}-1\right)}^{2}+{\left({\mathbf{n3}}-1\right)}^{2}\right)/3.0$.
ifail = 5${\mathbf{ifail}}=5$
Convergence was not achieved after itmax iterations.

Accuracy

The improvement in accuracy for each iteration depends on the size of the system and on the condition of the update matrix characterised by the seven-diagonal coefficient arrays. The ultimate accuracy obtainable depends on the above factors and on the machine precision. The rate of convergence obtained with the Strongly Implicit Procedure is not always smooth because of the cyclic use of nine acceleration parameters. The convergence may become slow with very large problems. The final accuracy obtained may be judged approximately from the rate of convergence determined from the sequence of values returned in the arrays resids and chngs and the magnitude of the maximum absolute value of the change vector on the last iteration stored in ${\mathbf{chngs}}\left({\mathbf{itused}}\right)$.

The time taken per iteration is approximately proportional to n1 × n2 × n3${\mathbf{n1}}×{\mathbf{n2}}×{\mathbf{n3}}$.
Convergence may not always be obtained when the problem is very large and/or the coefficients of the equations have widely disparate values. The latter case is often associated with a near ill-conditioned matrix.

Example

This example solves Laplace's equation in a rectangular box with a non-uniform grid spacing in the x$x$, y$y$, and z$z$ coordinate directions and with Dirichlet boundary conditions specifying the function on the surfaces of the box equal to
 e(1.0 + x) / y(n2) × cos(sqrt(2)y / y(n2)) × e( − 1.0 − z) / y(n2). $e(1.0+x)/y(n2)×cos(2y/y(n2))×e(-1.0-z)/y(n2).$
Note that this is the same problem as that solved in the example for nag_pde_3d_ellip_fd_iter (d03ub). The differences in the maximum residuals obtained at each iteration between the two test runs are explained by the fact that in nag_pde_3d_ellip_fd (d03ec) the residual at each node is normalized by dividing by the central coefficient, whereas this normalization has not been used in the example program for nag_pde_3d_ellip_fd_iter (d03ub).
```function nag_pde_3d_ellip_fd_example
fprintf('nag_pde_3d_ellip_fd example program results\n');

% We run through the calculation twice: once with a coarse mesh, when we
% output the results, and once with a finer mesh, when we collect the
% results for plotting.
for icalc = 1:2
if icalc == 1
n1 = 4;
n2 = 5;
n3 = 6;
else
n1 = 16;
n2 = 20;
n3 = 24;
end

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

% Set up coordinate arrays.
delta = 12.0/(n1*(n1-1));
for i = 1:n1
x(i) = delta*(i*(i-1))/2.0;
end
delta = 20.0/(n2*(n2-1));
for i = 1:n2
y(i) = delta*(i*(i-1))/2.0;
end
delta = 30.0/(n3*(n3-1));
for i = 1:n3
z(i) = delta*(i*(i-1))/2.0;
end

% Set up difference equation coefficients, source terms and initial
% approximation.
for k = 1:n3
for j = 1:n2
for i = 1:n1
if (i ~= 1 && i ~= n1 && j ~= 1 && j ~= n2 && ...
k ~= 1 && k ~= n3)
% Specification for internal nodes
a(i,j,k) = 2/((z(k)-z(k-1))*(z(k+1)-z(k-1)));
g(i,j,k) = 2/((z(k+1)-z(k))*(z(k+1)-z(k-1)));
b(i,j,k) = 2/((y(j)-y(j-1))*(y(j+1)-y(j-1)));
f(i,j,k) = 2/((y(j+1)-y(j))*(y(j+1)-y(j-1)));
c(i,j,k) = 2/((x(i)-x(i-1))*(x(i+1)-x(i-1)));
e(i,j,k) = 2/((x(i+1)-x(i))*(x(i+1)-x(i-1)));
d(i,j,k) = - a(i,j,k) - b(i,j,k) - c(i,j,k) - ...
e(i,j,k) - f(i,j,k) - g(i,j,k);
else
% Specification for boundary nodes
q(i,j,k) = exp((x(i)+1)/y(n2)) * ...
cos(sqrt(2)*y(j)/y(n2)) * exp((-z(k)-1)/y(n2));
end
end
end
end

[tOut, itcounOut, itused, resids, chngs, ifail] = nag_pde_3d_ellip_fd(int64(n1), ...
int64(n2), int64(n3), a, b, c, d, e, f, g, q, t, aparam, ...
int64(itmax), int64(itcoun), int64(ndir), int64(ixn), ...
int64(iyn), int64(izn), conres, conchn);
if ifail ~= 0
% Parameters out of range, or convergence problems.
% Print message and exit.
error('Warning: nag_pde_3d_ellip_fd returned with ifail = %1d ',ifail);
end

% Output results first time through (low resolution).
if icalc == 1
fprintf(['\nIteration      Maximum         Maximum\n', ...
' number        residual         change\n']);

for i = 1:int32(itused) % can't use int64 in loop range.
fprintf('    %1.0f         %1.4e     %1.4e\n', i, ...
resids(i), chngs(i));
end

fprintf('\nTable of calculated function values\n\n');
fprintf(['K  J   (I     T   )  (I     T   )  (I     T   )  ', ...
'(I     T   )\n']);

for  j = 1:n3
for k = 1:n2
fprintf('%1d  %1d   ', j, k)
for i = 1:n1
fprintf(' %1d    %1.3f   ', i, tOut(w(i), k, j))
end
fprintf(' \n');
end
end
end
end
% Plot results.
fig = figure('Number', 'off');
display_plot_iso(x, y, z, tOut, 0.2);
fig = figure('Number', 'off');
display_plot_mesh(x, y, z, tOut, 16);
fig = figure('Number', 'off');
display_plot(x, y, z, tOut, 1, [1 7 11 16]);
fig = figure('Number', 'off');
display_plot(x, y, z, tOut, 24, [1 7 11 16]);

function display_plot(x, y, z, tOut, zindex, xindex)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% Plot the function against y, with z fixed (using the given z array
% index), and for a range of x values (given by the array of x array
% indices).
plot(y, tOut(xindex,:,zindex));
% Prepare labels for the title and the legend, using the array of x array
% indices.
titleStr = '';
for ix = 1:length(xindex)
legendStr{ix} = ['X = ', num2str(x(xindex(ix)))];
titleStr = [titleStr, num2str(x(xindex(ix)))];
if ix < length(xindex)
titleStr = [titleStr, ', '];
end
end
% Set the title.
title(['Solution Profiles f(x=X,y,z=', num2str(z(zindex)), ...
') for X = ', titleStr], titleFmt{:});
% Label the axes
xlabel('y', labFmt{:});
ylabel(['f(x=X,y,z=', num2str(z(zindex)), ')'], labFmt{:});
legend(legendStr, 'Location', 'Best');

function display_plot_iso(x1d, y1d, z1d, fun, fval)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% x1d, y1d & z1d are the 1D vectors giving the coordinates in x y & z
% corresponding to the data values in fun (a 3D array).  The isosurface
% command requires 3D coordinates arrays, which must be produced from the
% meshgrid command.  Not sure why this requires the x1d and y1d arrays to
% be swapped round, but the dimensions of x3d etc won't match up with those
% of fun unless we do this.
[x3d, y3d, z3d] = meshgrid(y1d, x1d, z1d);
p = patch(isosurface(x3d, y3d, z3d, fun, fval));
isonormals(x3d, y3d, z3d, fun, p)
set(p, 'FaceColor', 'blue', 'EdgeColor', 'none');
% Set the relative scaling of the data units along the axes to be equal.
daspect([1 1 1]);
% Specify explicitly the limits of the axes.
axis([y1d(1) y1d(end) x1d(1) x1d(end) z1d(1) z1d(end)]);
camlight
lighting gouraud
% Set the title.
title(['Solution in xyz space for f =', num2str(fval)], titleFmt{:});
% Label the axes.  Note the comment above about swapping y and x.
xlabel('y', labFmt{:});
ylabel('x', labFmt{:});
zlabel('z', labFmt{:});
% Set the view to something nice (determined empirically).
view(25, 80);
function display_plot_mesh(x, y, z, tOut, ix)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% Plot a slice through the 3D array as a mesh and contours.  We use the
% given array index to select the X value for the slice (which makes a 1
% by M by N array), then use the squeeze function to convert this to a M
% by N array.
meshc(z, y, squeeze(tOut(ix,:,:)))
% Set the title.
title(['Solution in yz plane for x=', num2str(x(ix))], titleFmt{:});
% Label the axes.
xlabel('z', labFmt{:});
ylabel('y', labFmt{:});
zlabel(['f(x=', num2str(x(ix)), 'y,z)'], labFmt{:})
% Set the view to something nice (determined empirically).
view(60, 40);
```
```
nag_pde_3d_ellip_fd example program results

Iteration      Maximum         Maximum
number        residual         change
1         1.8221e+00     1.8221e+00
2         9.0252e-03     1.9701e-02
3         1.3576e-03     1.4961e-03
4         4.0131e-05     3.8483e-05
5         5.3208e-06     5.4813e-06
6         2.6951e-07     2.3331e-07

Table of calculated function values

K  J   (I     T   )  (I     T   )  (I     T   )  (I     T   )
1  1    1    1.000    2    1.105    3    1.350    4    1.822
1  2    1    0.990    2    1.094    3    1.336    4    1.804
1  3    1    0.911    2    1.007    3    1.230    4    1.661
1  4    1    0.661    2    0.731    3    0.892    4    1.205
1  5    1    0.156    2    0.172    3    0.211    4    0.284
2  1    1    0.905    2    1.000    3    1.221    4    1.649
2  2    1    0.896    2    0.990    3    1.210    4    1.632
2  3    1    0.825    2    0.912    3    1.114    4    1.503
2  4    1    0.598    2    0.662    3    0.809    4    1.090
2  5    1    0.141    2    0.156    3    0.190    4    0.257
3  1    1    0.741    2    0.819    3    1.000    4    1.350
3  2    1    0.733    2    0.811    3    0.991    4    1.336
3  3    1    0.675    2    0.747    3    0.913    4    1.230
3  4    1    0.490    2    0.543    3    0.664    4    0.892
3  5    1    0.116    2    0.128    3    0.156    4    0.211
4  1    1    0.549    2    0.607    3    0.741    4    1.000
4  2    1    0.543    2    0.601    3    0.734    4    0.990
4  3    1    0.500    2    0.554    3    0.677    4    0.911
4  4    1    0.363    2    0.402    3    0.492    4    0.661
4  5    1    0.086    2    0.095    3    0.116    4    0.156
5  1    1    0.368    2    0.407    3    0.497    4    0.670
5  2    1    0.364    2    0.403    3    0.492    4    0.664
5  3    1    0.335    2    0.371    3    0.454    4    0.611
5  4    1    0.243    2    0.270    3    0.330    4    0.443
5  5    1    0.057    2    0.063    3    0.077    4    0.105
6  1    1    0.223    2    0.247    3    0.301    4    0.407
6  2    1    0.221    2    0.244    3    0.298    4    0.403
6  3    1    0.203    2    0.225    3    0.274    4    0.371
6  4    1    0.148    2    0.163    3    0.199    4    0.269
6  5    1    0.035    2    0.038    3    0.047    4    0.063

```
```function d03ec_example
fprintf('d03ec example program results\n');

% We run through the calculation twice: once with a coarse mesh, when we
% output the results, and once with a finer mesh, when we collect the
% results for plotting.
for icalc = 1:2
if icalc == 1
n1 = 4;
n2 = 5;
n3 = 6;
else
n1 = 16;
n2 = 20;
n3 = 24;
end

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

% Set up coordinate arrays.
delta = 12.0/(n1*(n1-1));
for i = 1:n1
x(i) = delta*(i*(i-1))/2.0;
end
delta = 20.0/(n2*(n2-1));
for i = 1:n2
y(i) = delta*(i*(i-1))/2.0;
end
delta = 30.0/(n3*(n3-1));
for i = 1:n3
z(i) = delta*(i*(i-1))/2.0;
end

% Set up difference equation coefficients, source terms and initial
% approximation.
for k = 1:n3
for j = 1:n2
for i = 1:n1
if (i ~= 1 && i ~= n1 && j ~= 1 && j ~= n2 && ...
k ~= 1 && k ~= n3)
% Specification for internal nodes
a(i,j,k) = 2/((z(k)-z(k-1))*(z(k+1)-z(k-1)));
g(i,j,k) = 2/((z(k+1)-z(k))*(z(k+1)-z(k-1)));
b(i,j,k) = 2/((y(j)-y(j-1))*(y(j+1)-y(j-1)));
f(i,j,k) = 2/((y(j+1)-y(j))*(y(j+1)-y(j-1)));
c(i,j,k) = 2/((x(i)-x(i-1))*(x(i+1)-x(i-1)));
e(i,j,k) = 2/((x(i+1)-x(i))*(x(i+1)-x(i-1)));
d(i,j,k) = - a(i,j,k) - b(i,j,k) - c(i,j,k) - ...
e(i,j,k) - f(i,j,k) - g(i,j,k);
else
% Specification for boundary nodes
q(i,j,k) = exp((x(i)+1)/y(n2)) * ...
cos(sqrt(2)*y(j)/y(n2)) * exp((-z(k)-1)/y(n2));
end
end
end
end

[tOut, itcounOut, itused, resids, chngs, ifail] = d03ec(int64(n1), ...
int64(n2), int64(n3), a, b, c, d, e, f, g, q, t, aparam, ...
int64(itmax), int64(itcoun), int64(ndir), int64(ixn), ...
int64(iyn), int64(izn), conres, conchn);
if ifail ~= 0
% Parameters out of range, or convergence problems.
% Print message and exit.
error('Warning: d03ec returned with ifail = %1d ',ifail);
end

% Output results first time through (low resolution).
if icalc == 1
fprintf(['\nIteration      Maximum         Maximum\n', ...
' number        residual         change\n']);

for i = 1:int32(itused) % can't use int64 in loop range.
fprintf('    %1.0f         %1.4e     %1.4e\n', i, ...
resids(i), chngs(i));
end

fprintf('\nTable of calculated function values\n\n');
fprintf(['K  J   (I     T   )  (I     T   )  (I     T   )  ', ...
'(I     T   )\n']);

for  j = 1:n3
for k = 1:n2
fprintf('%1d  %1d   ', j, k)
for i = 1:n1
fprintf(' %1d    %1.3f   ', i, tOut(w(i), k, j))
end
fprintf(' \n');
end
end
end
end
% Plot results.
fig = figure('Number', 'off');
display_plot_iso(x, y, z, tOut, 0.2);
fig = figure('Number', 'off');
display_plot_mesh(x, y, z, tOut, 16);
fig = figure('Number', 'off');
display_plot(x, y, z, tOut, 1, [1 7 11 16]);
fig = figure('Number', 'off');
display_plot(x, y, z, tOut, 24, [1 7 11 16]);

function display_plot(x, y, z, tOut, zindex, xindex)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% Plot the function against y, with z fixed (using the given z array
% index), and for a range of x values (given by the array of x array
% indices).
plot(y, tOut(xindex,:,zindex));
% Prepare labels for the title and the legend, using the array of x array
% indices.
titleStr = '';
for ix = 1:length(xindex)
legendStr{ix} = ['X = ', num2str(x(xindex(ix)))];
titleStr = [titleStr, num2str(x(xindex(ix)))];
if ix < length(xindex)
titleStr = [titleStr, ', '];
end
end
% Set the title.
title(['Solution Profiles f(x=X,y,z=', num2str(z(zindex)), ...
') for X = ', titleStr], titleFmt{:});
% Label the axes
xlabel('y', labFmt{:});
ylabel(['f(x=X,y,z=', num2str(z(zindex)), ')'], labFmt{:});
legend(legendStr, 'Location', 'Best');

function display_plot_iso(x1d, y1d, z1d, fun, fval)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% x1d, y1d & z1d are the 1D vectors giving the coordinates in x y & z
% corresponding to the data values in fun (a 3D array).  The isosurface
% command requires 3D coordinates arrays, which must be produced from the
% meshgrid command.  Not sure why this requires the x1d and y1d arrays to
% be swapped round, but the dimensions of x3d etc won't match up with those
% of fun unless we do this.
[x3d, y3d, z3d] = meshgrid(y1d, x1d, z1d);
p = patch(isosurface(x3d, y3d, z3d, fun, fval));
isonormals(x3d, y3d, z3d, fun, p)
set(p, 'FaceColor', 'blue', 'EdgeColor', 'none');
% Set the relative scaling of the data units along the axes to be equal.
daspect([1 1 1]);
% Specify explicitly the limits of the axes.
axis([y1d(1) y1d(end) x1d(1) x1d(end) z1d(1) z1d(end)]);
camlight
lighting gouraud
% Set the title.
title(['Solution in xyz space for f =', num2str(fval)], titleFmt{:});
% Label the axes.  Note the comment above about swapping y and x.
xlabel('y', labFmt{:});
ylabel('x', labFmt{:});
zlabel('z', labFmt{:});
% Set the view to something nice (determined empirically).
view(25, 80);
function display_plot_mesh(x, y, z, tOut, ix)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% Plot a slice through the 3D array as a mesh and contours.  We use the
% given array index to select the X value for the slice (which makes a 1
% by M by N array), then use the squeeze function to convert this to a M
% by N array.
meshc(z, y, squeeze(tOut(ix,:,:)))
% Set the title.
title(['Solution in yz plane for x=', num2str(x(ix))], titleFmt{:});
% Label the axes.
xlabel('z', labFmt{:});
ylabel('y', labFmt{:});
zlabel(['f(x=', num2str(x(ix)), 'y,z)'], labFmt{:})
% Set the view to something nice (determined empirically).
view(60, 40);
```
```
d03ec example program results

Iteration      Maximum         Maximum
number        residual         change
1         1.8221e+00     1.8221e+00
2         9.0252e-03     1.9701e-02
3         1.3576e-03     1.4961e-03
4         4.0131e-05     3.8483e-05
5         5.3208e-06     5.4813e-06
6         2.6951e-07     2.3331e-07

Table of calculated function values

K  J   (I     T   )  (I     T   )  (I     T   )  (I     T   )
1  1    1    1.000    2    1.105    3    1.350    4    1.822
1  2    1    0.990    2    1.094    3    1.336    4    1.804
1  3    1    0.911    2    1.007    3    1.230    4    1.661
1  4    1    0.661    2    0.731    3    0.892    4    1.205
1  5    1    0.156    2    0.172    3    0.211    4    0.284
2  1    1    0.905    2    1.000    3    1.221    4    1.649
2  2    1    0.896    2    0.990    3    1.210    4    1.632
2  3    1    0.825    2    0.912    3    1.114    4    1.503
2  4    1    0.598    2    0.662    3    0.809    4    1.090
2  5    1    0.141    2    0.156    3    0.190    4    0.257
3  1    1    0.741    2    0.819    3    1.000    4    1.350
3  2    1    0.733    2    0.811    3    0.991    4    1.336
3  3    1    0.675    2    0.747    3    0.913    4    1.230
3  4    1    0.490    2    0.543    3    0.664    4    0.892
3  5    1    0.116    2    0.128    3    0.156    4    0.211
4  1    1    0.549    2    0.607    3    0.741    4    1.000
4  2    1    0.543    2    0.601    3    0.734    4    0.990
4  3    1    0.500    2    0.554    3    0.677    4    0.911
4  4    1    0.363    2    0.402    3    0.492    4    0.661
4  5    1    0.086    2    0.095    3    0.116    4    0.156
5  1    1    0.368    2    0.407    3    0.497    4    0.670
5  2    1    0.364    2    0.403    3    0.492    4    0.664
5  3    1    0.335    2    0.371    3    0.454    4    0.611
5  4    1    0.243    2    0.270    3    0.330    4    0.443
5  5    1    0.057    2    0.063    3    0.077    4    0.105
6  1    1    0.223    2    0.247    3    0.301    4    0.407
6  2    1    0.221    2    0.244    3    0.298    4    0.403
6  3    1    0.203    2    0.225    3    0.274    4    0.371
6  4    1    0.148    2    0.163    3    0.199    4    0.269
6  5    1    0.035    2    0.038    3    0.047    4    0.063

```