NAG FL Interface
d03eaf (dim2_laplace)
1
Purpose
d03eaf solves Laplace's equation in two dimensions for an arbitrary domain bounded internally or externally by one or more closed contours, given the value of either the unknown function or its normal derivative (into the domain) at each point of the boundary.
2
Specification
Fortran Interface
Subroutine d03eaf ( |
stage1, ext, dorm, n, p, q, x, y, n1p1, phi, phid, alpha, c, ldc, np4, icint, np1, ifail) |
Integer, Intent (In) |
:: |
n, n1p1, ldc, np4, np1 |
Integer, Intent (Inout) |
:: |
ifail |
Integer, Intent (Out) |
:: |
icint(n+1) |
Real (Kind=nag_wp), Intent (In) |
:: |
p, q, x(n1p1), y(n1p1) |
Real (Kind=nag_wp), Intent (Inout) |
:: |
phi(n), phid(n), alpha, c(ldc,n+4) |
Logical, Intent (In) |
:: |
stage1, ext, dorm |
|
C Header Interface
#include <nag.h>
void |
d03eaf_ (const logical *stage1, const logical *ext, const logical *dorm, const Integer *n, const double *p, const double *q, const double x[], const double y[], const Integer *n1p1, double phi[], double phid[], double *alpha, double c[], const Integer *ldc, const Integer *np4, Integer icint[], const Integer *np1, Integer *ifail) |
|
C++ Header Interface
#include <nag.h> extern "C" {
void |
d03eaf_ (const logical &stage1, const logical &ext, const logical &dorm, const Integer &n, const double &p, const double &q, const double x[], const double y[], const Integer &n1p1, double phi[], double phid[], double &alpha, double c[], const Integer &ldc, const Integer &np4, Integer icint[], const Integer &np1, Integer &ifail) |
}
|
The routine may be called by the names d03eaf or nagf_pde_dim2_laplace.
3
Description
d03eaf uses an integral equation method, based upon Green's formula, which yields the solution, , within the domain, given its value or that of its normal derivative at each point of the boundary (except possibly at a finite number of discrete points). The solution is obtained in two stages. The first, which is executed once only, determines the complementary boundary values, i.e., , where its normal derivative is known and vice versa. The second stage is entered once for each point at which the solution is required.
The boundary is divided into a number of intervals in each of which
and its normal derivative are approximated by constants. Of these half are evaluated by applying the given boundary conditions at one ‘nodal’ point within each interval while the remainder are determined (in stage
) by solving a set of simultaneous linear equations. Here this is achieved by means of auxiliary routine
f08aaf, which will yield the least squares solution of an overdetermined system of equations as well as the unique solution of a square nonsingular system.
In exterior domains the solution behaves as as tends to infinity, where is a constant, is the total integral of the normal derivative around the boundary and is the radial distance from the origin of coordinates. For the Neumann problem (when the normal derivative is given along the whole boundary) is fixed by the boundary conditions whilst is chosen by you. However, for a Dirichlet problem (when is given along the whole boundary) or for a mixed problem, stage produces a value of for which ; then as tends to infinity the solution tends to the constant .
4
References
Symm G T and Pitfield R A (1974) Solution of Laplace's equation in two dimensions NPL Report NAC 44 National Physical Laboratory
5
Arguments
-
1:
– Logical
Input
-
On entry: indicates whether the routine call is for stage
of the computation as defined in
Section 3.
- The call is for stage .
- The call is for stage .
-
2:
– Logical
Input
-
On entry: the form of the domain. If , the domain is unbounded. Otherwise the domain is an interior one.
-
3:
– Logical
Input
-
On entry: the form of the boundary conditions. If , the problem is a Dirichlet or mixed boundary value problem. Otherwise it is a Neumann problem.
-
4:
– Integer
Input
-
On entry: the number of intervals into which the boundary is divided (see
Sections 7 and
9).
-
5:
– Real (Kind=nag_wp)
Input
-
6:
– Real (Kind=nag_wp)
Input
-
On entry: to stage
,
p and
q must specify the
and
coordinates respectively of a point at which the solution is required.
When
stage1 is .TRUE.,
p and
q are ignored.
-
7:
– Real (Kind=nag_wp) array
Input
-
8:
– Real (Kind=nag_wp) array
Input
-
On entry: the
and
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.
-
9:
– Integer
Input
-
On entry: the value , where denotes the number of closed contours which make up the boundary.
-
10:
– Real (Kind=nag_wp) array
Input/Output
-
On entry: for stage
,
phi must contain the nodal values of
or its normal derivative (into the domain) as prescribed in each interval. For stage
it must retain its output values from stage
.
On exit: from stage , it contains the constants which approximate in each interval. It remains unchanged on exit from stage .
-
11:
– Real (Kind=nag_wp) array
Input/Output
-
On entry: for stage , must hold the value or accordingly as contains a value of or its normal derivative, for . For stage it must retain its output values from stage .
On exit: from stage
,
phid contains the constants which approximate the normal derivative of
in each interval. It remains unchanged on exit from stage
.
-
12:
– Real (Kind=nag_wp)
Input/Output
-
On entry: for stage
, the use of
alpha depends on the nature of the problem:
- if , alpha need not be set;
- if and , alpha must contain the prescribed constant (see Section 3);
- if and , alpha must contain an appropriate value (often zero) for the integral of around the outermost boundary.
For stage
, on every call
alpha must contain the value returned at stage
.
On exit: from stage
:
- if , alpha contains ;
- if and alpha is unchanged;
- if and alpha contains a computed estimate for .
From stage
:
- alpha contains the computed value of at the point (p,q).
-
13:
– Real (Kind=nag_wp) array
Workspace
-
14:
– Integer
Input
-
On entry: the first dimension of the array
c as declared in the (sub)program from which
d03eaf is called.
Constraint:
.
-
15:
– Integer
Input
-
On entry: this argument is no longer referenced, but is included for backwards compatability.
-
16:
– Integer array
Workspace
-
-
17:
– Integer
Input
-
On entry: this argument is no longer referenced, but is included for backwards compatability.
-
18:
– Integer
Input/Output
-
On entry:
ifail must be set to
,
or
to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of means that an error message is printed while a value of means that it is not.
If halting is not appropriate, the value
or
is recommended. If message printing is undesirable, then the value
is recommended. Otherwise, the value
is recommended.
When the value or is used it is essential to test the value of ifail on exit.
On exit:
unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
or
, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
-
Invalid tolerance used in an internal call to an auxiliary routine :
- Indicates too large a tolerance.
- 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 routine again.
-
Incorrect rank obtained by an auxiliary routine; is the computed rank.
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 7 in the Introduction to the NAG Library FL Interface for further information.
Your licence key may have expired or may not have been installed correctly.
See
Section 8 in the Introduction to the NAG Library FL Interface for further information.
Dynamic memory allocation failed.
See
Section 9 in the Introduction to the NAG Library FL Interface for further information.
7
Accuracy
The accuracy of the computed solution depends upon how closely and its normal derivative may be approximated by constants in each interval of the boundary and upon how well the boundary contours are represented by polygons with vertices at the selected points
, for
.
Consequently, in general, the accuracy increases as the boundary is subdivided into smaller and smaller intervals and by comparing solutions for successive subdivisions one may obtain an indication of the error in these solutions.
Alternatively, since the point of maximum error always lies on the boundary of the domain, an estimate of the error may be obtained by computing at a sufficient number of points on the boundary where the true solution is known. The latter method (not applicable to the Neumann problem) is most useful in the case where alone is prescribed on the boundary (the Dirichlet problem).
8
Parallelism and Performance
d03eaf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
d03eaf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the
X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the
Users' Note for your implementation for any additional implementation-specific information.
The time taken for stage 1, which is executed once only, is roughly proportional to
, being dominated by the time taken to compute the coefficients. The time for each stage
application is proportional to
n.
The intervals into which the boundary is divided need not be of equal lengths.
For many practical problems useful results may be obtained with to intervals per boundary contour.
10
Example
An interior Neumann problem to solve Laplace's equation in the domain bounded externally by the triangle with vertices
,
and
, and internally by the triangle with vertices
, (
) and
, given that the normal derivative of the solution
is zero on each side of each triangle and, for uniqueness that the total integral of
around the outer triangle is
(the length of the contour).
This trivial example has the obvious solution throughout the domain. However it provides a useful illustration of data input for a doubly-connected domain. The solution is computed at one corner of each triangle and at one point inside the domain.
The program is written to handle any of the different types of problem that the routine can solve. The array dimensions must be increased for larger problems.
10.1
Program Text
10.2
Program Data
10.3
Program Results