A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalised associated Legendre functions

Spherical harmonic expansions form partial sums of fully normalised associated Legendre functions (ALFs). However, when evaluated increasingly close to the poles, the ultra-high degree and order (e.g. 2700) ALFs range over thousands of orders of magnitude. This causes existing recursion techniques for computing values of individual ALFs and their derivatives to fail. A common solution in geodesy is to evaluate these expansions using Clenshaw's method, which does not compute individual ALFs or their derivatives. Straightforward numerical principles govern the stability of this technique. Elementary algebra is employed to illustrate how these principles are implemented in Clenshaw's method. It is also demonstrated how existing recursion algorithms for computing ALFs and their first derivatives are easily modified to incorporate these same numerical principles. These modified recursions yield scaled ALFs and first derivatives, which can then be combined using Horner's scheme to compute partial sums, complete to degree and order 2700, for all latitudes (except at the poles for first derivatives). This exceeds any previously published result. Numerical tests suggest that this new approach is at least as precise and efficient as Clenshaw's method. However, the principal strength of the new techniques lies in their simplicity of formulation and implementation, since this quality should simplify the task of extending the approach to other uses, such as spherical harmonic analysis

A modified equation analysis for immersed boundary methods based on volume penalization: applications to linear advection-diffusion and high-order discontinuous Galerkin schemes

The Immersed Boundary Method (IBM) is a popular numerical approach to impose boundary conditions without relying on body-fitted grids, thus reducing the costly effort of mesh generation. To obtain enhanced accuracy, IBM can be combined with high-order methods (e.g., discontinuous Galerkin). For this combination to be effective, an analysis of the numerical errors is essential. In this work, we apply, for the first time, a modified equation analysis to the combination of IBM (based on volume penalization) and high-order methods (based on nodal discontinuous Galerkin methods) to analyze a priori numerical errors and obtain practical guidelines on the selection of IBM parameters. The analysis is performed on a linear advection-diffusion equation with Dirichlet boundary conditions. Three ways to penalize the immerse boundary are considered, the first penalizes the solution inside the IBM region (classic approach), whilst the second and third penalize the first and second derivatives of the solution. We find optimal combinations of the penalization parameters, including the first and second penalizing derivatives, resulting in minimum errors. We validate the theoretical analysis with numerical experiments for one- and two-dimensional advection-diffusion equations

A projection hybrid high order finite volume/finite element method for incompressible turbulent flows

In this paper the projection hybrid FV/FE method presented in Busto et al. 2014 is extended to account for species transport equations. Furthermore, turbulent regimes are also considered thanks to the $k-\varepsilon$ model. Regarding the transport diffusion stage new schemes of high order of accuracy are developed. The CVC Kolgan-type scheme and ADER methodology are extended to 3D. The latter is modified in order to profit from the dual mesh employed by the projection algorithm and the derivatives involved in the diffusion term are discretized using a Galerkin approach. The accuracy and stability analysis of the new method are carried out for the advection-diffusion-reaction equation. Within the projection stage the pressure correction is computed by a piecewise linear finite element method. Numerical results are presented, aimed at verifying the formal order of accuracy of the scheme and to assess the performance of the method on several realistic test problems.Comment: arXiv admin note: text overlap with arXiv:1802.1058

Numerical simulation of reaction fronts in dissipative media

Fronts of reaction in certain systems (such as so-called solid flames and detonation fronts) can be simulated by a single-equation phenomenological model of Strunin (1999, 2009). This is a high-order nonlinear partial differential equation describing the shape of the front as a function of spatial coordinates and time. The equation is of active-dissipative type, with 6th-order spatial derivative. For one-dimensional case, the equation was previously solved using the Galerkin method, but only one numerical experiment with limited information on the dynamics was obtained. For two-dimensional case only two numerical ex- periments were reported so far, in which a low-accuracy infinite difference scheme was used. In this thesis, we use a more recent and sophisticated method, namely the one-dimensional integrated radial basis function networks (1D-IRBFN). The method had been developed by Tran-Cong and May-Duy (2001, 2003) and successfully applied to several problems such as structural analysis, viscoelastic flows and fluid-structure interaction. In contrast to commonly used approaches, where a function of interest is differentiated to give approximate derivatives, leading to a reduction in convergence rate for derivatives (and this reduction increases with derivative order, which magnifes errors), the 1D-IRBFN method uses the integral formulation. It utilizes spectral approximants to represent highest-order derivatives under consideration. They are then integrated analytically to yield approximate expressions for lower-order derivatives and the function itself. In this thesis the following main results are obtained. A numerical program implementing the 1D-IRBFN method is developed in Matlab to solve the equation of interest. The program is tested by (a) constructing a forced version of the equation, which allows analytical solution, and verifying the numerical solution against the analytical solution; (b) reproducing one-dimensional spinning waves obtained from the model previously. A modified version of the program is successfully applied to similar high-order equations modelling auto-pulses in fluid flows with elastic walls. We obtained numerically and analyzed a far richer variety of one-dimensional dynamics of the reaction fronts. Two kinds of boundary conditions were used: homogeneous conditions on the edges of the domain, and periodic conditions corresponding to periodicity of the front on a cylinder. The dependence of the dynamics on the size of the domain is explored showing how larger space accommodates multiple spinning waves. We determined the critical domain size (bifurcation point) at which non-trivial settled regimes become possible. We found a regime where the front is shaped as a pair of kinks separated by a relatively short distance. Interestingly, the pair moves in a stable joint formation far from the boundaries. A similar regime for three connected kinks is obtained. We demonstrated that the initial condition determines the direction of motion of the kinks, but not their size and velocity. This is typical for active-dissipative systems. The settled character of these regimes is demonstrated. We also applied the 1D-IRBFN method to two-dimensional topology corresponding to a solid cylinder. Stable spinning wave solutions are obtained for this case

Limiting techniques for the discontinuous Galerkin method on unstructured meshes

Many complex phenomena like shock-shock interactions, shock-vortex interactions, stratified flows, etc., are governed by nonlinear hyperbolic conservation laws. Higher order numerical schemes like the discontinuous Galerkin (DG) method are increasingly being used in computational fluid dynamics (CFD) codes to solve such conservation laws. Oscillations often develop in the numerical solutions obtained using higher-order methods. These spurious oscillations can lead to numerical instabilities and eventual degradation of the solution. Slope limiting is one of the mechanisms used to suppress such oscillations, and thus, stabilize the numerical solution. Slope limiters were originally introduced for finite volume (FV) methods, where the reconstructed slope in an element is modified to prevent oscillations while maintaining second-order accuracy and stability of the numerical solution. Limiters tailored specifically for the DG method have also been proposed. However, developing efficient higher order limiting techniques for the DG method on unstructured meshes remains an open problem. Many factors influence the design of a suitable limiter for the DG method, e.g., the domain discretization involved, the basis functions used to approximate the solution, the choice of limiting directions, and the neighborhood used to reconstruct the solution coefficients, to name a few. In this work, we propose high-order moment limiters for the DG method on unstructured two-dimensional (triangular, quadrilateral, curvilinear triangular) and three-dimensional (tetrahedral) meshes. The limiters work by relating the solution coefficients (moments) to directional derivatives along specified directions and limiting said directional derivatives independently using a one-dimensional slope limiter. The limiting is performed hierarchically starting at the highest moment and stops on reaching a set of moments/derivatives that remains unchanged, thereby preventing overlimiting. To use available computational resources efficiently, simulations often employ run-time adaptive mesh refinement strategies. In this thesis, we present a high-order limiter for the DG method on nonconforming triangular meshes that arise in such adaptive computations. Moreover, we also propose a simple algorithm to update the reconstruction stencil of elements in an adaptively refined triangular mesh. Our algorithm is implemented entirely on the graphics processing unit (GPU) and avoids race conditions. We provide numerical examples to show that limited solutions retain the theoretical rate of convergence and are robust in the presence of discontinuities. Finally, we perform wall clock studies to analyze the performance of the proposed limiter for the computational cost involved in executing the limiter

Fractional derivative models for the spread of diseases

This thesis considers the mathematical modelling of disease, using fractional differential equations in order to provide a tool for the description of memory effects. In Chapter 3 we illustrate a commensurate fractional order tumor model, and we find a critical value of the fractional derivative dependent on the parameter values of the model. For fractional derivatives of orders less than the critical value an unstable equilibrium point of the system becomes stable. In order to show changes in the observed areas of attraction of two stable points in the system, we then consider a fractional order SIR epidemic model and investigate the change from a monostable to a bistable system.;Chapter 4 considers a model for virus dynamics where the fractional orders for populations are different, called an incommensurate system. An approximate analytical solution for the characteristic equation of the incommensurate model is found when the different fractional orders are similar and close to the critical value of the fractional order of the commensurate system. In addition, the instability boundary is found as a function of both parameters. A comparison between analytical and numerical results shows the high accuracy of this approximation.;Chapter 5 consists of two parts, in the first part we generalise the integer Fisher's equation to be a space-time fractional differential equation and consider travelling wave solutions. In the second part we generalise an integer SIR model with spatial heterogeneity, which was studied by Murray [117], to a space-time fractional derivative model. We apply the (G0/G)-expansion method and find travelling wave solutions, although in this case we must consider the Jumarie's modified Riemann-Liouville fractional derivative. Finally, we consider the effect of changing the orders of time and space fractional derivatives on the location and speed of the travelling wave solution.This thesis considers the mathematical modelling of disease, using fractional differential equations in order to provide a tool for the description of memory effects. In Chapter 3 we illustrate a commensurate fractional order tumor model, and we find a critical value of the fractional derivative dependent on the parameter values of the model. For fractional derivatives of orders less than the critical value an unstable equilibrium point of the system becomes stable. In order to show changes in the observed areas of attraction of two stable points in the system, we then consider a fractional order SIR epidemic model and investigate the change from a monostable to a bistable system.;Chapter 4 considers a model for virus dynamics where the fractional orders for populations are different, called an incommensurate system. An approximate analytical solution for the characteristic equation of the incommensurate model is found when the different fractional orders are similar and close to the critical value of the fractional order of the commensurate system. In addition, the instability boundary is found as a function of both parameters. A comparison between analytical and numerical results shows the high accuracy of this approximation.;Chapter 5 consists of two parts, in the first part we generalise the integer Fisher's equation to be a space-time fractional differential equation and consider travelling wave solutions. In the second part we generalise an integer SIR model with spatial heterogeneity, which was studied by Murray [117], to a space-time fractional derivative model. We apply the (G0/G)-expansion method and find travelling wave solutions, although in this case we must consider the Jumarie's modified Riemann-Liouville fractional derivative. Finally, we consider the effect of changing the orders of time and space fractional derivatives on the location and speed of the travelling wave solution