Source-linked AI summary
A Finite Difference Method with Non-uniform Timesteps for Fractional Diffusion Equations
Santos B. Yuste, Joaquín Quintana-Murillo
TL;DR
Useful analytic solutions for fractional diffusion equations are often unavailable, and their finite difference computations can be costly because each step uses prior-time values. The paper develops an implicit Caputo-form method with non-uniform timesteps; it is unconditionally stable and convergent, and its adaptive use is illustrated on subdiffusive-particle dispersion.
Problem
Useful analytic solutions are often unavailable, while fractional finite difference computations scale as n^2 because they require solution values from previous times.
Source: §1 · p. 1Method
The method discretizes the Caputo derivative with a non-uniform-mesh L1 generalization and the Laplacian with a centered formula, allowing timesteps to adapt to solution behavior.
Source: §1 · p. 2, §6 · p. 10Results
The method is unconditionally stable regardless of the non-uniform temporal mesh, and the paper establishes convergence from consistency and stability under relatively weak conditions.
Source: §1 · p. 2Takeaways & Limitations
Adaptive timesteps sample regions of large solution variation more thoroughly and regions of mild variation more sparsely, reducing the timesteps needed in computation.
Source: §6 · p. 8The analyzed spatial mesh has constant spacing, and the alternative explicit method is unstable when timesteps are not small enough for the given parameters.
Source: §2 · pp. 3–4
Abstract
from arXiv · showhide
An implicit finite difference method with non-uniform timesteps for solving the fractional diffusion equation in the Caputo form is proposed. The method allows one to build adaptive methods where the size of the timesteps is adjusted to the behaviour of the solution in order to keep the numerical errors small without the penalty of a huge computational cost. The method is unconditionally stable and convergent. In fact, it is shown that consistency and stability implies convergence for a rather general class of fractional finite difference methods to which the present method belongs. The huge computational advantage of adaptive methods against fixed step methods for fractional diffusion equations is illustrated by solving the problem of the dispersion of a flux of subdiffusive particles stemming from a point source.
1. Introduction
Fractional diffusion equations model subdiffusion, but useful analytic solutions are often unavailable, making numerical methods important. Variable timesteps can support adaptive sampling, while the paper presents a Caputo-form finite difference scheme and establishes stability and convergence.
1. Introduction
- Analytic solution methods exist, but useful solutions are often unavailable, so numerical methods are needed for many problems.Finite difference methods are among the numerical approaches studied for fractional differential equations.Source: §1 · p. 1
- Variable timesteps support adaptive methods that sample rapidly changing solution regions more finely and use larger steps where change is mild.The paper presents its method as a basis for extending adaptive methods to fractional diffusion equations.Source: §1 · pp. 1–2
- Fractional finite difference calculations require values from all previous times, so operation counts scale as n^2 and large computations become slow and memory-demanding.The passage attributes this cost to the non-locality of fractional derivatives.Source: §1 · p. 2
- Unconditionally stable for any non-uniform temporal mesh, the method is also convergent under the paper’s consistency and stability argument.The convergence implication is shown in a rather general form for fractional finite difference methods.Source: §1 · p. 2
- The scheme discretizes the Caputo derivative on non-uniform meshes and reduces to the schemes of Liu et al. and Murio on uniform temporal meshes.The paper selects it for its suitability for variable timesteps and for handling discontinuous solutions.Source: §1 · pp. 2–3
2. Fractional difference algorithm with non-uniform timesteps
The algorithm combines a centered spatial Laplacian with a non-uniform-mesh L1 generalization for the Caputo derivative. Its implicit system uses prior-time solution values as memory, while an alternative explicit scheme is not pursued because it can be unstable.
2. Fractional difference algorithm with non-uniform timesteps
- The spatial Laplacian is discretized with the three-point centered formula on a constant spatial mesh.The passage introduces this as the spatial discretization used in the algorithm.Source: §2 · p. 3
- The Caputo derivative is discretized using a generalization of the L1 formula for non-uniform meshes.The method’s truncation error tends to zero as the time-space mesh spacing tends to zero, establishing consistency.Source: §2 · p. 3
- The implicit method obtains each new solution by solving a tridiagonal system, which can be handled efficiently with the Thomas algorithm.The passage states that the matrix is strictly diagonally dominant.Source: §2 · p. 4
- The memory operator at time n depends on solution values at all previous times.The method’s equation can be written in vector form, with this operator carrying the non-local memory of the fractional derivative.Source: §2 · p. 4
- The paper does not explore the explicit scheme because it is unstable when timesteps are not sufficiently small for the chosen parameters.The relevant parameters are γ, Δx, and K.Source: §2 · p. 4
3. Stability
The method’s stability is assessed by tracking perturbations through a von Neumann–Fourier analysis of individual modes. The resulting bound establishes unconditional stability for non-uniform temporal meshes.
3. Stability
- The perturbation evolves under the same linear difference equation as the solution when the source term is zero.It is defined as the difference between a perturbed solution and the reference solution.Source: §3 · p. 5
- The analysis expands the perturbation in discrete Fourier modes and checks the evolution of a generic mode.Stability of the full solution follows because it is a superposition of modes.Source: §3 · p. 5
- The method is unconditionally stable, with each Fourier-mode perturbation remaining no larger than its initial value.The bound applies for all time steps on the temporal mesh.Source: §3 · p. 5
- The stability condition bounds the norm of the full perturbation by its initial norm.The passage connects the mode-wise stability result to the Euclidean norm using Parseval’s relation.Source: §3 · p. 6
4. Consistency and stability implies convergence
For a broad class of fractional difference algorithms, consistency and stability together imply convergence. The paper establishes that its own method meets these conditions.
4. Consistency and stability implies convergence
- For a class of fractional difference algorithms, consistency and stability imply convergence as the spatiotemporal mesh is refined.The result extends a sufficient-condition argument associated with Lax’s equivalence theorem to fractional equations.Source: §4 · p. 6
- The paper’s method is convergent because it is consistent and A^-1 is bounded.A is a real symmetric tridiagonal Toeplitz matrix whose eigenvalues are all greater than 1.Source: §4 · p. 6
- The proof requires only generic properties of the operators M and A, including linearity, boundedness, and stability.The authors relate this proof strategy to the standard procedure for non-fractional diffusion problems.Source: §4 · p. 6
5. Dispersion of a flux of subdiffusive particles: Numerical results
The method is tested on subdiffusive particles released at a point source, with timestep size adjusted using solution curvature. The reported results show good agreement with exact solutions and substantially fewer steps than a fixed-step calculation.
5. Dispersion of a flux of subdiffusive particles: Numerical results
- The test models a unit flux of particles released at x = 0 at integer times in a one-dimensional infinite-medium problem.The computation uses K = 1, γ = 1/2, a domain length L = 20, and 100 spatial intervals of width 0.2.Source: §5 · pp. 6–7
- The algorithm handles each point-source discontinuity by updating the source node when a release occurs.The authors state that their scheme can handle discontinuities at discretization times, although the earlier analysis considers continuous solutions.Source: §5 · pp. 7–8
- The numerical and exact solutions agree well, though errors are relatively large shortly after each Dirac delta release.The authors attribute the short-time error to the crude approximation of the delta function.Source: §5 · p. 8
- 64 adaptive steps, versus 1,000 fixed steps, reach t = 1; the adaptive calculation’s CPU time is well above 100 times shorter.The comparison uses a fixed timestep of 0.001 and the adaptive rule in equation (41).Source: §5 · p. 8
- The adaptive rule varies the timestep with curvature at x = 0, using smaller steps as curvature increases.The rule sets minimum and maximum timestep values of 0.0001 and 0.02; it uses 5 steps from t = 0 to 0.01 and 59 more to reach t = 1.Source: §5 · p. 8
- The curvature estimate is used as a local proxy for the fractional time derivative because the latter is non-local in time.The authors describe the Laplacian estimate as an easier way to gauge the fractional solution’s temporal change.Source: §5 · p. 8
6. Conclusions and final remarks
The paper presents an implicit fractional-diffusion method with adaptive non-uniform timesteps, sampling finely where the solution changes rapidly and more sparsely where it varies mildly. It reports unconditional stability and convergence, and identifies reduced timestep counts as especially useful when computational cost grows quadratically with the number of timesteps.
6. Conclusions and final remarks
- Adaptive timesteps sample rapidly changing regions more thoroughly and mildly varying regions more sparsely, aiming to maintain accuracy while reducing timesteps.The figures compare exact and numerical solutions and show the error for fixed and variable timesteps, but the supplied passages do not state a numerical winner.Source: §6 · p. 8, Figure 1 · p. 9, Figure 2 · p. 9
- Fractional-problem computations scale quadratically with timestep count, making reductions in the number of timesteps especially useful.The paper contrasts this with the typical linear scaling for non-fractional diffusion problems.Source: §6 · pp. 8, 10
- The method is unconditionally stable, and consistency plus stability implies convergence under the stated general conditions.The consistency-and-stability implication is presented as a result for a broader class of fractional finite difference schemes, not only this method.Source: §6 · p. 10
- Possible extensions include higher-order non-uniform fractional-derivative discretizations, non-uniform spatial meshes, adaptive timestep-selection criteria, and other fractional equations.Examples named include diffusion-wave, Fokker–Planck, and Riemann–Liouville equations.Source: §6 · p. 10
Appendix A. Non-uniform L1 formula for the Caputo derivative
The appendix derives a non-uniform L1 formula for the Caputo derivative and discusses its temporal truncation error. It also gives a bound relating the accumulated error to the largest local error.
Appendix A. Non-uniform L1 formula for the Caputo derivative
- The derivation starts by rewriting the Caputo derivative and follows the L1 approach to approximate the ordinary derivative within the integral.The passage identifies the formula as a non-uniform-mesh generalization of the L1 method.Source: Appendix A · p. 11
- At mesh time t_m, the derivation distinguishes the left- and right-hand values of y(t).Source: Appendix A · p. 11
- The accumulated error R(t_n) is bounded by a constant times the largest local error among earlier steps.The constant is stated to be of order t^(1−γ).Source: Appendix A · p. 11
- The stated condition for the auxiliary function is that it is increasing while its derivative is decreasing; f(z)=z^(1−γ) is given as an example.Source: Appendix A · p. 11