The NAG Library is the largest commercially-available collection of numerical and statistical algorithms in the world. With over 1,600 tried and tested routines that are both flexible and portable it can be found at the core of thousands of programs and applications around the world. In this talk I will describe the work that goes into producing the Library, with a particular focus on the mathematical and computer science challenges that we face. These begin with how we identify new material to include in the Library, either internally or in conjunction with our numerous external collaborators. We then have to design the best way to make this material available to users of the Library. This requires providing interfaces to allow routines to be called from a variety of programming languages, and clear documentation that explains what the routines do and how to call them. Finally we need to ensure that our software remains accurate and reliable, and I will discuss the approaches we follow at NAG to ensure this.

As more diseases are targetted for eradication by the global community, the question arises of how best to target our efforts. Mathematical modelling can be an important tool in addressing these issues. I will present work on two case studies: household modelling in yaws eradication efforts; and the challenge of systematic non-adherence in mass drug administration. Yaws is a painful and disabling infectious disease, with a latent form of infection that complicates eradication efforts. We use mathematical modelling along with data from the Solomon Islands to analyse the efficacy of the new WHO guidelines for yaws eradication. 

It is well understood that the success or failure of a disease eradication campaign critically depends on the level of coverage achieved. However, the quality of the coverage achieved also affects the efficacy of the campaign: in particular, if the same people miss treatment every year, they can represent a reservoir of infection, posing a barrier to eradication efforts. We demonstrate that the modelling description used and the level of systematic non-adherence can have a profound effect on the outcome of a treatment campaign. 

Inverse problems are ubiquitous in many areas of Science and Engineering and, once discretised, they lead to ill-conditioned linear systems, often of huge dimensions: regularisation consists in replacing the original system by a nearby problem with better numerical properties, in order to find a meaningful approximation of its solution. After briefly surveying some classical regularisation methods, both iterative (such as many Krylov methods) and direct (such as Tikhonov method), this talk will introduce the recent class of the Krylov-Tikhonov methods, which merge an iterative and a direct approach to regularisation. In particular, strategies for choosing the regularization parameter and the regularization matrix will be emphasized.

That the linear systems arising upon the discretization of elliptic PDEs can be solved very efficiently is well-known, and many successful iterative solvers with linear complexity have been constructed (multigrid, Krylov methods, etc). Interestingly, it has recently been demonstrated that it is often possible to directly compute an approximate inverse to the coefficient matrix in linear (or close to linear) time, despite the fact that the inverse is almost always dense. The talk will survey some recent work in the field and will argue that direct solvers have several advantages, including improved stability and robustness, and dramatic improvements in speed in certain environments. Moreover, the direct solvers being proposed have low communication costs, and are very well suited to parallel implementations.

The widely used perfectly matched layer (PML) method employs a complex coordinate transformation to impose non-reflecting boundary conditions in the finite element solution of wave equations. The coordinate transformation is ap- plied in a layer surrounding the region of interest (the bulk) and in its exact form completely absorbs waves leaving the bulk, eliminating spurious reflections from the outer boundary. However, when the problem is discretised, if the solution in the PML region is under-resolved, spurious reflections are created, limiting the achievable accuracy of the computed solution in the bulk. Conversely, if the solution in the PML region is over-resolved then computational effort is wasted. The optimal balance between bulk and PML refinement is problem dependent and difficult to find a priori.

To address this problem we propose a PML which is optimal in the sense that the coordinate transformation ensures the solution varies linearly through the PML, making it trivial to discretise. We present an algebraic method for find- ing this optimal transformation which utilises information about the solution. While this makes the problem non-linear, we show that we can converge to the exact solution by iterating, using information from the previous solution. This iteration is natural if we are already performing spatial mesh adaptation in the bulk. We show that with this optimal PML, the numerical error is completely controlled by the refinement in the bulk.

Many PDEs have time-periodic solutions and it is often of interest to explore their dependence on problem parameters. The determination of such solutions by direct time-integration can be very inefficient as transients may take a long time to decay. Furthermore, this method cannot robustly find unstable solutions, which may be of interest. In this talk we present a finite-element based space-time approach that allows the direct computation of time-periodic solutions. We demonstrate the application of the methodology for the forced unsteady heat equation and the diffusive Van der Pol equation. Finally, we illustrate its application to the study of flow past an oscillating cylinder.

There are many areas of science where the uncertainties are critical for modeling scientific phenomena such as: climate models, weather models and nuclear reactor designs. My work aims to study the numerical solution of PDEs with random data using sparse grid stochastic collocation methods with three different choices of interpolation abscissas. These choices are Clenshaw-Curtis, Chebyshev-Gauss-Lobatto and Gauss-Patterson. We study the numerical solution of these two models: the first is the diffusion equation where the diffusion coefficient is a spatially dependent random field. The second model is Navier-Stokes equations where the viscosity is a random variable. In this talk, I will introduce the numerical method that I used to solve these models. Then, I will present the numerical results and error estimates for each case.

In the past couple of years, we have experienced a significant number of pivotal events which continue to influence our current lives, such as the Syrian refugee crisis, Brexit and the presidential election of the U.S.A. Central to all the challenges following these events is the question of how polarisations between different groups occur, and if similar future events can be predicted. This inflexibility in voting or political behaviour has critical consequences when applied to criminal activity; how likely is a criminal to be rehabilitated or, conversely, more deeply entrenched in crime?

In this talk I will introduce a mathematical model for \textit{recidivism}, the tendency of a convicted criminal to re-offend. Supported by empirical data, a key element of criminal activity and recidivism will be discussed, the degree of rehabilitation. The model is capable of describing both likely re-offenders as well as expected durations of low criminal activity. These results inform the best punitive measures to ensure a minimal number of repeat offenders in the penal system, taking into account the nature and location of the crime.

As well known, the Navier-Stokes equations are the governing equations of fluid motion. The set of these equations consists of the conservation of momentum equation and the mass conservation. The major difficulty arises when solving the equations in primitive variables. This results from the lack of an independent equation of the pressure and the boundary conditions. Most of the past and current research are using the definition of the stream function and vorticity to eliminate the pressure from these equations, so they can obtain a numerical solution. However, the stream function and vorticity formulation remains difficult to extend to three dimensions. In our work, we aim to implement a suitable numerical method in order to obtain the numerical solution to the Navier-Stokes in three dimensions which requires the use of the primitive variables. Indeed, we plan to use a hybrid approach that combines the classic finite difference scheme with a spectral method, particularly the Chebyshev collocation. Therefore, we exploit the sparsity pattern of the finite differences and the high spectral accuracy of the Chebyshev collocation. This method, after linearisation, leads to end with a sparse linear system which then can be solved using a direct solver.

Uranium reacts with oxygen in the air to form uranium oxide and/or uranium hydride, with the latter being formed only in the presence of moisture. The corrosion product ura- nium hydride, formed from the reaction of uranium with water vapour is undesirable as it is pyrophoric (i.e. self-ignites), hence prompting the current study. The kinetics of uranium oxidation in dry air to form uranium oxide is dominated by the diffusion of oxygen anions, formed when the adsorbed oxygen atoms on the surface accepts electrons from the underlying metal. There are several empirical models available in the literature for the oxidation of ura- nium in dry and moist air, but very few mathematical models exist.

In our current research, a suitable mathematical model, with unsteady diffusion to be the rate-determining step, has been formulated and numerically solved. A boundary-fixing transformation has been used to solve the governing equation in a fixed domain at the expense of a more-complicated nonlinear governing equation. The overall process is one of combined temperature and time-dependent chemical and ionic diffusion within a growing substrate. In this talk, I will discuss the results obtained for the dry-air oxidation problem.

In the hyperelastic modelling of fibrous soft tissues, such as skin, we construct strain energy functions to model the anisotropic and nonlinear stress-strain behaviour that they exhibit. Phenomenological models can fit experimental data well. Parameters in mi- crostructural models, however, are connected to the properties and arrangement of the tissue’s constituents that influence the macroscopic behaviour of the tissue. Microstruc- tural models thus enable us to understand how the tissue deforms and, potentially, to measure the values of the parameters directly in experiments [1].

We, therefore, introduce a new model that assumes that fibres of collagen, a strong fibrous protein, in the skin are crimped (i.e. wavy) according to a triangular distribution. Crimp is an important property because we assume that collagen fibres straighten as the skin is stretched, and only contribute mechanically once straightened. To test the new model, we compare its fit to four data sets of uniaxial loading on mammalian skin to those of a commonly used phenomenological model and a microstructural tendon model [1].

The new model achieves the closest relative fit to three data sets, and the phenomenolo- gical model to one. The new model achieves a closer fit than the tendon model for each data set. The tendon model also outperforms the phenomenological model for three data sets. The number of fit parameters in the new model is then reduced to that of the pheno- menological model. The new model still produces a better fit than the phenomenological model for three of the four data sets.

[1] T. Shearer, A new strain energy function for the hyperelastic modelling of ligaments and tendons based on fascicle microstructure. J. Biomech, Vol. 48, pp. 290–297, 2015.

Modelling the transmission and reflection of waves from periodic arrays has important applications in a variety of fields. We discuss scattering of a plane wave from an array of cylinders using the boundary element method, which requires an efficient means of evaluating the periodic Green's function. A new approach using an asymptotic correction is discussed and results are shown for cylinders with a variety of cross-sections using the. We also briefly discuss the alternative plane wave expansion approach and related problems.