An α-robust finite difference method for a time-fractional radially symmetric diffusion problem

https://doi.org/10.1016/j.camwa.2021.06.010Get rights and content

Abstract

A time-fractional diffusion problem is considered on a spherically symmetric domain in Rd for d=1,2,3. The solution of such a problem is shown in general to have a weak singularity near the initial time t=0, and bounds on certain derivatives of this solution are derived. To solve the problem numerically, spatial polar coordinates are used; a finite difference method is constructed on a mesh that is graded in time and spherical in space. The discretisation uses the L1 scheme in time and a polar-coordinate discretisation of the diffusion term. Its convergence is analysed and error bounds are derived that are robust in α, the order of the time derivative, as α1. Numerical experiments show that our results are sharp.

Introduction

Time-fractional diffusion equations arise in various physical processes and have been widely studied numerically; see [4] for a survey of recent papers on this topic. The case where the problem is posed on a spherical domain and has a radially-symmetric solution has received little attention from numerical analysts, although such problems are of interest in thermoelastics [9], [10], [11]; see also [3].

Let ΩRd (for d=1,2,3) be the ball of radius R>0 with centre at the origin. In this paper, we consider a time-fractional diffusion problem posed on Ω whose data and solution are radially symmetric. The numerical analysis of such a problem is considered in [6], but under the restrictive assumption that the solution is smooth in both time and space; in contrast, our numerical method and analysis deal with the weak singularity that appears at the initial time t=0 in typical solutions of this class of problem.

Since the problem is radially symmetric, its solution can be written as u=u(r,t), where r is the distance from the centre of Ω (so 0rR) and t is the time coordinate, with 0tT for some fixed T. We shall investigate the following time-fractional radially symmetric anomalous diffusion problem:Dtαu(r,t)=kαrd1r(rd1u(r,t)r)+f(r,t)for (r,t)(0,R)×(0,T],u(R,t)=0for t(0,T],u(r,0)=g(r)for r(0,R). Here α(0,1) is the anomalous diffusion exponent, kα>0 is the diffusion coefficient which depends on the physical properties of the system, and f(r,t) is the external source term, which we note is radially symmetric. The fractional time derivative Dtαu(r,t) is the Caputo derivative of order α, which is defined [2] byDtαu(r,t)=1Γ(1α)0tsu(r,s)(ts)αdsfor 0r<R and 0<tT.

The spatial differential operator ur(rd1u(r,t)r) of (1.1a) is the standard expression of the Cartesian-coordinate Laplacian Δu in polar coordinates (r,θ) if d=2 or spherical polar coordinates (r,θ,ϕ) if d=3, where the condition of radial symmetry is imposed, as then partial derivatives with respect to θ and ϕ become zero and can be removed.

Existence and uniqueness of the solution u(r,t) of (1.1a), (1.1b), (1.1c) will be discussed in Section 2, where we use Cartesian coordinates so that we can appeal to previous results for time-fractional problems posed on smooth domains. On the other hand, polar/spherical coordinates are more convenient for discretisation since they handle effortlessly the curved boundary of Ω, so in Section 3 we revert to the (r,t) formulation of (1.1a), (1.1b), (1.1c) to define a finite difference scheme in space. For the time discretisation we use the well-known L1 scheme on a graded mesh. The error analysis of our method is carried out in Section 4, where an α-robust error estimate is derived (this means that the generic constant in the estimate remains bounded if α1, which is in accord with the behaviour of the true solution of (1.1a), (1.1b), (1.1c); see [1] for further discussion of this desirable property). Finally, in Section 5 we give numerical results to demonstrate the sharpness of the error bounds for our finite difference scheme. Conclusions are given in Section 6.

Notation. In this paper C denotes a generic constant that depends on the data of the boundary value problem (1.1a) but is independent of any mesh used to solve (1.1a) numerically; note that C can take different values in different places.

The space of real-valued functions that are absolutely continuous on [0,T] is denoted by AC[0,T].

The L2(Ω) inner product is written as (,), and Wm,p(Ω) is the standard Sobolev space. We use H2 to denote the norm of W2,2(Ω)=H2(Ω); furthermore, denotes the norm of W0,2(Ω)=L2(Ω) and denotes the norm of L(Ω).

Section snippets

Regularity of the solution

To discuss existence, uniqueness and regularity of the solution of our time-fractional problem (1.1a), (1.1b), (1.1c), it is easier to work in Cartesian coordinates and treat Ω as a general smooth domain, since problems of this type have already been investigated by several authors. By standard transformations from polar coordinates (when d=2) or spherical polar coordinates (when d=3) to Cartesian coordinates, the problem (1.1a), (1.1b), (1.1c) can be written asDtαu˜(x,t)=kαΔu˜(x,t)+f˜(x,t)for 

The difference scheme

This section describes the construction of the finite difference scheme for our problem (1.1a), (1.1b), (1.1c).

Section 2 leads us to expect the unknown solution to have a weak singularity at t=0. To address this singularity we use a graded temporal mesh with N+1 points: tn:=T(n/N)ρ for n=0,1,,N, where ρ1 is the user-chosen mesh grading exponent and N is any positive integer. Similar meshes have been used by many authors when discretizing time-fractional problems. Set τn:=tntn1 for n=1,2,N.

Error analysis

In this section we prove an error bound for the solution of our difference scheme (3.21a), (3.21b), (3.21c).

Numerical experiments

In this section, we use our difference scheme to solve two problems of the form (1.1a), (1.1b), (1.1c) which are constructed such that their exact solutions are known and display typical singularities. These results show the scheme is 2nd-order accurate in space, while in time it is almost O(Nα) accurate on a uniform mesh and almost O(N(2α)) accurate on a suitable graded mesh.

Define the maximum error for t=t0,t1,,tN byError(N,J):=max1nNmax1jJ|ujnu(rj,tn)|, the temporal convergence

Conclusions

In this paper, we analysed the regularity of a time-fractional diffusion equation whose spatial domain is a circle or a ball in Rd with d=2 or 3. The solution is seen to have a weak singularity at t = 0. Using polar coordinates and a graded temporal mesh, a finite difference method for solving the problem is constructed and analysed. An α-robust error estimate is derived. Numerical experiments illustrate the sharpness of our theoretical error bounds.

References (15)

There are more references available in the full text version of this article.

Cited by (4)

1

The research of this author is supported in part by the National Natural Science Foundation of China under grant NSAF-U1930402.

View full text