PDF version (NAG web site
, 64bit version, 64bit version)
NAG Toolbox: nag_pde_2d_ellip_fd (d03eb)
Purpose
nag_pde_2d_ellip_fd (d03eb) uses the Strongly Implicit Procedure to calculate the solution to a system of simultaneous algebraic equations of fivepoint molecule form on a twodimensional topologicallyrectangular mesh. (‘Topological’ means that a polar grid, for example $\left(r,\theta \right)$, can be used, being equivalent to a rectangular box.)
Syntax
[
t,
itcoun,
itused,
resids,
chngs,
ifail] = d03eb(
n1,
a,
b,
c,
d,
e,
q,
t,
aparam,
itmax,
itcoun,
ndir,
ixn,
iyn,
conres,
conchn, 'n2',
n2)
[
t,
itcoun,
itused,
resids,
chngs,
ifail] = nag_pde_2d_ellip_fd(
n1,
a,
b,
c,
d,
e,
q,
t,
aparam,
itmax,
itcoun,
ndir,
ixn,
iyn,
conres,
conchn, 'n2',
n2)
Description
Given a set of simultaneous equations
(which could be nonlinear) derived, for example, from a finite difference representation of a twodimensional elliptic partial differential equation and its boundary conditions, the function determines the values of the dependent variable
$t$.
$q$ is a known vector of length
${n}_{1}\times {n}_{2}$ and
$M$ is a square
$\left({n}_{1}\times {n}_{2}\right)$ by
$\left({n}_{1}\times {n}_{2}\right)$ matrix.
The equations must be of fivediagonal form:
for
$i=1,2,\dots ,{n}_{1}$ and
$j=1,2,\dots ,{n}_{2}$, provided
${c}_{ij}\ne 0.0$. Indeed, if
${c}_{ij}=0.0$, then the equation is assumed to be
For example, if
${n}_{1}=3$ and
${n}_{2}=2$, the equations take the form:
The system is solved iteratively, from a starting approximation
${t}^{\left(1\right)}$, by the formulae
Thus
${r}^{\left(n\right)}$ is the residual of the
$n$th approximate solution
${t}^{\left(n\right)}$, and
${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 fivepoint molecule system of equations in the arrays
a,
b,
c,
d and
e, and the source terms in the array
q. The function derives the residual of the latest approximate solution and then uses the approximate
$LU$ factorization of the Strongly Implicit Procedure with the necessary acceleration argument adjustment by calling
nag_pde_2d_ellip_fd_iter (d03ua) at each iteration.
nag_pde_2d_ellip_fd (d03eb) combines the newly derived change with the old approximation to obtain the new approximate solution for
$t$. The new solution is checked for convergence against the usersupplied convergence criteria and if these have not been achieved 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$.
Problems in topologically nonrectangular regions can be solved using the function by surrounding the region by a circumscribing topological rectangle. The equations for the nodal values external to the region of interest are set to zero (i.e., ${c}_{ij}={t}_{ij}=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, an array of all zeros can be used as the initial approximation.
The function can be used to solve linear elliptic equations in which case the arrays
a,
b,
c,
d,
e and
q are constants 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 timedependent parabolic equation in two 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 positivedefiniteness, of the matrix
$M$ formed from the arrays
a,
b,
c,
d,
e 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 Laplace's equation with all Neumann boundary conditions, a argument is incorporated so that the solution can be rescaled by subtracting a specified nodal value from the whole solution $t$ after the completion of every iteration to keep rounding errors to a minimum for those cases when the convergence is slow.
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 multidimensional partial differential equations SIAM J. Numer. Anal. 5 530–558
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{a}\left(\mathit{lda},{\mathbf{n2}}\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}\right)$ must contain the coefficient of the ‘southerly’ term involving
${t}_{\mathit{i},\mathit{j}1}$ in the
$\left(\mathit{i},\mathit{j}\right)$th equation of the system
(1), for
$\mathit{i}=1,2,\dots ,{\mathbf{n1}}$ and
$\mathit{j}=1,2,\dots ,{\mathbf{n2}}$. The elements of
a, for
$j=1$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the rectangle.
 3:
$\mathrm{b}\left(\mathit{lda},{\mathbf{n2}}\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}\right)$ must contain the coefficient of the ‘westerly’ term involving
${t}_{i1,j}$ in the
$\left(i,j\right)$th equation of the system
(1), for
$\mathit{i}=1,2,\dots ,{\mathbf{n1}}$ and
$\mathit{j}=1,2,\dots ,{\mathbf{n2}}$. The elements of
b, for
$i=1$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the rectangle.
 4:
$\mathrm{c}\left(\mathit{lda},{\mathbf{n2}}\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}\right)$ must contain the coefficient of the ‘central’ term involving
${t}_{\mathit{i}\mathit{j}}$ in the
$\left(\mathit{i},\mathit{j}\right)$th equation of the system
(1), for
$\mathit{i}=1,2,\dots ,{\mathbf{n1}}$ and
$\mathit{j}=1,2,\dots ,{\mathbf{n2}}$. The elements of
c 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}_{\mathit{i}\mathit{j}}={q}_{\mathit{i}\mathit{j}}$. 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, by setting
${\mathbf{c}}\left(\mathit{i},\mathit{j}\right)=0.0$ at appropriate points, and the corresponding value of
${\mathbf{q}}\left(\mathit{i},\mathit{j}\right)$ to the appropriate value, namely the prescribed value of
${\mathbf{t}}\left(\mathit{i},\mathit{j}\right)$ in the Dirichlet case or zero at an external point.
 5:
$\mathrm{d}\left(\mathit{lda},{\mathbf{n2}}\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}\right)$ must contain the coefficient of the ‘easterly’ term involving
${t}_{\mathit{i}+1,\mathit{j}}$ in the
$\left(\mathit{i},\mathit{j}\right)$th equation of the system
(1), for
$\mathit{i}=1,2,\dots ,{\mathbf{n1}}$ and
$\mathit{j}=1,2,\dots ,{\mathbf{n2}}$. The elements of
d, for
$i={\mathbf{n1}}$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the rectangle.
 6:
$\mathrm{e}\left(\mathit{lda},{\mathbf{n2}}\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}\right)$ must contain the coefficient of the ‘northerly’ term involving
${t}_{\mathit{i},\mathit{j}+1}$ in the
$\left(\mathit{i},\mathit{j}\right)$th equation of the system
(1), for
$\mathit{i}=1,2,\dots ,{\mathbf{n1}}$ and
$\mathit{j}=1,2,\dots ,{\mathbf{n2}}$. The elements of
e, for
$j={\mathbf{n2}}$ must be zero after incorporating the boundary conditions, since they involve nodal values from outside the rectangle.
 7:
$\mathrm{q}\left(\mathit{lda},{\mathbf{n2}}\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}\right)$ must contain
${q}_{\mathit{i}\mathit{j}}$, for
$\mathit{i}=1,2,\dots ,{\mathbf{n1}}$ and
$\mathit{j}=1,2,\dots ,{\mathbf{n2}}$, i.e., the source term values at the nodal points for the system
(1).
 8:
$\mathrm{t}\left(\mathit{lda},{\mathbf{n2}}\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}\right)$ must contain the element
${t}_{\mathit{i}\mathit{j}}$ of the approximate solution to the equations supplied by the calling program as an initial starting value, for
$\mathit{i}=1,2,\dots ,{\mathbf{n1}}$ and
$\mathit{j}=1,2,\dots ,{\mathbf{n2}}$.
If no better approximation is known, an array of zeros can be used.
 9:
$\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}\right)/2.0$.
 10:
$\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 $30$ if ${\mathbf{n1}}={\mathbf{n2}}=10$ or $100$ if ${\mathbf{n1}}={\mathbf{n2}}=50$.
 11:
$\mathrm{itcoun}$ – int64int32nag_int scalar

On the first call of
nag_pde_2d_ellip_fd (d03eb),
itcoun must be set to
$0$. On subsequent entries, its value must be unchanged from the previous call.
 12:
$\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, for example, Laplace's equation with all Neumann boundary conditions, i.e., problems in which an arbitrary constant can be added to the solution,
ndir should be set to
$0$ and the values of the next two arguments must be specified. For such problems the function subtracts the value of the function derived at the node (
ixn,
iyn) from the whole solution after each iteration to reduce the possibility of large rounding errors. You must also ensure that for such problems the appropriate consistency condition on the source terms
q is satisfied.
 13:
$\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.
 14:
$\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.
 15:
$\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.
Clearly
conres should not be less than a reasonable multiple of the
machine precision.
 16:
$\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. Clearly
conchn should not be less than a reasonable multiple of the
machine precision 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{n2}$ – int64int32nag_int scalar

Default:
the second dimension of the arrays
a,
b,
c,
d,
e,
q,
t. (An error is raised if these dimensions are not equal.)
The number of nodes in the second coordinate direction, ${n}_{2}$.
Constraint:
${\mathbf{n2}}>1$.
Output Parameters
 1:
$\mathrm{t}\left(\mathit{lda},{\mathbf{n2}}\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 and
n2 but possibly different coefficients and/or source terms, as occur with nonlinear systems or with timedependent systems,
itcoun is the accumulated number of iterations. This applies to the second and subsequent calls to
nag_pde_2d_ellip_fd (d03eb). In this way a suitable cycling of the sequence of iteration arguments is obtained in the calls to
nag_pde_2d_ellip_fd_iter (d03ua).
 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 you want to know the maximum absolute residual of the solution which is returned you must calculate this in the calling program. 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
Errors or warnings detected by the function:
 ${\mathbf{ifail}}=1$

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

On entry,  $\mathit{lda}<{\mathbf{n1}}$. 
 ${\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}\right)/2.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
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 fivediagonal 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 arguments. The convergence may become slow with very large problems, for example when
${\mathbf{n1}}={\mathbf{n2}}=60$. The final accuracy may be judged approximately from the rate of convergence determined from the sequence of values returned in
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)$.
Further Comments
The time taken per iteration is approximately proportional to ${\mathbf{n1}}\times {\mathbf{n2}}$.
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 illconditioned matrix.
Example
This example solves Laplace's equation in a rectangle with a nonuniform grid spacing in the
$x$ and
$y$ coordinate directions and with Dirichlet boundary conditions specifying the function on the perimeter of the rectangle equal to
Open in the MATLAB editor:
d03eb_example
function d03eb_example
fprintf('d03eb example results\n\n');
aparam = 1;
itmax = int64(100);
itcoun = int64(0);
ndir = int64(1);
ixn = int64(0);
iyn = int64(0);
conres = 1e06;
conchn = 1e06;
n1 = int64(31);
n2 = 46;
x = 0:0.5:15;
y = 0:1:45;
a = zeros(n1, n2);
e = a; b = a; d = a; q = a; t = a;
for j = 2:n21
a(2:n11,j) = 2/((y(j)y(j1))*(y(j+1)y(j1)));
e(2:n11,j) = 2/((y(j+1)y(j))*(y(j+1)y(j1)));
end
for i = 2:n11
b(i,2:n21) = 2/((x(i)x(i1))*(x(i+1)x(i1)));
d(i,2:n21) = 2/((x(i+1)x(i))*(x(i+1)x(i1)));
end
c = a  b  d  e;
q(1,1:n2) = exp(6*(x(1) +1.0)/y(n2))*cos(7*y(1:n2)/y(n2));
q(n1,1:n2) = exp(6*(x(n1) +1.0)/y(n2))*cos(7*y(1:n2)/y(n2));
q(1:n1,1) = exp(6*(x(1:n1)+1.0)/y(n2))*cos(7*y(1) /y(n2));
q(1:n1,n2) = exp(6*(x(1:n1)+1.0)/y(n2))*cos(7*y(n2) /y(n2));
[u, itcoun, itused, resids, chngs, ifail] = ...
d03eb( ...
n1, a, b, c, d, e, q, t, aparam, itmax, ...
itcoun, ndir, ixn, iyn, conres, conchn);
fprintf('Number of iterations %4d\n',itused);
fprintf('Maximum residual %8.3e\n',resids(itused));
fig1 = figure;
display_plot(x, y, u);
function display_plot(x, y, tOut)
meshc(y, x, tOut);
title({'Solution to Laplace''s Equation', ...
'using the Strongly Implicit Procedure', ...
'on a Fivepoint Molecule Discretisation'});
xlabel('y');
ylabel('x');
zlabel('U(x,y)');
view(35, 20);
d03eb example results
Number of iterations 29
Maximum residual 3.750e08
PDF version (NAG web site
, 64bit version, 64bit version)
© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015