naginterfaces.library.pde.dim2_​ellip_​fd_​iter

naginterfaces.library.pde.dim2_ellip_fd_iter(a, b, c, d, e, aparam, it, r)

dim2_ellip_fd_iter 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,\theta )$, for example, can be used as it is equivalent to a rectangular box.)

For full information please refer to the NAG Library document for d03ua

https://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/d03/d03uaf.html

Parameters

afloat, array-like, shape $(\text{n1},\text{n2})$

$a[\text{i}-1,\text{j}-1]$ must contain the coefficient of the ‘southerly’ term involving $s_{i,j-1}$ in the $(i,j)$th equation of the system (2), for $\text{j}=1,2,\dots ,\text{n2}$, for $\text{i}=1,2,\dots ,\text{n1}$. The elements of $a$, for $j=1$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the rectangle.

bfloat, array-like, shape $(\text{n1},\text{n2})$

$b[\text{i}-1,\text{j}-1]$ must contain the coefficient of the ‘westerly’ term involving $s_{i-1,j}$ in the $(i,j)$th equation of the system (2), for $\text{j}=1,2,\dots ,\text{n2}$, for $\text{i}=1,2,\dots ,\text{n1}$. The elements of $b$, for $i=1$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the rectangle.

cfloat, array-like, shape $(\text{n1},\text{n2})$

$c[\text{i}-1,\text{j}-1]$ must contain the coefficient of the ‘central’ term involving $s_{\text{i}\text{j}}$ in the $(\text{i},\text{j})$th equation of the system (2), for $\text{j}=1,2,\dots ,\text{n2}$, for $\text{i}=1,2,\dots ,\text{n1}$. 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 $s_{\text{i}\text{j}}=r_{\text{i}\text{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 $c[\text{i}-1,\text{j}-1]=0.0$ at appropriate points. The corresponding value of $r[\text{i}-1,\text{j}-1]$ is set equal to the appropriate value, namely the difference between the prescribed value of $t_{\text{i}\text{j}}$ and the current value of $t_{\text{i}\text{j}}$ in the Dirichlet case, or zero at an external point.

dfloat, array-like, shape $(\text{n1},\text{n2})$

$d[\text{i}-1,\text{j}-1]$ must contain the coefficient of the ‘easterly’ term involving $s_{\text{i}+1,\text{j}}$ in the $(\text{i},\text{j})$th equation of the system (2), for $\text{j}=1,2,\dots ,\text{n2}$, for $\text{i}=1,2,\dots ,\text{n1}$. The elements of $d$, for $i=\text{n1}$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the rectangle.

efloat, array-like, shape $(\text{n1},\text{n2})$

$e[\text{i}-1,\text{j}-1]$ must contain the coefficient of the ‘northerly’ term involving $s_{\text{i},\text{j}+1}$ in the $(\text{i},\text{j})$th equation of the system (2), for $\text{j}=1,2,\dots ,\text{n2}$, for $\text{i}=1,2,\dots ,\text{n1}$. The elements of $e$, for $j=\text{n2}$, must be zero after incorporating the boundary conditions, since they involve nodal values from outside the rectangle.

aparamfloat

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$.

itint

The iteration number. It must be initialized, but not necessarily to $1$, before the first call, and must be incremented by one in the calling program for each subsequent call. dim2_ellip_fd_iter uses the counter to select the appropriate acceleration argument from a sequence of nine, each one being used twice in succession. (Note that the acceleration argument depends on the value of $aparam$.)

rfloat, array-like, shape $(\text{n1},\text{n2})$

$r[\text{i}-1,\text{j}-1]$ must contain the current residual $r_{\text{i}\text{j}}$ on the right-hand side of the $(\text{i},\text{j})$th equation of the system (2), for $\text{j}=1,2,\dots ,\text{n2}$, for $\text{i}=1,2,\dots ,\text{n1}$.

Returns

rfloat, ndarray, shape $(\text{n1},\text{n2})$

These residuals are overwritten by the corresponding components of solution $s$ to the system (2), i.e., the changes to be made to the vector $t$ to reduce the residuals supplied.

Raises

NagValueError

(errno $1$)

On entry, $\text{n2}=⟨value⟩$.

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

(errno $1$)

On entry, $\text{n1}=⟨value⟩$.

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

(errno $3$)

On entry, $aparam=⟨value⟩$.

Constraint: $aparam>0.0$.

(errno $4$)

On entry, $aparam=⟨value⟩$, $\text{n1}=⟨value⟩$ and $\text{n2}=⟨value⟩$.

Constraint: $aparam\le ({(\text{n1}-1)}^2+{(\text{n2}-1)}^2)/2$.

Notes

No equivalent traditional C interface for this routine exists in the NAG Library.

Given a set of simultaneous equations

$$Mt=q$$

(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^{\left(1\right)}$ by the formulae

$$\begin{array}{c}\begin{array}{ccc}{r}^{\left(n\right)}& =& q-M{t}^{\left(n\right)}\\ M{s}^{\left(n\right)}& =& r^{\left(n\right)}\\ t^{(n+1)}& =& t^{\left(n\right)}+s^{\left(n\right)}\text{.}\end{array}\end{array}$$

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.

dim2_ellip_fd_iter determines the approximate change vector $s$ corresponding to a given residual $r$, i.e., it determines an approximate solution to a set of equations

$$Ms=r$$

where $M$ is a square $(n_1\times n_2)\times (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_{ij}{s}_{i,j-1}+b_{ij}{s}_{i-1,j}+c_{ij}{s}_{ij}+d_{ij}{s}_{i+1,j}+e_{ij}{s}_{i,j+1}=r_{ij}\text{,}$$

for $i=1,2,\dots ,n_1$ and $j=1,2,\dots ,n_2$, provided that $c_{ij}\ne 0.0$. Indeed, if $c_{ij}=0.0$, then the equation is assumed to be

$$s_{ij}=r_{ij}\text{.}$$

For example, if $n_1=3$ and $n_2=2$, the equations take the form

$$\begin{array}{c}\left[\begin{array}{cccccc}{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}\right]\left[\begin{array}{c}{s}_{11}\\ s_{21}\\ s_{31}\\ s_{12}\\ s_{22}\\ s_{32}\end{array}\right]=\left[\begin{array}{c}{r}_{11}\\ r_{21}\\ r_{31}\\ r_{12}\\ r_{22}\\ r_{32}\end{array}\right]\text{.}\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 $LU$ 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_{ij}$ 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_{ij}=r_{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, one can use an array of all zeros as the initial approximation from which the first set of residuals are determined.

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)).

The function can also be used to solve each step of a time-dependent 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 positive-definiteness, of the matrix $M$ or of the update matrix formed from the arrays $a$, $b$, $c$, $d$ and $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), 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

Ames, W F, 1977, Nonlinear Partial Differential Equations in Engineering, (2nd Edition), Academic Press

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