Convergence in Positive Time for a Finite Difference Method Applied to a Fractional Convection-Diffusion Problem

José Luis Gracia; Eugene O’Riordan; Martin Stynes

doi:10.1515/cmam-2017-0019

Article Publicly Available

Convergence in Positive Time for a Finite Difference Method Applied to a Fractional Convection-Diffusion Problem

José Luis Gracia , Eugene O’Riordan and Martin Stynes

Published/Copyright: July 6, 2017

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Computational Methods in Applied Mathematics Volume 18 Issue 1

Abstract

A standard finite difference method on a uniform mesh is used to solve a time-fractional convection-diffusion initial-boundary value problem. Such problems typically exhibit a mild singularity at the initial time $t=0$ . It is proved that the rate of convergence of the maximum nodal error on any subdomain that is bounded away from $t=0$ is higher than the rate obtained when the maximum nodal error is measured over the entire space-time domain. Numerical results are provided to illustrate the theoretical error bounds.

Keywords: Caputo Fractional Derivative; Initial Boundary Value Problem; Weak Singularity; L1 Scheme

MSC 2010: 65M06; 65M12; 65M15

1 Introduction

In this paper, we examine the convergence rate of numerical approximations to a time-fractional convection-diffusion problem using a standard finite difference method on a uniform mesh. Initial-boundary value problems of this type, where the time derivative is fractional, have solutions that are mildly singular at the initial time $t=0$ ; that is, their temporal derivatives are unbounded on the closed space-time domain, but are bounded on any subdomain that is bounded away from $t=0$ . It is shown in [9, 10] that the case of solutions with bounded temporal derivatives on the closed space-time domain is very special and that the weakly singular solutions examined here are much more typical of how solutions to this class of problems behave. As one would expect, the rate of convergence of the computed numerical approximations is affected adversely by the presence of large temporal derivatives at $t=0$ .

This paper is a companion paper to [10], where it was shown that the convergence rate of the same finite difference scheme on a uniform mesh was $O⁢(N-α)$ , where $α∈(0,1)$ is the order of the fractional derivative and the mesh spacing in time is $O⁢(N-1)$ . Results related to the main result in [10] are available in [7, 8], using a finite element framework. In contrast to [10], we shall prove here for the same scheme on a uniform mesh that the convergence rate of the numerical solution is $O⁢(N-1)$ on any subdomain that is bounded away from $t=0$ .

Our analysis is carried out in the discrete $L∞$ norm; an analogous convergence result in the $L2$ norm was derived in [4]. Using an alternative formulation of the continuous problem, the phenomenon of higher-order convergence at some fixed distance away from the initial singularity is examined in [6] for a homogeneous version of (2.1) in the case of non-smooth initial data.

Notation.

In this paper C denotes a generic constant that depends on the data of the boundary value problem (2.1) but is independent of T and of any mesh used to solve (2.1) numerically. Note that C can take different values in different places. For all $x∈ℝ$ the ceiling function $⌈x⌉$ is the smallest integer greater than or equal to x. For any continuous function $z:Q→ℝ$ with $Q⊂ℝ2$ and any mesh function $zmn$ with $n=0,1,…,N$ and $m=0,1,…,M$ , we set

$∥z∥:=max(x,t)∈Q¯⁡|z⁢(x,t)| and ∥zn∥:=max0≤m≤M⁡|zmn|.$

2 The Continuous Problem

Consider the initial-boundary value problem

(2.1)

(2.1a)

$L⁢u:=Dtα⁢u-p⁢(x)⁢∂2⁡u∂⁡x2+q⁢(x)⁢∂⁡u∂⁡x+r⁢(x)⁢u=f⁢(x,t)$

$for (x,t)∈Q:=(0,l)×(0,T], with initial and boundary conditions$

(2.1b)

$u⁢(0,t)=u⁢(l,t)=0 for ⁢t∈(0,T],$

(2.1c)

$u⁢(x,0)=ϕ⁢(x) for ⁢x∈[0,l].$

$Here 0<α<1, p⁢(x)≥p0>0 on [0,l], the functions p⁢(x),q⁢(x) and r⁢(x) are smooth on [0,l] and are assumed to satisfy the constraint$

(2.1d)

$r⁢(x)-q′⁢(x)+p′′⁢(x)2≥0 for ⁢(x,t)∈Q.$

The initial condition ϕ is also smooth on $[0,l]$ and the function f is smooth on $Q¯$ . Furthermore, in (2.1a) $Dtα$ denotes the Caputo fractional derivative which is defined [1] by

$Dtα⁢g⁢(x,t):=[J1-α⁢(∂⁡g∂⁡t)]⁢(x,t) for ⁢0≤x≤l, 0<t≤T,$

where

$(J1-α⁢g)⁢(x,t):=[1Γ⁢(1-α)⁢∫s=0t(t-s)-α⁢g⁢(x,s)⁢𝑑s] for ⁢(x,t)∈Q¯$

is the Riemann–Liouville fractional integral operator of order $1-α$ .

There is no loss of generality in assuming homogeneous boundary conditions in (2.1b), because inhomogeneous boundary conditions are easily made homogeneous by a simple change of variable.

Under the transformation

$y⁢(x,t):=u⁢(x,t)⁢p⁢(0)p⁢(x)⁢e-12⁢∫0xq⁢(s)p⁢(s)⁢𝑑s,$

problem (2.1) becomes

(2.2)

(2.2a)

$Dtα⁢y-p⁢(x)⁢∂2⁡y∂⁡x2+r1⁢(x)⁢y=f1⁢(x,t)$

$for ⁢(x,t)∈Q:=(0,l)×(0,T],$

(2.2b)

$y⁢(0,t)=y⁢(l,t)=0$

$for ⁢t∈(0,T],$

(2.2c)

$y⁢(x,0)=ϕ1⁢(x):=ϕ⁢(x)⁢p⁢(0)p⁢(x)⁢e-12⁢∫0xq⁢(s)p⁢(s)⁢𝑑s$

$for ⁢x∈[0,l],$

(2.2d)

$where r1⁢(x):=r⁢(x)-(q′⁢(x)+p′′⁢(x))2+(q⁢(x)+p′⁢(x))24⁢p,$

(2.2e)

$f1⁢(x,t):=f⁢(x,t)⁢p⁢(0)p⁢(x)⁢e-12⁢∫0xq⁢(s)p⁢(s)⁢𝑑s.$

Note that no first-order derivative in space appears in (2.2a), and (2.1d) implies that $r1≥0$ . Consequently, (2.2) belongs to the class of problems analysed in [10]. In [10] it was assumed that $p⁢(x)$ was a positive constant, but the analysis of [10] can be extended to the case of a smooth variable positive coefficient $p⁢(x)$ . Thus after placing suitable regularity and compatibility conditions on the data of (2.2), one can invoke [10, Theorem 2.1] to conclude that (2.2) has a unique solution y whose derivatives satisfy certain bounds. Transforming back to the original problem (2.1), under certain conditions on its data one obtains the following bounds on the derivatives of u:

(2.3)

(2.3a)

$|∂k⁡u∂⁡xk⁢(x,t)|≤C$

$⁢for ⁢k=0,1,2,3,4,$

(2.3b)

$|∂ℓ⁡u∂⁡tℓ⁢(x,t)|≤C⁢(1+tα-ℓ)$

$⁢for ⁢ℓ=1,2,$

for all $(x,t)∈[0,l]×(0,T]$ .

In [10, Theorem 2.1] the estimates in (2.3) are proved assuming that $ϕ1∈D⁢(ℒ5/2)$ , $(f1)⁢(⋅,t)∈D⁢(ℒ5/2)$ , $(f1)t⁢(⋅,t)$ and $(f1)t⁢t⁢(⋅,t)$ are in $D⁢(ℒ1/2)$ for each $t∈(0,T]$ and

$∥(f1)⁢(⋅,t)∥ℒ5/2+∥(f1)t⁢(⋅,t)∥ℒ1/2+tρ⁢∥(f1)t⁢t⁢(⋅,t)∥ℒ1/2≤C1$

for all $t∈(0,T]$ and some constant $ρ<1$ , where $C1$ is a constant independent of t and $∥⋅∥ℒγ$ is the norm associated with the vector space $D⁢(ℒγ)$ . This space is defined by

$D⁢(ℒγ):={g∈L2⁢(0,l):∑i=1∞λi2⁢γ⁢|(g,ψi)|2<∞},γ≥0,$

where $(⋅,⋅)$ is the inner product in the Hilbert space $L2⁢(0,l)$ and ${(λi,ψi):i=1,2,…}$ are the eigenvalues and normalised eigenfunctions of the Sturm–Liouville two-point boundary value problem

$ℒ⁢ψi:=-p⁢ψi′′+c⁢ψi=λi⁢ψi on ⁢(0,l),ψi⁢(0)=ψi⁢(l)=0.$

3 The Discrete Problem

The solution of problem (2.1) is approximated by the solution of a finite difference scheme on a mesh ${(xm,tn):m=0,1,…,M,n=0,1,…,N}$ , that is uniform in both space and time. Let M and N be positive integers. Set $h=lM$ and $xm:=m⁢h$ for $m=0,1,…,M$ . Set $tn=n⁢τ=n⁢TN$ for $n=0,1,…,N$ . The nodal approximation to the solution u computed at the mesh point $(xm,tn)$ is denoted by $umn$ .

The first and second-order spatial derivatives are discretised using standard approximations:

$∂⁡u∂⁡x⁢(xm,tn)≈Dx0⁢umn:=um+1n-um-1n2⁢h,$

$∂2⁡u∂⁡x2⁢(xm,tn)≈δx2⁢unm:=um+1n-2⁢umn+um-1nh2.$

The Caputo fractional derivative $Dtα⁢u$ , which can be written as

$Dtα⁢u⁢(xm,tn)=1Γ⁢(1-α)⁢∑k=0n-1∫s=tktk+1(tn-s)-α⁢∂⁡u⁢(xm,s)∂⁡s⁢𝑑s,$

is approximated by the classical L1 approximation

(3.1)

$DNα⁢umn:=1Γ⁢(1-α)⁢∑k=0n-1umk+1-umkτ⁢∫s=tktk+1(tn-s)-α⁢𝑑s$

(3.1a)

$=1Γ⁢(2-α)⁢∑k=0n-1(umk+1-umk)⁢dn-k,$

$where$

(3.1b)

$dk:=k1-α-(k-1)1-α,k≥1.$

Thus, we approximate (2.1) by the discrete problem

(3.2)

(3.2a)

$LM,N⁢umn:=DNα⁢umn-p⁢(xm)⁢δx2⁢umn+q⁢(xm)⁢Dx0⁢umn+r⁢(xm)⁢umn=f⁢(xm,tn)$

$for ⁢1≤m≤M-1, 1≤n≤N,$

(3.2b)

$u0n=uMn=0$

$for ⁢0<n≤N,$

(3.2c)

$um0=ϕ⁢(xm)$

$for ⁢0≤m≤M.$

This discretisation of (2.1) is standard.

To ensure the stability of the discrete operator $LM,N$ by imposing the correct sign pattern in the associated matrix, we make the nonrestrictive assumption that N satisfies

$l⁢∥q∥2⁢p0<N.$

After some minor modifications in the proof of [10, Theorem 5.2] to handle the term $q⁢(xm)⁢Dx0⁢umn$ , it follows that the solution $umn$ of scheme (3.2) satisfies the error bound

(3.3)

$max(xm,tn)∈Q¯⁡|u⁢(xm,tn)-umn|≤C⁢(h2+N-α)$

for some constant C. In particular, the method has the low order of convergence $O⁢(N-α)$ in time when α is small. In the present paper we shall consider the rate of convergence in a subdomain $[0,l]×[κ,T]$ , where κ is a fixed positive value.

4 Error Analysis

The structure of our error analysis is the standard finite difference technique of estimating the truncation error at each mesh point, then invoking a stability argument to derive an error bound for the computed solution $umn$ . In this analysis the truncation error bound (4.2) indicates that the truncation error decreases as one moves further away from the initial time $t=0$ . The stability bound (4.5) shows that the error at any discrete time level depends on a weighted sum of the truncation errors at all the previous time levels.

The estimate of the truncation error in space is standard: using (2.3a), one gets

(4.1)

$∂⁡u∂⁡x⁢(xm,tn)=Dx0⁢u⁢(xm,tn)+O⁢(h2),$

$∂2⁡u∂⁡x2⁢(xm,tn)=δx2⁢u⁢(xm,tn)+O⁢(h2).$

The truncation error in time is more tricky to estimate and this is done in the next lemma.

Lemma 1.

Assume that u satisfies (2.3). Then there exists a positive constant C such that for each mesh point $(xm,tn)∈Q$ one has

(4.2)

$|(DNα-Dtα)⁢u⁢(xm,tn)|≤C⁢n-min⁡{2-α,α+1}.$

Proof.

We modify the argument of [10, Lemma 5.1]. By (3.1a) and the definition of $Dtα$ , for each mesh point $(xm,tn)∈Q$ one has

$(DNα-Dtα)⁢u⁢(xm,tn)=∑k=0n-1Tn⁢k,$

where for $n=1,2,…,N$ and $k=0,1,…,n-1$ we define the truncation error in the kth time cell $[tk,tk+1]$ to be

(4.3)

$Tn⁢k:=1Γ⁢(1-α)⁢∫s=tktk+1(tn-s)-α⁢[u⁢(xm,tk+1)-u⁢(xm,tk)τ-∂⁡u∂⁡s⁢(xm,s)]⁢𝑑s.$

The following four bounds are established in [10, equations (5.9), (5.10), (5.11) and (5.14)]:

(4.4)

(4.4a)

$∑k=1⌈n2⌉-1|Tn⁢k|≤C⁢n-(α+1) for ⁢1≤k<n-1,$

(4.4b)

$∑k=⌈n2⌉n-2|Tn⁢k|≤C⁢n-(2-α) for ⁢1≤k<n-1,$

(4.4c)

$|T10|≤C,$

(4.4d)

$|Tn,n-1|≤C⁢n-(2-α).$

It remains to bound $|Tn⁢0|$ for $n>1$ . We sharpen the bound [10, (5.12)] of this term. An integration by parts in (4.3) yields

$Tn⁢0=-αΓ⁢(1-α)⁢∫s=0t1(tn-s)-α-1⁢(ϕ-ψ)⁢(xm,s)⁢𝑑s,$

where

$ϕ⁢(xm,s):=s⁢[u⁢(xm,t1)-u⁢(xm,0)τ] and ψ⁢(xm,s):=u⁢(xm,s)-u⁢(xm,0).$

$|ϕ⁢(xm,s)|+|ψ⁢(xm,s)|≤2⁢∫0τ|ut⁢(xm,t)|⁢𝑑t≤C⁢∫0τ(1+tα-1)⁢𝑑t≤C⁢τα.$

Hence

$|Tn⁢0|≤C⁢τα⁢∫s=0t1(tn-s)-α-1⁢𝑑s≤C⁢τα⁢[(tn-t1)-α-tn-α]=C⁢[(n-1)-α-n-α]≤C⁢n-(α+1),$

by the Mean Value Theorem. Combine this bound with (4.4) to complete the proof. ∎

Observe that $min⁡{2-α,α+1}>1$ in (4.2) for all values of $α∈(0,1)$ ; thus this bound is sharper than the truncation error bound of $O⁢(n-α)$ proved in [10, Lemma 5.1]. This improvement is critical in establishing our main result later.

Next, we derive some new information about the stability constants that appear in [10, Section 4]. It follows from [10, Lemma 4.2] that the computed solution $umn$ of (3.2) satisfies

(4.5)

$∥u(xm,tn)-umn∥≤ταΓ(2-α)∑j=1nσn-j∥(LM,N(u(xm,tj))j-uj∥$

for $n=1,2,…,N,$ where the positive weights $σi$ are defined for $i=0,1,2,…,n-1$ by the recurrence relation

(4.6)

$σ0:=1,σi:=∑k=1i(dk-dk+1)⁢σi-k for ⁢i=1,2,….$

Note that when the mesh is uniform, the weights $θn,j$ in [10, Lemma 4.2] are the same as the weights $σn-j$ defined in (4.6).

Lemma 2.

The coefficients $σi$ satisfy $σi<(i+1)α-1$ for $i=1,2,…$ .

Proof.

First, $σ1=(d1-d2)⁢σ0=2-21-α<2-1+α$ as $0.5⁢w+2⁢w-1-2>0$ for all $w∈(1,2)$ , so the lemma is true when $i=1$ . The proof is completed by induction. Assume that $σj<(j+1)α-1$ for $j=1,2,…,i-1$ . We want to prove that $σi<(i+1)α-1$ . It is easy to check that $dk-dk+1>0$ for all k. Using this inequality and the inductive hypothesis, we require the inequality

$∑k=1i(dk-dk+1)⁢(i+1-k)α-1<(i+1)α-1,$

which is established in [5, Lemma 3.2]. ∎

The next result, which is a variant of [10, Lemma 4.3], bounds a weighted sum of the $σn-j$ that will be used in the proof of Theorem 4.

Lemma 3.

Let the parameter β satisfy $β>1$ . Then, for $n=1,2,…,N$ , one has

$τα⁢∑j=1nj-β⁢σn-j≤C⁢N-1⁢T⁢tnα-1.$

Proof.

By Lemma 2, we have

(4.7)

$τα⁢∑j=1nj-β⁢σn-j≤τα⁢∑j=1nj-β⁢(n+1-j)α-1≤τα⁢[(n2)α-1⁢∑j=1⌈n2⌉j-β+(n2)-(β-α)⁢∑j=⌈n2⌉+1nj-α⁢(n+1-j)α-1].$

But for $s≥j-1$ and $j≤n$ , one has

$(n+1-s)α-1≥(n+2-j)α-1≥2α-1⁢(n+1-j)α-1.$

Hence,

$∑j=⌈n2⌉+1nj-α⁢(n+1-j)α-1≤∑j=⌈n2⌉+1n(n+1-j)α-1⁢∫s=j-1js-α⁢𝑑s$

$≤∑j=⌈n2⌉+1n21-α⁢∫s=j-1j(n+1-s)α-1⁢s-α⁢𝑑s$

$≤21-α⁢∫s=0n+1(n+1-s)α-1⁢s-α⁢𝑑s$

$=21-α⁢Γ⁢(α)⁢Γ⁢(1-α),$

by [1, Theorem D.6]. Substituting this inequality into (4.7) and using $tn=n⁢τ$ and $β>1$ , we get

$τα⁢∑j=1nj-β⁢σn-j≤C⁢τα⁢nα-1⁢∑j=1⌈n2⌉j-β+C⁢τα⁢n-(β-α)≤C⁢τα⁢nα-1+C⁢τα⁢n-(β-α)=C⁢tnα-1⁢(τ+τ⁢n-(β-1))≤C⁢tnα-1⁢N-1⁢T.$

This completes the proof. ∎

We can now prove our main result.

Theorem 4.

Assume that u satisfies (2.3). Then, for $n=1,2,3,…,N$ , the solution $umn$ of scheme (3.2) satisfies

(4.8)

$max0≤m≤M⁡|u⁢(xm,tn)-umn|≤C⁢(Tα⁢h2+T⁢N-1⁢tnα-1)$

for some constant C.

Proof.

Fix $(xm,tn)∈Q$ . By (4.1) and Lemma 1, the truncation error at $(xm,tn)$ satisfies

$∥(LM,N⁢(u⁢(xm,tn))n-un)∥≤C⁢(h2+n-min⁡{2-α,1+α}).$

By (4.5) we then obtain

$max0≤m≤M⁡|u⁢(xm,tn)-umn|≤C⁢τα⁢∑j=1nh2⁢σn-j+C⁢τα⁢∑j=1nj-min⁡{2-α,1+α}⁢σn-j.$

Invoking Lemma 3 (with $β=min⁡{2-α,1+α}$ ) for the $j-min⁡{2-α,1+α}$ term and [10, Lemma 4.3] (with $β=0$ ) for the term involving $h2$ , we obtain (4.8). ∎

The bound in (4.8) implies that for any fixed $κ>0$ one has

(4.9)

$max(xm,tn)∈Q¯∩{tn≥κ>0}⁡|(xm,tn)-umn|≤C⁢Tα⁢(h2+N-1).$

That is, on any subdomain that is bounded away from $t=0$ , we observe an improved rate of convergence in time compared with the rate of convergence (in time) of $N-α$ on $Q¯$ that is given by (3.3).

5 Numerical Results

In this section we give numerical results for the numerical method (3.2) applied to two particular examples from the problem class (2.1). In the first example the exact solution of the problem is known; in the second example it is unknown, so we estimate the order of convergence using the double-mesh principle [2]. In these numerical experiments we always take $N=M$ . Hence the bounds in (3.3) and (4.8) imply that the spatial error term $C⁢h2$ will be dominated by the temporal error term $C⁢N-α$ or $C⁢N-1$ .

Example 5.1.

Consider the constant coefficient homogeneous problem

$Dtα⁢u-∂2⁡u∂⁡x2+∂⁡u∂⁡x=0 for ⁢(x,t)∈Q=(0,π)×(0,1],$

with initial condition $u⁢(x,0)=ex2⁢sin⁡x$ , $0<x<π$ , and boundary conditions $u⁢(0,t)=u⁢(π,t)=0$ , $0≤t≤1$ . The exact solution of this problem is

$u⁢(x,t)=Eα⁢(-1.25⁢tα)⁢ex2⁢sin⁡x,$

where $Eα$ is the Mittag-Leffler function which is defined [1] by

$Eα⁢(z):=∑k=0∞zkΓ⁢(α⁢k+1).$

In Figure 1 we display the computed solutions with scheme (3.2) for $α=0.4,0.8$ and $N=M=32$ and we observe that the solution has an initial layer at $t=0$ , which becomes sharper as the parameter α decreases.

Figure 1

Example 5.1: Computed solutions with scheme (3.2) for $N=M=32$ .

$(a) α=0.4{\alpha=0.4}$

(a)

$α=0.4$

$(b) α=0.8{\alpha=0.8}$

(b)

$α=0.8$

For Example 5.1 we computed the maximum errors

$eM,N:=max(xm,tn)∈Q′⁡|u⁢(xm,tn)-umn|$

and the orders of convergence

$pM,N:=log2⁡(eM,Ne2⁢M,2⁢N),$

where $Q′$ can be the entire domain $Q¯$ or the subdomain $Q¯*:=[0,π]×[0.1,1]$ . The numerical results in $Q¯$ (see Table 1) show that scheme (3.2) is $O⁢(N-α)$ convergent there (which agrees with [10, Theorem 5.2]), while it is $O⁢(N-1)$ convergent in the subdomain $Q¯*$ (see Table 2), which indicates that the error bound (4.9) is sharp.

Considering the convergence in time, identified by the factor $tnα-1$ in the error bound (4.8), the error $maxm⁡|u⁢(xm,tn)-umn|$ is compared with the error bound $N-1⁢tnα-1$ in Figure 2, for $α=0.2$ and $α=0.8$ . These plots indicate that the exponent $α-1$ in the error bound is sharp for small values of α. However, in the case of larger values of α close to one, the maximum error decreases at a faster rate than $α-1$ , as $tn$ increases.

Table 1

Example 5.1: Maximum errors and orders of convergence for scheme (3.2) in the domain $Q¯$ .

α	$N=M=128$	$N=M=256$	$N=M=512$	$N=M=1024$	$N=M=2048$
0.4	8.438E-2	6.714E-2	5.282E-2	4.120E-2	3.191E-2
	0.330	0.346	0.359	0.368
0.6	3.759E-2	2.512E-2	1.672E-2	1.109E-2	7.342E-3
	0.581	0.588	0.592	0.595
0.8	1.121E-2	6.401E-3	3.666E-3	2.102E-3	1.206E-3
	0.809	0.804	0.803	0.802

Table 2

Example 5.1: Maximum errors and orders of convergence for scheme (3.2) in the subdomain $Q¯*$ .

α	$N=M=128$	$N=M=256$	$N=M=512$	$N=M=1024$	$N=M=2048$
0.4	1.024E-2	4.966E-3	2.436E-3	1.214E-3	6.050E-4
	1.044	1.027	1.005	1.005
0.6	1.300E-2	6.432E-3	3.190E-3	1.595E-3	7.965E-4
	1.015	1.012	1.000	1.002
0.8	9.844E-3	5.123E-3	2.644E-3	1.361E-3	6.963E-4
	0.942	0.954	0.959	0.966

Figure 2

Example 5.1: Log-log plot of the error bound $N-1⁢tnα-1$ ( $⋄$ ) and the maximum errors $max0≤m≤M⁡|u⁢(xm,tn)-umn|$ ( $∘$ ) at each time level $t=tn,n=1,2,…,N$ , generated by scheme (3.2) for $N=M=100$ .

$(a) α=0.2{\alpha=0.2}$

(a)

$α=0.2$

$(b) α=0.8{\alpha=0.8}$

(b)

$α=0.8$

Example 5.2.

Consider the variable coefficient inhomogeneous problem

(5.1)

(5.1a)

$Dtα⁢u-∂2⁡u∂⁡x2+(1+x2)⁢∂⁡u∂⁡x+(1+x)⁢u=4π2⁢x⁢(π-x) for ⁢(x,t)∈Q,$

$with Q=(0,π)×(0,1] and$

(5.1b)

$u⁢(x,0)=0for ⁢0<x<π,u⁢(0,t)=u⁢(π,t)=0for ⁢0≤t≤1.$

Figure 3 displays the computed solution for $α=0.4,0.8$ and $N=M=64$ and we observe that the solution again exhibits an initial layer at $t=0$ .

Figure 3

Example 5.2: Computed solutions with scheme (3.2) for $N=M=64$ .

$(a) α=0.4{\alpha=0.4}$

(a)

$α=0.4$

$(b) α=0.8{\alpha=0.8}$

(b)

$α=0.8$

The exact solution of Example 5.2 is unknown and we shall estimate the order of convergence using the two-mesh principle [2]. Let $umn$ be the computed solution with scheme (3.2) on the mesh ${(xm,tn)}$ for $m=0,1,…,M$ , $n=0,1,…,N$ . To estimate the order of convergence, we compute a new approximation $zm/2n2$ using the same scheme defined on the finer mesh ${(xm/2,tn/2)}$ for $m=0,1,…,2⁢M$ , $n=0,1,…,2⁢N$ , where $xm+1/2:=12⁢(xm+1+xm)$ and $tn+1/2:=12⁢(tn+1+tn)/2$ . We then compute the two-mesh differences

$dM,N:=max(xm,tn)∈Q′⁡|umn-zmn|$

and hence the estimated orders of convergence

$qM,N:=log2⁡(dM,Nd2⁢M,2⁢N).$

Tables 3 and 4 give the maximum two-mesh differences and their corresponding orders of convergence for Example 5.2 in the domain $Q¯$ and the subdomain $Q¯*$ . The numerical results in both cases are again in agreement with Theorem 4: the order of convergence improves from $O⁢(N-α)$ on $Q¯$ to $O⁢(N-1)$ on $Q¯*$ .

Table 3

Example 5.2: Maximum two-mesh differences and orders of convergence for scheme (3.2) in the domain $Q¯$ .

α	$N=M=128$	$N=M=256$	$N=M=512$	$N=M=1024$	$N=M=2048$
0.4	1.031E-2	8.673E-3	7.123E-3	5.740E-3	4.558E-3
	0.250	0.284	0.311	0.333
0.6	4.935E-3	3.338E-3	2.234E-3	1.486E-3	9.857E-4
	0.564	0.579	0.588	0.593
0.8	1.661E-3	9.441E-4	5.368E-4	3.060E-4	1.748E-4
	0.815	0.815	0.811	0.808

Table 4

Example 5.2: Maximum two-mesh differences and orders of convergence for scheme (3.2) in the subdomain $Q¯*$ .

α	$N=M=128$	$N=M=256$	$N=M=512$	$N=M=1024$	$N=M=2048$
0.4	5.849E-4	2.783E-4	1.351E-4	6.711E-5	3.337E-5
	1.072	1.042	1.010	1.008
0.6	1.148E-3	5.457E-4	2.628E-4	1.291E-4	6.356E-5
	1.073	1.054	1.025	1.023
0.8	1.335E-3	6.752E-4	3.387E-4	1.703E-4	8.531E-5
	0.984	0.995	0.992	0.997

In [10], numerical results were given for the particular case of a fractional reaction-diffusion equation (i.e., with $q≡0$ in (2.1)) showing that the scheme also converges with order α when the whole domain is considered. Additional numerical results that illustrate the improved rate of convergence away from $t=0$ are given in [3].

Funding source: National Natural Science Foundation of China

Award Identifier / Grant number: 91430216

Award Identifier / Grant number: NSAF U1530401

Funding source: Ministerio de Economía y Competitividad

Award Identifier / Grant number: MTM2016-75139-R

Funding statement: The research of José Luis Gracia was partly supported by the Institute of Mathematics and Applications (IUMA), the project MTM2016-75139-R funded by the Spanish Ministry of Economy and Competitiveness (MINECO) and the Diputación General de Aragón. The research of Martin Stynes was supported in part by the National Natural Science Foundation of China under grants 91430216 and NSAF U1530401.

References

[1] K. Diethelm, The Analysis of Fractional Differential Equations, Lecture Notes in Math. 2004, Springer, Berlin, 2010. 10.1007/978-3-642-14574-2Search in Google Scholar

[2] P. A. Farrell, A. F. Hegarty, J. J. H. Miller, E. O’Riordan and G. I. Shishkin, Robust Computational Techniques for Boundary Layers, Appl. Math. Math. Comput. 16, Chapman & Hall/CRC, Boca Raton, 2000. 10.1201/9781482285727Search in Google Scholar

[3] J. L. Gracia, E. O’Riordan and M. Stynes, Convergence outside the initial layer for a numerical method for the time-fractional heat equation, Numerical Analysis and Its Applications, Lecture Notes in Comput. Sci. 10187, Springer, Cham (2017), 82–94. 10.1007/978-3-319-57099-0_8Search in Google Scholar

[4] B. Jin, R. Lazarov and Z. Zhou, An analysis of the L1 scheme for the subdiffusion equation with nonsmooth data, IMA J. Numer. Anal. 36 (2016), no. 1, 197–221. 10.1093/imanum/dru063Search in Google Scholar

[5] B. Jin and Z. Zhou, An analysis of Galerkin proper orthogonal decomposition for subdiffusion, ESAIM Math. Model. Numer. Anal. 51 (2017), no. 1, 89–113. 10.1051/m2an/2016017Search in Google Scholar

[6] W. McLean and K. Mustapha, Time-stepping error bounds for fractional diffusion problems withnon-smooth initial data, J. Comput. Phys. 293 (2015), 201–217. 10.1016/j.jcp.2014.08.050Search in Google Scholar

[7] K. Mustapha, An implicit finite-difference time-stepping method for a sub-diffusion equation, with spatial discretization by finite elements, IMA J. Numer. Anal. 31 (2011), no. 2, 719–739. 10.1093/imanum/drp057Search in Google Scholar

[8] K. Mustapha and W. McLean, Piecewise-linear, discontinuous Galerkin method for a fractional diffusion equation, Numer. Algorithms 56 (2011), no. 2, 159–184. 10.1007/s11075-010-9379-8Search in Google Scholar

[9] M. Stynes, Too much regularity may force too much uniqueness, Fract. Calc. Appl. Anal. 19 (2016), no. 6, 1554–1562. 10.1515/fca-2016-0080Search in Google Scholar

[10] M. Stynes, E. O’Riordan and J. L. Gracia, Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation, SIAM J. Numer. Anal. 55 (2017), no. 2, 1057–1079. 10.1137/16M1082329Search in Google Scholar

Received: 2017-1-27

Revised: 2017-5-27

Accepted: 2017-6-6

Published Online: 2017-7-6

Published in Print: 2018-1-1

Articles in the same Issue

DOI: doi.org/10.1515/cmam-2017-0019

Keywords for this article

Caputo Fractional Derivative; Initial Boundary Value Problem; Weak Singularity; L1 Scheme