Search

Electricity Markets are paramount to the economy of any country and therefore lie at the foundation of society. These markets are usually very complex when compared to traditional financial markets due to the physical nature of the asset traded, difficulty of storage and distribution. Supply and demand balancing, network constraints, renewable energy influx, desire of participants to maximise their own profits, their risk aversion, fuel costs required for production, and the market clearing mechanism, do not exhaust the list of factors that influence the market outcome, but are perhaps the most important ones. The market outcome includes the prices and quantities allocated whose understanding is of great importance for market regulators and market participants alike. This work is focused on improving the fundamental understanding of the Spot Electricity Markets which lie at the core of Electricity Trading and represent the majority of trading volume. In particular we focus on Day Ahead Auctions which represent the vast majority of Spot Electricity trading, due to the difficulty of storing electricity, coupled with uncertainty in demand and production patterns. Every market has its own peculiarities and, we consider as a case study the Central Western European (CWE) Day Ahead Auction (DAA) for the period of 01/01/2019−31/12/2019. The Central Western European market includes Austria, Belgium, Germany, France and the Netherlands. Following an introductory chapter, due to the lack of full and complete network constraints historical data, Chapter 2 is focused on the reconstruction of historical network constraints, that are a key input in any fundamental model considering the electricity grid's structure. Further, network constraints may impose price decoupling between countries, which is of great interest to practitioners, academics and regulators. The task is very challenging due to the input data being very large and sparse. We reconstruct the network constraints data, known as the Power Transmission Distribution Factors (PTDFs) and Remaining Available Margins (RAMs), by first recovering the underlying time dependent signals known as the Generation Shift Keys (GSKs) and Phase Angles (PAs), and the electricity grid characteristics, via a mathematical optimisation problem that is solved by exploiting problem structure. The second step maps GSKs and PAs to the network constraints via the electricity grid characteristics. Our reconstruction achieves good in- and out-of- sample relative errors when compared to a naive approach. With the network constraints available, and by obtaining the cost and capacity characteristics of suppliers, in Chapter 3 we formulate a fundamental model to resolve for the zonal prices and quantities assigned to each participant, based on a simplification of the CWE market clearing mechanism, known as the Social Welfare Maximisation Problem (SWM). We also propose a novel market clearing mechanism which we call the Total Cost Minimisation Problem, because as we show, the prices are generally much lower and the total cost of supplying electricity is guaranteed to be lower by definition. Solving the Cost Minimisation problem is more challenging than solving the SWM, due to zonal prices appearing as explicit variables, but we obtain strong (but partial) analytical results to aid numerical solutions and provide efficient numerical algorithms. For both clearing mechanisms, demand is considered as inflexible, and producers bid marginal prices linear in the quantities produced. Since as expected, the SWM better represents the target data, this is further used as the market clearing mechanism. We consider the case when producers truthfully bid (TB) their cost curves, which is usually assumed by the regulators and many academic papers. However, we also consider the case of players trying to maximise their own profit, tand thus forming a Game Theoretic (GT) framework. For our CWE case study, both GT and TB models generally give an accurate representation of the actual price outcomes. Further, we define a method based on simple hypothesis testing to identify if and when the players are bidding strategically (gaming). We show that in fact, while for a large majority of times we cannot say with high confidence whether the players are bidding strategically or bidding truthfully, the number of times we are confident that players are gaming is much higher than the number of times we are confident that the players are bidding truthfully. Our results show that, in general, strategic bidding appears to be especially pronounced at peak demand hours: between 8am − 1pm and 5pm − 9pm. A different regime, where renewable energy influx appears to dominate the price action is also identified and this effect is generally more pronounced at off-peak morning hours (0am − 6am). Our model can thus be used as a fundamental approach for market outcome and price prediction, understanding the impact various fundamental factors on the price, understanding of producers' behaviour, and identifying if and when players are not bidding truthfully. These are all of interest to academics, market participants, and regulators.

Discovering hidden partial differential equations (PDEs) and operators from data is an important topic at the frontier between machine learning and numerical analysis. Theoretical results and deep learning algorithms are introduced to learn Green's functions associated with linear partial differential equations and rigorously justify PDE learning techniques. A theoretically rigorous algorithm is derived to obtain a learning rate, which characterizes the amount of training data needed to approximately learn Green's functions associated with elliptic PDEs. The construction connects the fields of PDE learning and numerical linear algebra by extending the randomized singular value decomposition to non-standard Gaussian vectors and Hilbert--Schmidt operators, and exploiting the low-rank hierarchical structure of Green's functions using hierarchical matrices. Rational neural networks (NNs) are introduced and consist of neural networks with trainable rational activation functions. The highly compositional structure of these networks, combined with rational approximation theory, implies that rational functions have higher approximation power than standard activation functions. In addition, rational NNs may have poles and take arbitrarily large values, which is ideal for approximating functions with singularities such as Green's functions. Finally, theoretical results on Green's functions and rational NNs are combined to design a human-understandable deep learning method for discovering Green's functions from data. This approach complements state-of-the-art PDE learning techniques, as a wide range of physics can be captured from the learned Green's functions such as dominant modes, symmetries, and singularity locations.

The 'forward problem' of ground-state density functional theory (DFT) constitutes finding the ground-state density n(x) that minimises a Kohn-Sham total energy functional defined using some exchange-correlation (xc) functional Exc[n]. Towards this end, the associated Euler-Lagrange equations, i.e. the Kohn-Sham equations, are often solved in practice, which demand a procedure that iterates an initial guess density to a selfconsistent density (the solution). A new framework is presented for evaluating the performance of self-consistent field methods in Kohn–Sham DFT. The aims of this work are two-fold. First, we explore the properties of Kohn–Sham DFT as it pertains to the convergence of self-consistent field iterations. Sources of inefficiencies and instabilities are identified, and methods to mitigate these difficulties are discussed. Second, we introduce a framework to assess the relative utility of algorithms, comprising a representative benchmark suite of over fifty Kohn–Sham simulation inputs, the scf-xn suite. This provides a new tool to develop, evaluate and compare new algorithms in a fair, well-defined and transparent manner. The 'inverse problem' of time-dependent (ground-state) DFT constitutes finding the time-(in)dependent Kohn-Sham potential vKS(x,t) that yields a given reference density n(x,t) upon solution of the time-(in)dependent Kohn-Sham equations. This inverse map can be unstable, particularly in the presence of low-density regions, and thus methods are designed to alleviate numerical difficulties in the present context. On the other hand, linear response time-dependent DFT centres around the first-order response of the xc potential due to perturbing densities – the so-called xc kernel fxc(x,x',w). Computing exact xc kernels represents a linearised version of the previous inverse problem: this state of affairs, whilst still challenging, is more manageable. Methods to ensure the robustness of exact numerical fxc computations are set out. In the context of inhomogenous onedimensional finite systems, these developments permit an improved understanding of fxc in itself, and in relation to various applications, such as the optical spectrum and ground-state correlation energies using the adiabatic connection fluctuation-dissipation theorem. We expect that certain key insights derived from this work will assist in the informed development of improved functional approximations.

Nonlinear matrix recovery is an emerging paradigm in which specific classes of high-rank matrices can be recovered from an underdetermined linear system of measurements. In particular, we consider matrices whose columns, seen as data points, belong to an algebraic variety, namely, a set defined by finitely many polynomial equations. The nonlinear structure is exploited using an appropriately chosen feature map, which induces a low-rank structure in the feature space, thus allowing recovery of the original (high-rank) matrix. The same framework is applicable to matrices whose columns are divided in clusters, so that we can address clustering tasks with missing entries. We formulate the rank minimization of the feature matrix as the minimization of a smooth cost function on a Riemannian manifold, which then allows us to employ the theoretically rigorous Riemannian optimization framework, that calls on differential geometry tools to construct feasible iterative algorithms on specific Riemannian manifolds. In particular, we show the potential benefits of using second-order algorithmic variants --- instead of first-order ones --- as they are able to achieve recovery of the original data with high accuracy. We additionally develop an alternating minimization method to solve the recovery problem, for which we give convergence and complexity guarantees. Considering algebraic varieties also leads us to propose a new approach for the denoising of point clouds that approximately fit such structures, as well as for their registration, namely, the task of estimating a transformation which overlaps two point clouds representing the same shape in different coordinate systems. Finally, we consider general optimization problems with smooth equality constraints. Due to the challenge of achieving feasibility for generic nonlinear equality constraints, we must depart from the Riemannian framework and we propose an infeasible algorithm which minimizes a smooth penalty function, known as Fletcher's augmented Lagrangian. We provide optimal worst-case bounds on the number of iterations that are required to find an approximate first- and second-order critical point.

Magnetic Resonance Imaging (MRI) is a non-invasive, non-ionising imaging modality with unrivalled soft tissue contrast. A key consideration for MRI is data acquisition time, which is limited by inherent technological and physiological constraints. Compressed sensing is a relatively recent framework that can reduce the MRI acquisition time by undersampling randomly and exploiting presumed redundancies in the data. The Approximate Message Passing (AMP) algorithm is an iterative compressed sensing method that efficiently reconstructs signals that have been sampled with i.i.d. sub-Gaussian sensing matrices. However, when Fourier coefficients of a signal with non-uniform spectral density are sampled, such as in MRI, AMP performs poorly in practice. In response, this thesis proposes the Variable Density Approximate Message Passing (VDAMP) algorithm for undersampled MRI data. We present three versions of VDAMP: single-coil VDAMP, where receiver coil sensitivities are ignored, Parallel-VDAMP (P-VDAMP), which includes coil sensitivities, and Denoising-P-VDAMP (D-P-VDAMP), which incorporates the statistical modelling capabilities of neural networks. Central to VDAMP is a property that we term "coloured state evolution", where the difference between the intermediate image estimate at a given iteration and the ground truth is distributed according to a zero-mean Gaussian with known covariance. We demonstrate that coloured state evolution can be leveraged to yield an algorithm that converges rapidly, and to a competitive reconstruction quality, without the need to hand-tune model parameters.

In this thesis, we study the mathematical solution of wave scattering problems which describe the behaviour of waves incident on obstacles and are highly relevant to a raft of applications in the aerospace industry. The techniques considered in the present work can be broadly classed into two categories: analytically based methods which use special transforms and functions to provide a near-complete mathematical description of the scattering process, and numerical techniques which select an approximate solution from a general finite-dimensional space of possible candidates. The first part of this thesis addresses an analytical approach to the scattering of acoustic and vortical waves on an infinite periodic arrangement of finite-length flat blades in parallel mean flow. This geometry serves as an unwrapped model of the fan components in turbo-machinery. Our contributions include a novel semi-analytical solution based on the Wiener-Hopf technique that extends previous work by lifting the restriction that adjacent blades overlap, and a comprehensive study of the composition of the outgoing energy flux for acoustic wave scattering on this array of blades. These results provide an insight into the importance of energy conversion between the unsteady vorticity shed from the trailing edges of the cascade blades and the acoustic field. Furthermore, we show that the balance of incoming and outgoing energy fluxes of the unsteady field provides a convenient tool for understanding several interesting scattering symmetries on this geometry. In the second part of the thesis, we focus on numerical techniques based on the boundary integral method which allows us to write the governing equations for zero mean flow in the form of Fredholm integral equations. We study the solution of these integral equations using collocation methods for two-dimensional scatterers with smooth and Lipschitz boundaries. Our contributions are as follows: Firstly, we explore the extent to which least-squares oversampling can improve collocation. We provide rigorous analysis that proves guaranteed convergence for small amounts of oversampling and shows that superlinear oversampling can ensure faster asymptotic convergence rates of the method. Secondly, we examine the computation of the entries in the discrete linear system representing the continuous integral equation in collocation methods for hybrid numerical-asymptotic basis spaces on simple geometric shapes in the context of high-frequency wave scattering. This requires the computation of singular highly-oscillatory integrals and we develop efficient numerical methods that can compute these integrals at frequency-independent cost. Finally, we provide a general result that allows the construction of recurrences for the efficient computation of quadrature moments in a broad class of Filon quadrature methods, and we show how this framework can also be used to accelerate certain Levin quadrature methods.

In order to meet new stringent emissions and noise reduction targets, the turbomachinery community is increasingly concerned with addressing inefficiencies in aeroengines, and of particular interest is the issue of sound generated through aeroacoustic processes within the components of the turbomachine. At the same time, precipitous increases in computing power over the last thirty years are driving larger and more detailed simulations. Increasing data output from these requires the development of new techniques that are able to identify the underlying processes and facilitate future efforts for optimisation. This thesis seeks to address this by creating a direct-adjoint framework that can be applied in two sub-categories of flow problems commonly found in turbomachinery aeroacoustics, and lays the groundwork for future efforts in the efficient analysis and optimisation of these flows. For this, the in-house, high-order, time-domain pseudo-wave-based aeroacoustics code, previously used for isolated airfoil self-noise studies, has been extensively modified and adapted for application in turbomachinery geometries. Using this code as a starting point, self-excited instabilities and the resulting acoustic fields are simulated within a representative compressor blade row, at off-design conditions. Prior experiments and numerical calculations of the flow at these conditions have confirmed the presence of an unsteady reattaching laminar separation bubble on the suction side of the blade, shedding vortices that are subsequently scattered by the blunt trailing edge. Nonlinear flow analyses are carried out first within a single blade passage, followed by a linear impulse response analysis. A detailed linear direct-adjoint mean flow global stability analysis is performed, and the direct spectrum obtained and correlated with the nonlinear pressure signals, showing excellent agreement. The direct global modes are categorised and corresponding adjoint modes leveraged to compute structural sensitivities with wavemaker regions. The sensitivity analysis confirms the presence of aeroacoustic feedback loops between the pressure and suction sides of the blade, and is consistent with earlier observations in similar studies of isolated airfoils. The theory of block-circulant matrices is then detailed and exploited to achieve computationally efficient decompositions of the governing system matrices, and the previous global stability analysis reproduced with ten passages per period window. An increase in linear growth rates is observed, new low-frequency structures identified and analysed in the context of new feedback loop pathways which propagate disturbances throughout the domain. The effects of aerodynamic mistuning are then briefly addressed. Rotor-stator interaction is then considered in the context of an adjoint sensitivity problem. For this, a new high-order, time-domain adjoint sliding-plane interface treatment, which is capable of modelling and propagating self-excited instabilities and tonal blade row interaction noise, is developed and extensively validated. This work concludes with the application of this theory to a simplified rotor-stator sensitivity analysis within the context of a time-domain nonlinear-adjoint looping algorithm. Sensitivities identified here are briefly discussed. Applicability of these methods to different codes and numerical treatments is highlighted throughout this thesis.

Parallel-in-time (PinT) methods have become an increasingly popular tool in the numerical solution of time-dependent Partial Differential Equations (PDEs), in light of their potential for extracting additional concurrency when spatial parallelism saturates. Particularly complex applications, such as the ones involving plasma simulations for fusion research, would benefit greatly from the implementation of such methods, given the large computational time usually necessary for the numerical solution of their governing equations. In the first part of this thesis, we develop the concept of Space-Time Block Preconditioning (STBP), a technique for simplifying the task of time-parallelising the solution of a whole system of PDEs by reducing it to that of inverting space-time operators acting on single variables only, which are more amenable to be tackled by PinT algorithms. This technique borrows from block preconditioning strategies developed for spatial operators in the framework of time-stepping algorithms, and extends them to the whole space-time setting. After illustrating the guiding principles behind STBP, we measure its effectiveness with applications in fluid dynamics and magnetohydrodynamics. In both cases, the resulting pre-conditioners show excellent scaling properties with respect to mesh refinement, and their applications come at little-to-no overhead cost in terms of iterations to convergence with respect to their time-stepping counterpart, thus providing evidence for their potential parallel efficiency. Ultimately, however, the parallel performance of a space-time block preconditioner depends on that of its internal components, which still rely on already-existing PinT methods. This might become an issue when dealing with advection-dominated or hyperbolic equations, as classical PinT algorithms have shown to struggle in this case. For this reason, the second part of this thesis is dedicated to investigating the causes behind this degradation in performance, as well as to study more recent PinT approaches for overcoming this limitation.

Topology optimization finds the optimal material distribution of a continuum in a domain, subject to PDE and volume constraints. Density-based models often result in a PDE, volume and inequality constrained, nonconvex, infinite-dimensional optimization problem. These problems can exhibit many local minima. In practice, heuristics are used to aid the search for better minima, but these can fail even in the simplest of cases. In this thesis we address two core issues related to the nonconvexity of topology optimization problems: the convergence of the discretization and the computation of the solutions. First, we consider the convergence of a finite element discretization of a fluid topology optimization problem. Results available in literature show that there exists a sequence of finite element solutions that weakly(-*) converges to a solution of the infinite-dimensional problem. We improve on these classical results. In particular, by fixing any isolated minimizer, we show that there exists a sequence of finite element solutions that \emph{strongly} converges to that minimizer. Moreover, these results hold for both traditional conforming finite element methods and more sophisticated divergence-free discontinuous Galerkin finite element methods. We then focus on developing a solver that can systematically compute multiple minimizers of a general density-based topology optimization problem. This leads to the successful computation of 42 distinct solutions of a two-dimensional fluid topology optimization problem. Finally, by developing preconditioners for the linear systems that arise during the optimization process, we are able to apply the solver to three-dimensional fluid topology optimization problems. This culminates in an example where we compute 11 distinct three-dimensional solutions.

In this thesis we focus on the theoretical and computational aspects of nonlinear eigenvalue problems (NEPs), which arise in several fields of computational science and engineering, such as fluid dynamics, optics, and structural engineering. In the last twenty years several researchers devoted their time in studying efficient and precise ways to solve NEPs, which cemented their importance in numerical linear algebra. The most successful algorithms developed towards this goal are either based on contour integrals, or on rational approximations and linearizations. The first part of the thesis is dedicated to contour integral algorithms. In this framework, one computes specific integrals of a holomorphic function G(z) over the contour of a target region and exploits results of complex analysis to retrieve the eigenvalues of G(z) inside it. Our main contribution consists in having expanded the theory to include meromorphic functions, i.e., functions with poles inside the target region. We showed that under some circumstances, these algorithms can mistake a pole for an eigenvalue, but these situations are easily recognised and the main results from the holomorphic case can be extended. Furthermore, we proposed several heuristics to automatically choose the parameters that are needed to precisely retrieve the eigenpairs. In the second part of the thesis, we focus on rational approximations. Our goal was developing algorithms that construct robust, i.e., reliable for a given tolerance and scaling independent, rational approximations for functions given in split form or in black-box form. In the first case, we proposed a variant of the set-valued AAA, named weighted AAA, which guarantees to return an approximation with the chosen accuracy. In the second one, we built a two-phase approach, where an initial step performed by the surrogate AAA is followed by a cyclic Leja-Bagby refinement. We concluded the section with numerous numerical experiments based on the NLEVP library, comparing contour integral and rational approximation algorithms. The third and final part of the thesis is about tropical linear algebra. Our research on this topic started has a way to set the parameters of the aforementioned contour integral algorithms: in order to do that, we extended the theory of tropical roots from tropical polynomials to tropical Laurent series. Unlike in the polynomial case, a tropical Laurent series may have infinite tropical roots, but they are still in bijection with the slopes of the associated Newton polygon and they still provide annuli of exclusion for the eigenvalues of the Laurent series.

Finite Markov chains are probabilistic network models that are commonly used as representations of dynamical processes in the physical sciences, biological sciences, economics, and elsewhere. Markov chains that appear in realistic modelling tasks are frequently observed to be nearly reducible, incorporating a mixture of fast and slow processes that leads to ill-conditioning of the underlying matrix of probabilities for transitions between states. Hence, the wealth of established theoretical results that makes Markov chains attractive and convenient models often cannot be used straightforwardly in practice, owing to numerical instability associated with the standard computational procedures to evaluate the expressions. This work is concerned with the development of theory, algorithms, and simulation methods for the efficient and numerically stable analysis of finite Markov chains, with a primary focus on exact approaches that are robust and therefore applicable to nearly reducible networks. New methodologies are presented to determine representative paths, identify the dominant transition mechanisms for a particular process of interest, and analyze the local states that have a strong influence on the characteristics of the global dynamics. The novel approaches yield new insights into the behaviour of Markovian networks, addressing and overcoming numerical challenges. The methodology is applied to example models that are relevant to current problems in chemical physics, including Markov chains representing a protein folding transition, and a configurational transition in an atomic cluster. Relevant classical theory of finite Markov chains and a description of existing robust algorithms for their numerical analysis is given in Chapter 1. The remainder of this thesis considers the problem of investigating a transition from an initial set of states in a Markovian network to an absorbing (target) macrostate. A formal approach to determine a finite set of representative transition paths is proposed in Chapter 2, based on exact pathwise decomposition of the total productive flux. This analysis allows for the importance of competing dynamical processes to be rigorously quantified. A robust state reduction algorithm to compute the expectation of any path property for a transition between two endpoint states is also described in Chapter 2. Chapter 3 reports further numerically stable state reduction algorithms to compute quantities that characterize the features of a transition at a statewise level of detail, allowing for identification of the local states that play a key role in modulating the slow dynamics. An expression is derived for the probability that a state is visited on a path that proceeds directly to the absorbing state without revisiting the initial state, which characterizes the dynamical relevance of an individual state to the overall transition process. In Chapter 4, an unsupervised strategy is proposed to utilize a highly efficient simulation algorithm for sampling paths on a Markov chain. The framework employs a scalable community detection algorithm to obtain an initial clustering of the network into metastable sets of states, which is subsequently refined by a variational optimization procedure. The optimized clustering is then used as the basis for simulating trajectory segments that necessarily escape from the metastable macrostates. The thesis is concluded with an overview of recent related advances that are beyond the scope of the current work (Chapter 5), and a discussion of potential applications where the novel methodology reported herein may be applied to perform insightful analyses that were previously intractable.

Scattering problems are ubiquitous in scientific and engineering applications. These problems can often be straightforwardly modeled using relatively simple linear partial differential equations, such as the Helmholtz equation. The numerical solution of these equations, however, can be challenging due to several features of the problem, such as ill-conditioning, unbounded domains, corner singularities, and oscillatory solutions. In this thesis we present novel numerical methods for the efficient handling of some of these challenges. We first consider the case of homogeneous scattering on domains with corners and introduce a new solver in this regime, based on results in the rational approximation theory literature. The solver uses an exponentially clustered sum of dipoles to resolve the singularities in the solution and a simple expansion to resolve the smooth part. We show that the convergence is root-exponential with respect to the number of degrees of freedom and illustrate the behavior of the solver in several numerical experiments. We next turn our attention to the case of inhomogeneous scattering and present a new solver for the Lippmann-Schwinger equation, an integral equation reformulation of the problem. The new solver is based on the hierarchically block separable matrix format and exploits the geometry of the discretization and problem for accelerated inversion. The solver is shown to be effective both as a direct solver and preconditioner in a number of numerical experiments. For our last contribution we consider broadband applications in homogeneous scattering. We present a new technique for accelerating existing direct solvers in this regime. The technique works by computing basis matrices for all of the low-rank blocks in a rank-structured matrix format and leverages the observation that the resulting ranks are often comparable to those at just the highest frequency. Again, the viability of the method is illustrated in several numerical experiments.

Expressing a matrix as the sum of a low-rank matrix plus a sparse matrix is a flexible model capturing global and local features in data, and is the foundation of robust principal component analysis (Candès et al., 2011). This thesis is concerned with low-rank plus sparse matrix sensing - the problem of recovering a matrix that is formed as the sum of a low-rank and a sparse matrix, and the two components, from a number of measurements far smaller than the dimensionality of the matrix. It is well-known, that inverse problems over low-rank matrices, such as robust principal component analysis and matrix completion, require a low coherence between the low-rank matrix and the canonical basis. However, in this thesis, we demonstrate that the well-posedness issue is even more fundamental; in some cases, both robust principal component analysis and matrix completion can fail to have any solutions due to the fact that the set of low-rank plus sparse matrices is not closed. As a consequence, the lower restricted isometry constants (RICs) cannot be upper bounded for some low-rank plus sparse matrices unless further restrictions are imposed on the constituents. We close the set of low-rank plus sparse matrices by posing an additional bound on the Frobenius norm of the low-rank component, and ensure the optimisation is well-posed and that the RICs can be bounded. We show that constraining the incoherence of the low-rank component also closes the set provided μ < √mn/ (r √s) and satisfies a certain additivity property necessary for the analysis of recovery algorithms. Compressed sensing, matrix completion, and their variants have established that data satisfying low complexity models can be efficiently measured and recovered from a number of measurements proportional to the model complexity rather than the ambient dimension (Foucart and Rauhut, 2013). This thesis develops similar guarantees showing that m × n matrices that can be expressed as the sum of a rank-r matrix and a s-sparse matrix can be recovered by computationally tractable methods.

We present a modified Euler-Maruyama scheme using approximate random variables, produced by the inverse transform method, using cheap approximations to the inverse Gaussian cumulative distribution function. We analyse the error for two approximations: a piecewise constant approximation on equally spaced intervals, and a piecewise linear approximation using geometric intervals dense at the singularities. High speed implementations faster than Intel's MKL are provided, suitable for modern vector hardware. The error between the approximations from the exact and modified Euler-Maruyama schemes is bounded by the error from the approximate random variables. We incorporate this scheme into a multilevel Monte Carlo framework producing a nested scheme, and show that the discretisation error couples to the random variables' approximation error. The result directly extends to Lipschitz and differentiable payoff functions. For Lipschitz and non-differentiable payoff functions simulated using a time step delta, there is a transition from a variance decay of O(delta) to order(delta1/2) as the discretisation error becomes dominant. These variance bounds are demonstrated numerically for geometric Brownian motion and a variety of payoff functions of varying smoothness. For approximate random variables computed in low-precision, a model for the accumulated rounding error is developed and assessed. Half-precision is viable for a range of coarse path simulations, and can be extended further by incorporating a Kahan compensated summation. We empirically demonstrate these ideas are transferable to the Milstein scheme, and the more difficult Cox-Ingersoll-Ross process and its non-central chi-squared distribution. We estimate that under the Black-Scholes model, options can be priced using path simulations with approximate Gaussian random variables, obtaining a five-times or more speed improvement without losing accuracy.

In this thesis, we consider problems related to the numerical simulation of stochastic differential equations (SDEs). In particular, we are interested in methods that use both increments and areas of the Brownian path. For example, time integrals of Brownian motion carry useful information whilst being normally distributed and thus straightforward to generate. We will make precise the "path information" contained by such integrals and present a polynomial approximation theorem for Brownian motion. Since the Gaussian coefficients given by this expansion are independent, we can express Brownian motion as a polynomial with additional noise. We shall then use this simple observation to develop high order methods. The majority of these methods approximate SDEs through the use of ODEs, which can be accurately discretized using state-of-the-art solvers. To numerically demonstrate these approaches, we will simulate four SDEs: Inhomogeneous Geometric Brownian Motion, Cox-Ingersoll-Ross model, underdamped Langevin equation and the Lévy area of Brownian motion. The code for these examples can be found at github.com/james-m-foster. We also study problems related to variable step size methods for SDEs. Most notably, we investigate which numerical methods and variable step size strategies can ensure pathwise convergence to the true SDE solution. Our analysis is based on rough path theory and requires methods to be locally close to an ODE driven by a piecewise "Brownian polynomial". This result can therefore be applied to the Heun and midpoint methods. In order to prove convergence, we will require variable step size strategies to produce a nested sequence of partitions with mesh size tending to zero.

We consider the numerical approximation of incompressible non-Newtonian flow by means of the finite element method, where the constitutive law is defined through an implicit relation G(S ,D(u)) = 0. The setting considered in this work captures models widely used in applications, such as the Bingham and the Carreau–Yasuda constitutive relations. Since in general it is not possible to solve for the shear stress S in the constitutive relation, the emphasis is placed on formulations treating the shear stress as a variable. Under the assumption that the constitutive relation defines a monotone graph with r-growth, and that the finite element spaces satisfy appropriate inf-sup stability conditions, the first part of the thesis extends earlier results in the literature to provide a convergence result that guarantees that a subsequence of the numerical approximations converges weakly to a solution of the system, in the optimal range r > 2d/(d+2), where d is the spatial dimension. The qualitative nature of this convergence result is a consequence of the generality of the framework of implicitly constituted fluids, for which e.g. higher regularity estimates are not available. Computational examples show, nevertheless, that the numerical scheme considered exhibits the expected convergence rates, in the situations where these are available. In the second part of the thesis we develop an augmented Lagrangian preconditioner for a stress-velocity-pressure formulation of the steady system. The preconditioner involves a specialised multigrid algorithm that makes use of a space decomposition that captures the kernel of the divergence operator, and non-standard intergrid transfer operators. Although the current theory for robust multigrid works only for symmetric and positive-definite systems (and thus does not apply to the systems considered in this thesis), the resulting preconditioner exhibits remarkable robustness properties. In the final chapter of the thesis, the extension to the anisothermal case is carried out. We employ an implicit constitutive relation that allows for a temperature dependence of rheological parameters such as the viscosity and the yield stress, and provide convergence results for the unsteady forced convection system that takes the viscous dissipation term S : D (u) into account, and for the steady Oberbeck–Boussinesq approximation. For the latter an augmented Lagrangian preconditioner is also introduced; this preconditioner exhibits robust convergence behaviour when applied to the Navier–Stokes and power-law systems, including temperature-dependent viscosity, heat conductivity, and viscous dissipation.

Wrinkling is a universal instability occurring in a wide variety of engineering and biological materials. Despite extensive study across different systems, a full description is still lacking, particularly in the case of growing materials. Most studies consist of some combination of a substantially simplified elastic model, an analysis that does not go beyond first order, and numerical simulations in commercial software packages that do not necessarily correspond to the mathematical system studied. This thesis addresses all three of these shortcomings by providing a systematic analysis of a fully hyperelastic bilayer past the linear stability threshold into the weakly nonlinear regime, along with a carefully discretised numerical bifurcation analysis of the system. For comparison, we assume that wrinkling is generated either by the isotropic growth of the film or by the lateral compression of the entire film-substrate system (both in plane strain). We adopt a stream-function-based formulation and perform an exhaustive linear analysis of the wrinkling problem for all stiffness ratios and under a variety of additional boundary and material effects. Namely, we consider the effect of added pressure, surface tension, an upper substrate, and fibres. We obtain analytical estimates of the instability in the two asymptotic regimes of long and short wavelengths. We then carry out a weakly-nonlinear analysis to derive an amplitude equation that describes the evolution of the wrinkling amplitude beyond the bifurcation point, followed by a comprehensive numerical bifurcation analysis of the problem using the finite element method. We demonstrate excellent agreement between the weakly-nonlinear analysis and the numerical experiments and are also able to directly solve for the bifurcation point in our discretised system and characterise the effect of implementation details such as the aspect ratio of the computational domain on the observed bifurcation point. We then explore solutions of the amplitude equation in the case that the wrinkling amplitude is allowed to vary over long spatial and/or temporal scales. Finally, we demonstrate that our numerical experiments are able to identify secondary bifurcations in the system similar to those observed in experiments for which analytical methods have thus far been unable to provide a complete explanation.

This thesis consists of two parts. In the first, we develop two new strategies for spatial white noise and Gaussian-Mat\'ern field sampling that work within a non-nested multilevel (quasi) Monte Carlo (ML(Q)MC) hierarchy. In the second, we apply the techniques developed to quantify the level of uncertainty in a new stochastic model for tracer transport in the brain. The new sampling techniques are based on the stochastic partial differential equation (SPDE) approach, which recasts the sampling problem as the solution of an elliptic equation driven by spatial white noise. We present a new proof of an a priori error estimate for the finite element (FEM) solution of the white noise SPDE. The proof does not require the approximation of white noise in practice, and includes higher order elliptic operators and p-refinement. Within the SPDE approach, the efficient sampling of white noise realisations can be computationally expensive. In this thesis, we present two new sampling techniques that can be used to efficiently compute white noise samples in a FEM-MLMC and FEM-MLQMC setting. The key idea is to exploit the finite element matrix assembly procedure and factorise each local mass matrix independently, hence avoiding the factorisation of a large matrix. In a multilevel framework, the white noise samples must be coupled between subsequent levels. We show how our technique can be used to enforce this coupling even in the case of non-nested mesh hierarchies. In the MLQMC case, the QMC integrand variables must also be ordered in order of decaying importance to achieve fast convergence with respect to the number of samples. We express white noise as a Haar wavelet series whose hierarchical structure naturally exposes the leading order dimensions. We split this series in two terms which we sample via a hybrid standard Monte Carlo/QMC approach. We demonstrate the efficacy of our sampling methods with numerical experiments. In a multilevel setting, a good coupling is enforced and the telescoping sum is respected. In the MLQMC case, the asymptotic convergence rate is the same as standard Monte Carlo, but significant computational gains are obtained in practice thanks to a pre-asymptotic QMC-like regime. In the final part of the thesis, we employ a combination of the methods presented to solve a PDE with random coefficients describing tracer transport within the interstitial fluid of the brain. Numerical simulations support the claim that diffusion alone cannot explain the penetration of tracers within deep brain regions as observed in clinical experiments, even when uncertainties in the diffusivity have been accounted for. A convective velocity field may however increase tracer transport, provided that a directional structure is present in the interstitial fluid circulation.

High-order methods have become of increasing interest in recent years in computational physics. This is in part due to their perceived ability to, in some cases, reduce the computational overhead of complex problems through both an efficient use of computational resources and a reduction in the required degrees of freedom. One such high-order method in particular – Flux Reconstruction – is the focus of this thesis. This body of work relies and expands on the theoretical methods that are used to understand the behaviour of numerical methods – particularly related to their real-world application to industrial problems. The thesis begins by challenging some of the existing dogma surrounding computational fluid dynamics by evaluating the performance of high-order flux reconstruction. First, the use of the primitive variables as an intermediary step in the construction of flux terms is investigated. It is found that reducing the order of the flux function by using the conserved rather than primitive variables has a substantial impact on the resolution of the method. Critically, this is supported by a theoretical analysis, which shows that this mechanism of error generation becomes increasing important to consider as the order of accuracy increases. Next, the analysis of Flux Reconstruction was extended by analytically and numerically exploring the impact of higher dimensionality and grid deformation. It is found that both expanding and contracting grids – essential components of real-world domain decomposition – can cause dispersion overshoot in two dimensions, but that FR appears to suffer less that comparable Finite Difference approaches. Fully discrete analysis is then used to show that, depending on the correction function, small perturbations in incidence angle can cause large changes in group velocity. The same analysis is also used to theoretically demonstrate that Discontinuous Galerkin suffers less from dispersion error than Huynh’s FR scheme – a phenomenon that has previously been observed experimentally, but not explained theoretically. This thesis concludes with the presentation of a robust theoretical underpinning for determining stable correction functions for FR. Three new families of correction functions are presented, and their properties extensively explored. An important theoretical finding is introduced – that stable correction functions are not defined uniquely be a norm. As a result, a generalised approach is presented, which is able to recover all previously defined correction functions, but in some instances via a different norm to their original derivation. This new super-family of correction functions shows considerable promise in increasing temporal stability limits, reducing dispersion when fully discretised, and increasing the rate of convergence. Taken altogether, this thesis represents a considerable advance in the theoretical characterisation and understanding of a numerical method – one which, it has been shown, has enormous potential for forming the heart of future computational physics codes.

Multiphase flow through porous media can be modelled as a complex system of partial differential equations. Such models can be used to optimize the recovery of oil and gas from subsurface reservoirs. In the case of highly viscous oils, thermal recovery techniques are typically used to enhance their extraction. To simulate this, models describing the flow of fluids (typically oil, water, and gas) are coupled with a model for heat flow. Thermal reservoir simulation entails solving these highly coupled systems. Their complexity and the computational effort needed to solve them motivate the need for highly efficient solvers. In reservoir simulation, most of the computational time is spent on solving linearized systems with a preconditioned Krylov subspace iterative method. Industry-standard preconditioning techniques are based on the approach introduced by Wallis in 1983, the Constrained Pressure Residual method (CPR). This preconditioner is a two-stage process involving the solution of a restricted pressure system. While initially designed for isothermal reservoir simulation, CPR is also the standard for thermal cases. However, its treatment of the conservation of energy equation does not incorporate heat diffusion, which is often dominant in thermal cases. We are interested in preconditioners specifically designed for thermal reservoir simulation. In this thesis, we present an extension of CPR: the Constrained Pressure-Temperature Residual (CPTR) method, where a restricted pressure-temperature system is solved in the first stage. To study the effects of both pressure and temperature on fluid and heat flow, we first consider a model of non-isothermal single-phase flow through porous media. For this model, we develop a block preconditioner with an efficient Schur complement approximation. Then, we extend this method for multiphase flow as a solver for the first stage of CPTR. We present a comparison of the algorithmic performance of the different preconditioning approaches under mesh refinement and parallelization.

When designing a new car or a plane, engineers need to solve the Navier-Stokes equations to understand how air flows around the object. Based on experience and intuition, they modify the design slightly, then solve the equations again, and inspect the changes. This process is iterated many times until a final design that minimises or maximises some quantity of interest, such as drag or lift, is found. The goal of shape optimisation is to automate this type of process. In this thesis we address several issues related to shape optimisation. Focussing on the case when the shape is discretised using a mesh and when PDE constraints are solved using the finite element method, we describe a reformulation of the shape derivative as the derivative of the pushforward from the reference element. This viewpoint allows for automated calculation of shape derivatives in finite element software. When shape optimisation is performed by deforming an initial mesh, the choice of deformation is important. We propose a new Hilbert space structure on the space of deformations that results in high mesh quality of the deformed domains. We then focus on the solution of a particular PDE constraint given by the steady, incompressible Navier-Stokes equations that govern laminar flow. The solution of these equations becomes challenging for large Reynolds number. We develop augmented Lagrangian based preconditioners that exhibit robust performance as the Reynolds number is increased. The effectiveness and scalability of the developed solvers is demonstrated for a range of test problems.

Highly supersonic, magnetized and differentially rotating, accretion flows provide an environment for the propagation of waves and the growth of instabilities unlike those normally encountered on Earth. In this dissertation, I investigate the properties of oscillations and instabilities in accretion discs, developing theoretical and numerical models to explore the physical nature of variability encountered in a variety of astrophysical contexts. Through linear theory and semi-analytical calculations, I first consider the physical properties of magnetically altered inertial waves. 'Trapped inertial waves' provide an attractive explanation of the fast variability observed in the emission from low-mass black hole binary systems, but such oscillations can be affected by magnetic tension provided by large-scale poloidal magnetic fields threading the accretion disc. Through local and global analyses, I constrain the modification of trapped inertial waves by poloidal, toroidal and helical magnetic fields. I then investigate the excitation of oscillations in deformed discs with eccentric, non-circular streamlines. Many processes can lead to the growth of eccentricity in accretion discs, and turbulence deriving from the excitation of inertial waves by a local parametric instability provides one mechanism for curbing this growth. However, eccentric discs also provide an environment for the excitation of additional, inherently global oscillations, and I present a framework facilitating a semi-analytical investigation of these modes. I finally employ numerical simulations to explore the dynamics of accretion disc oscillations in the non-linear regime. I first follow the non-linear saturation of inertial waves driven by parametric resonance in non-relativistic discs, and confirm the growth of a second family of large-scale, low-frequency oscillations. Using a pseudo-Newtonian framework to approximate relativistic effects, I then demonstrate the excitation of trapped inertial waves through non-linear coupling with accretion disc deformations in a black hole accretion disc, providing preliminary evidence that trapped inertial waves can be excited even in the presence of MHD turbulence.

Minimising a nonlinear function is a ubiquitous problem in applications, and standard algorithms need accurate evaluations of the function and its derivatives. However, if the function is black-box, expensive to evaluate, or noisy, it may be impractical or even impossible to obtain accurate derivatives, and we require derivative-free optimisation (DFO). Model-based DFO methods perform well in practice, taking many features from derivative-based methods, but can have lower performance for noisy problems and a high linear algebra cost. In this thesis, we aim to improve the flexibility, robustness and scalability of state-of-the-art methods for model-based DFO by designing, analysing, implementing and comprehensively testing three new algorithms for nonlinear optimisation, with a special focus on nonlinear least-squares problems (such as parameter fitting). The first, DFO-LS, is designed for nonlinear least-squares problems, and is simpler than existing methods, constructing linear residual models with a flexible approach. DFO-LS can benefit from a reduced initialisation cost, and has improved scalability and runtime over existing methods without a performance penalty. DFO-LS also has features for noisy problems, including a novel multiple restarts mechanism, which are low cost in evaluations and linear algebra, giving DFO-LS better performance compared to existing techniques for handling noise. We then introduce Py-BOBYQA, an extension of BOBYQA (Powell, 2009) with a similar multiple restarts mechanism to DFO-LS, amongst other new enhancements. It matches or outperforms BOBYQA and other state-of-the-art solvers on both smooth and noisy problems. We also propose a simple extension of Py-BOBYQA that may enable it to escape local minima and progress towards the global optima. Lastly, we introduce, for large-scale nonlinear least-squares problems, DFBGN, the first model-based DFO method to optimise in low-dimensional subspaces at each iteration. DFBGN has lower linear algebra costs, improved runtime, and hence improved performance on large-scale problems than DFO-LS, as it can perform many more iterations within a reasonable runtime limit. It can also make progress under very small evaluation budgets, with a low runtime cost.

The subject of this thesis is a class of fast inpainting methods (image or video) based on the idea of filling the inpainting domain in successive shells from its boundary inwards. Image pixels (or video voxels) are filled by assigning them a color equal to a weighted average of either their already filled neighbors (the ``direct'' form of the method) or those neighbors plus additional neighbors within the current shell (the ``semi-implicit'' form). In the direct form, pixels (voxels) in the current shell may be filled independently, but in the semi-implicit form they are filled simultaneously by solving a linear system. We focus in this thesis mainly on the image inpainting case, where the literature contains several methods corresponding to the {\em direct} form of the method - the semi-implicit form is introduced for the first time here. These methods effectively differ only in the order in which pixels (voxels) are filled, the weights used for averaging, and the neighborhood that is averaged over. All of them are very fast, but at the same time all of them leave undesirable artifacts such as ``kinking'' (bending) or blurring of extrapolated isophotes. This thesis has two main goals. First, we introduce new algorithms within this class, which are aimed at reducing or eliminating these artifacts, and also target a specific application - the 3D conversion of images and film. The first part of this thesis will be concerned with introducing 3D conversion as well as Guidefill, a method in the above class adapted to the inpainting problems arising in 3D conversion. However, the second and more significant goal of this thesis is to study these algorithms as a class. In particular, we develop a mathematical theory aimed at understanding the origins of artifacts mentioned. Through this, we seek is to understand which artifacts can be eliminated (and how), and which artifacts are inevitable (and why). Most of the thesis is occupied with this second goal. Our theory is based on two separate limits - the first is a {\em continuum} limit, in which the pixel width →0, and in which the algorithm converges to a partial differential equation. The second is an asymptotic limit in which h is very small but non-zero. This latter limit, which is based on a connection to random walks, relates the inpainted solution to a type of discrete convolution. The former is useful for studying kinking artifacts, while the latter is useful for studying blur. Although all the theoretical work has been done in the context of image inpainting, experimental evidence is presented suggesting a simple generalization to video. Finally, in the last part of the thesis we explore shell-based video inpainting. In particular, we introduce spacetime transport, which is a natural generalization of the ideas of Guidefill and its predecessor, coherence transport, to three dimensions (two spatial dimensions plus one time dimension). Spacetime transport is shown to have much in common with shell-based image inpainting methods. In particular, kinking and blur artifacts persist, and the former of these may be alleviated in exactly the same way as in two dimensions. At the same time, spacetime transport is shown to be related to optical flow based video inpainting. In particular, a connection is derived between spacetime transport and a generalized Lucas-Kanade optical flow that does not distinguish between time and space.

This thesis is about computational theory in the setting of the mathematics of information. The first goal is to demonstrate that many commonly considered problems in optimisation theory cannot be solved with an algorithm if the input data is only known up to an arbitrarily small error (modelling the fact that most real numbers are not expressible to infinite precision with a floating point based computational device). This includes computing the minimisers to basis pursuit, linear programming, lasso and image deblurring as well as finding an optimal neural network given training data. These results are somewhat paradoxical given the success that existing algorithms exhibit when tackling these problems with real world datasets and a substantial portion of this thesis is dedicated to explaining the apparent disparity, particularly in the context of compressed sensing. To do so requires the introduction of a variety of new concepts, including that of a breakdown epsilon, which may have broader applicability to computational problems outside of the ones central to this thesis. We conclude with a discussion on future research directions opened up by this work.

Existing reinforced concrete bridges may be deemed inadequate to carry the ever-increasing traffic loads according to the current codes and standards before they reach the end of their design life. It may therefore be required to either strengthen or replace these structures, which can be costly and causes disruptions to the infrastructure. This work investigates experimentally the possibilities to extend the useful life of existing reinforced concrete slab-on-beam structures deficient in shear by means of structural strengthening with fibre-reinforced polymers (FRP). The experimental campaign involved mechanical testing of ten full-scale T-beam specimens, representative of typical existing slab-on-beam bridges. Two sizes of test specimen were used to investigate the effect of size on the ultimate shear capacity of the beams. The investigated shear-strengthening configurations included externally bonded carbon fibre reinforced polymer (CFRP) sheets in a U-wrap configuration with and without end-anchorage and deep embedded CFRP bars. Unstrengthened control specimens were also tested to provide baseline for comparison. The results from the experimental programme revealed that while the deep embedment strengthening solution provides an increase in shear capacity of up to 50%, the strengthening with CFRP U-wraps results in reduced capacity compared with the underlying control beam. This presents a major implication in terms of safe design predictions of shear capacity of reinforced concrete T-beams strengthened with CFRP sheets as this is the most commonly used shear-strengthening scheme in practice. The study also demonstrated that greater contribution from the externally bonded CFRP U-wraps can be achieved using end-anchorage systems, which delay the debonding of the CFRP. The applicability of current codes of standards and guidelines was studied as well as appropriateness of using advanced numerical methods for assessment of existing reinforced concrete structures. It was found that while the standards used for assessment greatly under-predict the shear capacity, the guidelines for FRP-strengthened beams either under- or over-predict the shear capacity of the tested beams. More accurate predictions are possible using advanced fracture mechanics-based methods for both the unstrengthened as well as the strengthened beams.

In this thesis, the experimental programmes and numerical investigations are described that have been conducted to partially cover the knowledge gap in the field of strengthening two-way reinforced concrete (RC) flat slabs. The conducted studies demonstrate that the most common method to strengthen two-way RC slabs is by applying fibre reinforced polymers (FRP) on the tension surface of the slabs. Applying prestressed FRP to strengthen two-way flat slabs combines the advantages of both FRP strengthening and prestressing to enhance the efficiency of the strengthening methods. Hence, two previous studies on strengthening two-way flat slabs with non-prestressed and prestressed FRP are analysed to clarify the effect of different strengthening methods on the behaviour of slabs. Both studies demonstrate the benefits of applying non-prestressed FRP to enhance the structures' capacities. However, for the case of strengthening RC slabs with prestressed FRP, the results seem to be controversial and more studies are necessary to arrive at a conclusion on whether it is feasible to strengthen RC slabs with prestressed FRP. Further analysis indicates that there is an optimum percentage of prestressing for the FRPs applied to RC slabs. Increasing the prestressing ratio of FRP to the optimal percentage increases the ultimate load capacity of the RC slabs. However, increasing the prestressing ratio of FRPs beyond the optimum value can cause de-bonding and loss of composite action, which prevents the RC slab from reaching its expected ultimate load capacity. The optimum prestressing ratio of FRP depends on the prestress load as well as the concrete tensile strength and slab depth. Eventually, a formula is proposed to estimate the optimum FRP-prestress ratio considering the effective parameters in both concrete and FRP. Moreover, this thesis elaborates on an investigation that was conducted to make a comparison between the effects of different strengthening methods such as FRP strengthening, applying vertical (shear) reinforcement, and their combination, on the behaviour of flat slabs with different conditions (tensile reinforcement ratios). To conduct the investigation, eight slab specimens were cast, which were classified into two categories: low and high tensile reinforcement ratios. The strengthening methods included applying FRP sheets to the tension surface of the RC slabs externally, applying vertical (shear) reinforcement, and a combination of both methods. The experimental and validated numerical results demonstrate that the most efficient strengthening strategy is a combination of strengthening methods in both categories. Strengthening with FRP sheets improves the slabs load capacity in both categories. However, applying vertical (shear) reinforcement does not significantly affect the behaviour of RC slabs with a low tensile reinforcement ratio. From the resulting analyses, it was concluded that the strut and tie model of the FRP strengthened structure changes compared with the control specimen. This enables researchers and designers to justify how FRP strengthening enhances the punching strength of the slab, an aspect that has not been explained in previous studies. The results also show that applying vertical (shear) reinforcement in the critical punching area strengthens the critical compressive strut of the RC slab. This shifts the critical punching area from the column vicinity to the outside of the shear reinforced zone and enhances the RC slabs load capacity. A comprehensive parametric study using calibrated finite element models has also been conducted to analyse the effect of varying different parameters such as tensile reinforcement ratio and compressive reinforcement as well as the pattern, number and thickness of FRP strips (applied to strengthen RC slabs) on the behaviour of flat slabs. The results demonstrated that enhancing the tensile reinforcement ratios (including both steel reinforcements and FRP strips) can enhance the ultimate load capacity of the strengthened slabs, but reduces the ductility of the structure. Based on the results, flat slabs with compressive reinforcements could reach more load capacity and deflection (which resulted in having more ductility) as compared with the samples that do not include compressive reinforcements.

Two-phase flows arise in many areas of application such as in the study of coastal and hydraulic processes. Often the fluids involved can be modelled by incompressible phases which have disparate physical properties such as density and viscosity. We utilise a two-phase flow model of immiscible Newtonian fluids. A key part of this model is a set of variable coefficient Navier-Stokes equations. This thesis focuses on the numerical solution of the linear systems which arise after linearisation of these equations. Solving these systems often dominates the computation time when running simulations. One of the challenges in solving the model is that the density and viscosity coefficients are discontinuous and can have large jumps between the two phases. In this work we consider preconditioned iterative Krylov methods to solve the large and sparse linear systems and pay particular attention to incorporating the highly varying coefficients into the block preconditioners that we propose. We will see that such considerations can be essential in order to obtain good performance. An important issue is the scalability of the solution methodology. Here, we will study how the convergence of the iterative solver depends on a grid parameter which controls the refinement of the computational mesh. We will see that the novel preconditioners we propose can lead to convergence which is effectively independent of the grid parameter. We also investigate dependence on other model parameters such as the Reynolds number as well as the density and viscosity ratios between the two fluids. Another topic we examine is the use of a multipreconditioned iterative method allowing more than one preconditioner to be used simultaneously. Our results using this approach show some promising features. Finally, we consider an implementation within a more realistic model used in practice for simulating complex air-water flows. In particular, we will provide results for a problem modelling the breaking of a dam.

We use a continuum phase field model to study pattern formation in Diblock Copolymers (DBCs), with an emphasis on terracing in thin films. Phase separation in DBCs is an important process in nanotechnology that yields an array of highly regular structures on very fine length scales that are of enormous utility. Terraces are ubiquitous in thin film DBC melts. It was long believed that modelling these structures was an intractable problem, but recent attempts to do so using the Self Consistent (Mean) Field Theory (SCFT) were successful. In this thesis, we use the simpler and more flexible Ohta-Kawasaki (non-local Cahn-Hilliard) Density Functional Theory (DFT) model, derived from the SCFT, with obstacle homogeneous, and long-range energy components. Our main result is that such models can be used to capture terraces in the simplified 2D scenarios we examine, and that these structures form as an inevitable consequence of a simple energy minimisation process. We consider a 2-phase model for DBC micro-structure evolution, which we analyse rigorously and comprehensively. We study stationary and dynamic (gradient flow) variants of the problem, using the latter as the basis for an efficient but provably convergent numerical method for simulating DBC evolution. Our theoretical analysis and numerical method are based on a common Moreau-Yosida regularisation of the homogeneous energy component. We use the Semi-Smooth Newton (SSN) method to resolve the non-linearities in our numerical method. Comprehensive numerical experiments confirm our theoretical predictions, including a Method of Manufactured Solutions (MMS) test for the order of convergence. We propose a complete specification of the equally spaced global energy minimisers of the symmetric, 2-phase problem in 1D, and offer a practical algorithm for finding them. We explicitly characterise the Order-Disorder Transition (ODT) in parameter space. We extend our 2-phase model and methods to a 3-phase situation involving a copolymer/homopolymer blend, in which the free surface of the copolymer can adjust to form terraces. We develop energy estimates for energy-minimising 1D morphologies in the sharp interface limit, which we then combine via a double-tangent construction to analyse low energy, 2D structures. We argue that thin films with terraces have a lower free energy than their flat counterparts. Time-dependent numerical simulations in 2D confirm these predictions and provide insight into the evolution of these structures. We replicate and extend recent SCFT results, by performing a numerical experiment in which a terrace forms from a flat initial state. The model, tools and methods developed in this thesis should serve as a solid foundation for researching more complex, multi-step 3D terraces.

This thesis is concerned with the evaluation, in-situ, of the elastic strain energy release rate of cracks. This can define the criteria for crack propagation, and it is usually necessary to obtain this via calculation from the geometry and applied load. A new method is proposed, based on the conjoint use of digital image correlation to measure full-field displacements and finite element to extract the strain energy release rate of surface cracks. It has been extended to 3-D datasets with the use of digital volume correlation and tomographic imaging. A finite element model with imported full-field displacements measured by DIC/DVC acting as boundary conditions is solved and the J-integral is calculated. For linear elastic materials, modal contributions can be separated via the interaction integral. The method has been benchmarked using synthetic datasets to assess its sensitivity to noise and experimental uncertainties. It is very robust to experimental noise and can be used without knowledge of the specimen geometry and applied loads. The application of the method in 2-D is demonstrated in an analysis of experimental data for a mode I fatigue crack, introduced to an aluminium alloy compact tension specimen. Analysis of mixed-mode cracks in 2-D is shown on a PMMA sample with the Arcan geometry. In 3-D, static loading of a fatigue crack in nodular graphite cast iron is studied and the results from the method are compared with those obtained via a field-fitting approach. Diffraction analysis of polycrystalline materials can determine the full tensor of the elastic strains within them. Maps of elastic strains can thus be obtained typically using synchrotron X-rays or neutrons. A method is presented to calculate the elastic strain energy release rate of a crack from 2-D diffraction strain maps. The diffraction data is processed via a finite element approach to obtain the parameters required to calculate the $J$-integral. A validation is presented using a synthetic dataset from a finite element model. Its experimental application is demonstrated in an analysis of synchrotron X-ray diffraction strain maps of a propagating fatigue crack in a bainitic steel, before and after an overload. Finally, a complex case study of stable fracture propagation in polygranular isotropic nuclear graphite is presented. Synchrotron X-ray tomography and strain mapping by diffraction were combined with DVC and image analysis to extract the full-field displacements and elastic crystal strains. The displacement fields have been analysed using the developed methods to extract the critical strain energy release rate for crack propagation. Non-linear properties described the effect of microcracking on the elastic modulus in the fracture process zone. The analysis was verified by comparison of the predicted and measured elastic strain fields.

This thesis presents a series of algorithmic innovations in Combinatorial Compressed Sensing and Persistent Homology. The unifying strategy across these contributions is in translating structural patterns in the underlying data into specific algorithmic designs in order to achieve: better guarantees in computational complexity, the ability to operate on more complex data, highly efficient parallelisations, or any combination of these.

Soft materials have a wide range of applications, which include the production of masks for nano–lithography, the separation of membranes with nano–pores, and the preparation of nano–size structures for electronic devices. Self–organization in soft matter is a primary mechanism for the formation of structure. Block copolymers are long chain molecules composed of several different polymer blocks covalently bonded into a single macromolecule, which belong to an important class of soft materials which can self–assemble into different nano–structures due to their natural ability to microphase separate. Experimental and theoretical studies of block copolymers are quite challenging and, without computer simulations, it is difficult and problematic to analyse modern experiments. The Cell Dynamics Simulation (CDS) technique is a fast and accurate computational technique, which has been used to investigate block copolymers. The stability has been analysed by making use of different discrete Laplacian operators using well–chosen time steps in CDS. This analysis offers stability conditions for phase–field, based on the Cahn–Hilliard Cook (CHC) equations of which CDS is the finite difference approximation. To overcome grid related artefacts (discretization errors) in the computational grid, the study has been done for employing an isotropic Laplacian operator in the CDS framework. Several 2D and 3D discrete Laplacians have been quantitatively compared for their isotropy. The novel 2D 9–point BV(D2Q9) isotropic stencil operators have been derived from the B.A.C. van Vlimmeren method and their isotropy measure has been determined optimally better than other exiting 2D 9–point discrete Laplacian operators. Overall, the stencils in 9–point family Laplacians in 2D and the 19–point stencil operators in 3D have been found to be optimal in terms of isotropy and time step stability. Considerable implementation of Laplacians with good isotropy has played an important role in achieving a proper structure factor in modelling methods of block copolymers. The novel models have been developed by implementing CDS via more stable implicit methods, including backward Euler, Crank–Nicolson (CN) and Alternating Direction Implicit (ADI) methods. The CN scheme were implemented for both one order and two order parameter systems in CDS and successful results were obtained compared to forward Euler method. Due to the implementation of implicit methods, the CDS has achieved second–order accuracy both in time and space and it has become stronger, robust and more stable technique for simulation of the phase–separation phenomena in soft materials.

This dissertation presents the design and development of a novel hydrodynamic thrust bearing test rig featuring a new (patent pending) pressure-feedback control system for maintaining static bearing alignment. This research aims to provide an enhanced understanding of how the critical operational characteristics of fixed-geometry hydrodynamic thrust bearings including minimum oil film thickness (MOFT), hydrodynamic pressure distribution, and bearing temperature are affected by variability in bearing pad taper geometry under different speed, load, and oil conditions. Further, a new (patent pending) additively manufactured (AM) thrust bearing fabricated using direct metal laser sintering (DMLS) is experimentally evaluated to determine in-service viability. To support the experimental data obtained in the variable taper experiment, a Matlab simulation code is developed using the Reynolds equation to generate numerically predicted performance data for direct comparison. The AM thrust bearing is experimentally compared to a traditionally manufactured cast alloy bearing with identical surface geometry. For the variable taper study, trends in performance established by the numerical analysis show mutually agreeable results compared to experimental data. The average percent deviation of the experimentally gathered change in MOFT as load is increased with respect to the numerically predicted values is 24%. Comparison of experimental to numerical pressure distribution data shows an overall average percent deviation of 32%. For the AM vs. traditionally manufactured bearing experiment, the AM bearing showed an average increase in minimum oil film thickness of 53%, an average increase in trailing edge hydrodynamic pressure of 116%, while exhibiting an average decrease in bearing temperature of 1%.

Subsea pipelines form an integral part of the infrastructure associated with offshore oil and gas developments. These pipelines fulfill a range of functions from linking extraction wells to other subsea infrastructure to transporting products onshore, or to a central processing facility. Ancillary pipelines may also be present for gas or water injection to the reservoir or transporting additives. Pipelines are typically installed directly onto the seabed and, in the absence of significant drivers to undertake burial operations, they may remain on the seabed for the remainder of their design life. This is typically the case for deepwater developments. Subsea pipelines are subjected to a wide range of load cases including, self weight, installation loads, thermal and pressure driven expansion and hydrodynamic loading. Design of pipeline systems to accommodate these load cases requires an understanding of pipe-soil interaction. This thesis reports the results of a research study investigating pipe-soil interaction on a clay seabed, as relevant to the design of subsea pipeline systems. This study has utilised numerical analysis techniques based on the finite difference code FLAC to investigate a range of problem definitions. These problem definitions include pipelines subject to both vertical loading (V) and combined vertical and horizontal (V-H) loading. Factors such as variation in interface conditions, large strain and large displacement effects, soil unit weight effects and variation in shear strength conditions were considered in these problem definitions. Reliability based analysis techniques have also been used to investigate both V and V-H loading problem definitions. The analyses and investigations undertaken as part of this study generally achieved the following; reproduction and validation of earlier research with additional interpretation, extension of problem definitions to deeper pipeline embedment depths and investigation of pipe-soil interaction problem definitions that have not previously been considered. Reliability based analysis techniques have also provided some interesting insights into the impact of soil shear strength variation as well as providing a fundamental link between safety factors and probability of failure. Application to design practice of this, and similar studies, has been considered as part of this thesis and potential areas for future research have also been suggested.

In this thesis we study methods for solving the neutron transport equation (or linear Boltzmann equation). This is an integro-differential equation that describes the behaviour of neutrons during a nuclear fission reaction. Applications of this equation include modelling behaviour within nuclear reactors and the design of shielding around x-ray facilities in hospitals. Improvements in existing modelling techniques are an important way to address environmental and safety concerns of nuclear reactors, and also the safety of people working with or near radiation. The neutron transport equation typically has seven independent variables, however to facilitate rigorous mathematical analysis we consider the monoenergetic, steady-state equation without fission, and with isotropic interactions and isotropic source. Due to its high dimension, the equation is usually solved iteratively and we begin by considering a fundamental iterative method known as source iteration. We prove that the method converges assuming piecewise smooth material data, a result that is not present in the literature. We also improve upon known bounds on the rate of convergence assuming constant material data. We conclude by numerically verifying this new theory. We move on to consider the use of a specific, well-known diffusion equation to approximate the solution to the neutron transport equation. We provide a thorough presentation of its derivation (along with suitable boundary conditions) using an asymptotic expansion and matching procedure, a method originally presented by Habetler and Matkowsky in 1975. Next we state the method of diffusion synthetic acceleration (DSA) for which the diffusion approximation is instrumental. From there we move on to explore a new method of seeing the link between the diffusion and transport equations through the use of a block operator argument. Finally we consider domain decomposition algorithms for solving the neutron transport equation. Such methods have great potential for parallelisation and for the local application of different solution methods. A motivation for this work was to build an algorithm applying DSA only to regions of the domain where it is required. We give two very different domain decomposed source iteration algorithms, and we prove the convergence of both of these algorithms. This work provides a rigorous mathematical foundation for further development and exploration in this area. We conclude with numerical results to illustrate the new convergence theory, but also solve a physically-motivated problem using hybrid source iteration/ DSA algorithms and see significant reductions in the required computation time.

This thesis discusses several topics related to interpolation and how it is used in numerical analysis. It begins with an overview of the aspects of interpolation theory that are relevant to the discussion at hand before presenting three new contributions to the field. The first new result is a detailed error analysis of the barycentric formula for trigonometric interpolation in equally-spaced points. We show that, unlike the barycentric formula for polynomial interpolation in Chebyshev points (and contrary to the main view in the literature), this formula is not always stable. We demonstrate how to correct this instability via a rewriting of the formula and establish the forward stability of the resulting algorithm. Second, we consider the problem of trigonometric interpolation in grids that are perturbations of equally-spaced grids in which each point is allowed to move by at most a fixed fraction of the grid spacing. We prove that the Lebesgue constant for these grids grows at a rate that is at most algebraic in the number of points, thus answering questions put forth by Trefethen and Weideman about the robustness of numerical methods based on trigonometric interpolation in points that are uniformly distributed but not equally-spaced. We use this bound to derive theorems about the convergence rate of trigonometric interpolation in these grids and also discuss the related question of quadrature. Specifically, we prove that if a function has V ≥ 1 derivatives, the Vth of which is Hölder continuous (with a Hölder exponent that depends on the size of the maximum allowable perturbation), then the interpolants converge uniformly to the function at an algebraic rate; larger values of V lead to more rapid convergence. A similar statement holds for the corresponding quadrature rule. We also consider what analogue, if any, there is for trigonometric interpolation of the famous 1/4 theorem of Kadec from sampling theory that restricts the size of the perturbations one can make to the integers and still be guaranteed to have a set of stable sampling for the Paley-Wiener space. We present numerical evidence suggesting that in the discrete case, the 1/4 threshold takes the form of a threshold for the boundedness of a "2-norm Lebesgue constant" and does not appear to have much significance in practice. We believe that these are the first results regarding this problem to appear in the literature. While we do not believe the results we establish are the best possible quantitatively, they do (rigorously) capture the main features of trigonometric interpolation in perturbations of equally-spaced grids. We make several conjectures as to what the optimal results may be, backed by extensive numerical results. Finally, we consider a new application of interpolation to numerical linear algebra. We show that recently developed methods for computing the eigenvalues of a matrix by dis- cretizing contour integrals of its resolvent are equivalent to computing a rational interpolant to the resolvent and finding its poles. Using this observation as the foundation, we develop a method for computing the eigenvalues of real symmetric matrices that enjoys the same advantages as contour integral methods with respect to parallelism but employs only real arithmetic, thereby cutting the computational cost and storage requirements in half.

We develop expansion methods as a new computational approach towards high-dimensional partial differential equations (PDEs), particularly of such type as arising in the valuation of financial derivatives. The proposed methods are extended from [41] and use principal component analysis (PCA) of the underlying process in combination with a Taylor expansion of the value function into solutions to low-dimensional PDEs. They enable calculation of highly accurate approximate solutions with computational complexity polynomial in the number of dimensions for PDEs with a low number of dominant principal components. For the case of PDEs with constant coefficients, we show existence of expansion solutions and prove theoretical error bounds. We give a precise characterisation of when our methods can be applied and construct specific examples of a first and second order version. We provide numerical results showing that the empirically observed convergence speeds are in agreement with the theoretical predictions. For the case of PDEs with varying coefficients, we give a heuristic motivation using the Parametrix approach and empirically test the methods' accuracy for a range of variable parameter stock models. We demonstrate the applicability of our expansion methods to real-world securities pricing problems by considering path-dependent and early-exercise options in the LIBOR market model. Using the example of Bermudan swaptions and Ratchet floors, which are considered difficult benchmark problems, we give a careful analysis of the numerical accuracy and computational complexity. We are able to demonstrate that for problems with medium to high dimensionality, around 60-100, and moderate time horizons, the presented PDE methods deliver results comparable in accuracy to benchmark state-of-the-art Monte Carlo methods in similar or (significantly) faster run time.

This thesis concerns the analytical and practical aspects of applying the Closest Point Method to solve elliptic partial differential equations (PDEs) on smooth surfaces and domains with smooth boundaries. A new numerical scheme is proposed to solve surface elliptic PDEs and a novel geometric multigrid solver is constructed to solve the resulting linear system. The method is also applied to coupled bulk-surface problems. A new embedding equation in a narrow band surrounding the surface is formulated so that it agrees with the original surface PDE on the surface and has a unique solution which is constant along the normals to the surface. The embedding equation is then discretized using standard finite difference scheme and barycentric Lagrange interpolation. The resulting scheme has 2nd-order accuracy in practice and is provably 2nd-order convergent for curves without boundary embedded in &Ropf;2. To apply the method to solve elliptic equations on surfaces and domains with boundaries, the "ghost" point approach is adopted to handle Dirichlet, Neumann and Robin boundary conditions. A systematic method is proposed to represent values of ghost points by values of interior points according to boundary conditions. A novel geometric multigrid method based on the closest point representation of the surface is constructed to solve the resulting large sparse linear systems. Multigrid solvers are designed for surfaces with or without boundaries and domains with smooth boundaries. Numerical results indicate that the convergence rate of the multigrid solver stays roughly the same as we refine the mesh, as is desired of a multigrid algorithm. Finally the above methods are combined to solve coupled bulk-surface PDEs with some applications to biology.

We propose a discontinuous Galerkin finite element method (DGFEM) for fully nonlinear elliptic Hamilton--Jacobi--Bellman (HJB) partial differential equations (PDE) of second order with Cordes coefficients. Our analysis shows that the method is both consistent and stable, with arbitrarily high-order convergence rates for sufficiently regular solutions. Error bounds for solutions with minimal regularity show that the method is generally convergent under suitable choices of meshes and polynomial degrees. The method allows for a broad range of hp-refinement strategies on unstructured meshes with varying element sizes and orders of approximation, thus permitting up to exponential convergence rates, even for nonsmooth solutions. Numerical experiments on problems with nonsmooth solutions and strongly anisotropic diffusion coefficients demonstrate the significant gains in accuracy and computational efficiency over existing methods. We then extend the DGFEM for elliptic HJB equations to a space-time DGFEM for parabolic HJB equations. The resulting method is consistent and unconditionally stable for varying time-steps, and we obtain error bounds for both rough and regular solutions, which show that the method is arbitrarily high-order with optimal convergence rates with respect to the mesh size, time-step size, and temporal polynomial degree, and possibly suboptimal by an order and a half in the spatial polynomial degree. Exponential convergence rates under combined hp- and τq-refinement are obtained in numerical experiments on problems with strongly anisotropic diffusion coefficients and early-time singularities. Finally, we show that the combination of a semismooth Newton method with nonoverlapping domain decomposition preconditioners leads to efficient solvers for the discrete nonlinear problems. The semismooth Newton method has a superlinear convergence rate, and performs very effectively in computations. We analyse the spectral bounds of nonoverlapping domain decomposition preconditioners for a model problem, where we establish sharp bounds that are explicit in both the mesh sizes and polynomial degrees. We then go beyond the model problem and show computationally that these algorithms lead to efficient and competitive solvers in practical applications to fully nonlinear HJB equations.

This thesis presents an investigation into sources of variability in activating and processing numerical information. Chapter 1 provides an overview of research literature exploring the ways in which magnitude information can be represented, and how models relating to number information have developed. These theoretical models are addressed in relation to the neural representation of number, and the range of behavioural markers which suggest an association between spatial and numerical processing. Chapter 2 using a dual-task paradigm investigated whether magnitude information is accessed on perceiving numbers, or if this information is linked to response selection or execution. Previous research studies investigating this question produced inconsistent findings (Oriet, Tombu & Jolicoeur, 2005; Sigman & Dehaene, 2005) with regard to the locus of magnitude processing; the findings of Experiments 1-3 reliably support access to magnitude information during response selection. Chapter 3 explored the activation of spatial-numerical response associations, where response-irrelevant magnitude information was not represented by a single stimulus (i.e. an Arabic digit) but by a numerosity representation. Experiments 4-7 found a strong association between spatial-orientation processing and numerical magnitude, but no association with perceptual-colour processing, extending previous work by Fias, Lammertyn and Lauwereyns (2001) regarding the neural overlap between the attended and irrelevant stimulus dimensions. However the strength of this association was found to be inconsistent across the number range. Chapter 4 investigated the impact of healthy aging on the presence of neural-overlap in processing spatial-numerical information, further developing the paradigm used in Chapter 3, and addressed direct predictions from the literature as to how age should influence these associations (Wood, Willmes, Nuerk & Fischer, 2008). Experiments 8-11 found evidence for spatial-numerical associations across the lifespan, but that the strength of these effects were moderated by 5 task instruction. Chapter 5 was designed to assess aging differences in numerical and spatial processing with a battery of tests and the extent to which other sources of individual difference (sex, embodiment) have a measureable impact. A range of standardised measures were used to assess verbal ability, mathematical processing, and spatial working memory alongside behavioural measures of spatial numerical associations. Experiment 12 provided evidence of aging and sex differences in different cognitive tasks and a marginal impact of embodiment on spatial-numerical processing; however the effect of embodiment was not supported in a larger more homogenous sample in Experiment 13. Chapter 6 reflects on the current findings and provides contextual information on how they align with previous research, outlining how evidence from the thesis extend current research paradigms and provides new evidence regarding the maintenance of spatial-numerical associations in healthy aging. Methodologies developed in the thesis are considered with relation to how they may be applied to assess individual differences in early number acquisition in children. Finally the discussion outlines methods and controversies within the field of numerical cognition, with consideration of new methods for measuring the strength of spatial-numerical associations (Pinhas, Tzelgov, & Ganor-Stern, 2012), alongside the potential application of modelling techniques to investigate individual differences in task performance.

Estimation of the fetal heart rate (FHR) has gained interest in the last century; low heart rate variability has been studied to identify intrauterine growth restricted fetuses (prepartum), and abnormal FHR patterns have been associated with fetal distress during delivery (intrapartum). Several monitoring techniques have been proposed for FHR estimation, including auscultation and Doppler ultrasound. This thesis focuses on the extraction of the non-invasive fetal electrocardiogram (NI-FECG) recorded from a limited set of abdominal sensors. The main challenge with NI-FECG extraction techniques is the low signal-to-noise ratio of the FECG signal on the abdominal mixture signal which consists of a dominant maternal ECG component, FECG and noise. However the NI-FECG offers many advantages over the alternative fetal monitoring techniques, the most important one being the opportunity to enable morphological analysis of the FECG which is vital for determining whether an observed FHR event is normal or pathological. In order to advance the field of NI-FECG signal processing, the development of standardised public databases and benchmarking of a number of published and novel algorithms was necessary. Databases were created depending on the application: FHR estimation with or without maternal chest lead reference or directed toward FECG morphology analysis. Moreover, a FECG simulator was developed in order to account for pathological cases or rare events which are often under-represented (or completely missing) in the existing databases. This simulator also serves as a tool for studying NI-FECG signal processing algorithms aimed at morphological analysis (which require underlying ground truth annotations). An accurate technique for the automatic estimation of the signal quality level was also developed, optimised and thoroughly tested on pathological cases. Such a technique is mandatory for any clinical applications of FECG analysis as an external confidence index of both the input signals and the analysis outputs. Finally, a Bayesian filtering approach was implemented in order to address the NI-FECG morphology analysis problem. It was shown, for the first time, that the NI-FECG can allow accurate estimation of the fetal QT interval, which opens the way for new clinical studies on the development of the fetus during the pregnancy.

Membrane protein mechanotransduction is the altered function of an integral membrane protein in response to mechanical force. Such mechanosensors are found in all kingdoms of life, and increasing numbers of membrane proteins have been found to exhibit mechanosensitivity. How they mechanotransduce is an active research area and the topic of this thesis. The methodology employed is classical molecular dynamics (MD) simulations. MD systems are complex, and two programs were developed to reduce this apparent complexity in terms of both visual abstraction and statistical analysis. Bendix detects and visualises helices as cylinders that follow the helix axis, and quantifies helix distortion. The functionality of Bendix is demonstrated on the symporter Mhp1, where a state is identified that had hitherto only been proposed. InterQuant tracks, categorises and orders proximity between parts of an MD system. Results from multiple systems are statistically interrogated for reproducibility and significant differences at the resolution of protein chains, residues or atoms. Using these tools, the interaction between membrane and the Escherichia coli mechanosensitive channel of small conductance, MscS, is investigated. Results are presented for crystal structures captured in different states, one of which features electron density proposed to be lipid. MD results supports this hypothesis, and identify differential lipid interaction between closed and open states. It is concluded that propensity for lipid to leave for membrane bulk drives MscS state stability. In a subsequent study, MscS is opened by membrane surface tension for the first time in an MD setup. The gating mechanism of MscS is explored in terms of both membrane and protein deformation in response to membrane stretch. Using novel tension methodology and the longest MD simulations of MscS performed to date, a molecular basis for the Dashpot gating mechanism is proposed. Lipid emerges as an active structural element with the capacity to augment protein structure in the protein structure-function paradigm.

The first contribution of this thesis is a new regularity theorem for time harmonic Maxwell's equations with less than Lipschitz complex anisotropic coefficients. By using the Lp theory for elliptic equations, it is possible to prove H1 and Hölder regularity results, provided that the coefficients are W1,p for some p = 3. This improves previous regularity results, where the assumption W1,∞ for the coefficients was believed to be optimal. The method can be easily extended to the case of bi-anisotropic materials, for which a separate approach turns out to be unnecessary. The second focus of this work is the boundary control of the Helmholtz and Maxwell equations to enforce local constraints inside the domain. More precisely, we look for suitable boundary conditions such that the corresponding solutions and their derivatives satisfy certain local non-zero constraints. Complex geometric optics solutions can be used to construct such illuminations, but are impractical for several reasons. We propose a constructive approach to this problem based on the use of multiple frequencies. The suitable boundary conditions are explicitly constructed and give the desired constraints, provided that a finite number of frequencies, given a priori, are chosen in a fixed range. This method is based on the holomorphicity of the solutions with respect to the frequency and on the regularity theory for the PDE under consideration. This theory finds applications to several hybrid imaging inverse problems, where the unknown coefficients have to be imaged from internal measurements. In order to perform the reconstruction, we often need to find suitable boundary conditions such that the corresponding solutions satisfy certain non-zero constraints, depending on the particular problem under consideration. The multiple frequency approach introduced in this thesis represents a valid alternative to the use of complex geometric optics solutions to construct such boundary conditions. Several examples are discussed.

Multiphysics systems with interface coupling are used to model a variety of physical phenomena, such as arterial blood flow, air flow around aeroplane wings, or interactions between surface and ground water flows. Numerical methods enable the practical application of these models through computer simulations. Specifically a high level of detail and accuracy is achieved in finite element methods by discretisations which use extremely large numbers of degrees of freedom, rendering the solution process challenging from the computational perspective. In this thesis we address this challenge by developing a twofold strategy for improving the efficiency of standard finite element coupled solvers. First, we propose to solve a monolithic coupled problem using block-preconditioned GMRES with a new Schur complement approximation. This results in a modular and robust method which significantly reduces the computational cost of solving the system. In particular, numerical tests show mesh-independent convergence of the solver for all the considered problems, suggesting that the method is well-suited to solving large-scale coupled systems. Second, we derive an adjoint-based formula for goal-oriented a posteriori error estimation, which leads to a time-space mesh refinement strategy. The strategy produces a mesh tailored to a given problem and quantity of interest. The monolithic formulation of the coupled problem allows us to obtain expressions for the error in the Lagrange multiplier, which often represents a physically relevant quantity, such as the normal stress on the interface between the problem components. This adaptive refinement technique provides an effective tool for controlling the error in the quantity of interest and/or the size of the discrete system, which may be limited by the available computational resources. The solver and the mesh refinement strategy are both successfully employed to solve a coupled Stokes-Darcy-Stokes problem modelling flow through a cartridge filter.

In mathematical finance, the sensitivities of option prices to various market parameters, also known as the “Greeks”, reflect the exposure to different sources of risk. Computing these is essential to predict the impact of market moves on portfolios and to hedge them adequately. This is commonly done using Monte Carlo simulations. However, obtaining accurate estimates of the Greeks can be computationally costly. Multilevel Monte Carlo offers complexity improvements over standard Monte Carlo techniques. However the idea has never been used for the computation of Greeks. In this work we answer the following questions: can multilevel Monte Carlo be useful in this setting? If so, how can we construct efficient estimators? Finally, what computational savings can we expect from these new estimators? We develop multilevel Monte Carlo estimators for the Greeks of a range of options: European options with Lipschitz payoffs (e.g. call options), European options with discontinuous payoffs (e.g. digital options), Asian options, barrier options and lookback options. Special care is taken to construct efficient estimators for non-smooth and exotic payoffs. We obtain numerical results that demonstrate the computational benefits of our algorithms. We discuss the issues of convergence of pathwise sensitivities estimators. We show rigorously that the differentiation of common discretisation schemes for Ito processes does result in satisfactory estimators of the the exact solutions’ sensitivities. We also prove that pathwise sensitivities estimators can be used under some regularity conditions to compute the Greeks of options whose underlying asset’s price is modelled as an Ito process. We present several important results on the moments of the solutions of stochastic differential equations and their discretisations as well as the principles of the so-called “extreme path analysis”. We use these to develop a rigorous analysis of the complexity of the multilevel Monte Carlo Greeks estimators constructed earlier. The resulting complexity bounds appear to be sharp and prove that our multilevel algorithms are more efficient than those derived from standard Monte Carlo.