Sinc-Approximations of Fractional Operators : A Computing Approach

We discuss a new approach to represent fractional operators by Sinc approximation using convolution integrals. A spin off of the convolution representation is an effective inverse Laplace transform. Several examples demonstrate the application of the method to different practical problems.


Introduction
During the past 300 years a large number of fractional operators were used by researchers to overcome the Leibniz's problem of fractional differentiation [1].One of the contributors in this field was Abel with his famous integral equation [2].The large numbers of novel applications of fractional differential equations in physics, chemistry, engineering, finance, and other sciences that have been developed in the last few decades created the demand to have numerical procedures to support the solution of these problems.In fact, the problem is that in all cases of fractional operators we must be able to deal with singular operators, and more specifically, with weak singular operators.In this paper we deal with only a few of these operators connected with their inventors Riemann, Liouville, Weyl, and Caputo [3].Based on these integral operators we define fractional integral and differential equations which turn out again to be singular integral equations using a singular convolution kernel of Abel-type [3].This type of singular kernels is the main obstacle of numerical methods in the field.
There are different approaches to overcome the problems with singularities.However, it turns out that most of the methods used to circumvent the problems are not effective.Standard methods for the numerical calculation of fractional integrals and derivatives are typically slow and memory consuming due to the non locality of the operators [4][5][6][7][8].For example Lubich proposes a discrete convolution procedure based on a linear multistep method [9].Yuan and Agrawal have proposed a more efficient approach based on Laguerre integral formulas for operators whose order is between 0 and 10.Their approach differs substantially from the traditional concepts [10,11].However, the accuracy of the results can be poor.Diethelm modified the algorithm, adapting it to the properties of the problem, and demonstrated that this leads to a significantly improved quality [12].However, both approaches use the same idea to convert the original problem to a system of first order initial value problems which seems to be an approach with a low convergent rate.In fact all the methods discussed in literature have polynomial convergence compare [8,13,14] and references therein.
In this paper we follow a new direction of approximation theory introduced by [15] in connection with convolution integrals [16].The Sinc methods, e.g., [15,16] have among other features a very fast convergence of exponential order exp ( −CN 1/2 ) with C a positive constant and in addition converges even when the explicit nature of singularities is unknown.Popular numeric solvers, i.e., finite difference and finite elements, have severe limitations with regard to speed and accuracy, as well as intrinsic restrictions making them ill fit for realistic, less simplified problems.Spectral methods also have substantial difficulties regarding complexity of implementation and accuracy of solution in the presence of singularities.However, in the case when the nature of the singularity is known, then suitable boundary element bases can be selected to produce an approximation that converges rapidly.On the other hand, it is often difficult to deduce a priori the nature of the singularity, and in such cases methods based on polynomial or finite element approximations converge very slowly.
We will restrict our discussions to the practical application of the Sinc method to different fractional differential equations (FDE).We note however that, apart from the approach that we will use, several other numerical methods for treating integrals, derivatives, and solutions of differential equations of fractional order have been developed.The most prominent one-the decomposition method attributed to Adomian [17]-has rather poor convergence properties in general [18,19].
The variational iteration method is based on the use of restricted variations combined with correction functional and the Lagrange multiplier technique.The method rely on the knowledge of the Lagrange multiplier which for non-linear problems is not directly accessible [20,21].The homotopy analysis method [22] and the homotopy perturbation method [20] are slow in convergence in connection with FDEs.The generalized differential transform method [23] and a few others deliver poor solutions and moreover, the convergence proofs given apply in very restrictive conditions [24].There are many other methods known, which also do not converge satisfactorily [25].This situation creates a demand for an effective method which is able to deal with singularities, which is highly accurate, and which converge rapidly.
The integrals we shall discuss in this paper are Riemann integrals which we assume to be locally and absolutely integrable on R + .Let B = [a, b] with (−∞ < a < b < ∞) be a finite interval on the real axis R. The Riemann-Liouville fractional integrals D −ν a,x and W −ν x,b of order ν ∈ C (Re(ν) > 0) are two typical examples of integrals we shall deal with and with where Γ(ν) is the gamma function.In literature [3], these two integrals are commonly called left-and right-sided fractional integrals.The common feature of both types of integrals is that they represent convolution integrals using a singular kernel of the Able-type; i.e., t −µ with 0 < µ ≤ 1 , as a typical example.Our aim in this paper will be the approximation of integral equations (IE) with such singular convolution operators.
The two integral operators D −ν a,x and W −ν x,b can be generalized to a non-linear version known as non-linear Volterra-Hammerstein integral equations of convolution type [7] L(u(x)) = ∫ x 0 G(x − s)K(s, u(s))ds + g(x) and x ∈ (a, b) ( where L, K, g, and G are explicitly known continuous generic functions and where u is unknown.Note L might be linear or non-linear while K is a generic non-linear function.We shall in addition to Equations ( 1) and ( 2) also consider such equations of this paper.

Approximation Method
This section introduced the basic ideas of Sinc methods.We shall introduce some ideas of Sinc approximation relevant to our paper [26] .We omit most of the proofs of the different important theorems and lemmas because these proofs are available in literature [16,[29][30][31].The following subsections collect information on the basic mathematical functions used in Sinc approximation.We introduce Sinc methods to represent indefinite integrations and convolution integrals.These types of integrals are essential for representing the fractional operators of differentiation and integration [29].

Sinc Bases
To start with we first introduce some definitions and theorems allowing us to specify the space of functions, domains, and arcs for a Sinc approximation.In the calculations which follow we use the conformal maps φ(z) = log((z − a)/(b − z)) and its inverse ψ(z) = (b + a exp(z))/(1 + exp(z)) for finite intervals (a, b).Other conformal maps for different arcs Γ can be found in [31].Note Burchard and Höllig [32] prove that there is no method that converges faster than Sinc, for approximating functions that are analytic on the interior of an arc, but may have singularities of unknown type at the end points.This argument is supported by the recent findings on Lebesgue measures for Sinc approximations [33].
Now let the positive numbers α and β belong to (0, 1], and let M M M α,β (D) denote the family of all functions g ∈ Hol Hol Hol(D), such that g(a) and g(b) are finite numbers, where g(a) = lim z→a g(z) and g(b) = lim z→b g(z), and such that u ∈ L L L α,β (D) where The two definitions allow us to formulate the following theorem for Sinc approximations.
Theorem 1. Sinc Approximation [16].Let u ∈ L L L α,β (D) for α > 0 and β > 0, take M = [βN/α], where [x] denotes the greatest integer in x, and then set m , and if h = (πd/(βN )) 1/2 then there exists a positive constant c 2 independent of N , such that where w j denote the basis functions defined by with S(j, h) the shifted Sinc based on the Sinc function The proof of this theorem is given in [16].Note the choice h = (πd/(βN )) 1/2 is close to optimal for an approximation in the space M M M α,β (D) in the sense that the error bound in Theorem 1 cannot be appreciably improved regardless of the basis [16].It is also optimal in the sense of the Lebesgue measure achieving an optimal value less than Chebyshev approximations [33].
The discrete shifting in Equation ( 8) allows us to cover the approximation interval (a, b) in a dense way while the conformal map is used to map the interval of approximation from an infinite range of values to a finite or a semi-infinite one, or even to an analytic arc.
A Sinc approximation is then defined as follows.Select positive integers N and M = [βN/α] so that m = M + N + 1.The step length is determined by h = (πd/(βN )) 1/2 where α, β, and d are real parameters.In addition assume there is a conformal map φ and its inverse ψ such that we can define Sinc points z j = ψ(jh), j ∈ Z [16].This allows us to define a row vector w w w T m of basis functions with w j defined as in Equation ( 7).For a given vector The superscript "T" here and in the following denotes the transpose.We now introduce the dot product as an approximation of the function u by If it is obvious from the context, we also will use the notation where the subscripts m are dropped, i.e., dot products, when there is no misunderstanding.Based on this notation, we will introduce the different integrals we shall require in the next few subsections [29].

Indefinite Integral
As a basis for the representation of fractional derivatives, we need indefinite integrals.This subsection describes how indefinite integrals can be defined [16] and how these definitions are related to the definition of definite integrals.Let us now discuss the collocation of the integrals J + and J − defined by Let Z denote the set of all integers, and let C denote the complex plane.Let Sinc(z) be given by Equation (9) and let e k be defined as with and Si(x) = ∫ x 0 t −1 sin(t)dt the sine integral.If now φ denotes a one-to-one transformation of the interval (a, b) onto the real line R as defined above, then let h denote a fixed positive step length, and let the Sinc points be defined on (a, b) by z k = φ −1 (kh), k ∈ Z, where φ −1 denotes the inverse function of the conformal map φ.Let M and N be positive integers, set m = M + N + 1, and for a given function u defined on (a, b), define a diagonal matrix D(u) by D(u) = diag [u (z −M ) , . . ., u (z N )].Let I (−1) be a square Töplitz matrix of order m having e i−j , as its (i, j) th entry, and i, j = −M, . . ., N , then ( Define square matrices A + m and A − m by The collocated representation of the indefinite integrals are thus given by These are collocated representations of the indefinite integrals defined in Equations ( 12) and ( 13), respectively [29].Quite recently it was proved that the eigenvalues of the matrix I (−1) lie in the open right half-plane [37].Note the matrices A ± m arises in several applications of Sinc methods and have positive eigenvalues on the open right half-plane [29].The matrices (±A ± m ) −1 yield accurate approximation of derivatives, i.e., their inverses exist, since the real parts of their eigenvalues are on the open right half plane, and so Comparing our results with the target equations for example Equations ( 1) and (2), we observe that we need one additional step to extend the indefinite integrals to convolution integrals.

Convolution Integrals
Indefinite convolution integrals can also be effectively collocated via Sinc methods [16].This section discusses the core procedure of this paper, for collocating the convolution integrals and for obtaining explicit approximations of the functions p and q defined by and where x ∈ Γ.In presenting these convolution results, we shall assume that Γ = (a, b) ⊆ R, unless otherwise indicated.Note also, that being able to collocate p and q enables us to collocate definite convolutions like ∫ b a f (x − t)g(t)dt.Recently such convolution integrals were determined by accurate tensor approximation [27,28].
Before we start to present the collocation of the Equations ( 21) and ( 22), we mention that there is a special approach to evaluate the convolution integrals by using a Laplace transform.Lubich [6,7] introduced this way of calculation by the following idea e st g(x − t)dt ds (23) for which the inner integral solves the initial value problem y ′ = sy + g with g(0) = 0.The Laplace transform with Furthermore, if the Laplace transform Equation ( 24) exists for all ℜ(s) > 0, then F (A ± m ) are well defined.
In view of Equation ( 19), we can expect that the approximations are also accurate, which is, in fact, the case [29].The procedure to calculate the convolution integrals is now as follows.The collocated integrals take the form where and with Σ = diag [s −M , . . ., s N ] the diagonal matrix for each of the matrices A ± m .Then the Laplace transform Equation (24) delivers the square matrices F (A ± m ) defined via the equations We can then get the approximation of Equations ( 25) and ( 26) from the expression These two formulas enable us to obtain vectors whose entries approximate p and q at the Sinc points.The convergence of the method is exponential as was proved in [16].
For the non-linear case of convolution with the target to represent Equation (3), we state the following result, taken from [31]: then P ∈ L α ,β .In the above notation, with F denoting the Laplace transform of G, we write i.e., Let M = N and h = (πd/(αN )) 1/2 , then there exists a constant c 3 which is independent of N such that and w w w T as defined in Equation (10).
A similar formulation can be given for The proof of this theorem can be found in [31].The use of Theorem 2 allows us to write for the non-linear convolution integral an approximation at Sinc points x k as follows where F kj is the (k, j) entry of matrix F (A + m ) and the error at a certain Sinc point 34) will be useful later on to estimate the error in a non-linear problem (see Theorem 8).

Sinc Representation of Fractional Operators
This section presents the main operators which are frequently used in practical problems.We present the basic notations first and then use them in collocated approximations.The collocation formulas have a common structure for linear operators.The non-linear version only deviates in the generic function K but also allows a separation of spectral and discrete properties of the operators.

Fractional Integrals
The Riemann-Liouville fractional integrals Equations ( 1) and ( 2) defined on the finite interval (a, b) of the real line R, are naturally extended to the half-line R + .The two integrals for the positive real line have corresponding to the representation of Miller and Ross [38] the forms: and Both integrals are of convolution type.We shall describe below the approximation of these integrals by use of our convolution Sinc methods.
Here N denotes the set of natural numbers.The Riemann-Liouville fractional derivatives on the real line D µ a,x and W µ x,b of order µ with µ ∈ C and ℜ(µ) ≥ 0 are defined by and These two definitions use Riemann-Liouville integrals to represent the derivatives.We again omit the straight forward extensions of the integrals Equations ( 37) and ( 38) to R + -see [3].
In some of the applications it is important to incorporate the initial conditions into the fractional derivative.For such cases another type of fractional derivatives was introduced by Caputo [34].Kilbas et al. define these types of derivatives as follows and The difference between the Riemann-Liouville and the Caputo derivative is that the integer order derivative is applied to the integrand of the integral.The Caputo derivative has proven to be an important tool in mathematical modeling of many phenomena in physics and engineering, see, e.g., [13,14,35] and the references therein.

Fractional Integral Equations
In this section we consider converting an fractional integral equation (FDI) to an integral equation (IE).Consider, for example, the simple Riemann-Liouville equation, Using Equation (35), this equation can also be written in the form which is just a linear Volterra equation of the form with G(x, t) = (x − t) ν−1 / Γ(ν) .A generalization to non-linear FIE is straight forward if we replace in Equation ( 41) the term D −ν 0,x u(x) with a generic non-linear function which is just The integral Equation ( 45) is a Volterra integral equation which is also known as Hammerstein's equation.

Fractional Differential Equations
Next let us discuss the formulation of linear fractional differential equation using the Riemann-Liouville fractional derivative D µ a,x u(x) of order µ > 0 in a Banach space of continuous functions [3].The Cauchy type problem for the fractional differential equation of order µ > 0 is given by and with the initial conditions where ⌈µ⌉ denotes the smallest integer greater than ℜ(µ) , and where we assume that denoting the weighted space of continuous functions and 0 ≤ γ < 1 is equivalent to the solution to the initial value problem Equations ( 46)-( 47) is thus given by The proof that Equation ( 48) is solution of Equations ( 46) and ( 47) is given by Kilbas et al. [3].
The formulation of the Cauchy type initial value problem using the Caputo fractional derivative c D ν a,x u(x) for µ > 0 and µ / ∈ N can be formulated in a similar way for such a simple equation as with initial conditions This problem again is equivalent to the solution of a Volterra integral equation of the second kind as shown in [3] The two types of fractional differential equations have the nice property that their solution is represented by convolution integrals.Such integrals can be easily dealt with our indefinite convolution and Sinc approximation procedures.Although we know the solutions Equations ( 48) and ( 51) of the FDEs Equations ( 46) and ( 49), we still have to solve a Volterra integral equation of the second kind.In the next section we shall describe our method of obtaining this solution.
Problems exist in applications which contain both fractional and ordinary derivatives [37].Such problems can also be readily dealt with using our methods.
So far we discussed only linear fractional differential equations.However, in some of the applications non-linear models are used [8].Another example of non-linear fractional differential equation was examined in [42,43].Diethelm discusses the numerical algorithm for the following fractional equation combined with initial conditions This initial value problem is equivalent to the Volterra integral Equation , [43] A generalization of this equation to a system of equations is given by Daftardar [39].Again we end up with a Volterra integral equation which uses an integral of convolution type.As we already discussed above, the convolution property is very useful in dealing with numerical solutions of these equations.All the FIEs and FDEs discussed so far are treated in literature with "standard" methods.The results gained by these approaches are limited to a certain accuracy combined with a high demand of computational efforts.In the following we will show how the accuracy can be improved and the computational tasks reduced to a minimum.
Let us introduce the discrete Sinc representation of integral equations of the different types discussed above.Using Equations ( 21) and (22) as basis, we know how to collocate convolution integrals in Volterra integral equations of second kind by applying Sinc methods.

Collocation of Fractional Integral Equations
Concerning the discrete representation of the Volterra integral Equations ( 41) and (43) we have the approximation The non-linear version of the FIE Equations ( 44) and ( 45) has the collocated representation The solutions in both cases follows by In the first case we have to solve a linear system of equations while the second case needs the solution of a non-linear system of equations.Which in both cases can be numerically achieved by Newton's method and Newton-Kantorovich iteration.

Collocation of Fractional Differential Equations
For fractional differential equations there exists a solution in terms of a second kind Volterra integral equation incorporating the initial conditions of the Cauchy problem as given in Equations ( 46) and (47).The discrete version of these equations ensue by collocating Equation (48) as follows: N collects initial conditions of the problem.The discrete version follows straight forward from Equation (58) to be Caputo's fractional differential Equation ( 51) allows the solution in a collocated Sinc representation as N collects initial conditions of the problem.The discrete version follows from Equation (60) to be Solving this linear system in both cases with respect to V V V m (u) allows us to approximate the solution by For the non-linear singular Volterra equation we can derive the determining equations based on Theorem 2. Taking into account the exponential decay of the error E(G, K) we can write Equation (3) in its discrete form as This non-linear system of Equation (63) delivers the unknowns V V V (u) by applying Newton's method and Newton-Kantorovich iteration method.The resulting approximation then enables us to represent the solution via Sinc approximation in the form Let us next proof the convergence of the non-linear integral Equation (3).Let D and L α,d be defined as above, corresponding to an interval (a, b).Let ν ∈ (0 , 1] , and let B ν denote the Banach space We wish to consider proving the existence of a solution to the IE Assumption 3. Let us assume that corresponding to some interval (C , D) to be made more specific below, that K(• , c) is analytic and uniformly bounded for all c ∈ (C , D) , that K(x , •) is analytic and informally bounded in D C,D , and also that g ∈ L ν ,d .
Let J + be defined on (a, b) as above, so that we can write the IE Equation (65) in the operator form where we have used the simpler notation J for J + .
Proof.Note to begin with, that since ν ∈ (0 , 1] , and since the eigenvalues of J + have non negative real parts and are located either on the real line or occur in conjugate pairs on the right half plane, the same is true for the eigenvalues of the operator J ν .The proof then follows by inspection. Lemma 5. Let the conditions of Lemma 4 be satisfied.66) has a unique solution.
Remark 1.For a finite interval (a , b) , we have for u ∈ B ν , that It thus follows from Remark 1 that the operator J ν K(• , •) is a contraction mapping for all |b − a| sufficiently small, implying the statement of Lemma 5.
and let us assume that Theorem 7. If Assumption 6 is satisfied, then the sequence {u n } ∞ n=0 starting with u 0 as in Equation ( 69) and with converges to the unique solution on (a , b) × (C , D) of Equation (66).
Next, we have, with Since P is decreasing on (a, b) , we have Q 0 − Q 1 < 0 , and since u 0 > 0 , we have u 2 − u 1 < 0 on (a , b) .It then follows by induction, that u n − u n−1 < 0 , i.e., that the sequence {u n } ∞ n=0 is a positive, and decreasing sequence.This sequence must therefore have a unique limit.
The above argument also holds in the case when u 0 is negative, or has sign changes on (a , b); in this case we simply replace the equal sign in Equation ( 72) by an inequality.Part 2. Consider now the case when P (x , •) is decreasing on (C , D) .Once again assuming that u 0 (x) > 0 on (a , b) , we have u 1 = (T 0 ) −1 u 0 < u 0 .Next, with the same definitions as in the above Part 1 proof, we now have Q 0 − Q 1 > 0 , implying that u 2 > u 1 .It thus follows from Equation (72) that we now get a nested sequence That is we again get convergence to a unique solution.This same argument again holds in the case when u 0 is negative, or when it changes sign on (a , b) ; by replacing the equal sign in Equation ( 72) with an inequality we then just get a monotonically convergent sequence of positive numbers, as we did in Part 1 of the proof.
Remark 2. We expect Theorem 7 to hold only under the assumption that ∂K(x , •)/∂y is strictly negative on (a , b) .One can imagine in such a case, the occurrence of a periodic, non-convergent sequence {u n } ∞ n=0 , but this is doubtful under our analytic assumptions on K .This is particularly the case when we are using the approximating operators J ± m , since the eigenvalues of the indefinite integration matrices of these approximating operators have strictly positive real parts.
To analyze the error of the approximation Equation (64), we assume on the involved functions L, G, and K to be analytic.Then the following Theorem delivers an estimate of the error.

Theorem 8. Non-linear Volterra Equation.
Let us assume that the functions g, L, G, and K are analytic and satisfy the assumptions of Theorem 2, and let u k , k = −N , ..., N denote the solutions to the approximating Equation (63).Then there exists a constant c 4 that is independent of N , such that Proof.To start with the proof for the error we substitute Equation (34) in Equation ( 3) to obtain and inserting the approximations given by Equation (64) into Equation (75), we get Using now continuity of L, we write for the local error at Sinc points Introducing the maximal error by ϵ = max This means The final step of the proof is based on the spectral properties of F (A + m ).The proof of the error formula becomes valid if we show that the spectral radius of any matrix A satisfies r σ (A) = max For the general interval (a, b) the following holds , showing that the norm of A + is proportional to the length of the interval b − a , which allows the estimation We have ample evidence from numerical computations that the spectral radius condition is satisfied for A + m and F (A + m ) on a standard interval (a, b) = (0, 1).Numerical calculations and Equation (81) suggest the decay r σ (A + m ) = λ 1 / √ N with N = 1, 2, . . ., where λ 1 is the largest eigenvalue of σ (A + m ).The value λ 1 can be determined by direct computation using the smallest matrix A + m ; i.e., N = 1.The value for the spectral radius is bound by r σ (A + m ) < 1/3 and the analytic computation using the smallest matrix delivers r σ (A + m ) ≈ 0.33227 . . .which is close to the given upper numerical bound.
The decay of the spectral radius r σ (F (A + m )) ≃ λ 1,ν N −ν/2 and λ 1,ν = exp(−ν) for ν > 0. Again the relations hold for a standard interval (a, b) = (0, 1).The important point here is that the spectral radius for both matrices A + m and F (A + m ) are smaller than 1.In this case the Neumann Lemma can be applied with and thus the error formula holds.Thus finally we get with c 4 a positive constant independent of N .Even more, since the spectral radii of both matrices are less than 1, successive iterations of the non-linear equation will converge to a fixed point according to the fixed point Theorem.This concludes the proof.
The result is that even for non-linear Volterra-Hammerstein integral equations the error in a Sinc approximation decays exponentially.
The different formulas presented in this section show that a Sinc approximation is basically derived from a collocation of the integral equations.The numeric steps of all of the calculations are straight forward and do not need any specialized numerical procedure.The core of the method is contained in the use of Sinc functions as a basis in connection with conformal maps.The formulas also show that an implementation uses standard objects like vectors and matrices which are simple to handle in any kind of programming language.However, since the method involves some symbolic steps we decided to do the implementation in Mathematica.A purely numeric implementation of Sinc methods in Matlab is provided by [29].

Examples
The examples presented in this section use implemented functions identifying the integral type used; e.g., RiemannLiouville and the type of the equation.I stands for integral equation FIE for fractional integral equation, NonlinearFIE[] for the non-linear version, etc.We implemented a set of functions using the discrete representation discussed in Section 3. In the following we will apply these functions to different problems and derive the related approximation of the solution.We start with simple applications like differentiation and integration, then solve linear integral equations, and finally deal with non-linear problems.

Generalized Gamma Function
In applications of fractional calculus it turned out that the generalization of the gamma function takes a prominent role [3,44,46].The frequent occurrence of this function asked for an effective way to generate its numerical values.The difficulty of the generalized gamma function is that for arguments x → 0 the functions possesses a singularity which causes numerical problems.One possible approach to deal with such a problem is to generate the generalized gamma function applying Sinc methods to the corresponding Riemann-Liouville integral of order ν The practical calculation is carried out by first defining in Mathematica the exponential function as f (x_, a_):=e ax which is related to the generalized gamma function after integrating by the following relation The actual numerical integration of the exponential is performed by our Mathematica function RiemannLiouvilleI[] using f , the order ν = 1/3, and N = 32 the number of Sinc points on an interval x ∈ (0, 1).The results of this integration for a specific choice of a = 1 is shown in Figure 1.In addition to the graph of the resulting function, we also plotted the location of the Sinc points.The distribution of the Sinc points gives us an idea how the discrete arrangement of the mesh along the x-axis is generated.The lower part of Figure 1 shows the local error of the approximation.We compare here our numeric results with the symbolic implementation of E x (a, ν) in Mathematica.The results displayed in Figure 1 show that the function possesses a vertical tangent if x → 0 causing numerical problems which need special treatments in "standard" approaches.In our calculations we do not need to care about such singularities at the end points because Sinc methods are able to deal with such behavior as an inherent property.

Extended Fractional Relaxation
In a paper examining primarily protein dynamics Glöckle and Nonnenmacher [44] derived a fractional relaxation equation possessing an analytic solution.It is well known today that for the generalization of a relaxation process which shows an exponential decay using the fractional counterpart is delivering as solution the generalization of the exponential function, a Mittag-Leffler function.In the same paper the two authors conjectured that the relaxation equation can be generalized to a driven fractional equation.At the time when the paper was written there was no method available to deal either analytically nor numerically with such kind of equation.In the next few lines we will show how such kind of driven fractional differential equation can be solved by Sinc methods using instead of a Rieman-Liouvill operator a Caputo fractional derivative.
The fractional differential equation which was solved by Glöckle and Nonnenmacher is given by which is equivalent to the fractional integral equation Equation (84) allows the general analytic solution where E β (t) is the Mittag-Leffler function.Selecting the initial value and the relaxation constant as ϕ 0 = 1 and τ 0 = 1/10, we can straight forward approximate the solution by using the function RiemannLiouvilleFIE[] delivering for N = 32 Sinc points an accurate representation of the solution The approximation and the local error of the approximation is shown in Figure 2. The local error (right panel) in Figure 2 is below an acceptable value of 10 −4 .If the number of Sinc points is increased this error decays exponentially.
Glöckle's representation of the fractional differential equation includes the right hand side of the solution which is the sum over initial conditions.The term ϕ 0 t −β / Γ(1 − β) is generated in the implemented solution; i.e., related to the initial value and thus not part of the Cauchy problem.Knowing this fact we are able to generate the approximation in two different ways using a function based on non-linear Caputo fractional differential equations.The first approach uses the linear term τ −β 0 ϕ(t) as argument and treats the exogenous term separately.
The second approach uses the exogenous term as part of the equation.If we use the part of the initial values we can suppress the initial value and use the given representation as The results of both approximations are shown in Figure 3.The top panel shows the two approximation approaches in addition to the exact solution.The bottom panel represents the local errors of the two different approximations.It becomes obvious that the second approximation is by 4 decades more accurate than the first approach.However, in both approximations the error is sufficient for an accurate representation of the exact result.
In Figure 4 we display approximations of the solution with varying β and fixed relaxation time τ 0 and initial condition ϕ 0 .The approximation is performed by a Caputo fractional operator.
The generalization Glöckle and Nonnenmacher suggested is a fractional differential equation of the form where E(t) is a harmonic function.Let us assume that E(t) = sin(π t) representing for example an external electric field applied to the relaxation process.The solution for such kind of equation is again generated by using the Caputo representation of the fractional derivative.An example for the selection of the amplitude χ 0 = 1/10 is shown in Figure 5.The graphs shown in Figure 5 are the driven solution (solid line) and as a reference the solution of Equation ( 88) with χ 0 = 0 using the same relaxation parameters.It is obvious that the driven solution follows the solution of the relaxation equation.The oscillations around the relaxation solution decay for larger times t.Note Figure 5 uses a logarithmic scale for ϕ(t).The initial condition and the relaxation time τ 0 , the amplitude of the driving force χ 0 , and the fractional order β was set to ϕ 0 = 0, τ 0 = 1/10, χ 0 = 1/10, and β = 1/2.The dased line represents the solution of the relaxation Equation (84) using the same parameters.The Sinc approximation of the Caputo representation used N = 64 Sinc points.
Another way to generalize the relaxation equation is to include a nonlinearity into the equation.One interesting nonlinearity is a periodic function like in a mathematical pendulum.The equation without exogenous term may be written as Figure 6 shows the approximation of Equation (89) in connection with the original linear model Equation (84).The two curves in Figure 6 represent the original linear model Equation (84) (far right curve) and the non-linear model Equation (89) (far left curve), respectively.It is obvious that the difference between the two models is a shift along the t-axis possessing the same asymptotic scaling behavior.This example demonstrates that fractional differential equations can be solved effectively with high accuracy if we use Sinc approximations.The variation of a model by extending it with exogenous or non-linear terms is in fact no problem for the used approximation method.In addition the error of the calculations decays exponential allowing a high precision representation of the approximation.

Fractal Filters
Modern design of electronic filters, such as low-pass, high-pass, and bandpass filters includes the construction of a rational function which satisfies the desired specifications for cutoff frequencies, pass-band gain, transition bandwidth, and stop-band attenuation.Recently these rational approximations were extended to fractal filters improving the design methodology in some direction.Following Podlubny [35] a fractional-order control system can be described by a fractional differential equation of the form: or by an equivalent continuous transfer function The use of the continuous transfer function is nowadays the standard approach to design filters.The transfer function Equation (91) uses however fractional exponents α i and β j to represent the fractional orders of the differential equations.One idea originally introduced by Roy is that the irrational transfer function G(s) can be replaced by a rational approximation G(s) [48].Such kind of rational approximation can be performed by a Sinc point based Thiele algorithm [31].Thiele's algorithm basically converts a rational expression to a continued fraction representation which is related to a rational expression possessing exponents of natural numbers.Such kind of approximation will change of course the representation of the rational character of the transfer function but allows us to use a fraction with integer orders if we need to invert the Laplace transformation to the original domain.The advantage here is that for rational functions with integer powers there exists a well known method by Heaviside to invert a Laplace transform.
To demonstrate the procedure let us consider a simple fractional relaxation equation.The equation for a fractional relaxation process is given by where 0 < q ≤ 1 is a positive number and u 0 is the initial condition.The Laplace transform of the above equation delivers the algebraic equation the Laplace transform of u is denoted by L s t u(t).The solution of this equation in Laplace space follows by solving the above equation with respect to the Laplace representation of u: Inverting this transfer function using inverse Laplace transforms and Mellin transforms [47] will deliver the exact solution The derived function G(s) using a Laplace transform from is in terms of the more general systems theory the transfer function of the system with q ∈ R. The definition of G allows us to perform the steps needed to approximate the irrational fraction by a rational one.First let us define G by The number of Sinc points used in the calculation are which also defines the step length for the inverse conformal map The Sinc points are generated next and used to map them to the basis points {(s i , G (s i ))} n i=−n which are used in Thiele's algorithm to generate the rational approximation.The parameters {u 0 , τ, q} used in the transfer function are given for this example by

}
The data points of the transfer function G are then generated by the mapping A rewriting of this continued fraction to a rational function is performed next.In addition we remove small quantities from this expression so that finally the rational expression follows cf2 = Simplify[cf2]//Chop  The result is an approximation of the transfer function G(s) ≈ G(s) = p m (s)/q n (s) where p m and q n are polynomials of order m and n satisfying m ≤ n , and m, n ∈ N which is a prerequisite of Heaviside's theorem.The approximated fraction G(s) and the original expression for the transfer function G(s) is shown in Figure 7.The left panel shows the original irrational fractional transfer function G(s) and its approximation using Thiele's algorithm.It is obvious that between the original function and its approximation there is no visible difference.However, the local error between the exact and approximated function is shown in the right panel of Figure 7.We observe in this Figure that the error first is a small quantity and that for larger values of s the local error decays more than 4 decades.If the approximation of the transfer function is known as a symbolic rational expression it becomes straight forward using Heaviside's expansion theorem for inverse Laplace transforms to find the symbolic solution for the fractional relaxation equation.Heaviside's theorem tells us that the fractional approximation can be represented in terms of exponential functions.The application of the inverse Laplace transform to the approximation of the transfer function G(s) which in fact represents the solution u(t) as an approximation in terms of exponential functions.This solution represents a continuous function which is shown in Figure 8.The results shown in Figure 8 make it obvious that an approximation of the solution of a fractional initial value problem is a direct method and delivers an accurate result with a minimal amount of calculation steps.In fact the steps needed are reduced to the collocation of the transfer function which is approximated by Thiele's algorithm as a continued fraction.This continued fraction is converted to a rational expression which is used in a Laplace inversion to get the solution.Such kind of inversion is always possible if the conditions for Heaviside's expansion theorem are satisfied.Thus any rational fraction of polynomials can be represented by a finite sum of exponential functions if the inverse Laplace transform is applied to this rational fraction in Laplace space.92) generated by the inverse Laplace transform of G(s).The left panel also includes the exact solution of the fractional relaxation equation given by u(t) = u 0 E q (−(t/τ ) q ) where E q (t) is the Mittag-Leffler function.The right panel shows the local error between the exact solution and the approximation.The number of approximation points used in the calculation was m = 2n + 1 = 35.The parameters to generate the plot are α = 1, τ = 1/2, q = 2/3.
In another experiment we examine the convergence of the solution u(t) in case if the number of approximation points n are increased.Since there is only an estimation formula for Thiele's algorithm which tells us that the error of approximation decays exponentially [31], we have to examine the error propagation after we have applied the inverse Laplace transform to the rational fraction G(s).However, since the inverse Laplace transform is carried out analytically we can assume that the error will not increase in this analytic step.In fact the numerical experiments show that the error decays exponentially after the application of the inverse Laplace transform.The decline of the global error based on the maximum norm follows a relation given by ϵ(n) = a exp (−bn c ) where a, b, and c are independent of n.The decay of ϵ(n) is shown in Figure 9 for a specific example.The dots in this figure represent the maximum norm of the error while the solid line represents a least square fit to these data using the predicted error formula.A model using exponential approximations was discussed in a quite different, biological, field examining the binding of oxygen in blood cells.Glöckle and Nonnenmacher used a transfer function in connection with protein dynamics in blood cells [44].The claim of the paper was that protein dynamics described in a standard approach as a superposition of many relaxation processes by a series of exponential functions as can be unified to a fractional relaxation equation of type Equation (92) where the weights w n and relaxation times τ n satisfy a certain scaling relation [44].In view of the method applied here such kind of scaling relation is naturally inherited by the Haeviside expansion theorem.According to Heaviside's theorem the inverse relaxation times are the roots of the denominator polynomial of G(s) and the weights are calculated from To demonstrate that such scaling behavior naturally is included in Heaviside's theorem, we analyzed an approximation with 40 exponential terms.Our observation is that there exist two scaling regimes with the same scaling exponent but with opposite sign.The two regimes are glued together at τ ≈ 1.The slopes of the two scaling regimes are ω = ±0.67512for the parameters selected above.We note that q = 2/3 for this calculation.The scaling behavior is shown in Figure 10 where the weights w n as a function of the inverse of the relaxation times τ n are shown in a double logarithmic plot.This log − log plot directly shows the scaling regions indicated by straight lines.The slopes determined by a least square fit for each of these regions was ω as given above.The scaling is present for at least 6 decades for both regions.In fact the scaling relation in [44] was based on experimental data for rebinding of CO to Mb after photo dissociation [45].These experimental data were represented by a sum of exponential functions as stated in Equation (96).However, it is quite interesting that Nonnenmacher and Nonnenmacher [46] derived scaling relations which are close to our estimations but based on a quite different approach using Markov chains.
To exemplify the scaling behavior for our calculations we also performed calculations with different fractional order q.It turned out that the slope of the scaling exponent determined in a least square fit is approximately representing the fractional order q.These relations are shown in Figure 11.This figure summarizes the calculations for different fractional orders q in a single plot.It becomes obvious that the transition from positive to negative scaling slopes occurs around τ ≈ 1.While for smaller q's there is a smooth transition for larger q values this change of slope occurs more rapidly.
Figure 11.Scaling regions of the weights w n as a function of 1 /τ n 's for a solution with 40 exponentials and parameter values u 0 = 1, τ = 1/2, q = {1/3, 1/2, 2/3, 3/4, 8/9}, and n = 51.The upper lines correspond to q = 1/3 while the lower line is assigned to q = 8/9.The conclusion of this example is that a rational approximation based on Thiele's algorithm is very effective if we have to solve initial value problems approximated or represented by rational transfer functions.For such cases the inverse Laplace transform is delivering the solution if the prerequisites of the Heaviside expansion theorem are satisfied.Knowing that the final solution of a fractional relaxation equation can be represented as a sum of exponential functions it becomes consistent with experimental findings that a rational approximation of a fractal transfer function is related to the fractal behavior of the physical system.In other words whenever a physical system allows a sum of exponential functions there is a chance that the system shows fractal behavior.

Inverse Laplace Transform
Knowing the facts about Sinc methods of convolution representations of integral equations and the property of linear systems allowing a fractional algebraic representation by a transfer function, we can take another route to find an approximation of the solution by using the inverse Laplace transform.The inverse Laplace transform is a spin off of the definition of the Laplace transform in Equation (24).The inverse transform of Laplace can be generated the following way.
Let the "Laplace transform" F be defined as in Equation (24).If J denotes the indefinite integral operator defined on (0, a) ⊂ (0, c) ⊂ (0, ∞), then [29,31] where " 1 1 1" denotes the function that is identically 1 on (0, a).Hence, with J m ≈ w w w we can proceed as follows to approximate f : • Assume that the spectrum; i.e., the matrix X + and vector s = (s −M , . . ., s N ) T have already been stored for some interval (α, β), corresponding to matrix A + m , make the replacement s → α/(β − α)s, and compute the column vector x = g * v with * a Hadamard product, (99) • All operations of this evaluation take a trivial amount of time, except for the last matrix vector equation.However, the size of these matrices is nearly always much smaller than the size of the DFT matrices for solving similar problems via FFT.
Then we have the approximation f (x j ) ≈ f j , which we can combine, if necessary, with our interpolation Formula (11) to get a continuous approximation of f on (0, a).
An example we already discussed above is related to the relaxation equation.It's transfer function is defined in general as G(s_, u0_, τ _, q_):= u0s q−1 s q + τ −q Specifying the related parameters parameters } ; we get the specific representation of the transfer function for any q and relaxation time τ as follows g1 = G(s, u0, τ, q)/.parameters 1 (s 2/7 + 2 2/7 ) s 5/7   The inverse Laplace transform of this transfer function delivers the time dependent relaxation function as shown in Figure 12. sol1 = stengerInverseLaplaceTransform[g1, {s}, {t, 0, 100}, 1, 0, 32]; Now it becomes straight forward to determine the time domain behavior of a system if the transfer function G(s) is known.A practical example of such an application is the modeling of ultracapacitors [49].Ultracapacitors or electrochemical double layer capacitors are handled as energy storage devices which may replace batteries due to their high capacity and long standing high energy density [50,51].As a result, ultracapacitors are becoming more widely used in energy storage systems especially in mobility applications [50][51][52].Such an ultracapacitor can be modeled as a ladder of resistors and capacitors building up an RC transmission line; for details see [49].The impedance transfer function of such a model has the structure textZs = 664 Figure 12.Solution of the fractional relaxation Equation (87) generated by the direct inverse Laplace transform of G(s).The left panel also includes the exact solution of the fractional relaxation equation given by u(t) = u 0 E q (−(t/τ ) q ) where E q (t) is the Mittag-Leffler function.The right panel shows the local error between the exact solution and the approximation.The number of approximation points used in the calculation was m = 2n + 1 = 65.The parameters to generate the plot are u 0 = 1, τ = 1/2, q = 2/7.
The numeric terms are due to the scaling given in time steps of months.The time dependent impedance of such a model now directly follows via the inverse Laplace transform for a large range approximation by It is obvious from the long term and short term behavior of the time dependent impedance that two plateaus exist (see Figure 13).One for the short time behavior which lasts about 10 −4 months which corresponds approximately to 45 min.The second plateau stabilizes on a level of approximately 220 after 10 3 month corresponding to approx.83 years.The drop from a level of about 640 to 300 lasts about 8 month.With such a characteristic the device is a highly promising energy storage device to replace the electrochemical based batteries.It becomes also straight forward to examine the impedance model for different fractional orders q if we generalize the model in [49] to the following relation zS4 = 448 2.06s q + s + 215 s ; The inverse Laplace transform for a sequence of fractional exponents is generated in the following line based on the sequence q ∈ { i 10 } 9 i=1 in a parallel calculation for each of the exponents q sol4 = ParallelTable[ { q, stengerInverseLaplaceTransform [ zS4, {s}, { t, 0, 10 6 } , 1, 0, 128 ]} , {q, 1/10, 9/10, 1/10}]; The results of these calculations are shown in Figure 14.It becomes apparent that the drop in the impedance is more rapidly for smaller fractal exponents than for lager ones 0 < q < 1 which may be a design indication.Figure 14.Time dependent impedance Z(t) derived from the inverse Laplace transform of the ultracapacitor with different fractal exponents q ∈ {i/10} 9 i=1 .The steepest transition corresponds to q = 1/10.This example demonstrates that Sinc methods are effective especially for practical engineering applications.The inverse Laplace transform of a linear fractional differential/integral equation is accessible with a minimum amount of calculation steps and a high accuracy.

Levinson's Integral Equation
Finally we discuss a model which originally was not formulated as a fractional model but in retrospect can be interpreted as such.Levinson in 1960 [40] based on a paper by Lin on the hydrodynamics of liquid Helium [41]  The solution u(t), u(t) + (6/5) sin(t), and (6/5) sin(t) are shown on the interval t ∈ [0, 2π] in the following Figure 15.In this Figure we compare the solution u(t) with the exogenous driving function and the addition of u and the driver function.It is apparent that the solution of the non-linear integral equation is not related any more to the exogenous function as the linear theory delivers.The structure of the solution is distorted but shows an oscillating structure nearly with the same period as the exogenous function.The results of such variation of µ is shown in Figure 17.

Conclusions
The different models examined here demonstrate that Sinc methods are well equipped to deal with fractional differential or integral equations.The Sinc approach to approximate convolution integrals is one of the keys to deal with linear and non-linear fractional equations delivering with a minimal amount of calculations a highly accurate result.If the equations are linear the approach of the transfer function becomes effective in connection with Thiele's algorithm and with a direct inversion of the Laplace transform.Both methods based on transfer functions are using Sinc points as discrete points which are a nearly optimal choice of the discrete mesh if the Lebesgue measure is used.The discussed models also demonstrate that the misnomer "fractal" equation is only on spot if the solution shows an asymptotic scaling law.From the variety of solutions generated with these few examples it seems to be more precise to call the "fractal" equations fractional equations which reflects more the diversity of solutions of such equations.Since all the fractional equations are singular integral or differential equations the Sinc method proves once more that it is self-sufficient of the singular behavior.

Definition 1 .
Domain and Conditions.Let D be a simply connected domain in the complex plane and z ∈ C having a boundary ∂D.Let a and b denote two distinct points of ∂D and let φ denote a conformal map of D onto D d , where D d = {z ∈ C : |Im(z)| < d}, such that φ(a) = −∞ and φ(b) = ∞.Let ψ = φ −1 denote the inverse conformal map, and let Γ be an arc defined by Γ = {z ∈ C : z = ψ(x), x ∈ R}.Given φ, ψ, and a positive number h, let us set z k = ψ(kh), k ∈ Z to be the Sinc points, and let us also define ρ(z) = e φ(z) .

Definition 2 .
Function Space.Let d ∈ (0, π), and let the domains D and D d be given as in Definition 1.If d ′ is a number such that d ′ > d, and if the function φ provides a conformal map of D ′ onto D d ′ , then D ⊂ D ′ .Let L L L α,β (D) be the set of all analytic functions, for which there exists a constant c 1 , such that

Theorem 2 .
Non-linear Convolution.Let D = D(a, b) be defined as in Definition 1, and let L α,β (D) be defined as in Definition 2. Assume that K(• , c) is analytic and uniformly bounded in D(a, b) for each c ∈ [c 1 , c 2 ] and also, that K(x , •) is analytic and uniformly bounded in D(c 1 , c 2 ) for each fixed x ∈ (a, b) .If u ∈ L α,β (D), and if P is defined by

Figure 1 .
Figure 1.Riemann-Liouville integral of order ν = −1/3 of the function f (x) = e ax where a = 1 (top panel).Shown are also the locations of the Sinc points along the solution.The bottom panel shows the local error for N = 32 Sinc points.

Figure 3 .
Figure 3. Approximations of the fractional relaxation Equation (84) using the fractional order β = 1/2 in a Caputo fractional representation.The initial condition and the relaxation time τ 0 was set to ϕ 0 = 1 and τ 0 = 1/10.The solid line represents the analytic solution E 1 (−10t).The Sinc approximation of the Caputo representation used N = 64 Sinc points.

Figure 4 .
Figure 4. Approximations of the fractional relaxation Equation (84) using different fractional orders β.The initial condition and the relaxation time τ 0 was set to ϕ 0 = 1 and τ 0 = 1/10.The solid line represents the analytic solution E 1 (−10t).The Sinc approximation of the Caputo representation used N = 64 Sinc points.

Figure 7 .
Figure 7. Approximation of the transfer function G(s) for a fractional relaxation equation given by Equation (94).The left panel shows the exact transfer function G(s) with u 0 = 1, τ = 1/2, and q = 2/3 in connection with the approximation using Thiele's algorithm.The right panel shows the local error between the exact function and its approximation.The total number of point used in the approximation is m = 2n + 1 = 35.

Figure 8 .
Figure 8. Solution of the fractional relaxation Equation (92) generated by the inverse Laplace transform of G(s).The left panel also includes the exact solution of the fractional relaxation equation given by u(t) = u 0 E q (−(t/τ ) q ) where E q (t) is the Mittag-Leffler function.The right panel shows the local error between the exact solution and the approximation.The number of approximation points used in the calculation was m = 2n + 1 = 35.The parameters to generate the plot are α = 1, τ = 1/2, q = 2/3.

Figure 13 .
Figure 13.Time dependent impedance Z(t) derived from the inverse Laplace transform of the ultracapacitor.
derived in a paper about super-fluidity the following non-linear Volterra equation of second kind Applying to Equation (101) the notation of fractional calculus we realize that this equation is equivalent to a non-linear fractional differential equation of the following form = 0 and where E α,β (t) denotes the generalized Mittag-Leffler function.Levinson proved in his paper that a solution of the nonlinerar integral equation exists and may be derived by successive iterations as an approximation.We take the position here that the existence of the solution is guaranteed by the proof of Levinson and approximate such kind of solution by applying non-iterative Sinc methods to Equation (102).The steps of the calculation are implemented in a function RiemannLiouvilleNonlinearFIE[] which sets us into position to directly solve the equation if we specify the non-linearity of the integral and the exogenous function g(t) = −c sin(t).The following line performs the calculation with a fractional order α = 1/2, with λ = −1 / √ π , and N = 64 Sinc points 0 (t − τ ) −1/2 u(τ ) 3 dτ = −c sin(t).(101)