Using the NAG Toolbox for MATLAB® - Part 4

The NAG Toolbox makes the full functionality of the NAG Library available from within MATLAB. This article is the latest in a series which illustrates the Toolbox’s functionality (earlier articles are here, here and here). Previously, we have presented separate examples of the use of specific routines from the Toolbox; by contrast, this article discusses a single larger-scale application which makes use of several routines, drawn from a variety of different Toolbox chapters. The application has been made available as part of the latest release of the Toolbox – see below.

Introduction and background

We describe an application that employs a number of NAG Toolbox routines to solve a partial differential equation (PDE) over a two-dimensional region.  PDEs are used to build mathematical models of many different scientific phenomena – for example, the propagation of heat or sound, the behaviour of electromagnetic fields, the flow of a fluid and the deformation of a solid.  Currently our application is able to solve the so-called Poisson PDE, which can be used to find the electric potential corresponding to a given charge distribution; in our application the user can interactively define the functional form of the right-hand side.  It can also solve the Laplace equation (which is the Poisson PDE with a right-hand side of zero); this governs – amongst other things – the way heat is distributed in a given medium.

We determine a numerical solution for the PDE using a simple finite element method, in which the region for which the PDE is to be solved is subdivided into a mesh, and the PDE is transformed into a set of simultaneous linear equations that can be solved using algorithms for large sparse systems.  An advantage of the finite element method is that it is well-suited for regions of arbitrary shape; we illustrate this aspect in our application, which uses MATLAB’s GUI tools to input and display the region.

Our application makes use of a number of routines from the NAG Toolbox for generating the mesh and solving the equations.  We discuss these in more detail in the following sections.

The D06 Chapter – Mesh Generation

The D06 chapter of the NAG library contains routines for automatic mesh generation over a user-defined domain.  It includes a routine for creating a mesh of the domain boundary, a choice of three routines for generating a mesh within the domain, and several other utility routines that can be used to improve the quality of the mesh.

Our application allows the user to create their own domain boundary by either specifying a polygon as a series of points, or selecting a regular shape such as a rectangle or circle.  The boundary can be of arbitrary shape, and the domain can also contain internal boundaries, or holes (see, for example, Figure 5 below).  Having created the boundary, the application sets up and calls d06ba from the NAG Toolbox to subdivide the boundary of the domain, thereby creating the boundary mesh.

Figure 1 is a screenshot of our application, showing the domain boundary which has been specified by the user as a circle (which is represented inside the application as an 11-sided polygon). 

Screenshot of our application
Figure 1: Screenshot of our application, showing (in the drawing area at left) the boundary of the domain to be meshed.  Controls for creating the boundary are available in the menu bar at the top, whilst the area on the right allows the user to edit the boundary conditions (see text), specify the meshing algorithm and solve the PDE.  The box at the bottom right displays context-sensitive help.

Having meshed the boundary, the next step is to create the mesh within the domain.  As noted above, the NAG Toolbox contains three different meshing algorithms, labelled Incremental (d06aa), Delaunay-Voronoi (d06ab), or Advancing Front (d06ac). The parameters of each routine can be varied to adjust the coarseness of the mesh (see figure 2).  It is also possible – though not implemented in our application – when using the Delaunay-Voronoi and Advancing Front methods to specify some internal points of the domain where a finer mesh is required; this is useful in problems such as the determination of the turbulent wake of an aircraft wing.  Figures 2-4 show the interior meshes generated by the three algorithms for our example. Examples of an interior mesh
Figure 2: Three examples of an interior mesh generated using the Incremental method.  One parameter for this algorithm controls the size of the internal triangles; its value is decreased going from the mesh on the left to that in the middle.  Another parameter controls the distance between boundary nodes; comparing the mesh in the middle to that on the right shows the effect of decreasing its value.

The interior mesh generated using the Delaunay-Voronoi method
Figure 3: The interior mesh generated using the Delaunay-Voronoi method.

The interior mesh generated using the Advancing front method
Figure 4: The interior mesh generated using the Advancing front method.

The F04 Chapter - Simultaneous linear equations

The interior meshing algorithms can handle domains of any shape, but the boundary must be properly specified.  In particular, the boundary edges are not allowed to cross.  We test for this in our application while the points are being input by solving a pair of simultaneous linear equations, using the f04aa routine from the NAG Toolbox.

Here is the code used:


function [cross] = simeqn(coor, grad, i, j, x, y)

% (x,y) is the coordinate on the new edge of the boundary.
% coor is the coordinates of all the previous boundary edges.
% grad contains the gradients of all the boundary edges.
% i and j are the indices of interest.

% Create the input arrays for the NAG routine.
% a is the jacobian for this system of equations, and
% b is the array containing the equations of the edges.
a = [-grad(i-1), 1; -grad(j), 1];
b = [y - grad(i-1)*x; coor(2,j) - grad(j)*coor(1,j)];

% Call the NAG routine to solve the equations.
[a, c, ifail] = f04aa(a, b);

% Pull out the coordinates of the crossing point of the two
% lines from the solution 'c'.
xx = c(1);
yy = c(2);

% The next bit of code determines whether the solution of the
% two equations above is within the bounds of the polygon
% edges that we are interested in - i.e. whether the polygon edges
% are crossing each other.

ab = norm(coor(:,j) - [xx;yy]) + norm(coor(:,j+1) - [xx;yy]);
distab = norm(coor(:,j) - coor(:,j+1));
AB = norm([x;y] - [xx;yy]) + norm(coor(:,i-1) - [xx;yy]);
distAB = norm([x;y] - coor(:,i-1));

if ab < distab + 1e-10 && AB < distAB + 1e-10
   % the lines of the polygon are crossing
    cross = 1;
    % the lines aren't crossing
    cross = 0;


This is a very simple example of the use of a routine from the F04 chapter, which includes powerful algorithms for solving systems of linear equations of various types.  More specifically, the matrix of coefficients may be real, complex, symmetric, Hermitian, positive-definite, positive-definite Toeplitz or banded.  The matrix may also be rectangular, in which case a least-squares solution is obtained.  F04 contains two types of routine, called Black Box (which are easier to use), and General Purpose (which provide more flexibility, for example in the solution of families of equations with multiple right-hand sides).

The F04 chapter can still be used when the coefficient matrix is sparse (that is, contains a high proportion of zero elements), but it is often better in this case to use a routine from the F11 chapter, which is discussed in the following section.

The F11 Chapter - Large scale linear systems

Having set up the mesh over the domain of interest, the solution of the PDE is transformed into a large system of simultaneous linear equations in which the number of coefficients reflects the size of the mesh – typically of the order of thousands.  However, the matrix of coefficients is also banded (i.e., only the diagonal and a few off-diagonal elements are non-zero).  The system of equations is therefore sparse, and routines from the F11 chapter – which exploit this feature by storing only the non-zero elements – can be used to solve it. 

The important prerequisites of the solution of a PDE are the boundary conditions, which are restrictions on the behaviour of the solution (or its derivative) on the boundary of the domain.  The user of our application is able to specify values for either the solution (a so-called Dirichlet boundary condition) or its derivative (a Neumann boundary condition).  The boundary conditions are set up in the code preceding the calls of the Toolbox F11 routines.  Our example here is solving the Laplace equation for the distribution of heat across the domain, and we use Dirichlet boundary conditions, specifying the temperature on some boundaries, and setting the initial temperature across the whole domain (see Figure 5).

The F11 chapter contains several routines for solving large sparse linear systems.  As in the case of the F04 chapter, these can be divided into Black Box (easier to use) and General Purpose (more flexible).  F11 also includes a number of utility routines to help with the determination of the solution.  In our application, we first call the utility routine f11zb to reorder the non-zero elements of the matrix, which increases the efficiency of the solver.  We then use the pair of routines called f11ja and f11jc. The former performs the incomplete Cholesky factorisation of the matrix; again, this preconditioning step increases the speed and efficiency of the solver. The Black Box routine f11jc calls five other F11 routines to solve a real sparse symmetric linear system using either a conjugate gradient method (as used in our application) or the Lanczos method.

Having computed the solution, our application allows its dynamic evolution to be played back in real time using MATLAB’s plotting tools.  Figure 5 shows the final (steady state) result for the distribution of heat around a set of irregularly shaped holes whose edges had a high initial temperature, whilst the rest of the domain was cold. 

Solution of the Laplace equation
Figure 5: Solution of the Laplace equation for a region containing (irregularly shaped) holesThe internal boundary has a temperature of 100, while the outer edge of the region has no boundary condition.


The application described in this article has been made available as part of the latest release (Mark 22) of The NAG Toolbox for MATLAB.  It can be found in the Demos pane of the MATLAB Help browser under the Toolboxes entry (look for the demo called Mesh Generation in the NAG demos: Applications folder).  Note that the examples described in the previous articles in this series (see, e.g. here) have also been made available in the NAG demos: MATLAB GUI folder. 

The NAG Toolbox for MATLAB is available for both 32 bit and 64 bit versions of Linux and Microsoft Windows.  More information about its availability can be found here.