Efficient computation of the Wright function and its applications to fractional diffusion-wave equations

In this article, we deal with the efficient computation of the Wright function in the cases of interest for the expression of solutions of some fractional differential equations. The proposed algorithm is based on the inversion of the Laplace transform of a particular expression of the Wright function for which we discuss in detail the error analysis. We also present a code package that implements the algorithm proposed here in different programming languages. The analysis and implementation are accompanied by an extensive set of numerical experiments that validate both the theoretical estimates of the error and the applicability of the proposed method for representing the solutions of fractional differential equations.


Introduction
The Wright function is a generalization of the exponential function defined by the convergent series in the whole complex plane where Γ is the Euler gamma function.It was initially introduced by E.M. Wright for  ≥ 0 in the framework of the asymptotic theory of partitions [1].The same author then studied the case  ∈ (−1, 0), which is now referred to in the literature as Wright function of the second kind (WF2K) [2].
Although several representations of the Wright function have been introduced such as integral representation and asymptotic expansions and many of its analytical properties have already been well-studied (see, e.g.[3][4][5][6]), its numerical evaluation is still an active research area.In this paper, we discuss the numerical evaluation of WF2K, since this is the most interesting case for applications.Indeed, this function plays an important role in describing non-Gaussian deterministic and stochastic processes and the transition from sub-diffusion processes to wave propagation.Noteworthy cases of WF2Ks are the functions Figure 1.Plots of the function  −,1− (1; ) for  = 1 /2, 3 /8, 1 /4, and 1 /8.related via the formula   () =   ().They were introduced by Mainardi in the 90s and are called auxiliary functions by virtue of their roles in solving the signalling problem and the Cauchy value problem, respectively, for the time-fractional diffusion-wave equation (see [7] and references therein).Except for some special cases where WF2K can be represented in terms of other elementary and special functions, such as  0 () = exp(−),  1/2 () = 1/ √  exp(− 2 /4),  1/3 () = 3 2/3 Ai(/3 1/3 ), for Ai the Airy function (Chapter 9 from [8]), most programming languages do not provide built-in functions for WF2K.In this paper we try to fill this gap by considering the numerical evaluation of the following function obtained from (1) by setting  = −||  .As an example, in Figure 1 we show the plots of some of these functions.
Starting from the integral representation of the Wright function several algorithms for its numerical computation have already been proposed in [9,10].The authors state that these could serve as the basis for the creation of a programming package for the numerical evaluation of the Wright function.However, the actual implementation is not straightforward as it seems.The proposed integrals have oscillatory behavior and selecting efficient quadrature formulas or a reliable truncation seems to be challenging.
In this paper, we consider a different approach based on the inversion of the Laplace transform in which a trapezoidal rule is applied on a parabolic contour.This approach solves the difficulties encountered in adopting the strategies discussed in [9,10].It is worth mentioning that methods of this type have already been successfully applied in [11,12] for the computation of Mittag-Leffler functions on the real line, in using a solve-then-discretize approach for certain fractional differential equations [13] and, more generally, for the computation of actions of semigroups with certified error bounds [14].
The paper is organized as follows.In Section 2 we first recall some generalities on the computation of Laplace-type integrals.Then we focus on our case of interest and delineate the relative error analysis leveraging corresponding results for the Mittag-Leffler function from [12].In Section 3 we look at several fractional differential problems whose solution can be expressed in terms of the Wright function to substantiate some possible use cases of the algorithm presented here.In Section 4 we present some numerical experiments which confirm our analysis from Section 2. Finally, Section 5 summarizes the main conclusions of this work.

Inversion of the Laplace transform
The aim of this section is to show that an effective way to compute  , (; ) given in ( 4) is based on the numerical inversion of the Laplace transform where (see, e.g.[Eq.( 11), 3]) Here the contour  is a suitable deformation of the Bromwich line that must be fixed in a region of the complex plane in which  , (; ) is analytic.Considering that the only singularity of  , (; ) lies in the branch point  = 0, we can choose the simplest contour so far discussed in the literature, namely the parabola [15].This contour is described by the equation where  is a parameter to be taken positive so that the parabola encloses the entire negative real axis; see Figure 2. Using this parabolic-shaped contour () we can write the integral (5) as where  , () :=  () (()) By approximating this integral by means of the finite trapezoidal rule with step-size ℎ we obtain We observe that for the symmetry of the contour with respect to the real axis, only the quadrature nodes in the upper (lower) half-plane can be considered; this leads to a halving of the computational cost.In addition, for the given  and , the three parameters , ℎ and  occurring in (9) must be set to minimize the error , (; ) Before carrying out a detailed error analysis, for the reader's convenience, we recall some basic facts about the well-known error estimates for the trapezoidal rule that will be useful for such analysis.

Error estimates for the trapezoidal rule
Consider the absolutely convergent integral and its infinite and finite trapezoidal approximations For the error we have where The quantities ℰ  and ℰ  are often referred to as the discretization error and the truncation error, respectively.
As for the discretization error, when () is a complex-valued function we can use the following theorem.
Theorem 2.1.( [15], Thm.2.1) Let  =  + , with  and  real.Suppose g(w) is analytic in the strip − <  < , for some  > 0,  > 0, with () → 0 uniformly as || → +∞ in that strip.Suppose further that for some  + > 0,  − > 0 the function () satisfies for all 0 <  < , 0 <  < .Then, where When the contribution of  + and  − is negligible then the estimates ℰ + ≈  −2/ℎ and ℰ − ≈  −2/ℎ are sufficient to obtain a satisfactory error analysis.Otherwise, these terms would have to be balanced with the decay of the exponential to get good estimates.For example, this is the situation encountered when dealing with the computation of Mittag-Leffler functions on the real line [12], and when dealing with the numerical approximation of strongly continuous semigroups on infinite-dimensional Hilbert spaces [14].
As for the truncation error, if () decays rapidly as  → ±∞, we have that is, the truncation error can be approximated by the last term retained of the trapezoidal sum.Following [16], to make a more complete error analysis, we also consider the roundoff errors.At this aim we introduce the term where   are the relative errors in the computed function values (ℎ) and they all satisfy |  | ≤ , the precision machine.The formula (11) is then extended to

Selecting the quadrature nodes
We use here the results from Section 2.1 to perform the error analysis for (9).Denoting by (see (8)) the infinite trapezoidal approximation to (7), we can bound the error in (10) by means of the discretization error and the truncation error which in this case are defined by Furthermore, based on the results of Theorem 2.1 we can use a non-symmetric strip of analyticity and therefore we also have that where the quantities  − (; ) depend on the behavior of  , on the upper and on the lower half-plane, respectively.In order to determine an estimate of these two quantities, we begin by studying the behavior of | , ( + )| for 0 <  <  < 1.Following the results reported in the proof of Theorem 2 given in [12] it is immediate to verify that Concerning the term |(( + )) − | we observe that (( + )) − =  log C ((+)) − and then where log(•) is the natural logarithm, and log C is the complex logarithm.Consequently, from (8) we have where  1 = 2 1−ℜ()  +2Im() .Using Lemma 1 of [12] we obtain where  = √  1 is a positive constant which does not depend on  and Ψ is the confluent hypergeometric function of the second kind; see, e.g., Section 13.2 from [8].Therefore, from Theorem 3 given in [12] we can deduce the following result.
where M is a suitable positive constant that depends on  but not on  and for log(•) the natural logarithm.
Using this result in Theorem 2.1 we obtain Hence, from ( 13)-( 14), an error  [ℎ, ] (; ) proportional to a given tolerance  > 0 is obtained after balancing [ℎ] + (; ) with the standard estimates provided in [15].However, as already mentioned at the end of Section 2.1, a more accurate analysis also takes into account the roundoff error.The following estimate applies to this error (substitute (8) in (12) and use the symmetry to half the domain of the integral) Assuming that ℎ = ( −2 ), we can compute the integral in (18) Therefore, the integral in relation ( 18) can be neglected together with the scaling term 2/ giving us Thus, the expression for the optimal parameters , ℎ and  is determined by requesting that asymptotically as ℎ → 0. First of all, in order to simplify the notation, we denote by (see (15)) In this way the balancing of the errors (20) leads to the following equations (see ( 16), (17), and ( 19)) where ℓ = − log .Solving them, we get the expressions for the optimal parameters: Therefore, an error  [ℎ, ] (; ) ≈  is obtained by requiring  log() ≈  −  /ℎ .Then, using the value of ℎ given in (22), we select a number of nodes or, more precisely, the greatest integer less than or equal to this quantity.
Let us observe now that  is a function of  which depends also on  as we may deduce from the expression of  in (21) and by taking into account (15).We distinguish now between the three different cases to uncover the determination for the discretization parameters  , ℎ, and .When ℜ() < 2 we have  = 2. Furthermore, as regards the value of  we can determine it with the aim of minimizing the number of nodes or, equivalently, the computational cost.Since  is a decreasing function in  and where ⌊•⌋ denotes the floor function.
When ℜ() ≥ 2, the value of  must be suitably fixed considering the effects of the term () (see (15)).We start by studying in more detail what happens when ℜ() > 2. Using the relation in the expression of  for this case Consequently, by substituting this quantity in (23) we have which is a positive function that admits at least a minimum since  () → +∞ for both  → 0 and  → 1.We refer to it as c.Then,  = ⌊ (c)⌋.Using c in (25) we can compute the corresponding value of .Finally, the so obtained values of ,  and  can be used in (22) to compute ℎ and .When ℜ() = 2 using arguments similar to those just considered in the case ℜ() > 2 and taking into account that This implies that Proceeding as in the case ℜ() > 2 we can determine the value of the free parameters , ℎ, and .
Remark 2.4.To determine the parameters ℎ,  in ( 22), we evaluate the minimum of the function  (), for the different ranges of , using Brent's method.This is implemented in MATLAB ® by the routine fminbnd for constrained minimization, see Figure 3.We observe that as we enlarge the value of  the value of  decreases causing an increase of the number of quadrature nodes.
We summarize the previous analysis in the following result.
Theorem 2.5.Let  be the working precision (e.g., single, double, quadruple), and let ℓ = − log().To obtain an absolute error of order  on the approximation of the Wright function We select instead for ℜ() > 2 )︂ ⌋︃ , and use the arg min  to compute While for ℜ() = 2 )︂ ⌋︃ and use the arg min  to compute Finally, for the case ℜ() ≥ 2 we select The method for evaluating the Wright function studied so far is applied to solve some fractional differential problems.

Time-fractional evolution equations on the real line
We consider here three time-fractional evolution problems whose solution can be expressed in terms of the Wright function.These model different types of anomalous diffusion, i.e., a diffusion process with a non-linear relationship between the mean squared displacement of the particles being diffused and time.They can be obtained through several generalizations of the classical diffusion process leading to the Brownian motion; see, e.g., [17] and the references therein.
Within this framework, we need to recall the definition of the Caputo fractional derivative for  () () = d   /d  .By of it, we focus first on the evolution equation [5,18,19] for (, ) a causal function of time, i.e., (, ) ≡ 0, ∀  < 0. Specifically (27) defines a fractional diffusion equation for  ∈ (0, 1 /2], and a fractional diffusion-wave equation for  ∈ ( 1 /2, 1].We complete such equation to a Cauchy problem by posing {︃ (, 0 + ; ) = (),  ∈ R, (±∞, ; ) = 0,  > 0; (28) whenever  ∈ ( 1 /2, 1] the initial values of the first time derivative   (, 0 + ; ) = () must be also imposed.The solution to this problem can be expressed by convolutions integrals of the functions in (28) with a characteristic function, usually called Green's function.This is expressed in terms of the Wright function [3,18] (see ( 2)) where we have denoted the Green's function by   (, ), i.e., the solution for () = () the Dirac delta.We then express the general solution of the Cauchy problem (Equations ( 17)-( 18) from [20]) as (, ; )( − ) for a primitive in time of   .Differential equation ( 27) could be also completed with the following initial and boundary conditions {︃ (, 0 + ; ) = 0,  > 0, (0 + , ; ) = ℎ(), (+∞, , ) = 0,  > 0. ( This gives us the signalling problem.To solve it, we introduce the Green's function, i.e., its solution for ℎ() = ( + ), As for the Cauchy case, to ensure the continuous dependence of the solution on the order parameter , we assume () = 0 for  ∈ ( 1 /2, 1], and express it as (Equation (2.6) from [18]) The last problem we focus on whose solution can be expressed in terms of the Wright function is the fractional heat conduction in nonhomogeneous media under perfect thermal contact from [21,22].Let us consider two infinite one-dimensional rods joined at  = 0 with perfect thermal contact, that is, with the same temperature and heat fluxes through the contact point.By using again the Caputo fractional derivative from (26) such problem is then stated as for  1 ,  2 the generalized thermal conductivities of the two semi infinite rods and the Riemann-Liouville derivative of fractional order , which is understood as the Riemann-Liouville fractional integral whenever its order turns out to be negative.The solution of (32) with initial condition  = 0,  1 =  0 ( − ),  > 0,  > 0, for the case  =  can be expressed in terms of the Wright function as )︂]︂ ,  ≥ 0, where  = 1 √ 2 /2 √ 1.

Codes and numerical examples
We use this section to numerically validate the results discussed in Sections 2 and 3, and to briefly discuss the implementation of the algorithms for computing the Wright function contained in the repository github.com/Cirdans-Home/mwright. For ease of use, the repository contains a MATLAB ® implementation of the core algorithm in double precision.For better performances and portability, implementations are available in a Fortran module for single, double, and quadruple precision.To test our proposal also in the cases for which alternative representations of the Wright function are not available, we have included for benchmark purposes a C implementation with arbitrary precision ball arithmetic [23] of the series representation (1).
We make use of the implementation in Section 4.1 to illustrate the error analysis for the inversion of the Laplace transform.Then, in Section 4.2 we analyze its usage for expressing the solutions of the problems discussed in Section 3.
All the experiments were performed on a Linux machine with an Intel ® Core i7-8750H CPU @ 2.20 GHz and 16 GB of memory.The m code are run in MATLAB ® version 9.6.0.1072779 (R2019a).Fortran and C codes are compiled with the gnu/9.4.0 suite.To use the single, double and quadruple precision in the Fortran environment we employ the iso fortran env intrinsic module that provides the needed constants, derived types, and intrinsic procedures.

Computing the Wright function
We consider first the Mainardi functions   (||), that are the cases of the form  , (−||),  = −,  = 1−, for which we have a closed form relation (3).For these cases, we report in Figure 4 the convergence with respect to the number of nodes  of the approximation in (9) to these close form expressions by looking at the relative error with respect to the exact representation.The last value of  is the one given in (24) from which we obtain the expected value.
For the two cases in which the function can be represented as an exponential, we can also test the routines for calculating the quadruple precision, see Figure 5.
To evaluate the accuracy in the general case, we make use of the extended precision ball arithmetic library [23].For this case we perform a comparison also for the case in which the parameter  ∈ C. We consider only the number of quadrature points determined by (23), and report the results in the heatmaps in Figure 6.
Each of the boxes in the plots considers the norm-wise relative error in 2-norm over the interval [−5, 0], and the color map is on a logarithmic scale.In all the cases we reach at least a relative error of the order of 10 −10 for the double precision routine.

Solving the time-fractional evolution equation
In this section, we deal with the approximation of the exact solution, represented by convolution with the Green's function, of the differential problems involving fractional derivatives in time discussed in Section 3.
We focus first on the analysis of the solution of the Cauchy problem in (30).To calculate the solution we need to perform a convolution between the values obtained with the computation of the Wright function and the initial condition through a fast Fourier transform.Figure 6.Heatmaps for the relative error (color map is in log-scale) between the values computed with the series computed to the 1000th entry in augmented 128 bits precision and the values computed with our algorithm for the cases in which a complex  is used.The number of quadrature points has been fixed as in the analysis in Section 2.2. to produce the solution for the case of  ∈ (0, 1 /2].The expected sub-diffusion behavior is qualitatively observed with respect to the case of the classical solution ( = 1 /2).That is, the convergence of the solution to the zero stationary state is slower, the more the order  approaches 0.
We observe that to obtain an approximate solution using a finite difference method, rather than a method based on a weak formulation, it would be necessary to set a domain truncation and study the decay to zero of the function to keep the error under control while also having to solve the associated linear system.Moreover, methods for time-stepping fractional differential equations have usually a low approximation order and need to way between the steps the time due to the treatment of the integration queue, which is instead an increasing quantity with the number of steps.
The other numerical experiment we perform is the heat conduction problem in (32).In Figure 9 we draw the solution (33) for the same physical parameters, and different values of the fractional order  computed by means of the MATLAB ® version of the Wright function.
We observe again that also, in this case, wanting to solve the problem with a discretization approach, we would face the same problem of having to discuss the truncation of the domain together with that of having to integrate the two fractional differential equations in time.It should also be noted that in this case, the coupling conditions are two fractional integrals, this makes the associated linear systems slightly more complex than those of the previous case.

Conclusions
In this work, we have proposed an algorithm for the computation of the Wright function of the second type based on the inversion of the Laplace transform.As far as we know, this is the first available and open-source code to effectively evaluate this function.The proposed method appears to be also a valid tool for solving fractional diffusion-wave equations thanks to its modest computational cost and high accuracy.

Figure 3 .
Figure 3. MATLAB ® listing of the optimization for  (left panel), values obtained through the optimization procedure together with the corresponding value of  (right panel) when double precision is used, and we pose  = 10 −15 and  = 2.2204 × 10 −16 .The dashed line represents the threshold value ℜ() = 2.