Second-order Convergence Research Articles

Efficient second-order accurate exponential time differencing for time-fractional advection–diffusion–reaction equations with variable coefficients

Time-fractional advection–diffusion–reaction type equations are useful for characterizing anomalous transport processes. In this paper, linearly implicit as well as explicit generalized exponential time differencing (GETD) schemes are proposed for solving a class of such equations having time–space dependent coefficients. The implicit scheme, being unconditionally stable, is robust in handling the numerical instabilities in problems where the advection term is dominant. Regarding the error analysis, uniformly optimal second-order convergence rates are derived using time-graded meshes to counter the effect of the inherent singularity of the continuous solution. Implementation of generalized exponential integrators requires computing the action of Mittag-Leffler function of matrices on a vector, or on a matrix in the case of the implicit scheme. For cost-effective implementation, using global Padé approximants these computation tasks get reduced to solving linear systems. A new approach based on Sylvester equation formulation of the resulting linear systems is developed in this paper. This technique leads to significantly faster algorithms for implementing the GETD schemes. Numerical experiments are provided to illustrate the theoretical findings and to assert the efficiency of the Sylvester equation based approach. Application of this approach to an existing GETD scheme for solving a nonlinear subdiffusion problem is also discussed.

A new linearized ADI compact difference method on graded meshes for a nonlinear 2D and 3D PIDE with a WSK

In this work, a new linearized alternating direction implicit (ADI) compact difference method (CDM) is proposed for solving nonlinear two-dimensional (2D) and three-dimensional (3D) partial integrodifferential equation (PIDE) with a weakly singular kernel (WSK). The time derivative is treated by Crank-Nicolson (CN) method and the Riemann-Liouville (R-L) integral by product integration (PI) rule on graded meshes. The linear interpolation combining with Taylor formula is applied in time to deal with nonlinear term v∇v in interval (t0,t1), and linear interpolation concerning two previous time points is employed to deal with v∇v in intervals (tn−1,tn),n≥2. A linearized semi-discrete scheme is obtained, which can achieve second-order convergence in time. Then via introducing two kinds of compact difference operators to discretize the spatial derivatives. To improve the computing efficiency, we construct a ADI compact difference method. It is the first time that the ADI compact difference method is applied for the nonlinear 2D and 3D PIDE with a WSK. The advantage of our proposed scheme is that it not only has second-order accuracy in time and fourth-order accuracy in space, but also fast computational speed, just by solving the linear coupled equations for tridiagonal matrices. In addition, we prove the existence, uniqueness and convergence of 2D scheme. Four numerical examples in 2D and 3D are present to demonstrate our proposed method.

A note on the numerical approximation of Greeks for American-style options

In this note, we consider the approximation of the Greeks Delta and Gamma of American-style options through the numerical solution of time-dependent partial differential complementarity problems (PDCPs). This approach is very attractive as it can yield accurate approximations to these Greeks at essentially no additional computational cost during the numerical solution of the PDCP for the pertinent option value function. For the temporal discretization, the Crank–Nicolson method is arguably the most popular method in computational finance. It is well-known, however, that this method can have an undesirable convergence behaviour in the approximation of the Greeks Delta and Gamma for American-style options, even when backward Euler damping (Rannacher smoothing) is employed.In this note, for the temporal discretization of the PDCP, we study an interesting family of diagonally implicit Runge–Kutta (DIRK) methods together with the two-stage Lobatto IIIC method. Through ample numerical experiments for one- and two-asset American-style options, it is shown that these methods can yield a regular second-order convergence behaviour for the option value as well as for the Greeks Delta and Gamma. A mutual comparison reveals that the DIRK method with suitably chosen parameter θ is preferable.

A convergence analysis for the approximation to the solution of an age-structured population model with infinite lifespan

Considering the numerical approximation of the density distribution for an age-structured population model with unbounded lifespan on a compact interval [0,T], we prove second order of convergence for a discretization that adaptively selects its truncated age-interval according to the exponential rate of decay with age of the solution of the model. It appears that the adaptive capacity of the length in the truncated age-interval of the discretization to the infinity lifespan is a very convenient approach for a long-time integration of the model to establish the asymptotic behavior of its dynamics numerically. The analysis of convergence uses an appropriate weighted maximum norm with exponential weights to cope with the unbounded age lifespan. We report experiments to exhibit numerically the theoretical results and the asymptotic behavior of the dynamics for an age-structured squirrel population model introduced by Sulsky.

An analysis of second-order sav-filtered time-stepping finite element method for unsteady natural convection problems

This paper presents an unconditionally stable time-filtering algorithm for natural convection equations. The algorithm is based on the scalar auxiliary variables in the exponential function and adopts a completely discrete Back-Euler combining time filter scheme. The proposed scheme requires minimal invasive modification of the existing program to improve the time accuracy from first-order to second-order without increasing the computational complexity, and we demonstrate the unconditional stability of the proposed algorithm and analyze its second-order convergence. In addition, due to the increasing demand for low-memory solvers, the application of a time-adaptive algorithm can improve the accuracy and efficiency of the proposed algorithm, so we extend the method to variable step sizes and construct an adaptive algorithm. Finally, the effectiveness of the proposed method and the accuracy of the theoretical results are verified by numerical experiments.

Analysis of efficient preconditioner for solving Poisson equation with Dirichlet boundary condition in irregular three-dimensional domains

This paper analyzes modified ILU (MILU)-type preconditioners for efficiently solving the Poisson equation with Dirichlet boundary conditions on irregular domains. In [1], the second-order accuracy of a finite difference scheme developed by Gibou et al. [2] and the effect of the MILU preconditioner were presented for two-dimensional problems. However, the analyses do not directly extend to three-dimensional problems. In this paper, we first demonstrate that the Gibou method attains second-order convergence for three-dimensional irregular domains, yet the discretized Laplacian exhibits a condition number of O(h−2) for a grid size h. We show that the MILU preconditioner reduces the order of the condition number to O(h−1) in three dimensions. Furthermore, we propose a novel sectored-MILU preconditioner, defined by a sectorized lexicographic ordering along each axis of the domain. We demonstrate that this preconditioner reaches a condition number of order O(h−1) as well. Sectored-MILU not only achieves a similar or better condition number than conventional MILU but also improves parallel computing efficiency, enabling very efficient calculation when the dimensionality or problem size increases significantly. Our findings extend the feasibility of solving large-scale problems across a range of scientific and engineering disciplines.

Positivity preserving and unconditionally stable numerical scheme for the three-dimensional modified Fisher–Kolmogorov–Petrovsky–Piskunov equation

This paper introduces a numerical approach for the practical solution of the modified Fisher–Kolmogorov–Petrovsky–Piskunov equation that describes population dynamics. The diffusion term and nonlinear term is based on the operator splitting method and interpolation method, respectively. The analytic proof of the discrete maximum principle and positivity preserving for the numerical algorithm is demonstrated. Numerical solution calculated using the proposed method remains stable without blowing up, which implies that the proposed method is unconditionally stable. Numerical studies show that the proposed method is second-order convergence in space and first-order convergence in time. The performance and applicability of the proposed scheme are studied through various computational tests that present the effects of model parameters and evolution dynamics.

Prediction–correction projection immersed interface method for interfacial stress-induced flows and applications to electrohydrodynamics

We propose a robust, high-resolution prediction–correction projection immersed interface method (IIM) for solving the unsteady incompressible Navier–Stokes equations with traction boundary conditions, which arise from free surface flows driven by capillary, electric, and elastic forces. This method combines the advantages of traditional body-fitted moving mesh methods and immersed boundary methods (IBM), allowing for the accurate imposition of boundary conditions on free surfaces using Cartesian grids and providing detailed interface information that is typically smoothed out in traditional IBM. The irregular liquid-phase domain is embedded within a square region and discretized using a dynamically adaptive Cartesian mesh. The free surface is captured using a narrow-band level set with hybrid reinitialization. A prediction–correction projection scheme is constructed, incorporating additional pressure predictions and corrections to enhance robustness. The resulting Helmholtz/Poisson equations are solved using the augmented IIM for boundary value problems. Grid refinement analysis demonstrates second-order convergence of the L∞ error, even in challenging cases such as the oscillating drop test with low viscosity. We further apply this method to electrohydrodynamic (EHD) problems, constructing an implicit augmented IIM to solve the equations governing the electric field and charge conservation. Numerical experiments demonstrate that this method accurately addresses highly intense EHD phenomena, such as the formation of Taylor cones, highlighting its robustness.

Efficient option pricing in the rough Heston model using weak simulation schemes

We provide an efficient and accurate simulation scheme for the rough Heston model in the standard (H>0) as well as the hyper-rough regime ( -1/2 $ ]]> H > − 1 / 2 ). The scheme is based on low-dimensional Markovian approximations of the rough Heston process derived in [Bayer and Breneis, arXiv:2309.07023], and provides a weak approximation to the rough Heston process. Numerical experiments show that the new scheme exhibits second-order weak convergence, while the computational cost increases linearly with respect to the number of time steps. In comparison, existing schemes based on discretization of the underlying stochastic Volterra integrals such as Gatheral's HQE scheme show a quadratic dependence of the computational cost. Extensive numerical tests for standard and path-dependent European options and Bermudan options show the method's accuracy and efficiency.

Complex Roots-Finding Method of Nonlinear Equations

Various researchers developed numerical methods to solve nonlinear equations. But still, there is a need to improve these methods to get accurate results. Some techniques exist where convergence is not assured. These methods increase the problem's cost using a second/higher-order derivative. So, cost efficiency is also considered an issue in many fast-converging ways. This paper proposes a new numerical method to find complex and real roots of nonlinear equations. By taking an initial guess w0 ∈ C, the new approach contributes to a series of iterations converging to a real or complex solution. Moreover, the method evaluates complex roots even if a real number is taken as an initial guess. This innovation presents that the proposed method achieves second-order convergence. The number of iterations is also determined to check the performance of the new iterative scheme. Using Python 3.10.9 package, the authors have tested the method's efficacy on several numerical problems, presented the results in the tables, and illustrated them with the help of graphs.

Second-order error analysis of a corrected average finite difference scheme for time-fractional Cable equations with nonsmooth solutions

In this paper, we first formulate an average finite difference scheme with appropriate correction terms for the time-fractional initial-value problem of y′(t)+0RLDt1−αy(t)=g(t), where 0RLDt1−α denotes the Riemann–Liouville fractional derivative of order 1−α(α∈(0,1)). It is shown that, under some suitable conditions on the data, the numerical solution converges to the exact solution with order O(τ2), where τ is the time step size. The method is then extended to solve the time-fractional Cable equation, combined with a standard discretisation of the spatial derivatives on a uniform mesh. The second-order time convergence rate is proved. Several numerical examples are given to illustrate the good agreement with the theoretical analysis of the presented methods.

A Posteriori Error Analysis of Defect Correction Method for Singular Perturbation Problems With Discontinuous Coefficient and Point Source

Abstract The paper presents a defect correction method to solve singularly perturbed problems with discontinuous coefficient and point source. The method combines an inexpensive, lower-order stable, upwind difference scheme and a higher-order, less stable central difference scheme over a layer-adapted mesh. The mesh is designed so that most mesh points remain in the regions with rapid transitions. A posteriori error analysis is presented. The proposed numerical method is analyzed for consistency, stability, and convergence. The error estimates of the proposed numerical method satisfy parameter-uniform second-order convergence on the layer-adapted grid. The convergence obtained is optimal because it is free from any logarithmic term. The numerical analysis confirms the theoretical error analysis.

Unconditional optimal [formula omitted]-norm error estimate and superconvergence analysis of a linearized nonconforming finite element variable-time-step BDF2 method for the nonlinear complex Ginzburg–Landau equation

In this paper, a linearized fully discrete scheme combining the variable-time-step two-step backward differentiation formula (VSBDF2) in time and the nonconforming finite element methods (FEMs) in space is constructed and analyzed for the nonlinear complex Ginzburg–Landau equation. A novel convergence analysis approach is proposed, which shows that the H1-norm error of this scheme can reach optimal order in space and sharp second-order convergence in time. The key ingredients of our analysis are temporal–spatial error splitting approach, a new positive definiteness property of discrete orthogonal convolution (DOC) kernels with respect to complex sequences, an energy projection operator and a discrete Laplace operator Δh related to the nonconforming element space, and a relationship between Δh and classical L2 projection operator in complex space. Furthermore, through a specific treatment of the error equation and the combination technique of interpolating operator and projection operator, we prove that the proposed scheme has the global superconvergence result of order O(τ2+h2) in H1-norm for the first time, which further improves the efficiency of numerical computation. Here τ and h are the maximum time and spatial step sizes, respectively. It should be pointed out that all these error results are rigorously established under mild assumptions of the time step, and get rid of the usual grid ratio constraints between temporal and spatial step sizes in linearized scheme. At last, numerical experiments are provided to verify the theoretical analysis.

Automatic adjoint-based inversion schemes for geodynamics: reconstructing the evolution of Earth's mantle in space and time

Abstract. Reconstructing the thermo-chemical evolution of Earth's mantle and its diverse surface manifestations is a widely recognised grand challenge for the geosciences. It requires the creation of a digital twin: a digital representation of Earth's mantle across space and time that is compatible with available observational constraints on the mantle's structure, dynamics and evolution. This has led geodynamicists to explore adjoint-based approaches that reformulate mantle convection modelling as an inverse problem, in which unknown model parameters can be optimised to fit available observational data. Whilst there has been a notable increase in the use of adjoint-based methods in geodynamics, the theoretical and practical challenges of deriving, implementing and validating adjoint systems for large-scale, non-linear, time-dependent problems, such as global mantle flow, has hindered their broader use. Here, we present the Geoscientific ADjoint Optimisation PlaTform (G-ADOPT), an advanced computational modelling framework that overcomes these challenges for coupled, non-linear, time-dependent systems by integrating three main components: (i) Firedrake, an automated system for the solution of partial differential equations using the finite-element method; (ii) Dolfin-Adjoint, which automatically generates discrete adjoint models in a form compatible with Firedrake; and (iii) the Rapid Optimisation Library, ROL, an efficient large-scale optimisation toolkit; G-ADOPT enables the application of adjoint methods across geophysical continua, showcased herein for geodynamics. Through two sets of synthetic experiments, we demonstrate the application of this framework to the initial condition problem of mantle convection, in both square and annular geometries, for both isoviscous and non-linear rheologies. We confirm the validity of the gradient computations underpinning the adjoint approach, for all cases, through second-order Taylor remainder convergence tests and subsequently demonstrate excellent recovery of the unknown initial conditions. Moreover, we show that the framework achieves theoretical computational efficiency. Taken together, this confirms the suitability of G-ADOPT for reconstructing the evolution of Earth's mantle in space and time. The framework overcomes the significant theoretical and practical challenges of generating adjoint models and will allow the community to move from idealised forward models to data-driven simulations that rigorously account for observational constraints and their uncertainties using an inverse approach.

Macroscopic finite-difference scheme and modified equationsof the general propagation multiple-relaxation-time lattice Boltzmann model.

In this paper we first present the general propagation multiple-relaxation-time lattice Boltzmann (GPMRT-LB) model and obtain the corresponding macroscopic finite-difference (GPMFD) scheme on conservative moments. Then based on the Maxwell iteration method, we conduct the analysis on the truncation errors and modified equations(MEs) of the GPMRT-LB model and GPMFD scheme at both diffusive and acoustic scalings. For the nonlinear anisotropic convection-diffusion equation(NACDE) and Navier-Stokes equations(NSEs), we also derive the first- and second-order MEs of the GPMRT-LB model and GPMFD scheme. In particular, for the one-dimensional convection-diffusion equation(CDE) with the constant velocity and diffusion coefficient, we can develop a fourth-order GPMRT-LB (F-GPMRT-LB) model and the corresponding fourth-order GPMFD (F-GPMFD) scheme at the diffusive scaling. Finally, three benchmark problems, the Gauss hill problem, the CDE with nonlinear convection and diffusion terms, and the Taylor-Green vortex flow in two-dimensional space, are used to test the GPMRT-LB model and GPMFD scheme, and it is found that the numerical results not only are in good agreement with corresponding analytical solutions, but also have a second-order convergence rate in space. Additionally, a numerical study on one-dimensional CDE also demonstrates that the F-GPMRT-LB model and F-GPMFD scheme can achieve a fourth-order accuracy in space, which is consistent with our theoretical analysis.

A fitted parameter convergent finite difference scheme for two-parameter singularly perturbed parabolic differential equation

The objective of this paper is to develop a numerical scheme that is uniform in its parameters for a specific type of time-dependent parabolic problem with two perturbation parameters. The existence of these two parameters in the terms with the highest-order derivatives results in the formation of boundary layer(s) in the solution of such problems. Solving these model problems using classical methods does not yield satisfactory results due to the layer behavior. Therefore, nonstandard finite difference schemes have been developed as a means to obtain numerical solutions for these problems. To develop the scheme, we employ the Crank-Nicolsondiscretization on a uniform time mesh and apply a fitted operator method with a uniform spatial mesh. We have established the stability and convergence of the proposed scheme. The proposed scheme exhibits uniform convergence of second order in the temporal direction and first order in the spatial direction. However, temporal mesh refinements is employed to enhance the order to two in both directions. Model examples are provided to validate the practicality of the proposed numerical scheme.

Improvements on the discretisation of boundary conditions to the momentum balance for glacial ice

Abstract The flow of glacial ice is typically approximated as a nonNewtonian viscous fluid, with the momentum balance described by (an approximation to) the Stokes equations, and the nonlinear rheology described by a flow law. The most commonly used rheological law for glacial ice, Glen's flow law, yields infinite viscosity in the case of zero deformation, which can be the case at the ice surface. This poses a problem when solving the momentum balance numerically. We show that two commonly-used discretisation schemes for the boundary conditions at the ice surface and base, which yield proper numerical convergence when applied to simpler problems, produce poor numerical convergence and large errors, when used to solve the momentum balance with Glen's flow law. We show that a discretisation scheme based on the concept of ghost nodes, which substitutes the boundary conditions directly into the momentum balance equations, yields second-order numerical convergence and errors that can be up to four orders of magnitude smaller. These results are robust across different momentum balance approximations. We show that the improved boundary conditions are particularly useful for solving the 3-D higher-order Blatter-Pattyn Approximation (BPA). In general, this work underlines the importance of thoroughly verifying the numerical solvers used in ice-sheet models, before applying them to future projections of ice-sheet mass loss.

A fast solvable operator-splitting scheme for time-dependent advection diffusion equation

It is well known that the implicit central difference discretization for unsteady advection diffusion equation (ADE) suffers from being time-consuming to solve when the advection term dominates. In this paper, we propose an operator-splitting scheme for the unsteady ADE, in which the ADE is firstly discretized by Crank-Nicolson (CN) scheme in time and central difference scheme in space; and then the discrete advection-diffusion problem is split as advection sub-problem and diffusion sub-problem at each time-level. The significance of the new scheme is that these sub-problems can be fast and directly solved within a linearithmic complexity (a linear-times-logarithm complexity) by means of fast sine transforms (FSTs). In particular, the complexity is independent of the dominance of the advection term. Theoretically, we show that proposed scheme is unconditionally stable and of second-order convergence in time and space. Numerical results are reported to show the efficiency of the proposed scheme.

Analysis and simulation of the (2 + 1)-dimensional Fisher’s reaction-diffusion equation with higher order finite element method

In this article, a finite element analysis of Fisher’s reaction-diffusion equation is performed, by using quadratic shape elements. The higher-order elements in the finite element method (FEM) are known to yield better results in the approximation of solutions. However, such elements are less explored in the literature. Therefore, quadratic shape functions are taken into consideration for the finite element analysis of Fisher’s reaction-diffusion equation. A priori error estimates for the semi-discrete solution are derived by applying the Galerkin approximation. For the fully discrete solution, a priori error estimates are derived using the Crank-Nicolson method and the predictor-corrector scheme. Also, the stability of the time-discrete scheme is analyzed by the energy technique. Third-order convergence for space discretization in L 2 ( Ω ) norm and second-order convergence in H 1 ( Ω ) norm are observed. Furthermore, second-order convergence for time discretization is also found. Finally, we validate our theoretical results by considering a variety of numerical examples in one and two-dimensional spaces, which also illustrate the reliability of the developed finite element algorithm. Through comparison tables, 3D plots, and convergence plots, the results of the current study are compared with the results from the existing literature.

Temporal error analysis of an unconditionally energy stable second-order BDF scheme for the square phase-field crystal model

In this paper, we first propose and study the second-order time-discrete numerical scheme for the sixth-order nonlinear parabolic problem of the square phase-field crystal model. Then, we demonstrate the two-step backward differentiation formula (BDF-2) scheme with mass conservation and energy dissipation, where the higher order nonlinear term is treated implicitly. Moreover, a rigorous error analysis is presented and we prove the optimal second-order convergence rate O(τ2) in H1- norm, where τ is the time step. Finally, some numerical results are provided to confirm our theoretical analysis.