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_2d_ellip_mgrid (d03ed)

## Purpose

nag_pde_2d_ellip_mgrid (d03ed) solves seven-diagonal systems of linear equations which arise from the discretization of an elliptic partial differential equation on a rectangular region. This function uses a multigrid technique.

## Syntax

[a, rhs, ub, us, u, numit, ifail] = d03ed(ngx, ngy, a, rhs, ub, maxit, acc, iout)
[a, rhs, ub, us, u, numit, ifail] = nag_pde_2d_ellip_mgrid(ngx, ngy, a, rhs, ub, maxit, acc, iout)

## Description

nag_pde_2d_ellip_mgrid (d03ed) solves, by multigrid iteration, the seven-point scheme
 Ai,j6ui − 1,j + 1 + Ai,j7ui,j + 1 + Ai,j3ui − 1,j + Ai,j4uij + Ai,j5ui + 1,j + Ai,j1ui,j − 1 + Ai,j2ui + 1,j − 1 = fij,  i = 1,2, … ,nx​ and ​j = 1,2, … ,ny,
$Ai,j6ui-1,j+1 + Ai,j7ui,j+1 + Ai,j3ui-1,j + Ai,j4uij + Ai,j5ui+1,j + Ai,j1ui,j-1 + Ai,j2ui+1,j-1=fij, i=1,2,…,nx​ and ​j=1,2,…,ny,$
which arises from the discretization of an elliptic partial differential equation of the form
 α (x,y)Uxx + β (x,y)Uxy + γ (x,y)Uyy + δ (x,y)Ux + ε (x,y)Uy + φ (x,y)U = ψ (x,y) $α (x,y)Uxx+β (x,y)Uxy+γ (x,y)Uyy+δ (x,y)Ux+ε (x,y)Uy+ϕ (x,y)U=ψ (x,y)$
and its boundary conditions, defined on a rectangular region. This we write in matrix form as
 Au = f. $Au=f.$
The algorithm is described in separate reports by Wesseling (1982a), Wesseling (1982b) and McCarthy (1983).
Systems of linear equations, matching the seven-point stencil defined above, are solved by a multigrid iteration. An initial estimate of the solution must be provided by you. A zero guess may be supplied if no better approximation is available.
A ‘smoother’ based on incomplete Crout decomposition is used to eliminate the high frequency components of the error. A restriction operator is then used to map the system on to a sequence of coarser grids. The errors are then smoothed and prolongated (mapped onto successively finer grids). When the finest cycle is reached, the approximation to the solution is corrected. The cycle is repeated for maxit iterations or until the required accuracy, acc, is reached.
nag_pde_2d_ellip_mgrid (d03ed) will automatically determine the number l$l$ of possible coarse grids, ‘levels’ of the multigrid scheme, for a particular problem. In other words, nag_pde_2d_ellip_mgrid (d03ed) determines the maximum integer l$l$ so that nx${n}_{x}$ and ny${n}_{y}$ can be expressed in the form
 nx = m2l − 1 + 1 ,   ny = n2l − 1 + 1 ,   with   m ≥ 2 ​ and ​ n ≥ 2 . $nx = m2l-1 + 1 , ny = n2l-1 + 1 , with m≥2 ​ and ​ n≥2 .$
It should be noted that the rate of convergence improves significantly with the number of levels used (see McCarthy (1983)), so that nx${n}_{x}$ and ny${n}_{y}$ should be carefully chosen so that nx1${n}_{x}-1$ and ny1${n}_{y}-1$ have factors of the form 2l${2}^{l}$, with l$l$ as large as possible. For good convergence the integer l$l$ should be at least 2$2$.
nag_pde_2d_ellip_mgrid (d03ed) has been found to be robust in application, but being an iterative method the problem of divergence can arise. For a strictly diagonally dominant matrix A$A$
 |Aij4| > ∑k ≠ 4 |Aijk| ,   i = 1,2, … ,nx ​ and ​ j = 1,2, … ,ny $| Aij4 | > ∑k≠4 | Aijk | , i=1,2,…,nx ​ and ​ j=1,2,…,ny$
no such problem is foreseen. The diagonal dominance of A$A$ is not a necessary condition, but should this condition be strongly violated then divergence may occur. The quickest test is to try the function.

## References

McCarthy G J (1983) Investigation into the multigrid code MGD1 Report AERE-R 10889 Harwell
Wesseling P (1982a) MGD1 – a robust and efficient multigrid method Multigrid Methods. Lecture Notes in Mathematics 960 614–630 Springer–Verlag
Wesseling P (1982b) Theoretical aspects of a multigrid method SIAM J. Sci. Statist. Comput. 3 387–407

## Parameters

### Compulsory Input Parameters

1:     ngx – int64int32nag_int scalar
The number of interior grid points in the x$x$-direction, nx${n}_{x}$. ngx1${\mathbf{ngx}}-1$ should preferably be divisible by as high a power of 2$2$ as possible.
Constraint: ngx3${\mathbf{ngx}}\ge 3$.
2:     ngy – int64int32nag_int scalar
The number of interior grid points in the y$y$-direction, ny${n}_{y}$. ngy1${\mathbf{ngy}}-1$ should preferably be divisible by as high a power of 2$2$ as possible.
Constraint: ngy3${\mathbf{ngy}}\ge 3$.
3:     a(lda,7$7$) – double array
a(i + (j1) × ngx,k)${\mathbf{a}}\left(\mathit{i}+\left(\mathit{j}-1\right)×{\mathbf{ngx}},\mathit{k}\right)$ must be set to aijk${{\mathbf{a}}}_{\mathit{i}\mathit{j}}^{\mathit{k}}$, for i = 1,2,,ngx$\mathit{i}=1,2,\dots ,{\mathbf{ngx}}$, j = 1,2,,ngy$\mathit{j}=1,2,\dots ,{\mathbf{ngy}}$ and k = 1,2,,7$\mathit{k}=1,2,\dots ,7$.
4:     rhs(lda) – double array
rhs(i + (j1) × ngx)${\mathbf{rhs}}\left(\mathit{i}+\left(\mathit{j}-1\right)×{\mathbf{ngx}}\right)$ must be set to fij${f}_{\mathit{i}\mathit{j}}$, for i = 1,2,,ngx$\mathit{i}=1,2,\dots ,{\mathbf{ngx}}$ and j = 1,2,,ngy$\mathit{j}=1,2,\dots ,{\mathbf{ngy}}$.
5:     ub(${\mathbf{ngx}}×{\mathbf{ngy}}$) – double array
ub(i + (j1) × ngx)${\mathbf{ub}}\left(i+\left(j-1\right)×{\mathbf{ngx}}\right)$ must be set to the initial estimate for the solution uij${u}_{ij}$.
6:     maxit – int64int32nag_int scalar
The maximum permitted number of multigrid iterations. If maxit = 0${\mathbf{maxit}}=0$, no multigrid iterations are performed, but the coarse-grid approximations and incomplete Crout decompositions are computed, and may be output if iout is set accordingly.
Constraint: maxit0${\mathbf{maxit}}\ge 0$.
7:     acc – double scalar
The required tolerance for convergence of the residual 2$2$-norm:
 ‖r‖2 = sqrt( ∑ k = 1(rk)2) $‖r‖2=∑k=1 ngx×ngy (rk) 2$
where r = fAu$r=f-Au$ and u$u$ is the computed solution. Note that the norm is not scaled by the number of equations. The function will stop after fewer than maxit iterations if the residual 2$2$-norm is less than the specified tolerance. (If maxit > 0${\mathbf{maxit}}>0$, at least one iteration is always performed.)
If on entry acc = 0.0${\mathbf{acc}}=0.0$, then the machine precision is used as a default value for the tolerance; if acc > 0.0${\mathbf{acc}}>0.0$, but acc is less than the machine precision, then the function will stop when the residual 2$2$-norm is less than the machine precision and ifail will be set to 4$4$.
Constraint: acc0.0${\mathbf{acc}}\ge 0.0$.
8:     iout – int64int32nag_int scalar
Controls the output of printed information to the advisory message unit as returned by nag_file_set_unit_advisory (x04ab):
iout = 0${\mathbf{iout}}=0$
No output.
iout = 1${\mathbf{iout}}=1$
The solution uij${u}_{\mathit{i}\mathit{j}}$, for i = 1,2,,ngx$\mathit{i}=1,2,\dots ,{\mathbf{ngx}}$ and j = 1,2,,ngy$\mathit{j}=1,2,\dots ,{\mathbf{ngy}}$.
iout = 2${\mathbf{iout}}=2$
The residual 2$2$-norm after each iteration, with the reduction factor over the previous iteration.
iout = 3${\mathbf{iout}}=3$
As for iout = 1${\mathbf{iout}}=1$ and iout = 2${\mathbf{iout}}=2$.
iout = 4${\mathbf{iout}}=4$
As for iout = 3${\mathbf{iout}}=3$, plus the final residual (as returned in ub).
iout = 5${\mathbf{iout}}=5$
As for iout = 4${\mathbf{iout}}=4$, plus the initial elements of a and rhs.
iout = 6${\mathbf{iout}}=6$
As for iout = 5${\mathbf{iout}}=5$, plus the Galerkin coarse grid approximations.
iout = 7${\mathbf{iout}}=7$
As for iout = 6${\mathbf{iout}}=6$, plus the incomplete Crout decompositions.
iout = 8${\mathbf{iout}}=8$
As for iout = 7${\mathbf{iout}}=7$, plus the residual after each iteration.
The elements a(p,k)${\mathbf{a}}\left(p,k\right)$, the Galerkin coarse grid approximations and the incomplete Crout decompositions are output in the format:
• Y-index = j$\text{}=j$
• X-index = ia(p,1)a(p,2)a(p,3)a(p,4)a(p,5)a(p,6)a(p,7)$\text{}=i{\mathbf{a}}\left(p,1\right){\mathbf{a}}\left(p,2\right){\mathbf{a}}\left(p,3\right){\mathbf{a}}\left(p,4\right){\mathbf{a}}\left(p,5\right){\mathbf{a}}\left(p,6\right){\mathbf{a}}\left(p,7\right)$
• where p = i + (j1) × ngx$p=\mathit{i}+\left(\mathit{j}-1\right)×{\mathbf{ngx}}$, for i = 1,2,,ngx$\mathit{i}=1,2,\dots ,{\mathbf{ngx}}$ and j = 1,2,,ngy$\mathit{j}=1,2,\dots ,{\mathbf{ngy}}$.
The vectors u(p)${\mathbf{u}}\left(p\right)$, ub(p)${\mathbf{ub}}\left(p\right)$, rhs(p)${\mathbf{rhs}}\left(p\right)$ are output in matrix form with ngy rows and ngx columns. Where ngx > 10${\mathbf{ngx}}>10$, the ngx values for a given j$j$ value are produced in rows of 10$10$. Values of iout > 4${\mathbf{iout}}>4$ may yield considerable amounts of output.
Constraint: 0iout8$0\le {\mathbf{iout}}\le 8$.

None.

lda

### Output Parameters

1:     a(lda,7$7$) – double array
Is overwritten.
2:     rhs(lda) – double array
The first ${\mathbf{ngx}}×{\mathbf{ngy}}$ elements are unchanged and the rest of the array is used as workspace.
3:     ub(${\mathbf{ngx}}×{\mathbf{ngy}}$) – double array
The corresponding component of the residual r = fau$r=f-{\mathbf{a}}u$.
4:     us(lda) – double array
The residual 2$2$-norm, stored in element us(1)${\mathbf{us}}\left(1\right)$.
5:     u(lda) – double array
The computed solution uij${u}_{\mathit{i}\mathit{j}}$ is returned in u(i + (j1) × ngx)${\mathbf{u}}\left(\mathit{i}+\left(\mathit{j}-1\right)×{\mathbf{ngx}}\right)$, for i = 1,2,,ngx$\mathit{i}=1,2,\dots ,{\mathbf{ngx}}$ and j = 1,2,,ngy$\mathit{j}=1,2,\dots ,{\mathbf{ngy}}$.
6:     numit – int64int32nag_int scalar
The number of iterations performed.
7:     ifail – int64int32nag_int scalar
${\mathrm{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:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

ifail = 1${\mathbf{ifail}}=1$
 On entry, ngx < 3${\mathbf{ngx}}<3$, or ngy < 3${\mathbf{ngy}}<3$, or lda is too small, or acc < 0.0${\mathbf{acc}}<0.0$, or maxit < 0${\mathbf{maxit}}<0$, or iout < 0${\mathbf{iout}}<0$, or iout > 8${\mathbf{iout}}>8$.
ifail = 2${\mathbf{ifail}}=2$
maxit iterations have been performed with the residual 2$2$-norm decreasing at each iteration but the residual 2$2$-norm has not been reduced to less than the specified tolerance (see acc). Examine the progress of the iteration by setting iout2${\mathbf{iout}}\ge 2$.
ifail = 3${\mathbf{ifail}}=3$
As for ${\mathbf{ifail}}={\mathbf{2}}$, except that at one or more iterations the residual 2$2$-norm did not decrease. It is likely that the method fails to converge for the given matrix A$A$.
W ifail = 4${\mathbf{ifail}}=4$
On entry, acc is less than the machine precision. The function terminated because the residual norm is less than the machine precision.

## Accuracy

See acc (Section [Parameters]).

The rate of convergence of this function is strongly dependent upon the number of levels, l$l$, in the multigrid scheme, and thus the choice of ngx and ngy is very important. You are advised to experiment with different values of ngx and ngy to see the effect they have on the rate of convergence; for example, using a value such as ngx = 65${\mathbf{ngx}}=65$ ( = 26 + 1$\text{}={2}^{6}+1$) followed by ngx = 64${\mathbf{ngx}}=64$ (for which l = 1$l=1$).

## Example

```function nag_pde_2d_ellip_mgrid_example
ngx = 9;
ngy = 9;
a = zeros(133,7);
rhs = zeros(133, 1);
ub = zeros(81, 1);
maxit = int64(15);
acc = 0.0001;
iout = int64(0);

hx = 1/double(ngx+1);
hy = 1/double(ngy+1);
alpha = 1.7;

% Set up operator, right-hand side and initial guess for
% step-lengths hx and hy
for i = 1:ngy
for j = 1:ngx
k = (j-1)*ngx + i;
a(k,1) = 1.0 - 0.5*alpha;
a(k,2) = 0.5*alpha;
a(k,3) = 1.0 - 0.5*alpha;
a(k,4) = -4.0 + alpha;
a(k,5) = 1.0 - 0.5*alpha;
a(k,6) = 0.5*alpha;
a(k,7) = 1.0 - 0.5*alpha;
rhs(k) = -4.0*hx*hy;
end
end

% Correction for the boundary conditions
% Horizontal boundaries --
for i = 2:ngx-1
% Boundary condition on y=0 -- u=0
ix = i;
a(ix,1) = 0.0;
a(ix,2) = 0.0;
% Boundary condition on y=1 -- u=0
ix = i + (ngy-1)*ngx;
a(ix,6) = 0.0;
a(ix,7) = 0.0;
end
% Vertical boundaries --
for j = 2:ngy-1
% Boundary condition on x=0 -- u=0
iy = (j-1)*ngx + 1;
a(iy,3) = 0.0;
a(iy,6) = 0.0;
% Boundary condition on x=1 -- u=1
iy = j*ngx;
rhs(iy) = rhs(iy) - a(iy,5) - a(iy,2);
a(iy,2) = 0.0;
a(iy,5) = 0.0;
end

% Now the four corners --
% Bottom left corner
k = 1;
a(k,1) = 0.0;
a(k,2) = 0.0;
a(k,3) = 0.0;
a(k,6) = 0.0;
% Top left corner
k = 1 + (ngy-1)*ngx;
a(k,3) = 0.0;
a(k,6) = 0.0;
a(k,7) = 0.0;
% Bottom right corner
% Use average value at discontinuity ( = 0.5 )
k = ngx;
rhs(k) = rhs(k) - a(k,2)*0.5 - a(k,5);
a(k,1) = 0.0;
a(k,2) = 0.0;
a(k,5) = 0.0;
% Top right corner
k = ngx*ngy;
rhs(k) = rhs(k) - a(k,2) - a(k,5);
a(k,2) = 0.0;
a(k,5) = 0.0;
a(k,6) = 0.0;
a(k,7) = 0.0;

[aOut, rhsOut, ubOut, us, u, numit, ifail] = ...
nag_pde_2d_ellip_mgrid(int64(ngx), int64(ngy), a, rhs, ub, maxit, acc, iout);
aOut, numit, ifail
```
```

aOut =

0         0         0   -2.3000    0.1500         0    0.1500
0         0   -0.0652   -2.2902    0.1500    0.8598    0.1500
0         0   -0.0655   -2.2902    0.1500    0.8598    0.1500
0         0   -0.0655   -2.2902    0.1500    0.8598    0.1500
0         0   -0.0655   -2.2902    0.1500    0.8598    0.1500
0         0   -0.0655   -2.2902    0.1500    0.8598    0.1500
0         0   -0.0655   -2.2902    0.1500    0.8598    0.1500
0         0   -0.0655   -2.2902    0.1500    0.8598    0.1500
0         0   -0.0655   -2.2902         0    0.8598    0.1500
-0.0652   -0.3754         0   -1.9674    0.2063         0    0.1500
-0.0655   -0.3754   -0.1049   -1.9457    0.2063    0.8657    0.1500
-0.0655   -0.3754   -0.1060   -1.9455    0.2063    0.8659    0.1500
-0.0655   -0.3754   -0.1060   -1.9455    0.2063    0.8659    0.1500
-0.0655   -0.3754   -0.1060   -1.9455    0.2063    0.8659    0.1500
-0.0655   -0.3754   -0.1060   -1.9455    0.2063    0.8659    0.1500
-0.0655   -0.3754   -0.1060   -1.9455    0.2063    0.8659    0.1500
-0.0655   -0.3754   -0.1060   -1.9455    0.2063    0.8659    0.1500
-0.0655         0   -0.1060   -2.2683         0    0.8659    0.1500
-0.0762   -0.4449         0   -1.9034    0.2167         0    0.1500
-0.0771   -0.4451   -0.1139   -1.8784    0.2168    0.8671    0.1500
-0.0771   -0.4451   -0.1154   -1.8780    0.2168    0.8673    0.1500
-0.0771   -0.4451   -0.1154   -1.8780    0.2168    0.8673    0.1500
-0.0771   -0.4451   -0.1154   -1.8780    0.2168    0.8673    0.1500
-0.0771   -0.4451   -0.1154   -1.8780    0.2168    0.8673    0.1500
-0.0771   -0.4451   -0.1154   -1.8780    0.2168    0.8673    0.1500
-0.0771   -0.3817   -0.1154   -1.9329    0.2073    0.8673    0.1500
-0.0661         0   -0.1072   -2.2679         0    0.8661    0.1500
-0.0788   -0.4616         0   -1.8879    0.2192         0    0.1500
-0.0799   -0.4618   -0.1161   -1.8620    0.2193    0.8674    0.1500
-0.0799   -0.4618   -0.1178   -1.8617    0.2193    0.8677    0.1500
-0.0799   -0.4618   -0.1178   -1.8616    0.2193    0.8677    0.1500
-0.0799   -0.4618   -0.1178   -1.8616    0.2193    0.8677    0.1500
-0.0799   -0.4618   -0.1178   -1.8616    0.2193    0.8677    0.1500
-0.0799   -0.4487   -0.1178   -1.8730    0.2173    0.8677    0.1500
-0.0776   -0.3819   -0.1160   -1.9324    0.2073    0.8674    0.1500
-0.0661         0   -0.1073   -2.2678         0    0.8661    0.1500
-0.0795   -0.4658         0   -1.8840    0.2199         0    0.1500
-0.0806   -0.4661   -0.1167   -1.8579    0.2199    0.8675    0.1500
-0.0806   -0.4661   -0.1184   -1.8575    0.2199    0.8678    0.1500
-0.0806   -0.4661   -0.1184   -1.8575    0.2199    0.8678    0.1500
-0.0806   -0.4661   -0.1184   -1.8575    0.2199    0.8678    0.1500
-0.0806   -0.4632   -0.1184   -1.8599    0.2195    0.8678    0.1500
-0.0801   -0.4489   -0.1180   -1.8727    0.2173    0.8677    0.1500
-0.0776   -0.3819   -0.1161   -1.9324    0.2073    0.8674    0.1500
-0.0661         0   -0.1073   -2.2678         0    0.8661    0.1500
-0.0796   -0.4669         0   -1.8830    0.2200         0    0.1500
-0.0807   -0.4672   -0.1169   -1.8568    0.2201    0.8675    0.1500
-0.0808   -0.4672   -0.1185   -1.8564    0.2201    0.8678    0.1500
-0.0808   -0.4672   -0.1185   -1.8564    0.2201    0.8678    0.1500
-0.0808   -0.4666   -0.1185   -1.8569    0.2200    0.8678    0.1500
-0.0806   -0.4633   -0.1185   -1.8598    0.2195    0.8678    0.1500
-0.0801   -0.4489   -0.1180   -1.8727    0.2173    0.8677    0.1500
-0.0776   -0.3819   -0.1161   -1.9324    0.2073    0.8674    0.1500
-0.0661         0   -0.1073   -2.2678         0    0.8661    0.1500
-0.0797   -0.4672         0   -1.8827    0.2201         0    0.1500
-0.0808   -0.4674   -0.1169   -1.8565    0.2201    0.8675    0.1500
-0.0808   -0.4675   -0.1186   -1.8561    0.2201    0.8678    0.1500
-0.0808   -0.4673   -0.1186   -1.8562    0.2201    0.8678    0.1500
-0.0808   -0.4666   -0.1186   -1.8569    0.2200    0.8678    0.1500
-0.0807   -0.4633   -0.1185   -1.8598    0.2195    0.8678    0.1500
-0.0801   -0.4489   -0.1180   -1.8727    0.2173    0.8677    0.1500
-0.0776   -0.3819   -0.1161   -1.9324    0.2073    0.8674    0.1500
-0.0661         0   -0.1073   -2.2678         0    0.8661    0.1500
-0.0797   -0.4673         0   -1.8827    0.2201         0    0.1500
-0.0808   -0.4675   -0.1169   -1.8564    0.2201    0.8675    0.1500
-0.0808   -0.4675   -0.1186   -1.8561    0.2201    0.8678    0.1500
-0.0808   -0.4673   -0.1186   -1.8562    0.2201    0.8678    0.1500
-0.0808   -0.4666   -0.1186   -1.8569    0.2200    0.8678    0.1500
-0.0807   -0.4633   -0.1185   -1.8598    0.2195    0.8678    0.1500
-0.0801   -0.4489   -0.1180   -1.8727    0.2173    0.8677    0.1500
-0.0776   -0.3819   -0.1161   -1.9324    0.2073    0.8674    0.1500
-0.0661         0   -0.1073   -2.2678         0    0.8661    0.1500
-0.0797   -0.4673         0   -1.8826    0.2201         0         0
-0.0808   -0.4675   -0.1169   -1.8564    0.2201         0         0
-0.0808   -0.4675   -0.1186   -1.8561    0.2201         0         0
-0.0808   -0.4673   -0.1186   -1.8562    0.2201         0         0
-0.0808   -0.4666   -0.1186   -1.8569    0.2200         0         0
-0.0807   -0.4633   -0.1185   -1.8598    0.2195         0         0
-0.0801   -0.4489   -0.1180   -1.8727    0.2173         0         0
-0.0776   -0.3819   -0.1161   -1.9324    0.2073         0         0
-0.0661         0   -0.1073   -2.2678         0         0         0
0         0         0   -2.7250   -0.1375         0   -0.1375
0         0    0.0505   -2.7181   -0.1375    0.8569    0.1500
0         0    0.0506   -2.7180   -0.1375    0.8424    0.1500
0         0    0.0506   -2.7180   -0.1375    0.8424    0.1500
0         0    0.0506   -2.7180         0    0.8424   -0.1375
0.0505   -0.3153         0   -2.4479    0.1973         0   -0.1375
-0.0552   -0.3099   -0.0806   -2.0147    0.1965    0.8389    0.1500
-0.0552   -0.3099   -0.0975   -2.0115    0.1965    0.8646    0.1500
-0.0552   -0.3099   -0.0977   -2.0114    0.1074    0.8647    0.1500
0.0506         0   -0.0534   -2.7123         0    0.8580   -0.1375
0.0562   -0.4164         0   -2.3680    0.2125         0   -0.1375
-0.0745   -0.4298   -0.0897   -1.8981    0.2145    0.8377    0.1500
-0.0746   -0.4299   -0.1130   -1.8929    0.2145    0.8669    0.1500
-0.0746   -0.3163   -0.1133   -1.9931    0.1065    0.8670    0.1500
0.0507         0   -0.0534   -2.7123         0    0.8580   -0.1375
0.0581   -0.4413         0   -2.3473    0.2162         0   -0.1375
-0.0790   -0.4580   -0.0921   -1.8712    0.2187    0.8373    0.1500
-0.0792   -0.4350   -0.1169   -1.8854    0.2153    0.8675    0.1500
-0.0753   -0.3163   -0.1142   -1.9927    0.1065    0.8671    0.1500
0.0507         0   -0.0534   -2.7123         0    0.8580   -0.1375
0.0586   -0.4475         0   -2.3422   -0.0704         0         0
-0.0802   -0.4601    0.0300   -2.3117   -0.0685         0         0
-0.0796   -0.4351    0.0296   -2.3337   -0.0722         0         0
-0.0753   -0.3163    0.0309   -2.4400   -0.1810         0         0
0.0507         0    0.0742   -2.7046         0         0         0
0         0         0   -3.9375   -0.5312         0   -0.5312
0         0    0.1349   -3.8658   -0.5312    0.9217    0.1500
0         0    0.1374   -3.8645         0    0.8294   -0.5312
0.1349   -0.2384         0   -3.6461    0.1858         0   -0.5312
-0.0388   -0.2146   -0.0509   -2.1067    0.0360    0.8229    0.1500
0.1375         0   -0.0171   -3.8639         0    0.8526   -0.5312
0.1457   -0.3906         0   -3.5386   -0.4727         0         0
-0.0712   -0.2207    0.1336   -3.6756   -0.6485         0         0
0.1375         0    0.1764   -3.7500         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0

numit =

4

ifail =

0

```
```function d03ed_example
ngx = 9;
ngy = 9;
a = zeros(133,7);
rhs = zeros(133, 1);
ub = zeros(81, 1);
maxit = int64(15);
acc = 0.0001;
iout = int64(0);

hx = 1/double(ngx+1);
hy = 1/double(ngy+1);
alpha = 1.7;

% Set up operator, right-hand side and initial guess for
% step-lengths hx and hy
for i = 1:ngy
for j = 1:ngx
k = (j-1)*ngx + i;
a(k,1) = 1.0 - 0.5*alpha;
a(k,2) = 0.5*alpha;
a(k,3) = 1.0 - 0.5*alpha;
a(k,4) = -4.0 + alpha;
a(k,5) = 1.0 - 0.5*alpha;
a(k,6) = 0.5*alpha;
a(k,7) = 1.0 - 0.5*alpha;
rhs(k) = -4.0*hx*hy;
end
end

% Correction for the boundary conditions
% Horizontal boundaries --
for i = 2:ngx-1
% Boundary condition on y=0 -- u=0
ix = i;
a(ix,1) = 0.0;
a(ix,2) = 0.0;
% Boundary condition on y=1 -- u=0
ix = i + (ngy-1)*ngx;
a(ix,6) = 0.0;
a(ix,7) = 0.0;
end
% Vertical boundaries --
for j = 2:ngy-1
% Boundary condition on x=0 -- u=0
iy = (j-1)*ngx + 1;
a(iy,3) = 0.0;
a(iy,6) = 0.0;
% Boundary condition on x=1 -- u=1
iy = j*ngx;
rhs(iy) = rhs(iy) - a(iy,5) - a(iy,2);
a(iy,2) = 0.0;
a(iy,5) = 0.0;
end

% Now the four corners --
% Bottom left corner
k = 1;
a(k,1) = 0.0;
a(k,2) = 0.0;
a(k,3) = 0.0;
a(k,6) = 0.0;
% Top left corner
k = 1 + (ngy-1)*ngx;
a(k,3) = 0.0;
a(k,6) = 0.0;
a(k,7) = 0.0;
% Bottom right corner
% Use average value at discontinuity ( = 0.5 )
k = ngx;
rhs(k) = rhs(k) - a(k,2)*0.5 - a(k,5);
a(k,1) = 0.0;
a(k,2) = 0.0;
a(k,5) = 0.0;
% Top right corner
k = ngx*ngy;
rhs(k) = rhs(k) - a(k,2) - a(k,5);
a(k,2) = 0.0;
a(k,5) = 0.0;
a(k,6) = 0.0;
a(k,7) = 0.0;

[aOut, rhsOut, ubOut, us, u, numit, ifail] = ...
d03ed(int64(ngx), int64(ngy), a, rhs, ub, maxit, acc, iout);
aOut, numit, ifail
```
```

aOut =

0         0         0   -2.3000    0.1500         0    0.1500
0         0   -0.0652   -2.2902    0.1500    0.8598    0.1500
0         0   -0.0655   -2.2902    0.1500    0.8598    0.1500
0         0   -0.0655   -2.2902    0.1500    0.8598    0.1500
0         0   -0.0655   -2.2902    0.1500    0.8598    0.1500
0         0   -0.0655   -2.2902    0.1500    0.8598    0.1500
0         0   -0.0655   -2.2902    0.1500    0.8598    0.1500
0         0   -0.0655   -2.2902    0.1500    0.8598    0.1500
0         0   -0.0655   -2.2902         0    0.8598    0.1500
-0.0652   -0.3754         0   -1.9674    0.2063         0    0.1500
-0.0655   -0.3754   -0.1049   -1.9457    0.2063    0.8657    0.1500
-0.0655   -0.3754   -0.1060   -1.9455    0.2063    0.8659    0.1500
-0.0655   -0.3754   -0.1060   -1.9455    0.2063    0.8659    0.1500
-0.0655   -0.3754   -0.1060   -1.9455    0.2063    0.8659    0.1500
-0.0655   -0.3754   -0.1060   -1.9455    0.2063    0.8659    0.1500
-0.0655   -0.3754   -0.1060   -1.9455    0.2063    0.8659    0.1500
-0.0655   -0.3754   -0.1060   -1.9455    0.2063    0.8659    0.1500
-0.0655         0   -0.1060   -2.2683         0    0.8659    0.1500
-0.0762   -0.4449         0   -1.9034    0.2167         0    0.1500
-0.0771   -0.4451   -0.1139   -1.8784    0.2168    0.8671    0.1500
-0.0771   -0.4451   -0.1154   -1.8780    0.2168    0.8673    0.1500
-0.0771   -0.4451   -0.1154   -1.8780    0.2168    0.8673    0.1500
-0.0771   -0.4451   -0.1154   -1.8780    0.2168    0.8673    0.1500
-0.0771   -0.4451   -0.1154   -1.8780    0.2168    0.8673    0.1500
-0.0771   -0.4451   -0.1154   -1.8780    0.2168    0.8673    0.1500
-0.0771   -0.3817   -0.1154   -1.9329    0.2073    0.8673    0.1500
-0.0661         0   -0.1072   -2.2679         0    0.8661    0.1500
-0.0788   -0.4616         0   -1.8879    0.2192         0    0.1500
-0.0799   -0.4618   -0.1161   -1.8620    0.2193    0.8674    0.1500
-0.0799   -0.4618   -0.1178   -1.8617    0.2193    0.8677    0.1500
-0.0799   -0.4618   -0.1178   -1.8616    0.2193    0.8677    0.1500
-0.0799   -0.4618   -0.1178   -1.8616    0.2193    0.8677    0.1500
-0.0799   -0.4618   -0.1178   -1.8616    0.2193    0.8677    0.1500
-0.0799   -0.4487   -0.1178   -1.8730    0.2173    0.8677    0.1500
-0.0776   -0.3819   -0.1160   -1.9324    0.2073    0.8674    0.1500
-0.0661         0   -0.1073   -2.2678         0    0.8661    0.1500
-0.0795   -0.4658         0   -1.8840    0.2199         0    0.1500
-0.0806   -0.4661   -0.1167   -1.8579    0.2199    0.8675    0.1500
-0.0806   -0.4661   -0.1184   -1.8575    0.2199    0.8678    0.1500
-0.0806   -0.4661   -0.1184   -1.8575    0.2199    0.8678    0.1500
-0.0806   -0.4661   -0.1184   -1.8575    0.2199    0.8678    0.1500
-0.0806   -0.4632   -0.1184   -1.8599    0.2195    0.8678    0.1500
-0.0801   -0.4489   -0.1180   -1.8727    0.2173    0.8677    0.1500
-0.0776   -0.3819   -0.1161   -1.9324    0.2073    0.8674    0.1500
-0.0661         0   -0.1073   -2.2678         0    0.8661    0.1500
-0.0796   -0.4669         0   -1.8830    0.2200         0    0.1500
-0.0807   -0.4672   -0.1169   -1.8568    0.2201    0.8675    0.1500
-0.0808   -0.4672   -0.1185   -1.8564    0.2201    0.8678    0.1500
-0.0808   -0.4672   -0.1185   -1.8564    0.2201    0.8678    0.1500
-0.0808   -0.4666   -0.1185   -1.8569    0.2200    0.8678    0.1500
-0.0806   -0.4633   -0.1185   -1.8598    0.2195    0.8678    0.1500
-0.0801   -0.4489   -0.1180   -1.8727    0.2173    0.8677    0.1500
-0.0776   -0.3819   -0.1161   -1.9324    0.2073    0.8674    0.1500
-0.0661         0   -0.1073   -2.2678         0    0.8661    0.1500
-0.0797   -0.4672         0   -1.8827    0.2201         0    0.1500
-0.0808   -0.4674   -0.1169   -1.8565    0.2201    0.8675    0.1500
-0.0808   -0.4675   -0.1186   -1.8561    0.2201    0.8678    0.1500
-0.0808   -0.4673   -0.1186   -1.8562    0.2201    0.8678    0.1500
-0.0808   -0.4666   -0.1186   -1.8569    0.2200    0.8678    0.1500
-0.0807   -0.4633   -0.1185   -1.8598    0.2195    0.8678    0.1500
-0.0801   -0.4489   -0.1180   -1.8727    0.2173    0.8677    0.1500
-0.0776   -0.3819   -0.1161   -1.9324    0.2073    0.8674    0.1500
-0.0661         0   -0.1073   -2.2678         0    0.8661    0.1500
-0.0797   -0.4673         0   -1.8827    0.2201         0    0.1500
-0.0808   -0.4675   -0.1169   -1.8564    0.2201    0.8675    0.1500
-0.0808   -0.4675   -0.1186   -1.8561    0.2201    0.8678    0.1500
-0.0808   -0.4673   -0.1186   -1.8562    0.2201    0.8678    0.1500
-0.0808   -0.4666   -0.1186   -1.8569    0.2200    0.8678    0.1500
-0.0807   -0.4633   -0.1185   -1.8598    0.2195    0.8678    0.1500
-0.0801   -0.4489   -0.1180   -1.8727    0.2173    0.8677    0.1500
-0.0776   -0.3819   -0.1161   -1.9324    0.2073    0.8674    0.1500
-0.0661         0   -0.1073   -2.2678         0    0.8661    0.1500
-0.0797   -0.4673         0   -1.8826    0.2201         0         0
-0.0808   -0.4675   -0.1169   -1.8564    0.2201         0         0
-0.0808   -0.4675   -0.1186   -1.8561    0.2201         0         0
-0.0808   -0.4673   -0.1186   -1.8562    0.2201         0         0
-0.0808   -0.4666   -0.1186   -1.8569    0.2200         0         0
-0.0807   -0.4633   -0.1185   -1.8598    0.2195         0         0
-0.0801   -0.4489   -0.1180   -1.8727    0.2173         0         0
-0.0776   -0.3819   -0.1161   -1.9324    0.2073         0         0
-0.0661         0   -0.1073   -2.2678         0         0         0
0         0         0   -2.7250   -0.1375         0   -0.1375
0         0    0.0505   -2.7181   -0.1375    0.8569    0.1500
0         0    0.0506   -2.7180   -0.1375    0.8424    0.1500
0         0    0.0506   -2.7180   -0.1375    0.8424    0.1500
0         0    0.0506   -2.7180         0    0.8424   -0.1375
0.0505   -0.3153         0   -2.4479    0.1973         0   -0.1375
-0.0552   -0.3099   -0.0806   -2.0147    0.1965    0.8389    0.1500
-0.0552   -0.3099   -0.0975   -2.0115    0.1965    0.8646    0.1500
-0.0552   -0.3099   -0.0977   -2.0114    0.1074    0.8647    0.1500
0.0506         0   -0.0534   -2.7123         0    0.8580   -0.1375
0.0562   -0.4164         0   -2.3680    0.2125         0   -0.1375
-0.0745   -0.4298   -0.0897   -1.8981    0.2145    0.8377    0.1500
-0.0746   -0.4299   -0.1130   -1.8929    0.2145    0.8669    0.1500
-0.0746   -0.3163   -0.1133   -1.9931    0.1065    0.8670    0.1500
0.0507         0   -0.0534   -2.7123         0    0.8580   -0.1375
0.0581   -0.4413         0   -2.3473    0.2162         0   -0.1375
-0.0790   -0.4580   -0.0921   -1.8712    0.2187    0.8373    0.1500
-0.0792   -0.4350   -0.1169   -1.8854    0.2153    0.8675    0.1500
-0.0753   -0.3163   -0.1142   -1.9927    0.1065    0.8671    0.1500
0.0507         0   -0.0534   -2.7123         0    0.8580   -0.1375
0.0586   -0.4475         0   -2.3422   -0.0704         0         0
-0.0802   -0.4601    0.0300   -2.3117   -0.0685         0         0
-0.0796   -0.4351    0.0296   -2.3337   -0.0722         0         0
-0.0753   -0.3163    0.0309   -2.4400   -0.1810         0         0
0.0507         0    0.0742   -2.7046         0         0         0
0         0         0   -3.9375   -0.5312         0   -0.5312
0         0    0.1349   -3.8658   -0.5312    0.9217    0.1500
0         0    0.1374   -3.8645         0    0.8294   -0.5312
0.1349   -0.2384         0   -3.6461    0.1858         0   -0.5312
-0.0388   -0.2146   -0.0509   -2.1067    0.0360    0.8229    0.1500
0.1375         0   -0.0171   -3.8639         0    0.8526   -0.5312
0.1457   -0.3906         0   -3.5386   -0.4727         0         0
-0.0712   -0.2207    0.1336   -3.6756   -0.6485         0         0
0.1375         0    0.1764   -3.7500         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0
0         0         0         0         0         0         0

numit =

4

ifail =

0

```