Two Time-Stepping Schemes for Sub-Diffusion Equations with Singular Source Terms

The subdiffusion equation with a Caputo fractional derivative of order � 2 (0, 1) in time arises in a wide variety of practical applications, and it is often adopted to model anomalous subdiffusion processes in heterogeneous media. The L1 scheme is one of the most popular and successful numerical methods for dis- cretizing the Caputo fractional derivative in time. The scheme was analyzed earlier independently by Lin and Xu (2007) and Sun and Wu (2006), and an O(� 2−� ) convergence rate was established, under the assumption that the solution is twice continuously differentiable in time. However, in view of the smoothing property of the subdiffusion equation, this regularity condition is restrictive, since it does not hold even for the homogeneous problem with a smooth initial data. In this work, we revisit the error analysis of the scheme, and establish an O(�) convergence rate for both smooth and nonsmooth initial data. The analysis is valid for more gen- eral sectorial operators. In particular, the L1 scheme is applied to one-dimensional space-time fractional diffusion equations, which involves also a Riemann-Liouville derivative of order � 2 (3/2,2) in space, and error estimates are provided for the fully discrete scheme. Numerical experiments are provided to verify the sharpness of the error estimates, and robustness of the scheme with respect to data regularity.

On many occasions the Navier-Stokes equations are integrated over one or more space directions to remove small-scale features from the solution, such as the size of individual rotor blades in large wind farms, or the micro-pores in flow through porous media. The result of these manipulations is the appearance of additional terms on the right-hand side of the conservation laws. These additional terms may exhibit discontinuous or singular behavior, even if the underlying behavior they represent is regular. Naïve treatment of these additional terms results in non-physical features in the solution or inaccurate predictions of mass, momentum, or energy transfer. This paper proposes a rigorous and straightforward approach for accurately handling source terms with singularities or discontinuities for systems governed by conservation laws, discretized using finite volume schemes with Godunov-type fluxes. We distinguish between volume and surface source terms. Both types can be handled by carefully constructed upwind procedures. This discretization exhibits second-order accuracy, for formally first-order numerical models. The method is validated using the Euler equations with large source terms, with both preconditioned and non-preconditioned flux functions. We also present a novel, elegant method for introducing jump conditions in the solution via singular source terms. This ensures global compatibility relations are met without the need for internal boundaries, as commonly done in actuator disk models. The robustness of the method is demonstrated on a Riemann problem with abrupt area change, matching analytical solutions for 10 different wave patterns. Finally, the new methods are applied to a multi-dimensional, open rotor problem and compared to other results in the literature.

A reduced-dimension method for unknown Crank-Nicolson finite element solution coefficient vectors of elastic wave equation with singular source term

ABSTRACTIn this paper, a second-order backward differentiation formula (BDF) orthogonal spline collocation (OSC) method with the truncation error of order in time and in space is presented for a kind of partial differential equation (PDE) with memory. The stability and convergence of the BDF OSC method in a new norm are proved. Besides, we provide three numerical experiments to demonstrate the results of theoretical analysis and show the accuracy and effectiveness of the BDF OSC method. Numerical experiments also exhibit the optimal error estimates in the , and-norms.

In this paper, a backward Euler method combined with finite element discretization in spatial direction is discussed for the equations of motion arising in the two-dimensional Oldroyd model of viscoelastic fluids of order one with the forcing term independent of time or in $L^{\infty }$ in time. It is shown that the estimates of the discrete solution in Dirichlet norm is bounded uniformly in time. Optimal a priori error estimate in $\textbf {L}^2$-norm is derived for the discrete problem with nonsmooth initial data. This estimate is shown to be uniform in time, under the assumption of uniqueness condition. Finally, we present some numerical results to validate our theoretical results.

We study the space-time finite element discretizations for time fractional parabolic integro-differential equations in a bounded convex polygonal domain in $${\mathbb {R}}^d (d=1,2,3)$$ . Both spatially semidiscrete and fully discrete finite element approximations are considered and analyzed. We use piecewise linear and continuous finite elements to approximate the space variable whereas the time discretization uses two fully discrete schemes based on the convolution quadrature, namely the backward Euler and the second-order backward difference. For the spatially discrete scheme, optimal order a priori error estimates are derived for smooth initial data, i.e., when $$u_0\in H_0^1(\varOmega )\cap H^2(\varOmega )$$ . Moreover, for the homogeneous problem, almost optimal error estimates for positive time are established for nonsmooth initial data, i.e., when the initial function $$u_0$$ is only in $$ L^2(\varOmega )$$ . The error estimates for the fully discrete methods are shown to be optimal in time for both smooth and nonsmooth initial data under the specific choice of the kernel operator in the integral. Finally, we provide some numerical illustrations to verify our theoretical analysis.

Fully discrete approximation methods are studied for parabolic boundary value problems evolving from nonsmooth (and noncompatible) initial data. Galerkin methods are used for the spatial discretization, and single step or backwards differentiation multistep methods are used in time. The parabolic equation is assumed to be linear, but the coefficients are allowed to depend on space and time. An error estimate is given for the combined space-time discretization scheme which contains the expected optional order space and time error contributions. However, another term appears which limits the accuracy of the method to second order (or fractionally better) in time if the time stepping method is of third or higher order. This latter term arises from the time dependence of the equation’s highest order coefficients, and hence does not appear in analyses of problems with only spatially dependent coefficients for which optimal order estimates are known to hold. It is shown through more detailed analysis of a specific problem, and the formulation of a lower error bound, that the error estimate is indeed sharp. Thus, it follows that standard higher order time stepping schemes cannot be expected to perform well for parabolic problems with time dependent coefficients in situations when nonsmooth initial data are encountered.

Implicit methods are well known to have greater stability than explicit methods for stiff systems, but they often are not used in practice due to perceived computational complexity. This paper applies the backward Euler (BE) method and a second-order one-step two-stage composite backward differentiation formula (C-BDF2) for the monodomain equations arising from mathematically modeling the electrical activity of the heart. The C-BDF2 scheme is an L-stable implicit time integration method and easily implementable. It uses the simplest forward Euler and BE methods as fundamental building blocks. The nonlinear system resulting from application of the BE method for the monodomain equations is solved for the first time by a nonlinear elimination method, which eliminates local and nonsymmetric components by using a Jacobian-free Newton solver, called Newton--Krylov solver. Unlike other fully implicit methods proposed for the monodomain equations in the literature, the Jacobian of the global system after the nonlinear elimination has much smaller size, is symmetric and possibly positive definite, which can be solved efficiently by standard optimal solvers. Numerical results are presented demonstrating that the C-BDF2 scheme can yield accurate results with less CPU times than explicit methods for both a single patch and spatially extended domains.

Using the weighted and shifted Lubich difference (WSLD) operator, we propose an efficient and stable difference scheme for solving space tempered fractional diffusion equations. The scheme has fourth order accuracy in space and second order accuracy in time. To improve the computational efficiency further, a fourth order backward Euler scheme is used to achieve fourth order accuracy in time. The fourth order backward Euler method is shown to be stable under a minor requirement. For two-dimensional problems, the Crank--Nicolson alternating direction implicit (CN-ADI) scheme is used to reduce computational complexity. These proposed methods are applied to the Carr--Geman--Madan--Yor (CGMY) model in finance, and several numerical experiments are designed to verify the efficiency of the proposed methods.

In this paper, a finite volume element (FVE) method is considered for spatial approximations of time fractional diffusion equations involving a Riemann-Liouville fractional derivative of order α ∈ (0, 1) in time. Improving upon earlier results [Karaa et al., IMA J. Numer. Anal. 37 (2017) 945–964], error estimates in L2 (Ω)- and H1 (Ω)-norms for the semidiscrete problem with smooth and mildly smooth initial data, i.e., v ∈ H2(Ω) ∩ H01(Ω) and v ∈ H01(Ω) are established. For nonsmooth data, that is, v ∈ L2 (Ω), the optimal L2 (Ω)-error estimate is shown to hold only under an additional assumption on the triangulation, which is known to be satisfied for symmetric triangulations. Superconvergence result is also proved and as a consequence, a quasi-optimal error estimate is established in the L∞ (Ω)-norm. Further, two fully discrete schemes using convolution quadrature in time generated by the backward Euler and the second-order backward difference methods are analyzed, and error estimates are derived for both smooth and nonsmooth initial data. Based on a comparison of the standard Galerkin finite element solution with the FVE solution and exploiting tools for Laplace transforms with semigroup type properties of the FVE solution operator, our analysis is then extended in a unified manner to several time fractional order evolution problems. Finally, several numerical experiments are conducted to confirm our theoretical findings.

We propose a Legendre–Laguerre spectral approximation to price the European and double barrier options in the time-fractional framework. By choosing an appropriate basis function, the spectral discretization is used for the approximation of the spatial derivatives of the time-fractional Black–Scholes equation. For the time discretization, we consider the popular \(L1\) finite difference approximation, which converges with order \(\mathcal{O}((\Delta \tau)^{2-\alpha})\) for functions which are twice continuously differentiable. However, when using the \(L1\) scheme for problems with nonsmooth initial data, only the first-order accuracy in time is achieved. This low-order accuracy is also observed when solving the time-fractional Black–Scholes European and barrier option pricing problems for which the payoffs are all nonsmooth. To increase the temporal convergence rate, we therefore consider a Richardson extrapolation method, which when combined with the spectral approximation in space, exhibits higher order convergence such that high accuracies over the whole discretization grid are obtained. Compared with the traditional finite difference scheme, numerical examples clearly indicate that the spectral approximation converges exponentially over a small number of grid points. Also, as demonstrated, such high accuracies can be achieved in much fewer time steps using the extrapolation approach. doi:10.1017/S1446181121000286

We derive rigorous error bounds for exponential Runge–Kutta discretizations of parabolic equations with nonsmooth initial data. Our analysis is carried out in the framework of abstract semilinear evolution equations, allowing for operators with non-dense domains. The results provide a foundation for establishing error estimates for nonsmooth data in prototypical problems such as the Allen-Cahn and Burgers equations. Furthermore, we apply these estimates to the analysis of split exponential integrators, yielding convergence results expressed explicitly in terms of the prescribed data.

Physics-Informed Neural Networks (PINNs) have emerged as a promising approach for solving acoustic seismic wave equations in geophysical modeling. However, they often encounter challenges such as spectral bias, gradient instability, and slow convergence when applied to complex media, high-frequency wavefields, and singular source terms. To address these limitations, we propose Adaptive SIREN-PINNs, an enhanced framework that integrates learnable amplitude and frequency modulation parameters into the Sinusoidal Representation Networks (SIRENs), which employs periodic activation functions to naturally represent continuous signals and their derivatives, along with a principled weight initialization strategy. This enhancement improves the model’s frequency-aware representation capability, enabling it to better capture both low- and high-frequency components while enhancing gradient propagation stability during training. Furthermore, to mitigate the singularity associated with point sources, we adopt a Gaussian approximation of the Dirac function and reformulate the wave equation in terms of the scattered wavefield, following established regularization practices in seismic wave modeling. Additionally, we conduct a spectral analysis of the Neural Tangent Kernel (NTK) to gain insights into the model’s learning dynamics across different frequency components during training. Numerical experiments conducted on isotropic homogeneous medium, the Marmousi model, and anisotropic TTI/VTI medium demonstrate that the proposed method significantly outperforms conventional PINNs and baseline SIREN-PINNs in terms of wavefield reconstruction accuracy and convergence efficiency. The results validate the effectiveness of the adaptive SIREN mechanism in improving both the expressive power and training stability of PINNs, offering a mesh-free and solver-free alternative to traditional numerical methods.

The Poisson–Boltzmann equation (PBE) is a fundamental implicit solvent continuum model for calculating the electrostatic potential of large ionic solvated biomolecules. However, its numerical solution encounters severe challenges arising from its strong singularity and nonlinearity. In (P. Benner, V. Khoromskaia, B. Khoromskij, C. Kweyu, and M. Stein, “Regularization of Poisson-Boltzmann type equations with singular source terms using the range-separated tensor format,” SIAM J. Sci. Comput., vol. 43, no. 1, pp. A415–A445, 2021; C. Kweyu, V. Khoromskaia, B. Khoromskij, M. Stein, and P. Benner, “Solution decomposition for the nonlinear Poisson-Boltzmann equation using the range-separated tensor format,” arXiv:2109.14073, 2021), the effect of strong singularities was eliminated by applying the range-separated (RS) canonical tensor format (P. Benner, V. Khoromskaia, and B. N. Khoromskij, “Range-separated tensor format for many-particle modeling,” SIAM J. Sci. Comput., vol. 40, no. 2, pp. A1034–A1062, 2018; B. N. Khoromskij, “Range-separated tensor representation of the discretized multidimensional Dirac delta and elliptic operator inverse,” J. Comput. Phys., vol. 401, p. 108998, 2020) to construct a solution decomposition scheme for the PBE. The RS tensor format allows deriving a smooth approximation to the Dirac delta distribution in order to obtain a regularized PBE (RPBE) model. However, solving the RPBE is still computationally demanding due to its high dimension N $\mathcal{N}$ , where N $\mathcal{N}$ is always in the millions. In this study, we propose to apply the reduced basis method (RBM) and the (discrete) empirical interpolation method ((D)EIM) to the RPBE in order to construct a reduced order model (ROM) of low dimension N ≪ N $N\ll \mathcal{N}$ , whose solution accurately approximates the nonlinear RPBE. The long-range potential can be obtained by lifting the ROM solution back to the N $\mathcal{N}$ -space while the short-range potential is directly precomputed analytically, thanks to the RS tensor format. The sum of both provides the total electrostatic potential. The main computational benefit is the avoidance of computing the numerical approximation of the singular electrostatic potential. We demonstrate in the numerical experiments, the accuracy and efficacy of the reduced basis (RB) approximation to the nonlinear RPBE (NRPBE) solution and the corresponding computational savings over the classical nonlinear PBE (NPBE) as well as over the RBM being applied to the classical NPBE.

A diffuse interface–lattice Boltzmann model for conjugate heat transfer with imperfect interface