d03ee:: Partial Differential Equations (NAG Toolbox)

nag_pde_2d_ellip_discret (d03ee) discretizes a second-order linear elliptic partial differential equation of the form

α (x, y) \frac{\partial^{2} U}{\partial x^{2}} + β (x, y) \frac{\partial^{2} U}{\partial x \partial y} + γ (x, y) \frac{\partial^{2} U}{\partial y^{2}} + δ (x, y) \frac{\partial U}{\partial x} + ε (x, y) \frac{\partial U}{\partial y} + ϕ (x, y) U = ψ (x, y)

(1)

on a rectangular region

\begin{matrix} x_{A} \leq x \leq x_{B} \\ y_{A} \leq y \leq y_{B} \end{matrix}

subject to boundary conditions of the form

a (x, y) U + b (x, y) \frac{\partial U}{\partial n} = c (x, y)

where

\frac{\partial U}{\partial n}

denotes the outward pointing normal derivative on the boundary. Equation (1) is said to be elliptic if

4 α (x, y) γ (x, y) \geq {(β (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

n_{x}

grid points in the

x

-direction and

n_{y}

grid points in the

y

-direction. The grid spacing used is therefore

\begin{matrix} h_{x} = (x_{B} - x_{A}) / (n_{x} - 1) \\ h_{y} = (y_{B} - y_{A}) / (n_{y} - 1) \end{matrix}

and the coordinates of the grid points

(x_{i}, y_{j})

are

\begin{matrix} x_{i} = x_{A} + (i - 1) h_{x}, i = 1, 2, \dots, n_{x}, \\ y_{j} = y_{A} + (j - 1) h_{y}, j = 1, 2, \dots, n_{y} . \end{matrix}

The following approximations are used for the second derivatives:

\begin{array}{l} \frac{\partial^{2} U}{\partial x^{2}} & ≃ & \frac{1}{h_{x}^{2}} (u_{E} - 2 u_{O} + u_{W}) \\ \frac{\partial^{2} U}{\partial y^{2}} & ≃ & \frac{1}{h_{y}^{2}} (u_{N} - 2 u_{O} + u_{S}) \\ \frac{\partial^{2} U}{\partial x \partial y} & ≃ & \frac{1}{2 h_{x} h_{y}} (u_{N} - u_{NW} + u_{E} - 2 u_{O} + u_{W} - u_{SE} + u_{S}) . \end{array}

Upwind Differences

\begin{array}{l} \frac{\partial U}{\partial x} & ≃ & \frac{1}{h_{x}} (u_{O} - u_{W}) & if δ (x, y) > 0 \\ \frac{\partial U}{\partial x} & ≃ & \frac{1}{h_{x}} (u_{E} - u_{O}) & if δ (x, y) < 0 \\ \frac{\partial U}{\partial y} & ≃ & \frac{1}{h_{y}} (u_{N} - u_{O}) & if ε (x, y) > 0 \\ \frac{\partial U}{\partial y} & ≃ & \frac{1}{h_{y}} (u_{O} - u_{S}) & if ε (x, y) < 0 . \end{array}

The approximations used for the first derivatives may be written in a more compact form as follows:

\begin{array}{l} \frac{\partial U}{\partial x} & ≃ & \frac{1}{2 h_{x}} ((k_{x} - 1) u_{W} - 2 k_{x} u_{O} + (k_{x} + 1) u_{E}) \\ \frac{\partial U}{\partial y} & ≃ & \frac{1}{2 h_{y}} ((k_{y} - 1) u_{S} - 2 k_{y} u_{O} + (k_{y} + 1) u_{N}) \end{array}

where

k_{x} = sign δ

and

k_{y} = sign ε

for upwind differences, and

k_{x} = k_{y} = 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:

\begin{array}{l} A_{i j}^{6} u_{i - 1, j + 1} & + & A_{i j}^{7} u_{i, j + 1} \\ + & A_{i j}^{3} u_{i - 1, j} & + & A_{i j}^{4} u_{i j} & + & A_{i j}^{5} u_{i + 1, j} \\ + & A_{i j}^{1} u_{i, j - 1} & + & A_{i j}^{2} u_{i + 1, j - 1} = f_{i j}, i = 1, 2, \dots, n_{x} ​ and ​ j = 1, 2, \dots, n_{y}, \end{array}

where the coefficients are given by

\begin{array}{l} A_{i j}^{1} & = & β (x_{i}, y_{j}) \frac{1}{2 h_{x} h_{y}} + γ (x_{i}, y_{j}) \frac{1}{h_{y}^{2}} + ε (x_{i}, y_{j}) \frac{1}{2 h_{y}} (k_{y} - 1) \\ A_{i j}^{2} & = & - β (x_{i}, y_{j}) \frac{1}{2 h_{x} h_{y}} \\ A_{i j}^{3} & = & α (x_{i}, y_{j}) \frac{1}{h_{x}^{2}} + β (x_{i}, y_{j}) \frac{1}{2 h_{x} h_{y}} + δ (x_{i}, y_{j}) \frac{1}{2 h_{x}} (k_{x} - 1) \\ A_{i j}^{4} & = & - α (x_{i}, y_{j}) \frac{2}{h_{x}^{2}} - β (x_{i}, y_{j}) \frac{1}{h_{x} h_{y}} - γ (x_{i}, y_{j}) \frac{2}{h_{y}^{2}} - δ (x_{i}, y_{j}) \frac{k_{y}}{h_{x}} - ε (x_{i}, y_{j}) \frac{k_{y}}{h_{y}} - ϕ (x_{i}, y_{j}) \\ A_{i j}^{5} & = & α (x_{i}, y_{j}) \frac{1}{h_{x}^{2}} + β (x_{i}, y_{j}) \frac{1}{2 h_{x} h_{y}} + δ (x_{i}, y_{j}) \frac{1}{2 h_{x}} (k_{x} + 1) \\ A_{i j}^{6} & = & - β (x_{i}, y_{j}) \frac{1}{2 h_{x} h_{y}} \\ A_{i j}^{7} & = & β (x_{i}, y_{j}) \frac{1}{2 h_{x} h_{y}} + γ (x_{i}, y_{j}) \frac{1}{h_{y}^{2}} + ε (x_{i}, y_{j}) \frac{1}{2 h_{y}} (k_{y} + 1) \\ f_{i j} & = & ψ (x_{i}, y_{j}) \end{array}

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

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:

μ = \min_{i j} (- \{\frac{2}{h_{x}^{2}} + \frac{2}{h_{y}^{2}}\}, A_{i j}^{4}) .

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 = 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.

If this condition is not satisfied then the function returns with

ifail = 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'

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

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.

Accuracy

Further Comments

Example

The program solves the elliptic partial differential equation

\frac{\partial^{2} U}{\partial x^{2}} + \frac{\partial^{2} U}{\partial y^{2}} + 50 \{\frac{\partial U}{\partial x} + \frac{\partial U}{\partial y}\} = f (x, y)

on the unit square

0 \leq x

y \leq 1

, with boundary conditions

$\frac{\partial U}{\partial n}$ given on $x = 0$ and $y = 0$ ,
$U$ given on $x = 1$ and $y = 1$ .

function d03ee_example

% d03ee example: Solve elliptic equation by discretising using d02ee
%                and solving the disretised system using d03ed.


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

%  u= sin(kx)sin(ly)
global k l;
k = 4; l = 5;

% Discretise on unit square usind d03ee
xmin = 0; xmax = 1;
ymin = 0; ymax = 1;
ngx = int64(33); ngy = ngx;
scheme = 'Upwind';
[a, rhs, ifail] = ...
d03ee( ...
       xmin, xmax, ymin, ymax, @pdef, @bndy, ngx, ngy, scheme);

% Solve using d03ed
maxit = int64(100); acc = 0.0001; iout = int64(0);
ub = zeros(ngx*ngy, 1);
[a, rhs, ub, resid, u, numit, ifail] = ...
d03ed( ...
       ngx, ngy, a, rhs, ub, maxit, acc, iout);

fprintf('Number of iterations to solution = %2d.\n', numit);
fprintf('Maximum residual                   %8.1e\n',resid(numit));

% plot solution
u2 = u(1:ngx*ngy);
umat = reshape(u2,ngx,ngy);
hx = 1/double(ngx+1); 
hy = hx;
x(1:ngx) = hx:hx:1-hx;
y(1:ngy) = hy:hy:1-hy;
[xs,ys] = meshgrid(y,x);

fig1 = figure;
meshc(xs,ys,umat);

xlabel('x'); ylabel('y'); zlabel ('u(x,y)');
exact = sprintf('u = sin %dx sin %dy',k,l);
title({'d03ee example:  U_{xx}+U_{yy}+50(U_x+U_y) = f;',exact});
view(-24,28);

 
function [alpha, beta, gamma, delta, epsilon, phi, psi] = pdef(x,y)
  global k l;
  u = sin(k*x)*cos(l*y);
  ux = k*cos(k*x)*sin(l*y); uy = l*sin(k*x)*cos(l*y);
  uxx = -k*k*u; uyy = -l*l*u;
  uxy = k*l*cos(k*x)*cos(l*y);
  alpha = 1; beta = 0; gamma = 1; delta = 50; epsilon = 50;
  phi = 0;
  psi = alpha*uxx + beta*uxy + gamma*uyy + delta*ux + epsilon*uy + phi*u;


function [a, b, c] = bndy(x, y, ibnd)
  global k l;
  u = sin(k*x)*sin(l*y);
  if (ibnd == 2 || ibnd == 1)
    % Solution prescribed
    a = 1;
    b = 0;
    c = u;
  elseif (ibnd == 0)
    % Derivative prescribed: uy(y=0) = l*sin(k*x)
    a = 0;
    b = 1;
    c = -l*sin(k*x);
  elseif (ibnd == 3)
    % Derivative prescribed: ux(x=0) = k*sin(l*y)
    a = 0;
    b = 1;
    c = -k*sin(l*y);
  end

On entry,	$xmin \geq xmax$ ,
or	$ymin \geq ymax$ ,
or	$ngx < 3$ ,
or	$ngy < 3$ ,
or	$lda < ngx \times ngy$ ,
or	scheme is not one of 'C' or 'U'.

NAG Toolbox: nag_pde_2d_ellip_discret (d03ee)

▸▿ Contents

Purpose

Syntax

Description