d03ua:: Partial Differential Equations (NAG Toolbox)

nag_pde_2d_ellip_fd_iter (d03ua) performs at each call one iteration of the Strongly Implicit Procedure. It is used to calculate on successive calls a sequence of approximate corrections to the current estimate of the solution when solving a system of simultaneous algebraic equations for which the iterative update matrix is of five-point molecule form on a two-dimensional topologically-rectangular mesh. (‘Topological’ means that a polar grid

(r, θ)

, for example, can be used as it is equivalent to a rectangular box.)

Syntax

Description

Given a set of simultaneous equations

M t = q

(1)

(which could be nonlinear) derived, for example, from a finite difference representation of a two-dimensional elliptic partial differential equation and its boundary conditions, the solution

t

may be obtained iteratively from a starting approximation

t^{(1)}

by the formulae

\begin{array}{lcl} r^{(n)} & = & q - M t^{(n)} \\ M s^{(n)} & = & r^{(n)} \\ t^{(n + 1)} & = & t^{(n)} + s^{(n)} . \end{array}

Thus

r^{(n)}

is the residual of the

n

th approximate solution

t^{(n)}

, and

s^{(n)}

is the update change vector.

nag_pde_2d_ellip_fd_iter (d03ua) determines the approximate change vector

s

corresponding to a given residual

r

, i.e., it determines an approximate solution to a set of equations

M s = r

(2)

where

M

is a square

(n_{1} \times n_{2})

(n_{1} \times n_{2})

matrix and

r

is a known vector of length

n_{1} \times n_{2}

. The set of equations (2) must be of five-diagonal form

a_{i j} s_{i, j - 1} + b_{i j} s_{i - 1, j} + c_{i j} s_{i j} + d_{i j} s_{i + 1, j} + e_{i j} s_{i, j + 1} = r_{i j},

for

i = 1, 2, \dots, n_{1}

and

j = 1, 2, \dots, n_{2}

, provided that

c_{i j} \neq 0.0

. Indeed, if

c_{i j} = 0.0

, then the equation is assumed to be

s_{i j} = r_{i j} .

For example, if

n_{1} = 3

and

n_{2} = 2

, the equations take the form

[\begin{array}{l} c_{11} & d_{11} & e_{11} \\ b_{21} & c_{21} & d_{21} & e_{21} \\ b_{31} & c_{31} & e_{31} \\ a_{12} & c_{12} & d_{12} \\ a_{22} & b_{22} & c_{22} & d_{22} \\ a_{32} & b_{32} & c_{32} \end{array}] [\begin{array}{l} s_{11} \\ s_{21} \\ s_{31} \\ s_{12} \\ s_{22} \\ s_{32} \end{array}] = [\begin{array}{l} r_{11} \\ r_{21} \\ r_{31} \\ r_{12} \\ r_{22} \\ r_{32} \end{array}] .

The calling program supplies the current residual

r

at each iteration and the coefficients of the five-point molecule system of equations on which the update procedure is based. The function performs one iteration, using the approximate

L U

factorization of the Strongly Implicit Procedure with the necessary acceleration argument adjustment, to calculate the approximate solution

s

of the set of equations (2). The change

s

overwrites the residual array for return to the calling program. The calling program must combine this change stored in

r

with the old approximation to obtain the new approximate solution for

t

. It must then recalculate the residuals and, if the accuracy requirements have not been satisfied, commence the next iterative cycle.

Clearly there is no requirement that the iterative update matrix passed in the form of the five-diagonal element arrays a, b, c, d and e is the same as that used to calculate the residuals, and therefore the one governing the problem. However, the convergence may be impaired if they are not equal. Indeed, if the system of equations (1) is not precisely of the five-diagonal form illustrated above but has a few additional terms, then the methods of deferred or defect correction can be employed. The residual is calculated by the calling program using the full system of equations, but the update formula is based on a five-diagonal system (2) of the form given above. For example, the solution of a system of nine-diagonal equations each involving the combination of terms with

t_{i \pm 1, j \pm 1}, t_{i \pm 1, j}, t_{i, j \pm 1}

and

t_{i j}

could use the five-diagonal coefficients on which to base the update, provided these incorporate the major features of the equations.

Problems in topologically non-rectangular regions can be solved using the function by surrounding the region with a circumscribing topological rectangle. The equations for the nodal values external to the region of interest are set to zero (i.e.,

c_{i j} = r_{i j} = 0

) and the boundary conditions are incorporated into the equations for the appropriate nodes.

The function can be used to solve linear elliptic equations in which case the arrays a, b, c, d, e and the quantities

q

will be unchanged during the iterative cycles, or for solving nonlinear elliptic equations in which case some or all of these arrays may require updating as each new approximate solution is derived. Depending 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 (see Jacobs (1972)).

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), the calling program should subtract a typical nodal value from the whole solution

t

at every iteration to keep rounding errors to a minimum.

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

The improvement in accuracy for each iteration, i.e., on each call, depends on the size of the system and on the condition of the update matrix characterised by the five-diagonal coefficient arrays. The ultimate accuracy obtainable depends on the above factors and on the machine precision. However, since nag_pde_2d_ellip_fd_iter (d03ua) works with residuals and the update vector, the calling program can, in most cases where at each iteration all the residuals are usually of about the same size, calculate the residuals from extended precision values of the function, source term and equation coefficients if greater accuracy is required. 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. The final accuracy obtained can be judged approximately from the rate of convergence determined from the changes to the dependent variable

t

and in particular the change on the last iteration.

Further Comments

When used with deferred or defect correction, the residual is calculated in the calling program from a different system of equations to those represented by the five-point molecule coefficients used by nag_pde_2d_ellip_fd_iter (d03ua) as the basis of the iterative update procedure. When using deferred correction the overall rate of convergence depends not only on the items detailed in Accuracy but also on the difference between the two coefficient matrices used.

Example

This example solves Laplace's equation in a rectangle with a non-uniform 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

e^{(1.0 + x) / y (n_{2})} \times \cos (y / y (n_{2}))

function d03ua_example


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

n1 = int64(6);
n2 = int64(10);
nits = 10;
for i=1:n1
  x(i) = double((i*(i-1))/2);
end
for i=1:n2
  y(i) = double((i*(i-1))/2);
end

a = zeros(n1,n2);
b = a; d = a; e = a;
q = a; t = a;

for j=2:n2-1
  a(2:n1-1,j) = 2/((y(j)-y(j-1))*(y(j+1)-y(j-1)));
  e(2:n1-1,j) = 2/((y(j+1)-y(j))*(y(j+1)-y(j-1)));
end
for i=2:n1-1
  b(i,2:n2-1) = 2/((x(i)-x(i-1))*(x(i+1)-x(i-1)));
  d(i,2:n2-1) = 2/((x(i+1)-x(i))*(x(i+1)-x(i-1)));
end

c = -a - b - d - e;

% Specification for boundary nodes
q( 1,: ) = exp((x( 1)+1)/y(n2))*cos(y(: )/y(n2));
q(n1,: ) = exp((x(n1)+1)/y(n2))*cos(y(: )/y(n2));
q(: , 1) = exp((x(: )+1)/y(n2))*cos(y( 1)/y(n2));
q(: ,n2) = exp((x(: )+1)/y(n2))*cos(y(n2)/y(n2));

aparam = 1;

% Iterative loop
for it = int64(1:nits)
  [r] = resid(n1,n2,a,b,c,d,e,q,t);
  [r, ifail] = d03ua(n1, a, b, c, d, e, aparam, it, r);
  t = t + r;
end

fprintf('Final residual after %d iterations = %10.1e\n',nits,norm(r));
fprintf('\nApproximate solution is:\n\ny/x');
fprintf('%8d',1:n1);
for i=1:n2
  fprintf('\n%3d',i);
  fprintf('%8.3f',t(:,i));
end
fprintf('\n');



function [r] = resid(n1,n2,a,b,c,d,e,q,t)
  for j=1:n2
    for i=1:n1
      if (c(i,j)==0)
        r(i,j) = q(i,j) - t(i,j);
      else
        r(i,j) = q(i,j) - a(i,j)*t(i,j-1) - b(i,j)*t(i-1,j) - c(i,j)*t(i,j) ...
                        - d(i,j)*t(i+1,j) - e(i,j)*t(i,j+1);
      end
    end
  end

d03ua example results

Final residual after 10 iterations =    4.1e-10

Approximate solution is:

y/x       1       2       3       4       5       6
  1   1.022   1.045   1.093   1.168   1.277   1.427
  2   1.022   1.045   1.093   1.168   1.277   1.427
  3   1.020   1.043   1.091   1.166   1.274   1.424
  4   1.013   1.036   1.083   1.158   1.266   1.414
  5   0.997   1.020   1.066   1.140   1.246   1.392
  6   0.966   0.988   1.033   1.104   1.207   1.348
  7   0.913   0.934   0.976   1.044   1.141   1.274
  8   0.831   0.850   0.888   0.950   1.038   1.160
  9   0.712   0.728   0.762   0.814   0.890   0.994
 10   0.552   0.565   0.591   0.631   0.690   0.771

On entry,	$n1 < 2$ ,
or	$n2 < 2$ .

NAG Toolbox: nag_pde_2d_ellip_fd_iter (d03ua)

▸▿ Contents

Purpose