NAG Toolbox: nag_pde_2d_ellip_mgrid (d03ed)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{ngx}$ – int64int32nag_int scalar

The number of interior grid points in the $x$-direction, $n_x$. $\mathbf{ngx}-1$ should preferably be divisible by as high a power of $2$ as possible.

Constraint: $\mathbf{ngx}\ge 3$.

2:     $\mathrm{ngy}$ – int64int32nag_int scalar

The number of interior grid points in the $y$-direction, $n_y$. $\mathbf{ngy}-1$ should preferably be divisible by as high a power of $2$ as possible.

Constraint: $\mathbf{ngy}\ge 3$.

3:     $\mathrm{a}\left(\mathit{lda},7\right)$ – double array

$\mathbf{a}\left(\mathit{i}+\left(\mathit{j}-1\right)\times \mathbf{ngx},\mathit{k}\right)$ must be set to ${\mathbf{a}}_{\mathit{i}\mathit{j}}^{\mathit{k}}$, for $\mathit{i}=1,2,\dots ,\mathbf{ngx}$, $\mathit{j}=1,2,\dots ,\mathbf{ngy}$ and $\mathit{k}=1,2,\dots ,7$.

4:     $\mathrm{rhs}\left(\mathit{lda}\right)$ – double array

$\mathbf{rhs}\left(\mathit{i}+\left(\mathit{j}-1\right)\times \mathbf{ngx}\right)$ must be set to $f_{\mathit{i}\mathit{j}}$, for $\mathit{i}=1,2,\dots ,\mathbf{ngx}$ and $\mathit{j}=1,2,\dots ,\mathbf{ngy}$.

5:     $\mathrm{ub}\left(\mathbf{ngx}\times \mathbf{ngy}\right)$ – double array

$\mathbf{ub}\left(i+\left(j-1\right)\times \mathbf{ngx}\right)$ must be set to the initial estimate for the solution $u_{ij}$.

6:     $\mathrm{maxit}$ – int64int32nag_int scalar

The maximum permitted number of multigrid iterations. If $\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: $\mathbf{maxit}\ge 0$.

7:     $\mathrm{acc}$ – double scalar

The required tolerance for convergence of the residual $2$-norm:

$${‖r‖}_2=\sqrt{\sum _{k=1}^{\mathbf{ngx}\times \mathbf{ngy}}{\left(r_k\right)}^2}$$

where $r=f-Au$ and $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$-norm is less than the specified tolerance. (If $\mathbf{maxit}>0$, at least one iteration is always performed.)

If on entry $\mathbf{acc}=0.0$, then the machine precision is used as a default value for the tolerance; if $\mathbf{acc}>0.0$, but acc is less than the machine precision, then the function will stop when the residual $2$-norm is less than the machine precision and ifail will be set to $4$.

Constraint: $\mathbf{acc}\ge 0.0$.

8:     $\mathrm{iout}$ – int64int32nag_int scalar

Controls the output of printed information to the advisory message unit as returned by nag_file_set_unit_advisory (x04ab):

$\mathbf{iout}=0$

No output.

$\mathbf{iout}=1$

The solution $u_{\mathit{i}\mathit{j}}$, for $\mathit{i}=1,2,\dots ,\mathbf{ngx}$ and $\mathit{j}=1,2,\dots ,\mathbf{ngy}$.

$\mathbf{iout}=2$

The residual $2$-norm after each iteration, with the reduction factor over the previous iteration.

$\mathbf{iout}=3$

As for $\mathbf{iout}=1$ and $\mathbf{iout}=2$.

$\mathbf{iout}=4$

As for $\mathbf{iout}=3$, plus the final residual (as returned in ub).

$\mathbf{iout}=5$

As for $\mathbf{iout}=4$, plus the initial elements of a and rhs.

$\mathbf{iout}=6$

As for $\mathbf{iout}=5$, plus the Galerkin coarse grid approximations.

$\mathbf{iout}=7$

As for $\mathbf{iout}=6$, plus the incomplete Crout decompositions.

$\mathbf{iout}=8$

As for $\mathbf{iout}=7$, plus the residual after each iteration.

The elements $\mathbf{a}\left(p,k\right)$, the Galerkin coarse grid approximations and the incomplete Crout decompositions are output in the format:

• Y-index $=j$
• X-index $=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=\mathit{i}+\left(\mathit{j}-1\right)\times \mathbf{ngx}$, for $\mathit{i}=1,2,\dots ,\mathbf{ngx}$ and $\mathit{j}=1,2,\dots ,\mathbf{ngy}$.

The vectors $\mathbf{u}\left(p\right)$, $\mathbf{ub}\left(p\right)$, $\mathbf{rhs}\left(p\right)$ are output in matrix form with ngy rows and ngx columns. Where $\mathbf{ngx}>10$, the ngx values for a given $j$ value are produced in rows of $10$. Values of $\mathbf{iout}>4$ may yield considerable amounts of output.

Constraint: $0\le \mathbf{iout}\le 8$.

Optional Input Parameters

None.

Output Parameters

1:     $\mathrm{a}\left(\mathit{lda},7\right)$ – double array

Is overwritten.

2:     $\mathrm{rhs}\left(\mathit{lda}\right)$ – double array

The first $\mathbf{ngx}\times \mathbf{ngy}$ elements are unchanged and the rest of the array is used as workspace.

3:     $\mathrm{ub}\left(\mathbf{ngx}\times \mathbf{ngy}\right)$ – double array

The corresponding component of the residual $r=f-\mathbf{a}u$.

4:     $\mathrm{us}\left(\mathit{lda}\right)$ – double array

The residual $2$-norm, stored in element $\mathbf{us}\left(1\right)$.

5:     $\mathrm{u}\left(\mathit{lda}\right)$ – double array

The computed solution $u_{\mathit{i}\mathit{j}}$ is returned in $\mathbf{u}\left(\mathit{i}+\left(\mathit{j}-1\right)\times \mathbf{ngx}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{ngx}$ and $\mathit{j}=1,2,\dots ,\mathbf{ngy}$.

6:     $\mathrm{numit}$ – int64int32nag_int scalar

The number of iterations performed.

7:     $\mathrm{ifail}$ – int64int32nag_int scalar

$\mathbf{ifail}=\mathbf{0}$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

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

   $\mathbf{ifail}=1$

On entry,	$\mathbf{ngx}<3$,
or	$\mathbf{ngy}<3$,
or	lda is too small,
or	$\mathbf{acc}<0.0$,
or	$\mathbf{maxit}<0$,
or	$\mathbf{iout}<0$,
or	$\mathbf{iout}>8$.

   $\mathbf{ifail}=2$

maxit iterations have been performed with the residual $2$-norm decreasing at each iteration but the residual $2$-norm has not been reduced to less than the specified tolerance (see acc). Examine the progress of the iteration by setting $\mathbf{iout}\ge 2$.

   $\mathbf{ifail}=3$

As for $\mathbf{ifail}=\mathbf{2}$, except that at one or more iterations the residual $2$-norm did not decrease. It is likely that the method fails to converge for the given matrix $A$.

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

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

Further Comments

Example

Open in the MATLAB editor: d03ed_example

function d03ed_example


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

ngx = int64(65); 
ngy = ngx; 
nn  = ngx*ngy;
lda = 4*(ngx+1)*(ngy+1)/3;
a   = zeros(lda,7);
rhs = zeros(lda, 1);
ub  = zeros(nn, 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
a(1:nn,[1 3 5 7]) = 1 - 0.5*alpha;
a(1:nn,[2 6])     = 0.5*alpha;
a(1:nn,4)         = -4 + alpha;
rhs(1:nn)         = -4*hx*hy;

% Boundary conditions u(x,y=0,1) = 0
a(2:ngx-1,1:2)       = 0.0;
a(nn-ngx+2:nn-1,6:7) = 0.0;
for j = 2:ngy-1
  % Boundary condition u(x=0,y) = 0
  iy = (j-1)*ngx + 1;
  a(iy,[3 6]) = 0.0;
  % Boundary condition u(x=1,y) = 1
  iy = j*ngx;
  rhs(iy) = rhs(iy) - a(iy,5) - a(iy,2);
  a(iy,[2 5]) = 0.0;
end

% Bottom left corner
a(1,[1:3 6]) = 0.0;
% Top left corner
k = nn - ngx + 1;
a(k,[3 6: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:2 5]) = 0.0;
% Top right corner
k = ngx*ngy;
rhs(k) = rhs(k) - a(k,2) - a(k,5);
a(k,[2 5:7]) = 0.0;
u = zeros(ngx*ngy,1);

[a, rhs, ub, resid, u, numit, ifail] = ...
d03ed(...
       ngx, ngy, a, rhs, ub, maxit, acc, iout);

 % Output solution statistics
fprintf('Number of iterations  %4d\n',numit);
fprintf('Maximum residual       %8.3e\n',resid(numit));

% Plot the solution u
fig1 = figure;
u2 = u(1:nn);
umat = reshape(u2,ngx,ngy);
x(1:ngx) = hx:hx:1-hx; y(1:ngy) = hy:hy:1-hy;
[xs,ys] = meshgrid(y,x);
mesh(xs,ys,umat);
xlabel('x'); ylabel('y'); zlabel ('u(x,y)');
title('d03ed example:  U_{xx}-\alpha U_{xy}+U_{yy}=-4, \alpha=1.7');

d03ed example results

Number of iterations     4
Maximum residual       -1.060e-08