NAG Toolbox: nag_pde_2d_laplace (d03ea)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{stage1}$ – logical scalar

Indicates whether the function call is for stage $1$ of the computation as defined in Description.

$\mathbf{stage1}=\mathit{true}$

The call is for stage $1$.

$\mathbf{stage1}=\mathit{false}$

The call is for stage $2$.

2:     $\mathrm{ext}$ – logical scalar

The form of the domain. If $\mathbf{ext}=\mathit{true}$, the domain is unbounded. Otherwise the domain is an interior one.

3:     $\mathrm{dorm}$ – logical scalar

The form of the boundary conditions. If $\mathbf{dorm}=\mathit{true}$, then the problem is a Dirichlet or mixed boundary value problem. Otherwise it is a Neumann problem.

4:     $\mathrm{p}$ – double scalar

5:     $\mathrm{q}$ – double scalar

To stage $2$, p and q must specify the $x$ and $y$ coordinates respectively of a point at which the solution is required.

When stage1 is true, p and q are ignored.

6:     $\mathrm{x}\left(\mathbf{n1p1}\right)$ – double array

7:     $\mathrm{y}\left(\mathbf{n1p1}\right)$ – double array

The $x$ and $y$ coordinates respectively of points on the one or more closed contours which define the domain of the problem.

Note: each contour is described in such a manner that the subscripts of the coordinates increase when the domain is kept on the left. The final point on each contour coincides with the first and, if a further contour is to be described, the coordinates of this point are repeated in the arrays. In this way each interval is defined by three points, the second of which (the nodal point) always has an even subscript. In the case of the interior Neumann problem, the outermost boundary contour must be given first, if there is more than one.

8:     $\mathrm{phi}\left(\mathbf{n}\right)$ – double array

For stage $1$, phi must contain the nodal values of $\varphi$ or its normal derivative (into the domain) as prescribed in each interval. For stage $2$ it must retain its output values from stage $1$.

9:     $\mathrm{phid}\left(\mathbf{n}\right)$ – double array

For stage $1$, $\mathbf{phid}\left(\mathit{i}\right)$ must hold the value $0.0$ or $1.0$ accordingly as $\mathbf{phi}\left(\mathit{i}\right)$ contains a value of $\varphi$ or its normal derivative, for $\mathit{i}=1,2,\dots ,\mathbf{n}$. For stage $2$ it must retain its output values from stage $1$.

10:   $\mathrm{alpha}$ – double scalar

For stage $1$, the use of alpha depends on the nature of the problem:

• if $\mathbf{dorm}=\mathit{true}$, alpha need not be set;
• if $\mathbf{dorm}=\mathit{false}$ and $\mathbf{ext}=\mathit{true}$, alpha must contain the prescribed constant $c$ (see Description);
• if $\mathbf{dorm}=\mathit{false}$ and $\mathbf{ext}=\mathit{false}$, alpha must contain an appropriate value (often zero) for the integral of $\varphi$ around the outermost boundary.

For stage $2$, on every call alpha must contain the value returned at stage $1$.

Optional Input Parameters

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

Default: the dimension of the arrays phi, phid. (An error is raised if these dimensions are not equal.)

The number of intervals into which the boundary is divided (see Accuracy and Further Comments).

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

Default: the dimension of the arrays x, y. (An error is raised if these dimensions are not equal.)

The value $2\left(\mathbf{n}+M\right)-1$, where $M$ denotes the number of closed contours which make up the boundary.

Output Parameters

1:     $\mathrm{phi}\left(\mathbf{n}\right)$ – double array

From stage $1$, it contains the constants which approximate $\varphi$ in each interval. It remains unchanged on exit from stage $2$.

2:     $\mathrm{phid}\left(\mathbf{n}\right)$ – double array

From stage $1$, phid contains the constants which approximate the normal derivative of $\varphi$ in each interval. It remains unchanged on exit from stage $2$.

3:     $\mathrm{alpha}$ – double scalar

From stage $1$:

• if $\mathbf{ext}=\mathit{false}$, alpha contains $0.0$;
• if $\mathbf{ext}=\mathit{true}$ and $\mathbf{dorm}=\mathit{false}$ alpha is unchanged;
• if $\mathbf{ext}=\mathit{true}$ and $\mathbf{dorm}=\mathit{true}$ alpha contains a computed estimate for $c$.

From stage $2$:

• alpha contains the computed value of $\varphi$ at the point (p,q).

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

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

Error Indicators and Warnings

   $\mathbf{ifail}=1$

Invalid tolerance used in an internal call to an auxiliary function:

$\mathit{icint}\left(1\right)=0$

Indicates too large a tolerance.

$\mathit{icint}\left(1\right)>0$

Indicates too small a tolerance.

Note:  this error is only possible in stage 1, and the circumstances under which it may occur cannot be foreseen. In the event of a failure, it is suggested that you change the scale of the domain of the problem and apply the function again.

   $\mathbf{ifail}=2$

Incorrect rank obtained by an auxiliary function; $\mathit{icint}\left(1\right)$ contains the computed rank.

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

function d03ea_example


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


% This example illustrates using an interior Neumann problem to
% solve the Laplace equation in the domain bounded by two triangles.

% Stage 1: Set up phi and phid
stage1 = true;
ext = false;
dorm = false;
p = 0;
q = 0;
% Outer and inner triangles with mid-points
x = [ 3.0;  1.5;  0.0; -1.5; -3.0;  0.0;  3.0;
      3.0;
      2.0;  0.0; -2.0; -1.0;  0.0;  1.0;  2.0];
y = [ 0.0;  2.0;  4.0;  2.0;  0.0;  0.0;  0.0;
      0.0;
      1.0;  1.0;  1.0;  2.0;  3.0;  2.0;  1.0];
% Some interesting Neumann conditions
phi(1:6) = [3, 4, 4, 5, 4, 3];
phid(1:6) = [0, 0, 0, 0, 0, 0];
% Integral of solution around outer triangle
alpha = 12;
[phi, phid, alpha, ifail] = ...
    d03ea( ...
           stage1, ext, dorm, p, q, x, y, phi, phid, alpha);

% Stage 2 - solution at a set of points
stage1 = false;
phisol = zeros(401,601);
solmin = 1.e10;
solmax = -1.0e10;
i = 0;
for xx = -3:0.01:3
  i = i + 1;
  j = 0;
  xl = 3-abs(xx);
  y1 = (4/3)*xl;
  for yy = 0:0.01:4
    j = j + 1;
    if ( (yy<=y1) && ((yy<=1) || (yy>=xl)));
      % In Outer triangle and outside inner triangle
      [phi, phid, sol, ifail] = ...
      d03ea( ...
             stage1, ext, dorm, xx, yy, x, y, phi, phid, alpha);
      phisol(j,i) = sol;
      solmax = max(solmax,sol);
      solmin = min(solmin,sol);
    end
  end
end

fprintf('The maximum value of solution = %8.3f\n',solmax);
fprintf('The minimum value of solution = %8.3f\n',solmin);

% Plot Solution
fig1 = figure;
contourf(phisol,[solmin:-1,1,2:0.2:solmax]);
title('Laplace Equation Solution between two triangles');
set(gca,'XTick',[],'XTickLabel',[],'YTick',[],'YTickLabel',[]);

d03ea example results

The maximum value of solution =    5.430
The minimum value of solution =   -2.979