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_ellip_discret (d03ee)

Purpose

nag_pde_2d_ellip_discret (d03ee) discretizes a second-order elliptic partial differential equation (PDE) on a rectangular region.

Syntax

[a, rhs, ifail] = d03ee(xmin, xmax, ymin, ymax, pdef, bndy, ngx, ngy, scheme)
[a, rhs, ifail] = nag_pde_2d_ellip_discret(xmin, xmax, ymin, ymax, pdef, bndy, ngx, ngy, scheme)

Description

nag_pde_2d_ellip_discret (d03ee) discretizes a second-order linear elliptic partial differential equation of the form
α(x,y) (2U)/(x2) + β(x,y) (2U)/(x y) + γ(x,y) (2U)/(y2) + δ(x,y) (U)/(x) + ε(x,y)(U)/(y) + φ(x,y) U = ψ(x,y)
α(x,y) 2U x2 + β(x,y) 2U x y + γ(x,y) 2U y2 + δ(x,y) U x +ε(x,y) U y + ϕ(x,y) U = ψ(x,y)
(1)
on a rectangular region
xAxxB
yAyyB
xAxxB yAyyB
subject to boundary conditions of the form
a(x,y) U + b(x,y) (U)/(n) = c(x,y)
a(x,y) U+ b(x,y) U n = c(x,y)
where (U)/(n) U n  denotes the outward pointing normal derivative on the boundary. Equation (1) is said to be elliptic if
4 α(x,y) γ(x,y) (β(x,y))2
4 α(x,y) γ(x,y) ( β (x,y) ) 2
for all points in the rectangular region. The linear equations produced are in a form suitable for passing directly to the multigrid function nag_pde_2d_ellip_mgrid (d03ed).
The equation is discretized on a rectangular grid, with nxnx grid points in the xx-direction and nyny grid points in the yy-direction. The grid spacing used is therefore
hx = (xBxA) / (nx1)
hy = (yByA) / (ny1)
hx=(xB-xA)/(nx-1) hy=(yB-yA)/(ny-1)
and the coordinates of the grid points (xi,yj)(xi,yj) are
xi = xA + (i1)hx,  i = 1,2,,nx,
yj = yA + (j1)hy,  j = 1,2,,ny.
xi=xA+(i-1)hx,  i=1,2,,nx, yj=yA+(j-1)hy,  j=1,2,,ny.
At each grid point (xi,yj)(xi,yj) six neighbouring grid points are used to approximate the partial differential equation, so that the equation is discretized on the seven-point stencil shown in Figure 1.
Figure 1
Figure 1
For convenience the approximation uij uij  to the exact solution U(xi,yj) U(xi,yj)  is denoted by uOuO, and the neighbouring approximations are labelled according to points of the compass as shown. Where numerical labels for the seven points are required, these are also shown.
The following approximations are used for the second derivatives:
(2U)/(x2) 1/(hx2)(uE2uO + uW)
(2U)/(y2) 1/(hy2)(uN2uO + uS)
(2U)/(xy) 1/(2hxhy)(uNuNW + uE2uO + uWuSE + uS).
2U x2 1hx2 (uE-2uO+uW) 2U y2 1hy2 (uN-2uO+uS) 2U xy 12hxhy (uN-uNW+uE-2uO+uW-uSE+uS).
Two possible schemes may be used to approximate the first derivatives:
Central Differences
(U)/(x) 1/(2hx)(uEuW)
(U)/(y) 1/(2hy)(uNuS)
U x 12hx (uE-uW) U y 12hy (uN-uS)
Upwind Differences
(U)/(x) 1/(hx)(uOuW) if  δ(x,y) > 0
(U)/(x) 1/(hx)(uEuO) if  δ(x,y) < 0
(U)/(y) 1/(hy)(uNuO) if  ε(x,y) > 0
(U)/(y) 1/(hy)(uOuS) if  ε(x,y) < 0.
U x 1hx(uO-uW) if  δ(x,y)>0 U x 1hx(uE-uO) if  δ(x,y)<0 U y 1hy(uN-uO) if  ε(x,y)>0 U y 1hy(uO-uS) if  ε(x,y)<0.
Central differences are more accurate than upwind differences, but upwind differences may lead to a more diagonally dominant matrix for those problems where the coefficients of the first derivatives are significantly larger than the coefficients of the second derivatives.
The approximations used for the first derivatives may be written in a more compact form as follows:
(U)/(x) 1/(2hx) ((kx1)uW2kxuO + (kx + 1)uE)
(U)/(y) 1/(2hy) ((ky1)uS2kyuO + (ky + 1)uN)
U x 12hx ((kx- 1)uW- 2kxuO+(kx+ 1)uE) U y 12hy ((ky- 1)uS- 2kyuO+(ky+ 1)uN)
where kx = signδkx=signδ and ky = signεky=signε for upwind differences, and kx = ky = 0kx=ky=0 for central differences.
At all points in the rectangular domain, including the boundary, the coefficients in the partial differential equation are evaluated by calling pdef, and applying the approximations. This leads to a seven-diagonal system of linear equations of the form:
Aij6ui1,j + 1 + Aij7ui,j + 1
+ Aij3ui1,j + Aij4uij + Aij5ui + 1,j
+ Aij1ui,j1 + Aij2ui + 1,j1 = fij,  i = 1,2,,nx​ and ​j = 1,2,,ny,
Aij6ui-1,j+1 + Aij7ui,j+1 + Aij3ui-1,j + Aij4uij + Aij5ui+1,j + Aij1ui,j-1 + Aij2ui+1,j-1=fij,  i=1,2,,nx​ and ​j=1,2,,ny,
where the coefficients are given by
Aij1 = β (xi,yj)1/(2hxhy) + γ (xi,yj)1/(hy2) + ε (xi,yj)1/(2hy)(ky1)
Aij2 = β (xi,yj)1/(2hxhy)
Aij3 = α (xi,yj)1/(hx2) + β (xi,yj)1/(2hxhy) + δ (xi,yj)1/(2hx)(kx1)
Aij4 = α (xi,yj)2/(hx2)β (xi,yj)1/(hxhy)γ (xi,yj)2/(hy2)δ (xi,yj)(ky)/(hx)ε (xi,yj)(ky)/(hy)φ (xi,yj)
Aij5 = α (xi,yj)1/(hx2) + β (xi,yj)1/(2hxhy) + δ (xi,yj)1/(2hx)(kx + 1)
Aij6 = β (xi,yj)1/(2hxhy)
Aij7 = β (xi,yj)1/(2hxhy) + γ (xi,yj)1/(hy2) + ε (xi,yj)1/(2hy)(ky + 1)
fij = ψ (xi,yj)
Aij1 = β (xi,yj)12hxhy +γ (xi,yj)1hy2 +ε (xi,yj)12hy (ky- 1) Aij2 = -β (xi,yj)12hxhy Aij3 = α (xi,yj)1hx2 +β (xi,yj)12hxhy +δ (xi,yj)12hx (kx- 1) Aij4 = -α (xi,yj)2hx2 -β (xi,yj)1hxhy -γ (xi,yj)2hy2 -δ (xi,yj)kyhx-ε (xi,yj)kyhy-ϕ (xi,yj) Aij5 = α (xi,yj)1hx2 +β (xi,yj)12hxhy +δ (xi,yj)12hx (kx+ 1) Aij6 = -β (xi,yj)12hxhy Aij7 = β (xi,yj)12hxhy +γ (xi,yj)1hy2 +ε (xi,yj)12hy (ky+ 1) fij = ψ (xi,yj)
These equations then have to be modified to take account of the boundary conditions. These may be Dirichlet (where the solution is given), Neumann (where the derivative of the solution is given), or mixed (where a linear combination of solution and derivative is given).
If the boundary conditions are Dirichlet, there are an infinity of possible equations which may be applied:
μuij = μfij , μ0 .
μuij = μfij , μ0 .
(2)
If nag_pde_2d_ellip_mgrid (d03ed) is used to solve the discretized equations, it turns out that the choice of μμ can have a dramatic effect on the rate of convergence, and the obvious choice μ = 1μ=1 is not the best. Some choices may even cause the multigrid method to fail altogether. In practice it has been found that a value of the same order as the other diagonal elements of the matrix is best, and the following value has been found to work well in practice:
μ = minij ({2/(hx2) + 2/(hy2)},Aij4) .
μ=minij (-{2hx2 +2hy2 } ,Aij4) .
If the boundary conditions are either mixed or Neumann (i.e., B0B0 on return from bndy), then one of the points in the seven-point stencil lies outside the domain. In this case the normal derivative in the boundary conditions is used to eliminate the ‘fictitious’ point, uoutsideuoutside:
(U)/(n)1/(2h)(uoutsideuinside).
U n 12h (uoutside-uinside).
(3)
It should be noted that if the boundary conditions are Neumann and φ(x,y)0ϕ(x,y)0, then there is no unique solution. The function returns with ifail = 5ifail=5 in this case, and the seven-diagonal matrix is singular.
The four corners are treated separately. bndy is called twice, once along each of the edges meeting at the corner. If both boundary conditions at this point are Dirichlet and the prescribed solution values agree, then this value is used in an equation of the form (2). If the prescribed solution is discontinuous at the corner, then the average of the two values is used. If one boundary condition is Dirichlet and the other is mixed, then the value prescribed by the Dirichlet condition is used in an equation of the form given above. Finally, if both conditions are mixed or Neumann, then two ‘fictitious’ points are eliminated using two equations of the form (3).
It is possible that equations for which the solution is known at all points on the boundary, have coefficients which are not defined on the boundary. Since this function calls pdef at all points in the domain, including boundary points, arithmetic errors may occur in pdef which this function cannot trap. If you have an equation with Dirichlet boundary conditions (i.e., B = 0B=0 at all points on the boundary), but with PDE coefficients which are singular on the boundary, then nag_pde_2d_ellip_mgrid (d03ed) could be called directly only using interior grid points at your discretization.
After the equations have been set up as described above, they are checked for diagonal dominance. That is to say,
|Aij4| > k4 |Aijk| ,   i = 1,2,,nx ​ and ​ j = 1,2,,ny .
| Aij4 | > k4 | Aijk | ,   i=1,2,,nx ​ and ​ j=1,2,,ny .
If this condition is not satisfied then the function returns with ifail = 6ifail=6. The multigrid functionnag_pde_2d_ellip_mgrid (d03ed) may still converge in this case, but if the coefficients of the first derivatives in the partial differential equation are large compared with the coefficients of the second derivative, you should consider using upwind differences (scheme = 'U'scheme='U').
Since this function is designed primarily for use with nag_pde_2d_ellip_mgrid (d03ed), this document should be read in conjunction with the document for that function.

References

Wesseling P (1982) MGD1 – a robust and efficient multigrid method Multigrid Methods. Lecture Notes in Mathematics 960 614–630 Springer–Verlag

Parameters

Compulsory Input Parameters

1:     xmin – double scalar
2:     xmax – double scalar
The lower and upper xx coordinates of the rectangular region respectively, xAxA and xBxB.
Constraint: xmin < xmaxxmin<xmax.
3:     ymin – double scalar
4:     ymax – double scalar
The lower and upper yy coordinates of the rectangular region respectively, yAyA and yByB.
Constraint: ymin < ymaxymin<ymax.
5:     pdef – function handle or string containing name of m-file
pdef must evaluate the functions α(x,y)α(x,y), β(x,y)β(x,y), γ(x,y)γ(x,y), δ(x,y)δ(x,y), ε(x,y)ε(x,y), φ(x,y)ϕ(x,y) and ψ(x,y)ψ(x,y) which define the equation at a general point (x,y)(x,y).
[alpha, beta, gamma, delta, epslon, phi, psi] = pdef(x, y)

Input Parameters

1:     x – double scalar
2:     y – double scalar
The xx and yy coordinates of the point at which the coefficients of the partial differential equation are to be evaluated.

Output Parameters

1:     alpha – double scalar
2:     beta – double scalar
3:     gamma – double scalar
4:     delta – double scalar
5:     epslon – double scalar
6:     phi – double scalar
7:     psi – double scalar
alpha, beta, gamma, delta, epslon, phi and psi must be set to the values of α(x,y)α(x,y), β(x,y)β(x,y), γ(x,y)γ(x,y), δ(x,y)δ(x,y), ε(x,y)ε(x,y), φ(x,y)ϕ(x,y) and ψ(x,y)ψ(x,y) respectively at the point specified by x and y.
6:     bndy – function handle or string containing name of m-file
bndy must evaluate the functions a(x,y)a(x,y), b(x,y)b(x,y), and c(x,y)c(x,y) involved in the boundary conditions.
[a, b, c] = bndy(x, y, ibnd)

Input Parameters

1:     x – double scalar
2:     y – double scalar
The xx and yy coordinates of the point at which the boundary conditions are to be evaluated.
3:     ibnd – int64int32nag_int scalar
Specifies on which boundary the point (x,y) lies. ibnd = 0ibnd=0, 11, 22 or 33 according as the point lies on the bottom, right, top or left boundary.

Output Parameters

1:     a – double scalar
2:     b – double scalar
3:     c – double scalar
a, b and c must be set to the values of the functions appearing in the boundary conditions.
7:     ngx – int64int32nag_int scalar
8:     ngy – int64int32nag_int scalar
The number of interior grid points in the xx- and yy-directions respectively, nxnx and nyny. If the seven-diagonal equations are to be solved by nag_pde_2d_ellip_mgrid (d03ed), then ngx1ngx-1 and ngy1ngy-1 should preferably be divisible by as high a power of 22 as possible.
Constraints:
  • ngx3ngx3;
  • ngy3ngy3.
9:     scheme – string (length ≥ 1)
The type of approximation to be used for the first derivatives which occur in the partial differential equation.
scheme = 'C'scheme='C'
Central differences are used.
scheme = 'U'scheme='U'
Upwind differences are used.
Constraint: scheme = 'C'scheme='C' or 'U''U'.
Note: generally speaking, if at least one of the coefficients multiplying the first derivatives (delta or epslon as returned by pdef) are large compared with the coefficients multiplying the second derivatives, then upwind differences may be more appropriate. Upwind differences are less accurate than central differences, but may result in more rapid convergence for strongly convective equations. The easiest test is to try both schemes.

Optional Input Parameters

None.

Input Parameters Omitted from the MATLAB Interface

lda

Output Parameters

1:     a(lda,77) – double array
lda = (4 × (ngx + 1) × (ngy + 1)) / 3lda=(4×(ngx+1)×(ngy+1))/3.
a(i,j)aij, for i = 1,2,,ngx × ngyi=1,2,,ngx×ngy and j = 1,2,,7j=1,2,,7, contains the seven-diagonal linear equations produced by the discretization described above. If lda > ngx × ngylda>ngx×ngy, the remaining elements are not referenced by the function, but if lda(4 × (ngx + 1) × (ngy + 1)) / 3lda(4×(ngx+1)×(ngy+1))/3 then the array a can be passed directly to nag_pde_2d_ellip_mgrid (d03ed), where these elements are used as workspace.
2:     rhs(lda) – double array
lda = (4 × (ngx + 1) × (ngy + 1)) / 3lda=(4×(ngx+1)×(ngy+1))/3.
The first ngx × ngyngx×ngy elements contain the right-hand sides of the seven-diagonal linear equations produced by the discretization described above. If lda > ngx × ngylda>ngx×ngy, the remaining elements are not referenced by the function, but if lda(4 × (ngy + 1) × (ngy + 1)) / 3lda(4×(ngy+1)×(ngy+1))/3 then the array rhs can be passed directly to nag_pde_2d_ellip_mgrid (d03ed), where these elements are used as workspace.
3:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Note: nag_pde_2d_ellip_discret (d03ee) may return useful information for one or more of the following detected errors or 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 = 1ifail=1
On entry,xminxmaxxminxmax,
oryminymaxyminymax,
orngx < 3ngx<3,
orngy < 3ngy<3,
orlda < ngx × ngylda<ngx×ngy,
orscheme is not one of 'C' or 'U'.
  ifail = 2ifail=2
At some point on the boundary there is a derivative in the boundary conditions (b0b0 on return from bndy) and there is a nonzero coefficient of the mixed derivative (2U)/(xy) 2U xy  (beta0beta0 on return from pdef).
  ifail = 3ifail=3
A null boundary has been specified, i.e., at some point both a and b are zero on return from a call to bndy.
W ifail = 4ifail=4
The equation is not elliptic, i.e., 4 × alpha × gamma < beta24×alpha×gamma<beta2 after a call to pdef. The discretization has been completed, but the convergence of nag_pde_2d_ellip_mgrid (d03ed) cannot be guaranteed.
W ifail = 5ifail=5
The boundary conditions are purely Neumann (only the derivative is specified) and there is, in general, no unique solution.
W ifail = 6ifail=6
The equations were not diagonally dominant. (See Section [Description].)

Accuracy

Not applicable.

Further Comments

If this function is used as a preprocessor to the multigrid function nag_pde_2d_ellip_mgrid (d03ed) it should be noted that the rate of convergence of that function is strongly dependent upon the number of levels in the multigrid scheme, and thus the choice of ngx and ngy is very important.

Example

function nag_pde_2d_ellip_discret_example
xmin = 0;
xmax = 1;
ymin = 0;
ymax = 1;
ngx = int64(9);
ngy = int64(9);
scheme = 'Central';
[a, rhs, ifail] = ...
    nag_pde_2d_ellip_discret(xmin, xmax, ymin, ymax, @pdef, @bndy, ngx, ngy, scheme);
 a(1:81,:)
 rhs(1:81)

function [alpha, beta, gamma, delta, epsilon, phi, psi] = pdef(x,y)

      alpha = 1;
      beta = 0;
      gamma = 1;
      delta = 50;
      epsilon = 50;
      phi = 0;
      psi = (-alpha-gamma+phi)*sin(x)*sin(y) + beta*cos(x)*cos(y) + ...
                                   delta*cos(x)*sin(y) + epsilon*sin(x)*cos(y);

function [a, b, c] = bndy(x, y, ibnd)
      if (ibnd == 2 || ibnd == 1)
%        Solution prescribed
         a = 1.0d0;
         b = 0.0d0;
         c = sin(x)*sin(y);
      elseif (ibnd == 0)
%        Derivative prescribed
         a = 0.0d0;
         b = 1.0d0;
         c = -sin(x);
      elseif (ibnd == 3)
%        Derivative prescribed
         a = 0.0d0;
         b = 1.0d0;
         c = -sin(y);
      end
 
Warning: nag_pde_2d_ellip_discret (d03ee) returned a warning indicator (6) 

ans =

     0     0     0  -256   128     0   128
     0     0  -136  -256   264     0   128
     0     0  -136  -256   264     0   128
     0     0  -136  -256   264     0   128
     0     0  -136  -256   264     0   128
     0     0  -136  -256   264     0   128
     0     0  -136  -256   264     0   128
     0     0  -136  -256   264     0   128
     0     0     0  -256     0     0     0
  -136     0     0  -256   128     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
     0     0     0  -256     0     0     0
  -136     0     0  -256   128     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
     0     0     0  -256     0     0     0
  -136     0     0  -256   128     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
     0     0     0  -256     0     0     0
  -136     0     0  -256   128     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
     0     0     0  -256     0     0     0
  -136     0     0  -256   128     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
     0     0     0  -256     0     0     0
  -136     0     0  -256   128     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
     0     0     0  -256     0     0     0
  -136     0     0  -256   128     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
     0     0     0  -256     0     0     0
     0     0     0  -256     0     0     0
     0     0     0  -256     0     0     0
     0     0     0  -256     0     0     0
     0     0     0  -256     0     0     0
     0     0     0  -256     0     0     0
     0     0     0  -256     0     0     0
     0     0     0  -256     0     0     0
     0     0     0  -256     0     0     0
     0     0     0  -256     0     0     0


ans =

         0
    1.9948
    3.9585
    5.8604
    7.6708
    9.3616
   10.9062
   12.2807
         0
    1.9948
   12.3391
   18.2519
   23.8799
   29.1353
   33.9360
   38.2072
   41.8822
  -26.8570
    3.9585
   18.2519
   23.8489
   29.0736
   33.8447
   38.0877
   41.7363
   44.7336
  -53.2949
    5.8604
   23.8799
   29.0736
   33.8136
   38.0260
   41.6449
   44.6140
   46.8870
  -78.9012
    7.6708
   29.1353
   33.8447
   38.0260
   41.6139
   44.5524
   46.7956
   48.3087
 -103.2762
    9.3616
   33.9360
   38.0877
   41.6449
   44.5524
   46.7646
   48.2470
   48.9766
 -126.0396
   10.9062
   38.2072
   41.7363
   44.6140
   46.7956
   48.2470
   48.9455
   48.8802
 -146.8363
   12.2807
   41.8822
   44.7336
   46.8870
   48.3087
   48.9766
   48.8802
   48.0211
 -165.3416
         0
  -26.8570
  -53.2949
  -78.9012
 -103.2762
 -126.0396
 -146.8363
 -165.3416
 -181.2668


function d03ee_example
xmin = 0;
xmax = 1;
ymin = 0;
ymax = 1;
ngx = int64(9);
ngy = int64(9);
scheme = 'Central';
[a, rhs, ifail] = ...
    d03ee(xmin, xmax, ymin, ymax, @pdef, @bndy, ngx, ngy, scheme);
 a(1:81,:)
 rhs(1:81)

function [alpha, beta, gamma, delta, epsilon, phi, psi] = pdef(x,y)

      alpha = 1;
      beta = 0;
      gamma = 1;
      delta = 50;
      epsilon = 50;
      phi = 0;
      psi = (-alpha-gamma+phi)*sin(x)*sin(y) + beta*cos(x)*cos(y) + ...
                                   delta*cos(x)*sin(y) + epsilon*sin(x)*cos(y);

function [a, b, c] = bndy(x, y, ibnd)
      if (ibnd == 2 || ibnd == 1)
%        Solution prescribed
         a = 1.0d0;
         b = 0.0d0;
         c = sin(x)*sin(y);
      elseif (ibnd == 0)
%        Derivative prescribed
         a = 0.0d0;
         b = 1.0d0;
         c = -sin(x);
      elseif (ibnd == 3)
%        Derivative prescribed
         a = 0.0d0;
         b = 1.0d0;
         c = -sin(y);
      end
 
Warning: nag_pde_2d_ellip_discret (d03ee) returned a warning indicator (6) 

ans =

     0     0     0  -256   128     0   128
     0     0  -136  -256   264     0   128
     0     0  -136  -256   264     0   128
     0     0  -136  -256   264     0   128
     0     0  -136  -256   264     0   128
     0     0  -136  -256   264     0   128
     0     0  -136  -256   264     0   128
     0     0  -136  -256   264     0   128
     0     0     0  -256     0     0     0
  -136     0     0  -256   128     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
     0     0     0  -256     0     0     0
  -136     0     0  -256   128     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
     0     0     0  -256     0     0     0
  -136     0     0  -256   128     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
     0     0     0  -256     0     0     0
  -136     0     0  -256   128     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
     0     0     0  -256     0     0     0
  -136     0     0  -256   128     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
     0     0     0  -256     0     0     0
  -136     0     0  -256   128     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
     0     0     0  -256     0     0     0
  -136     0     0  -256   128     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
  -136     0  -136  -256   264     0   264
     0     0     0  -256     0     0     0
     0     0     0  -256     0     0     0
     0     0     0  -256     0     0     0
     0     0     0  -256     0     0     0
     0     0     0  -256     0     0     0
     0     0     0  -256     0     0     0
     0     0     0  -256     0     0     0
     0     0     0  -256     0     0     0
     0     0     0  -256     0     0     0
     0     0     0  -256     0     0     0


ans =

         0
    1.9948
    3.9585
    5.8604
    7.6708
    9.3616
   10.9062
   12.2807
         0
    1.9948
   12.3391
   18.2519
   23.8799
   29.1353
   33.9360
   38.2072
   41.8822
  -26.8570
    3.9585
   18.2519
   23.8489
   29.0736
   33.8447
   38.0877
   41.7363
   44.7336
  -53.2949
    5.8604
   23.8799
   29.0736
   33.8136
   38.0260
   41.6449
   44.6140
   46.8870
  -78.9012
    7.6708
   29.1353
   33.8447
   38.0260
   41.6139
   44.5524
   46.7956
   48.3087
 -103.2762
    9.3616
   33.9360
   38.0877
   41.6449
   44.5524
   46.7646
   48.2470
   48.9766
 -126.0396
   10.9062
   38.2072
   41.7363
   44.6140
   46.7956
   48.2470
   48.9455
   48.8802
 -146.8363
   12.2807
   41.8822
   44.7336
   46.8870
   48.3087
   48.9766
   48.8802
   48.0211
 -165.3416
         0
  -26.8570
  -53.2949
  -78.9012
 -103.2762
 -126.0396
 -146.8363
 -165.3416
 -181.2668



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