hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_pde_2d_laplace (d03ea)

Purpose

nag_pde_2d_laplace (d03ea) 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.

Syntax

[phi, phid, alpha, ifail] = d03ea(stage1, ext, dorm, p, q, x, y, phi, phid, alpha, 'n', n, 'n1p1', n1p1)
[phi, phid, alpha, ifail] = nag_pde_2d_laplace(stage1, ext, dorm, p, q, x, y, phi, phid, alpha, 'n', n, 'n1p1', n1p1)

Description

nag_pde_2d_laplace (d03ea) 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 11) by solving a set of simultaneous linear equations. Here this is achieved by means of auxiliary functions nag_eigen_real_gen_qu_svd (f02wd) and nag_linsys_real_gen_solve (f04jg), 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 c + s log(r) + O(1 / r)c+s log(r) + O(1/r) as rr tends to infinity, where cc is a constant, ss is the total integral of the normal derivative around the boundary and rr is the radial distance from the origin of coordinates. For the Neumann problem (when the normal derivative is given along the whole boundary) ss is fixed by the boundary conditions whilst cc is chosen by you. However, for a Dirichlet problem (when φϕ is given along the whole boundary) or for a mixed problem, stage 11 produces a value of cc for which s = 0s=0; then as rr tends to infinity the solution tends to the constant cc.

References

Symm G T and Pitfield R A (1974) Solution of Laplace's equation in two dimensions NPL Report NAC 44 National Physical Laboratory

Parameters

Compulsory Input Parameters

1:     stage1 – logical scalar
Indicates whether the function call is for stage 11 of the computation as defined in Section [Description].
stage1 = truestage1=true
The call is for stage 11.
stage1 = falsestage1=false
The call is for stage 22.
2:     ext – logical scalar
The form of the domain. If ext = trueext=true, the domain is unbounded. Otherwise the domain is an interior one.
3:     dorm – logical scalar
The form of the boundary conditions. If dorm = truedorm=true, then the problem is a Dirichlet or mixed boundary value problem. Otherwise it is a Neumann problem.
4:     p – double scalar
5:     q – double scalar
To stage 22, p and q must specify the xx and yy coordinates respectively of a point at which the solution is required.
When stage1 is true, p and q are ignored.
6:     x(n1p1) – double array
7:     y(n1p1) – double array
The xx and yy 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:     phi(n) – double array
For stage 11, phi must contain the nodal values of φϕ or its normal derivative (into the domain) as prescribed in each interval. For stage 22 it must retain its output values from stage 11.
9:     phid(n) – double array
For stage 11, phid(i)phidi must hold the value 0.00.0 or 1.01.0 accordingly as phi(i)phii contains a value of φϕ or its normal derivative, for i = 1,2,,ni=1,2,,n. For stage 22 it must retain its output values from stage 11.
10:   alpha – double scalar
For stage 11, the use of alpha depends on the nature of the problem:
  • if dorm = truedorm=true, alpha need not be set;
  • if dorm = falsedorm=false and ext = trueext=true, alpha must contain the prescribed constant cc (see Section [Description]);
  • if dorm = falsedorm=false and ext = falseext=false, alpha must contain an appropriate value (often zero) for the integral of φϕ around the outermost boundary.
For stage 22, on every call alpha must contain the value returned at stage 11.

Optional Input Parameters

1:     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 Sections [Accuracy] and [Further Comments]).
2:     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(n + M)12(n+M)-1, where MM denotes the number of closed contours which make up the boundary.

Input Parameters Omitted from the MATLAB Interface

c ldc np4 icint np1

Output Parameters

1:     phi(n) – double array
From stage 11, it contains the constants which approximate φϕ in each interval. It remains unchanged on exit from stage 22.
2:     phid(n) – double array
From stage 11, phid contains the constants which approximate the normal derivative of φϕ in each interval. It remains unchanged on exit from stage 22.
3:     alpha – double scalar
From stage 11:
  • if ext = falseext=false, alpha contains 0.00.0;
  • if ext = trueext=true and dorm = falsedorm=false alpha is unchanged;
  • if ext = trueext=true and dorm = truedorm=true alpha contains a computed estimate for cc.
From stage 22:
  • alpha contains the computed value of φϕ at the point (p,q).
4:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
Invalid tolerance used in an internal call to an auxiliary function:
icint(1) = 0icint1=0
Indicates too large a tolerance.
icint(1) > 0icint1>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.
  ifail = 2ifail=2
Incorrect rank obtained by an auxiliary function; icint(1)icint1 contains the computed rank.

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 (x(i),y(i)) (xi,yi) , for i = 1,2,,2(n + M)1i=1,2,,2(n+M)-1 .
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).

Further Comments

The time taken for stage 1, which is executed once only, is roughly proportional to n2n2, being dominated by the time taken to compute the coefficients. The time for each stage 22 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 2020 to 4040 intervals per boundary contour.

Example

function nag_pde_2d_laplace_example
stage1 = true;
ext = false;
dorm = false;
p = 0;
q = 0;
x = [3;
     1.5;
     0;
     -1.5;
     -3;
     0;
     3;
     3;
     2;
     0;
     -2;
     -1;
     0;
     1;
     2];
y = [0;
     2;
     4;
     2;
     0;
     0;
     0;
     0;
     1;
     1;
     1;
     2;
     3;
     2;
     1];
phi = [0;
     0;
     0;
     0;
     0;
     0];
phid = [1;
     1;
     1;
     1;
     1;
     1];
alpha = 16;
[phiOut, phidOut, alphaOut, ifail] = ...
    nag_pde_2d_laplace(stage1, ext, dorm, p, q, x, y, phi, phid, alpha)
 

phiOut =

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000


phidOut =

     0
     0
     0
     0
     0
     0


alphaOut =

     0


ifail =

                    0


function d03ea_example
stage1 = true;
ext = false;
dorm = false;
p = 0;
q = 0;
x = [3;
     1.5;
     0;
     -1.5;
     -3;
     0;
     3;
     3;
     2;
     0;
     -2;
     -1;
     0;
     1;
     2];
y = [0;
     2;
     4;
     2;
     0;
     0;
     0;
     0;
     1;
     1;
     1;
     2;
     3;
     2;
     1];
phi = [0;
     0;
     0;
     0;
     0;
     0];
phid = [1;
     1;
     1;
     1;
     1;
     1];
alpha = 16;
[phiOut, phidOut, alphaOut, ifail] = ...
    d03ea(stage1, ext, dorm, p, q, x, y, phi, phid, alpha)
 

phiOut =

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000


phidOut =

     0
     0
     0
     0
     0
     0


alphaOut =

     0


ifail =

                    0



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013