On the time fractional heat equation with obstacle
Carlo Alberini, Raffaela Capitanelli, Mirko D'Ovidio, Stefano Finzi Vita
OOn the time fractional heat equation with obstacle
C. Alberini , R. Capitanelli , M. D’Ovidio , and S. Finzi Vita Dipartimento SBAI, Sapienza Universit`a di Roma, [email protected] Dipartimento SBAI, Sapienza Universit`a di Roma, raff[email protected] Dipartimento SBAI, Sapienza Universit`a di Roma, [email protected] Dipartimento di Matematica, Sapienza Universit`a di Roma, stefano.fi[email protected]
December 23, 2020
Abstract.
We study a Caputo time fractional degenerate diffusion equation whichwe prove to be equivalent to the fractional parabolic obstacle problem, showing that itssolution evolves for any α ∈ (0 ,
1) to the same stationary state, the solution of the classicelliptic obstacle problem. The only thing which changes with α is the convergence speed.We also study the problem from the numerical point of view, comparing some finite differentapproaches, and showing the results of some tests. These results extend what recentlyproved in [1] for the case α = 1. In the works of the last decades, the use of fractional calculus has known an ever increasinguse in describing phenomena from the most disparate scientific fields, from biology tomechanics, to superslow diffusion in porous media till financial interests.Aim of this paper is to study the following problem, which generalizes the one addressedin [1], for 0 < α < ∂ αt u − H ( u − ψ ) ∆ u = 0 a.e. in Ω , for all t ∈ (0 , T ) ,u (0) = u in Ω ,u = 0 on ∂ Ω , for all t ∈ (0 , T ) , (1)with T >
0, where H is the extended Heaviside function such that H (0) = 0, that is H ( r ) = r >
00 for r ≤ a r X i v : . [ m a t h . A P ] D ec nd Ω is a bounded domain in R n , n ∈ N , with smooth boundary. For simplicity, weomitted the presence of a forcing term in the problem.Here ∂ αt u denotes the Caputo fractional derivative, that is, ∂ αt u ( x, t ) := 1Γ(1 − α ) (cid:90) t ( t − s ) − α ∂ s u ( x, s ) ds, (3)where Γ is the Gamma function.The Riemann-Liouville derivative is defined as R ∂ αt u ( x, t ) := ∂∂t − α ) (cid:90) t ( t − s ) − α u ( x, s ) ds. (4)We observe that R ∂ αt C = 1Γ(1 − α ) ∂ t (cid:90) t C ( t − s ) α ds = C t − α Γ(1 − α )and the following relation between derivatives holds ∂ αt u ( x, t ) = R ∂ αt (cid:0) u ( x, t ) − u ( x, (cid:1) = R ∂ αt u ( x, t ) − u ( x, t − α Γ(1 − α ) . (5)We refer to the book [5] for further details on fractional derivatives.Let us consider the complementarity system, for all t > w ( x, t ) ≥ ψ ( x ) in Ω ∂ αt w ≥ ∆ w a.e. in Ω( w − ψ )( ∂ αt w − ∆ w ) = 0 a.e. in Ω w ( x,
0) = u in Ω w ( x, t ) = 0 on ∂ Ω . (6)In Section 2 we will prove that - by analogy with the classical case, (1) is equivalent to the complementarity system (6); - asymptotically, the solution evolves for each α towards the same steady state of theclassical problem, that is¯ u ≥ ψ, − ∆¯ u ≥ , (¯ u − ψ )∆¯ u = 0 in Ω × (0 , T ) (7)with different convergence speed (for α = 1 exponential, for α ∈ (0 ,
1) polynomial).2n Section 3 we present three possible finite difference schemes for the numerical ap-proximation of problem (1) or its equivalent form given by the complementary system (6).The Caputo derivative has been discretized using standard methods found in the literature,the so-called L1 or Convolution Quadrature (CQ) approaches (see e. g. [9] or [7]). Thespace discretization has been carried out through the semi implicit finite differences schemeintroduced in [1] for problem (1), or through the implicit scheme of [3] for the evolutiveobstacle problem (6). Note that for α → − all these schemes give back the known resultsof the classic heat equation with obstacle.Our aim was not to find an optimal strategy of approximation for the problem, but onlyto derive working schemes in orders to confirm through explicit simulations the behaviorof the solution as characterized by the results of Section 2. For that, in Section 4, we havetested the previous schemes on a couple of one-dimensional examples. Such simulationsalso allow to compare their reliability and computational cost. The semi implicit approachrequires time step restrictions (strongly increasing as α → + ) in order to detect the correctcontact set with the obstacle, restrictions unnecessary for the implicit approach. On theother hand, in the first case each time iteration is much less expensive, since it reducesto a single linear system solution. So, after all, all the proposed schemes appear to becompetitive.For the following, we assume that the initial datum u and the (independent of time)obstacle ψ satisfy the following conditions u ∈ H (Ω) , ψ ∈ H (Ω) , ψ ≤ ∂ Ω . (8)We define as solution of problem (1) a function u ∈ L (cid:0) , T ; H (Ω) ∩ H (Ω) (cid:1) with ∂ αt u ( x, t ) ∈ L (cid:0) , T ; L (Ω) (cid:1) which solves problem (1). In the present section we analyze problem (1), showing that, under suitable conditions, itis equivalent to complementarity system (6). Let us introduce the following hypoteses: H : u > ψ a.e. in Ω; H : ∆ ψ ≤ Proposition 2.1.
Assume that conditions (8) , H and H hold. Then problems (1) and (6) are equivalent. Proof.
First we note that if w solves problem (6), then w solves (1). In fact, initial andboundary conditions are the same in (1) and (6). Since ( w − ψ )( ∂ αt w − ∆ w ) = 0 a.e. in Ω,where w > ψ , then we obtain ∂ αt w − ∆ w = 0 a.e in Ω whereas, where w ( x, t ) = ψ ( x ) then ∂ t w = 0 which implies ∂ αt w = 0. 3ow we prove that if u solves problem (1) then u coincides with solution w of (6).Initial and boundary conditions are the same in (1) and (6). Let us now prove that if u ( x, t ) is a solution of (1) then necessarily u ( ., t ) ≥ ψ ( . ) in Ω for any time t . Assumethat u < ψ . Then ∂ αt u = 0 from (1). From H , u > ψ . Moreover, due to the fact that u < ψ < u , we obtain from (5) ∂ αt u = R ∂ αt ( u − u ) < , which does not agree with ∂ αt u = 0. Then u ≥ ψ by contradiction and the first inequalityof (6) holds. The equation in the third line of (6) is trivially satisfied where u ( x, t ) = ψ ( x );where u ( x, t ) > ψ ( x ) , from (1) we obtain ∂ αt − ∆ u = 0, so it is always true. Concerningthe second inequality of (6), we have already seen that it is satisfied (with the equal sign)when u > ψ . But when u ( x, t ) = ψ ( x ), (1) and assumption H imply that ∂ αt u − ∆ u = − ∆ ψ ≥ . Then u solves problem (6).We recall that the same problem for α = 1 has been faced in [1], proving in particularthe equivalence of problem (1) with a parabolic obstacle problem. Let v be the uniquesolution of problem (6) with α = 1, that is v ( x, t ) ≥ ψ ( x ) in Ω ∂ t v ≥ ∆ v a.e. in Ω( v − ψ )( ∂ t v − ∆ v ) = 0 a.e. in Ω v ( x,
0) = u in Ω v ( x, t ) = 0 . on ∂ Ω (9)Let us introduce u ( x, t ) = (cid:90) ∞ v ( x, s ) (cid:96) ( s, t ) ds (10)where the probability density (cid:96) : [0 , ∞ ) × [0 , ∞ ) → [0 , ∞ ) satisfies (see [6]) ∂ αt (cid:96) = − ∂(cid:96)∂s , ( s, t ) ∈ (0 , ∞ ) × (0 , ∞ ) (cid:96) ( s,
0) = δ ( s ) , s > ,(cid:96) ( s, t ) = 0 , s < . and R ∂ αt (cid:96) = − ∂(cid:96)∂s , ( s, t ) ∈ (0 , ∞ ) × (0 , ∞ ) ,(cid:96) ( s,
0) = δ ( s ) , s > ,(cid:96) (0 , t ) = t − α Γ(1 − α ) , t > ,(cid:96) ( s, t ) = 0 , s < . roposition 2.2. The function u defined in (10) solves (6) . Proof. As v ≥ ψ , then (cid:90) ∞ v ( x, s ) (cid:96) ( s, t ) ds ≥ (cid:90) ∞ ψ ( x ) (cid:96) ( s, t ) ds = ψ ( x ) , and we obtain that u ≥ ψ. Now we prove that ( u − ψ )( ∂ αt u − ∆ u ) = 0 . For u > ψ , it holds the classical theoryabout fractional Cauchy problems (see [2], [10] and [4, Theorem 5.2]) and therefore we havethat, in L (Ω), ∂ αt u ( x, t ) = ∆ u ( x, t ) , u ( x,
0) = u ( x ) . Concerning the second inequality of (6), we observe that if ∂ t v − ∆ v ≥
0, then∆ (cid:90) ∞ v ( x, s ) (cid:96) ( s, t ) ds ≤ (cid:90) ∞ ∂ s v ( x, s ) (cid:96) ( s, t ) ds = − v ( x, t − α Γ(1 − α ) − (cid:90) ∞ v ( x, s ) ∂ s (cid:96) ( s, t ) ds = − v ( x, t − α Γ(1 − α ) + R ∂ αt u ( x, t ) . Since u ( x,
0) = v ( x,
0) from the relation (5) we obtain ∂ αt u − ∆ u ≥ u follow trivially fromthe initial condition and the boundary condition for the function v .In [1] the asymptotic solution of (9) has been characterized as the solution of thecorresponding stationary (elliptic) obstacle problem. In particular, it has been proved thatunder conditions (8), H and H , the solution v ( t ) to the parabolic obstacle problem (9)converges strongly in H (Ω) , for t → ∞ , to the unique solution u of the correspondingstationary obstacle problem u ∈ H (Ω) , u ≥ ψ, − ∆ u ≥ , ( u − ψ )(∆ u ) = 0 a.e. in Ω . (11)Moreover, it has been proved that there is a constant C > t ≥ (cid:107) v ( t ) − u (cid:107) H (Ω) ≤ e − Ct . (12)In the next Theorem we prove that a similar result holds for any α ∈ (0 , . Theorem 2.1.
Assume conditions (8) , H and H hold. Let u ( t ) be the function definedin (10) solution of (6) . Let u be the unique solution of the obstacle problem (11) . Then, u ( t ) converges to u for t → ∞ , and there is a constant C > such that, for every t ≥ , (cid:107) u ( t ) − u (cid:107) L (Ω) ≤ κ Ω E α ( − Ct α ) (13)5 or any α ∈ (0 , , where κ Ω = (cid:18)(cid:90) Ω dx (cid:19) / and E α ( − Ct α ) J ( t ) → , as t → ∞ , (14) with J ( t ) = 1 C t − α Γ(1 − α ) . (15) Proof.
Since (cid:107) v − ¯ u (cid:107) L (Ω) ≤ κ Ω (cid:107) v − ¯ u (cid:107) L (Ω) we have (cid:107) u ( t ) − ¯ u (cid:107) L (Ω) = (cid:90) Ω (cid:12)(cid:12)(cid:12)(cid:12) (cid:90) ∞ v ( x, s ) (cid:96) ( s, t ) ds − ¯ u ( x ) (cid:12)(cid:12)(cid:12)(cid:12) dx = (cid:90) Ω (cid:12)(cid:12)(cid:12)(cid:12) (cid:90) ∞ ( v ( x, s ) − ¯ u ( x )) (cid:96) ( s, t ) ds (cid:12)(cid:12)(cid:12)(cid:12) dx ≤ (cid:90) Ω (cid:90) ∞ (cid:12)(cid:12)(cid:12)(cid:12) v ( x, s ) − ¯ u ( x ) (cid:12)(cid:12)(cid:12)(cid:12) (cid:96) ( s, t ) dsdx = (cid:90) ∞ (cid:107) v ( x, s ) − ¯ u ( x ) (cid:107) L (Ω) (cid:96) ( s, t ) ds ≤ κ Ω (cid:90) ∞ e − Cs (cid:96) ( s, t ) ds = κ Ω E α ( − Ct α ) , ∀ α ∈ (0 , E α is the Laplace transform of the density (cid:96) .For the asymptotic behaviour of the Mittag-Leffler, consult the book [8, formula (4.4.17)]. The starting idea for a numerical approximation of the problem has been to combine classi-cal time-stepping discretization schemes for the Caputo derivative, such as the convolutionquadrature (CQ) or the finite difference (L1) scheme, see e.g. [9], with schemes usuallyworking for the parabolic obstacle problem (6), such as the one proposed in [3], but alsothe semi-implicit f.d. scheme tested in [1] for the equivalent Heaviside function formulation(1) of the problem.For sake of simplicity we work in the one-dimensional case, with Ω = ( a, b ). Let us call h > h = ( b − a ) /N , so that we have ( N −
1) internal nodes6n Ω, x i = a + ih , for i = 1 , ..., N −
1) and τ = T /M the time discretization step (with M time instants t m = mτ , for m = 1 , ..., M ); with γ α = τ α /h we denote the parabolicratio between the steps related to a specific α . Note that for fixed steps h and τ , if α decreases to zero then γ α quickly grows: in other words for small α , in order to keep γ α small on a fixed mesh, the step τ has to be considerably reduced, with a significant increaseof computational costs.Here we analyze three possible approaches:1. Scheme S1: solves problem (1) with L1 for the Caputo derivative and thesemi implicit f.d. scheme of [1] in space.
The time discretization of the Caputo derivative by the L1 scheme leeds to the formula(see [9]): ∂ αt u ( x, t m ) (cid:39) − α ) τ α (cid:40) u ( x, t m ) − m − (cid:88) k =0 C m,k u ( x, t k ) (cid:41) , where C m, := f ( m ) , C m,k := f ( m − k ) − f ( m − ( k − k = 1 , ..., m − , and f ( r ) := r − α − ( r − − α , for r ≥ . Then, after semidiscretization in time, we need to solve for any instant t m the equation1Γ(2 − α ) τ α (cid:40) u ( x, t m ) − m − (cid:88) k =0 C m,k u ( x, t k ) (cid:41) − H ( u − ψ ) u xx = 0 in Ω , (16)with the same boundary conditions of (1). Since Ω was splitted in N subintervalsthrough the nodes x = { x i } i , the initial data will be the vector u ∈ R N − , with u i = u ( x i ); applying the semi implicit finite difference scheme of [1] in space, thesolution u at the first discrete time t = τ will solve at any node the relation:1 gτ α (cid:0) u i − C , u i (cid:1) = H ( u i − ψ i ) δ u i := H ( u i − ψ i ) u i − − u i + u i +1 h , with g = Γ(2 − α ) (note that 0 . ≤ g ≤ , ∀ α ∈ [0 , γ α ). If we set v ki = u ki − ψ i , redistributing all the termsbetween the two members, it is equivalent to solve u i − gτ α H ( v i ) δ u i = C , u i for every i ;with vector notations it means that u solves the linear system: B u := ( I + gγ α H ( v ) ∗ A ) u = C , u , A is the usual tridiagonal matrix ( N − × ( N −
1) with values 2 on the maindiagonal and − { ( H ( v ) ∗ A ) u } i := H ( v i )( Au ) i = H ( v i ) N − (cid:88) j =1 A i,j u j . Since the discrete solution could overstep the obstacle at some nodes, in particularwhen a large value of γ α is used, the following correction is needed at any iteration: u i = max( u i , ψ ( x i )) . In the same way we see that u solves for any i gτ α (cid:0) u i − C , u i − C , u i (cid:1) = H ( u i − ψ i ) δ u i , that is the system B u = ( I + gγ α H ( v ) ∗ A ) u = C , u + C , u , with the same matrix A and the subsequent correction. In general, at any time step u m solves the linear system B m − u m = b m , (17)where we have set B m − = I + gγ α H ( v m − ) ∗ A, b m = m − (cid:88) k =0 C m,k u k , (18)followed by the correction u mi = max( u mi , ψ ( x i )) . (19)Note that all the matrices B m are symmetric positive definite (and M-matrices), sinceso it is A , while H ( v ) ≥
0. Then all the previous linear systems are well posed.With respect to the classic parabolic obstacle problem approach discussed in [1] thereis here an important difference. In that case (which corresponds to the case α = 1),when the solution at time t m − touches the obstacle at node x i , then v m − i = 0, H ( v m − i ) = 0, and (17) trivially yields: u mi = u m − i ;in other words, once touched the obstacle at a particular node the solution does notchange anymore there, but only at the remaining free nodes. In the general case of8 ∈ (0 ,
1) on the contrary, at the contact time system (17) immediately yields for the i -th component: u mi = m − (cid:88) k =0 C m,k u ki = C m, u i + ... + C m,m − ψ i > ψ i m − (cid:88) k =0 C m,k = ψ i , since at least u i > ψ i and (cid:80) m − k =0 C m,k = 1. It follows that the solution has a littlerebound at x i which detaches it again from the obstacle, and produces an (innatural)oscillating evolution from that time on. The width of such rebound depends on thesize of γ α . Then, even if in principle the semi implicit scheme does not not requirestability restrictions on the discretization steps, a small value of γ α will be necessaryto reduce the oscillations (they would vanish for τ → b m of (18) by the vector ˆ b m defined byˆ b m = max( H ( v m − ) ∗ b m , u m − ); (20)if H ( v m − i ) = 0, that is u m − i = ψ i , then ˆ b mi = u m − i and u mi = u m − i ; otherwiseˆ b mi = b mi , since b mi ≥ u m − i , and all remains as before. No rebound is still possibleafter a contact.2. Scheme S2: solves problem (1) with CQ for the Caputo derivative and thesemi implicit f.d. scheme of [1] in space.
In this case the Caputo derivative is approximated through the so-called convolu-tion quadrature (CQ) method, proposed by Lubich for the discretization of Volterraintegral equations. In particular, for a function φ ( t ) (with φ (0) = 0), the Riemann-Liouville derivative R ∂ αt φ defined in (4) can be approximated by the discrete convo-lution: R ∂ ατ φ m := 1 τ α m (cid:88) j =0 c j φ m − j , where φ m = φ ( t m ), and the coefficients { c j } are obtained from a suitable powerseries expansion, connected to a specific approximation method for the ODE (see[9]). In the case of the Euler backward method, it is known as the Grunwald-Letnikovapproximation, and provides the following recursive formula for the coefficients: c = 1 , c j = − α − j + 1 j c j − . Then, using the relation (5) between the Caputo and the Riemann-Liouville deriva-tives we can rewrite the initial problem as R ∂ αt ( u − u ) − H ( u − ψ )∆ u = 0 , τ α m (cid:88) j =0 c j ( u m − j − u ) + 1 h H ( u m − − ψ ) Au m = 0 , equivalent to the solution at any iteration of the linear system B m − u m = b m , (21)where this time we have set B m − = I + γ α H ( u m − − ψ ) ∗ A, b m = u − m − (cid:88) j =1 c j ( u m − j − u ) , (22)followed again by the correction u mi = max( u mi , ψ ( x i )) . (23)Even in this case the vector b m has to be modified inside the contact set in orderto remove the memory effect and prevent rebounds, as done in (20). In fact when u m − i = ψ i then again b mi = 0, and from (22) we get u mi = u i − m − (cid:88) j =1 c j ( u m − ji − u i ) > u i − ( ψ i − u i ) m − (cid:88) j =1 c j > ψ i , since u m − ji ≥ ψ i for any j , and (cid:80) m − j =1 c j ≥ − Scheme S3: solves problem (6) with L1 for the Caputo derivative and thescheme of [3] for the evolutive obstacle problem.
If we discretize the equation of system (6) through finite differences, using the L1scheme for the Caputo derivative, we get the equation( u m − ψ ) T (cid:32) gτ α ( u m − m − (cid:88) k =0 C m,k u k ) + 1 h Au m (cid:33) = 0 . Setting y m = u m − ψ , and remembering that (cid:80) m − k =0 C m,k = 1, it is equivalent to y m (cid:32) y m + gγ α A ( y m + ψ ) − m − (cid:88) k =0 C m,k y k (cid:33) = 0 . Then y m = max(0 , x m ) is solution of the previous equation if x m solves( I + gγ α AP ( x )) x = b m , (24)10here now b m = m − (cid:88) k =0 C m,k y k − gγ α Aψ = m − (cid:88) k =0 C m,k u k − ψ − gγ α Aψ , (25)while P ( x ) = diag ( p ( x i )) denotes the diagonal matrix with p ( x i ) = 1 if x i > p ( x i ) = 0 otherwise. As seen in [3], (24) can be solved by the so-called Picarditerations: P = O, ( I + gγ α AP n ) x n +1 = b m , P n +1 = diag ( p ( x n +1 )) per n = 0 , , ... until P n +1 = P n ( O is the null matrix); at that point x n +1 is the sought solution.Of course other schemes could be obtained by different combinations of specific numer-ical approaches, but for our purposes the three previous schemes were sufficient to performin the next section explicit simulations of the problem. We have applied the schemes described in Section 3 to some specific examples, for differentvalues of α . We have choosen a sufficiently large final time T , but also added a stoppingtime criterium in order to put in evidence the convergence towards the asymptotic solution.Since this convergence corresponds to the stabilization of the solution vector and to thesatisfaction of the asymptotic complementarity relation( u − ψ )∆ u = 0 , which means u harmonic (linear in 1D) outside the contact set, we have used the criterium ST OP when max( (cid:107) u m − u m − (cid:107) ∞ , (cid:107) ( u m − ψ ) Au m (cid:107) ∞ ) < tol , (26)for a given tolerance parameter tol . Example 1.
Ω = ( − , u ( x ) = 0 . − . x , ψ ( x ) = 0 . − x (see Fig.1).In the following Table 1 we reported some results obtained by the simulations withschemes S1, S2 and S3, for different values of α , N and γ α , with tol = 10 − . We adoptedthe following notations: • FC time = full contact time, that is the first time at which the solution has reachedthe whole contact set (no more changing in the successive iterations); • STOP time = the exit time according to criterium (26);11
Contact set: C=[-0.12,0.12]
Example 1 initial datumstationary solutionobstacle
Figure 1:
Data of Example 1. • • • • S = working schemes (the ones detecting the correct contact set).Looking at the table some remarks and comments are possible: • The stationary solution is the same for any α , as stated by Theorem 2.1, and cor-responds to the one of the stationary problem (7). Different is only the speed ofconvergence. The detected right extremum of the symmetric contact set C on theused meshes is 0.125 (the continuous value should be approximately 0 . • Schemes S1 and S2 have a very similar behaviour; in both cases for large values of γ α the contact set can be overestimated, a problem not present for S3, due to theimplicit nature of the quasi-Newton scheme of [3]. The semi implicit schemes S1and S2 pay for the delay with which the contact information with the obstacle isachieved, allowing an uncorrect evolution of the solution. To avoid that, strong τ -step restrictions are necessary: experiments show that all the schemes work correctlyfor example if τ α < .
1, a bound which becomes particularly heavy when α is small.In the Table we reported approximately for any number of nodes the largest valuesof γ α for which all the three schemes give the same correct solution. • Computational costs: the previous remark suggests that S3 is the more reliable andeven the less expensive of the three schemes, allowing larger time steps and hence less12 able 1 α N γ α τ FC time STOP time .
16 11739 12 140868 all32 75 0.016 0.11 93 .
15 5580 13 72540 S364 220 0.0059 0.49 1.33 226 21 4746 all128 850 0.005 0.2 1.66 315 42 13230 all0.5 32 25 0.009 0.16 8.38 880 8 7040 all32 50 0.038 0.19 8.39 221 11 2431 S364 100 0.009 0.39 2.13 225 15 3375 all128 400 0.009 0.8 2.67 282 28 7896 all0.7 32 50 0.097 0.29 5.81 61 11 671 all32 100 0.026 0.52 10.44 41 14 574 S364 200 0.097 0.38 7.08 74 21 1554 all128 400 0.036 0.5 4.82 135 29 3915 all1 32 15 0.058 0.27 0.99 18 7 126 all32 20 0.078 0.31 1.09 15 8 120 S364 60 0.058 0.35 1.05 19 13 247 all128 240 0.058 0.64 1.05 19 23 437 alliterations. Anyway, any single time iteration of S3 is much more expensive, manyPicard iterations (growing with N ) with respect to a single linear system necessaryto be solved for S1 and S2. Then, comparing the computational cost of the schemes,even these two schemes reveal competitive in terms of the total number of linearsystem solved. • For a fixed α the full contact time grows with the number of nodes, and does notseem to depend from γ α . On the contrary the stabilization time grows with γ α butdecreases with the number of nodes, since a finer mesh reduces the error at theboundary of the contact set. • All the schemes correctly work also for the case α = 1: it easy to see that in thatcase both the used approximations of the Caputo derivative reduce to the standardincremental ratio in time. • The case α = 0 is a sort of control test: since in that case ∂ αt u = u − u , in absence13f an obstacle the equation (1) would reduce to the stationary equation − ∆ u + u = u . (27)In presence of an obstacle we expect the discrete solution to satisfy (27) only outsidethe contact set, and from the first iteration. It is in fact clear from (18) and (22) thatsince H ( u − ψ ) = 1 and γ α = 1 /h , the first iteration of all the schemes becomes B u = (cid:18) I + 1 h A (cid:19) u = u , which is essentially the discrete version of (27). If u goes over the obstacle, suchidentity will be satisfied only where u > ψ . In Figure 2 it can be seen what happensin our example with N = 128 nodes and scheme S3: the solution u is plotted inblue, the obstacle in red, the initial datum u in black and the quantity − ∆ u + u through asterisks: the last two quantities coincide in the detachment set, with anatural discontinuity at the boundary of such a region. Since the time step τ has noeffect on the solution, the semi implicit schemes S S -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1-0.500.511.522.533.544.55 S3, alfa=0, N=128 -Du+uu0ucu
Figure 2:
Solution for α = 0 and 128 nodes, scheme S3. Example 2.
Ω = ( − , u ( x ) = 1 − x , ψ ( x ) = 0 . − (2 x − . (see Fig.3).On this example we tested numerically the estimate (13) of Theorem 2.1, using schemeS3. The L norm of the error at time t m on the given mesh was approximated by a naturalquadrature formula, that is (cid:107) u ( t m ) − u (cid:107) L (Ω) (cid:39) h (cid:88) i | u mi − ¯ u i | Contact set: C=[-0.6,-0.52] [0.52,0.6]
Example 2 initial datumstationary solutionobstacle
Figure 3:
Data of Example 2.(the vector ¯ u on the mesh was computed in advance with a sufficiently high precision).Such an error was computed for different values of α and of the time t m . In Table 2 wereported (in the third column) the discrete L error at T = 10 with N = 32 nodes andthe same γ α = 50 for any α . In the fourth column the corresponding values of quantity J ( T ) = C T − α Γ(1 − α ) of (15) are shown; for the constant C we adopted a computed estimate( C = 26). It is evident that the two quantities decay at the same rate. In the second columnit is also possible to see the number of instant times M needed for each α to keep the same γ α , and the consequent increasing complexity of the computation. In order to confirm theorder of decay of the error, polynomial for α ∈ (0 ,
1) and exponential (up to the machineprecision) for α = 1, we illustrate in Figure 4 such behavior plotting the continuous curves J ( t ) in the time interval (0 ,
20) for α = 0 . , . .
9, and the computed errors on tendifferent discrete times (marked by asterisks). These values correctly lie on the curves. time L1 e rr o r ||u(t)-ubar|| for alfa=0.3, 0.6, 0.9 Figure 4:
Example 2: polynomial speed of stabilization for α = 0 . , . . able 2: Example 2: N = 32, γ α = 50, scheme S α M L error at T = 10 T − α / ( C Γ(1 − α ))1 51 8 .
07 10 − (cid:39) .
27 10 − .
09 10 − .
37 10 − .
32 10 − .
62 10 − .
56 10 − .
41 10 − .
35 10 − .
87 10 − .
86 10 − .
01 10 − .
02 10 − .
44 10 − .
48 10 − .
98 10 − .
08 10 − . − (estimate) 2 .
85 10 − References [1] Alberini C., Capitanelli R. and Finzi Vita S.,
A numerical study of an Heavisidefunction driven degenerate diffusion equation , preprint 2020, submitted.[2] Bazhlekova E.G.,
Subordination principle for fractional evolution equations , FractionalCalculus and Applied Analysis. , No. 3, (2000), 213–230.[3] Brugnano L. and Sestini A., Iterative solution of piecewise linear systems for thenumerical solution of obstacle problems , J. Num. Anal. Ind. Appl. Math. (JNAIAM), , No. 3-4, (2011), 67–82.[4] Capitanelli R. and D’Ovidio M., Fractional equations via convergence of forms , Frac-tional Calculus and Applied Analysis, , (2019), 844–870.[5] Diethelm K., The analysis of fractional differential equations. An application-orientedexposition using differential operators of Caputo type , Lecture Notes in Mathematics,2004. Springer-Verlag, Berlin, (2010).[6] D’Ovidio M.,
On the fractional counterpart of the higher-order equations , Statisticsand Probability Letters, , (2011), 1929–1939.[7] Giga Y., Liu Q. and Mitake H., On a discrete scheme for time fractional fully nonlinearevolution equations , Asymptotic Analysis, , No. 1-2, (2020), 151–162.168] Gorenflo R., Kilbas A.A., Mainardi F. and Rogosin S.V.,
Mittag-Leffler Functions, Re-lated Topics and Applications , Springer Monographs in Mathematics, Springer-VerlagBerlin Heidelberg (2014).[9] Jin B, Lazarov R. and Zhou Z.,
Numerical methods for time-fractional evolution equa-tions with nonsmooth data: a concise overview , Comput. Methods Appl. Mech. Engrg., , (2019), 332–358.[10] Kochubei A.N.,
The Cauchy problem for evolution equations of fractional order , Dif-ferential Equations,25