0
Research Papers

An Explicit Difference Method for Solving Fractional Diffusion and Diffusion-Wave Equations in the Caputo Form

[+] Author and Article Information
Joaquín Quintana Murillo

Departamento de Física, Universidad de Extremadura, E-06071 Badajoz, Spain

Santos Bravo Yuste1

Departamento de Física, Universidad de Extremadura, E-06071 Badajoz, Spainsantos@unex.es

1

Corresponding author.

J. Comput. Nonlinear Dynam 6(2), 021014 (Nov 15, 2010) (6 pages) doi:10.1115/1.4002687 History: Received January 18, 2010; Revised September 17, 2010; Published November 15, 2010; Online November 15, 2010

An explicit difference method is considered for solving fractional diffusion and fractional diffusion-wave equations where the time derivative is a fractional derivative in the Caputo form. For the fractional diffusion equation, the L1 discretization formula of the fractional derivative is employed, whereas the L2 discretization formula is used for the fractional diffusion-wave equation. In both equations, the spatial derivative is approximated by means of the three-point centered formula. The accuracy of the present method is similar to other well-known explicit difference schemes, but its region of stability is larger. The stability analysis is carried out by means of a kind of fractional von Neumann (or Fourier) method. The stability bound so obtained, which is given in terms of the Riemann zeta function, is checked numerically.

FIGURES IN THIS ARTICLE
<>
Copyright © 2011 by American Society of Mechanical Engineers
Your Session has timed out. Please sign back in to continue.

References

Figures

Grahic Jump Location
Figure 1

Numerical solutions (symbols) and exact solutions u(π/2,t)=Eγ(−tγ) (lines) at the midpoint x=π/2 of the fractional diffusion problem described in the main text for γ=1 (squares), γ=0.75 (triangles), and γ=0.5 (circles). We have used Δx=π/20 in all cases and S¯=0.5 for γ=1, S¯=0.44 for γ=3/4, and S¯=0.38 for γ=1/2.

Grahic Jump Location
Figure 2

Numerical solutions (symbols) and exact solutions u(π/2,t)=Eγ(−tγ) (lines) at the midpoint x=π/2 of the fractional diffusion problem described in the main text for γ=1.25 (triangles), γ=1.5 (circles), γ=1.75 (squares), and γ=2 (stars) with Δt=0.01 and Δx=π/20

Grahic Jump Location
Figure 3

Numerical solution (circles) for the subdiffusion problem considered in Fig. 1 (γ=0.5, f(x)=sin(x), and Δx=π/20) after 1000 time steps when S=(Δt)γ/(Δx)2=0.44. Note that this value of S is larger than the stability bound S×=(1−21.5)ζ(−0.5)/Γ(1.5)≃0.429 provided by Eq. 21. The broken line is to guide the eye.

Grahic Jump Location
Figure 4

Numerical solution (symbols) for the subdiffusion problem considered in Fig. 1 (γ=0.5, f(x)=sin(x), and Δx=π/20) after 100 time steps (circles), 500 time steps (squares), 2000 time steps (triangles), 5000 time steps (stars) when S¯=0.48, that is, S=(Δt)γ/(Δx)2≃0.36. Note that this value is below the stability bound S×=(1−21.5)ζ(−0.5)/Γ(1.5)≃0.43 provided by Eq. 21 so we are inside the stability region. The solid lines are the corresponding exact solutions.

Grahic Jump Location
Figure 5

Numerical solution (circles) for the fractional diffusion-wave equation considered in Fig. 2 (γ=1.5, f(x)=sin(x), g(x)=0, and Δx=π/20) after 1000 time steps when S¯=0.77, that is, for S=(Δt)γ/(Δx)2≃0.87. Note that this value of S is larger than the stability bound S×=2−0.5(21.5−23)ζ(−0.5)/Γ(1.5)≃0.86 provided by Eq. 26. The broken line is to guide the eye.

Grahic Jump Location
Figure 6

Numerical solution (symbols) for the fractional diffusion-wave equation considered in Fig. 2 (γ=1.5, f(x)=sin(x), g(x)=0, and Δx=π/20). Hollow symbols are numerical solutions after 10 (circles), 20 (triangles), 40 (squares), and 50 (stars) time steps when Δt=0.075 so that S=(Δt)γ/(Δx)2≃0.83. This value of S is smaller than the stability bound S×=2−0.5(21.5−23)ζ(−0.5)/Γ(1.5)≃0.86 given by Eq. 26. Filled symbols are the numerical solution after 100 (circles), 200 (triangles), 400 (squares), and 500 (stars) time steps when Δt=0.0075 so that S=(Δt)γ/(Δx)2≃0.026. Solid lines are the corresponding exact solutions. Numerical solutions obtained for Δt=0.0075 are noticeably better than those obtained with Δt=0.075 because the (truncation) error is of order Δt (see Sec. 4).

Grahic Jump Location
Figure 7

Numerical values of Scrit corresponding to the onset of instability versus the subdiffusion exponent γ. The solid line is the prediction S× of the von Neumann analysis and the symbols denote the numerical results with the criterion in Eq. 35: squares for M=200 and circles for M=1000. The broken line is the stability bound corresponding to the explicit methods of Refs. 34,36,40.

Tables

Errata

Discussions

Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging and repositioning the boxes below.

Related Journal Articles
Related eBook Content
Topic Collections

Sorry! You do not have access to this content. For assistance or to subscribe, please contact us:

  • TELEPHONE: 1-800-843-2763 (Toll-free in the USA)
  • EMAIL: asmedigitalcollection@asme.org
Sign In