An introduction to fractional calculus: Numerical methods and application to HF dielectric response

The aim of this work is to introduce the main concepts of Fractional Calculus, followed by one of its application to classical electrodynamics, illustrating how non-locality can be interpreted naturally in a fractional scenario. In particular, a result relating fractional dynamics to high frequency dielectric response is used as motivation. In addition to the theoretical discussion, a comprehensive review of two numerical procedures for fractional integration is carried out, allowing one immediately to build numerical models applied to high frequency electromagnetics and correlated fields.


Introduction
Fractional calculus has a long history, at least as old as the usual calculus [1]. Concepts on derivatives and integrals of general order have remained in the realm of pure mathematics for a long time, but have been applied to a variety of different fields in the last decades, ranging from control theory to electrodynamics, economics and quantum mechanics [2]. Among several applications and interpretations, fractional calculus (FC) seems to constitute a well-suited apparatus to model systems with memory, i.e., systems whose response to a local stimulus depends on their whole "history" [2,3]. Such property makes FC an important tool in linear classical electrodynamics, since field analysis inside bulk matter takes into account non-local contributions of the stimulus, both in time and space [4]. In fact, there is a formal similarity between fractional integral development and Maxwell equations for non-local media, since in both scenarios convolution integrals play the main role [3,4,5], connecting stimulus and response fields. In addition, there are evidences that for sufficiently high frequencies, dielectrics of completely distinct kinds present a power-law response. This so-called dielectric universal response was discovered first by Jonscher [6,7] rises naturally in theory when one uses fractional integrals in Maxwell non-local equations, as shown by Tarasov in [8,9,10].
It is natural that, in face of the potential applicability of FC to classical electrodynamics, numerical methods are devised to approximate solutions of complex problems. More essentially, FC machinery is not trivial 1 , so that numerical methods might be useful even for simple estimates. In 1 Very often, integral or derivative of simple functions such as exp(t) this context, the present work intends to present the fundamentals of FC, as well on numerical methods for fractional integral computation. Such concepts are useful for those who seek for alternative formulations of high-frequency response or, in a wide sense, to any non-local response theory [3,9].
The article is organized as follows: • In Section 2 FC fundamentals are exposed, ranging from basic definitions up to identification criteria for fractional operators. The development focuses on two main definitions: Riemann-Liouville integral and Grünwald-Letnikov derivative, due to their interconnection.
• Section 3 carries out a review on non-local fields in classical linear electrodynamics and dielectric response, where Jonscher's dielectric universality is discussed.
• In Section 4 a constructive review on numerical methods for fractional integration is carried out. The analysis focuses on two different tools: linear multi-step and Newton-Cotes methods, discussed in both classical and fractional scenarios. Section 4.7 presents a few illustrations of the preceding methods.
A deep treatise on the mathematical aspects of fractional calculus is found in Samko et al. [11], whereas more applied approaches are given by Hermann [2], Oldham and Spanier [12] and Tarasov [9]. An introductory course on linear electrodynamics can be found in Greiner [13], whereas a much more detailed development on field response can be found in Oughstun [4]. Classical texts such as Dahlquist [14], Ralston and Rabinowitz [15] and Golub and Ortega [16] provide the whole background necessary to the understanding of the considered numerical methods, at least in the classical scenario. The generalization for fractional integrals through linear multi-step methods is strongly based on [17], supported by the formalism introduced in [18,19,20,21]. On the other hand, fractional Newton-Cotes formulae development can be found, for example, in Li and Zeng [22]. lead to results expressed by unusual special functions. Vide, e.g., Hermann [2].

Fractional calculus 2.1. Fundamental concepts
Usual differentiation and integration are well-known concepts in science and technology, and have a strong theoretical basis, built on the efforts of Newton, Leibniz, Weirstrass, among many other. However, as old as the usual differential and integral calculus, there is the socalled fractional calculus apparatus. According to Ross [1], fractional calculus origins goes back to a dialog between L'Hospital and Leibniz, in which the meaning of the operation d 1/2 /dt 1/2 was discussed. FC main idea consists in generalizing the order of a derivative,say n, allowing one to perform derivations -and also integrations -of α-order, with α ∈ C. Immediately, one may ask about the meaning of the operation and how it is operationalized. An intuitive definition of fractional derivative was given by Fourier [2,5]: it is well known that Fourier transform maps the derivative of its operand, say f (t), to a product between Fourier transform of the original function and a power of ω, with  = √ −1, i.e., Generalizing the derivative order to α ∈ C and performing an inversion, one obtains where F(ω) = f (t) exp(−ωt)dt is the Fourier transform of f . For example, if f (t) = sin(ω 0 t), one obtains D α [f ](t) = |ω 0 | α sin ω 0 t + π 2 α , which leads to the usual result when α ∈ Z. This is not a coincidence, but a fundamental requirement for fractional derivatives, as will be discussed later.
Despite its simplicity, Fourier definition of fractional derivative is not unique: there are many others. For instance, consulting [2], [5] and [23], one finds more than thirty definitions, not generally equivalent. Before introducing two of them, it is important to emphasize the fact that Fourier fractional derivative, Eq. (3), is an integral in the classical sense. This is a general characteristic of fractional operators: derivatives and integrals are deeply connected, similarly to what occurs in the usual scenario, where Fundamental Theorem of Calculus relates both operations. To visualize such connection in the fractional scenario, we first remind that an arbitrary number α can be decomposed on its integer and fractional parts, i.e., α = [α] + {α}, with [α] ∈ Z and α ∈ [0, 1[. Therefore, Eq. (4) uses the important semi-group property, which basically ensures the commutativity of D [·]. In order to connect Eq. (4) to a fractional integral [2], one expands D α into D n D α−n , with n = [α] + 1: that is, an α-order fractional derivative is related to an (1 − α)-order iterated fractional integral, denoted by I[·], followed by an integer n-order derivative, with n = [α] + 1.
This intricate connection has led Oldham and Spanier [12] to coin picturesque terms such "differintegral" and "differentegration", since if one defines an operation identified as a fractional integral, its derivative is readily provided, according Eq.(5). In fact, for most of the definitions of fractional operators, such as the Caputo and Riemann-Liouville approaches, one starts with the fractional integral, since the derivative counterpart is already defined. One remarkable exception is the Grünwald-Letnikov derivative, which is based on the generalization of finite differences.
In what follows we shall focus on two definitions: Riemann-Liouville and Grünwald-Letnikov. These two approaches were chosen due to their interesting interconnection, of great importance for numerical purposes.

Riemann-Liouville approach
The following result, due to Cauchy 2 [5,2,11] relates the iterated integral to a simpler form, given by Generalizing the order in Eq. (7) to an arbitrary α, one obtains left and right-sided Riemman-Liouville integrals of order α > 0, respectively: (8) and with t < a, and. (9) In both Eqs. (8) and (9), Γ(·) refers to the Gamma function. According to Eq. (5), to obtain the Riemann-Liouville fractional derivatives, it suffices to perform an integer differentiation of order on the expressions above.

Grünwald-Letnikov approach
It is well-known that the usual n-order derivative of f (t) can be written as the limit in which is the binomial coefficient. Suggestively, one might gener- Now, the binomial coefficient is given by .
The Grünwald-Letnikov fractional derivative is given by the following limit: If h > 0, one deals with backward differences, otherwise one has forward difference. Care must be taken when α < 0, once the series given by Eq. (12) can easily diverge. However, as pointed out by Samko et al. [11], the series shown in Eq. (12) eventually converges if f (t) is a non-periodic function which smoothly -and fastly -decays to zero at infinity. Within these conditions, Eq. (14) turns into Grünwald-Letnikov fractional integral [11,3].

On the plurality of definitions: identification criteria
As pointed out earlier, the number of different definitions for fractional integrals and derivatives is remarkable. Surprisingly, different definitions lead to different results when applied on the same operand. In light of such plurality, it is natural to ask how can one identify a "legitimate" fractional operator. In 1975, Ross [1] proposed five criteria in order to ensure the fractionality of a given operator: 1. The α-order derivative of an analytic function f (z) must be analytic on z and α; 2. If α ∈ Z + , the result must be the same obtained by usual -integer order -differentiation. If α ∈ Z − , the result must be the same of the α-fold iterated integral; 3. The operator must be linear; 4. If α = 0, the operand must remain unaltered. In other words, 0 is the neutral element; 5. Given α, β ∈ C, the semi-group property must hold: In addition to the criteria proposed by Ross, there is a proof given by Tarasov [24] which shows that for a legitimate fractional operator, the usual Leibniz rule, given by 16) must be broken. Ortigueira and Machado [25] used this result to improve Ross criteria, stating that 6. For a true fractional operator, the generalized Leibniz product rule must hold: In this same work, the authors show that Riemann-Liouville and Grünwald-Letnikov fractional derivatives -among several other -obey Eq. (17). The six criteria shown above work as a sort of filter, since there is a profusion of allegedly fractional operators in the literature.

Linear media and temporal non-locality
In an arbitrary material media, Maxwell equations are given by in which ρ and j are the free charge and current densities, E, D, B and H are the electric field, electric displacement, magnetic induction and magnetic field, respectively. The fields E, P and D are related by in which P is the polarization of the medium and 0 is the vacuum permittivity. The simplest interpretation for P states that this field represents the dipole moment per unit (of a non-infinitesimal) volume. For linear media, there is a direct relation between E and P, expressed as in which χ is the electric susceptibility. Combining Eqs. (18) and (19), one obtains the so-called constitutive relation for linear media: where the permittivity is given by ≡ 0 (1 + χ). One should notice that we are assuming isotropic media, since is a scalar function 3 . The simplicity of Eq. (20) is due to two strong implicit assumptions: • a stimulus provided by E at a given spatial coordinate r affects D only at r ; and • a stimulus provided by E at a given time instant t is immediately felt by D.
In fact, in the general case, linearity between D and E is expressed by [4], [13] where (·) is the integral kernel. If the medium is local [4], assumption 1 above is valid, since and then Spatial locality roughly means that the medium molecular constituents are independent from each other [4]. Such assumption is, in some measure, admissible. On the other hand, assumption 2 is not realistic: it assumes that information propagates at infinity speed in the medium. In other words, it violates the principle of causality. However, in static and quasi-static scenarios, Eq. (20) is a good approximation [13]. If we make a further assumption, assuming that the medium is temporally homogeneous 4 , we have which leads to a temporal convolution between D and E: (25) Eq. (25) makes explicit the fundamental characteristic of temporal non-locality, since the response D depends, in a certain manner dictated by , on the whole history of the stimulus E. A last improvement of Eq. (25) requires us to take into account the principle of causality: since response can not precede stimulus, we must have (r; t − t ) = 0 when t − t < 0. Therefore, one deals now with a causal temporal convolution: A medium in which Eq. (26) is valid is denominated as linear temporally dispersive [4]. Applying Fourier transform on both sides of Eq. (26) leads us to in which k and ω are the wave vector and angular frequency, respectively. The convolution in Eq. (27) is due to the fact that in order to ensure causality, we have performed a product in time-space between (r, t)E(r, t) and the Heaviside function, θ(·), whose Fourier transform pair is given by Under adequate assumptions of symmetry and analyticity of stimulus and response functions, one obtains from Eq. (27) the important Kramers-Kronig (KK) relations [4,13], which connect real and imaginary parts of the system's optical response. KK relations are built on the theory of Hilbert transformation, and the interested reader may consult the works of Oughstun [4], Greiner [13] and Witthaker et al. [26] for further details.

Empirical models for polarization
Once the role played by temporal dispersion was made explicit, one might ask about the behavior of the field P when subjected to a stimulus from E . There are many empirical models which relate these fields [4,27]. In general lines, one tries to find a phenomenological model for an individual dipole moment p and then computes its average structure, related to the macroscopic polarization P: in which the brackets indicate ensemble average and f i are dipolar statistical weights. For example, a working model for non-dense media [13] starts from the hypothesis that each electron of the medium is subject to classical damping and elastic forces. Therefore, one has the movement equationẍ in which γ i , ω i and m are the i-th damping constant, characteristic frequency and electronic mass, respectively. The i-th state is characterized by γ i and ω i , i.e., f i equals the number of electrons with damping constant and angular frequency γ i and ω i . If the field is harmonic, i.e., E(t) = E 0 exp(−ωt), Eq. (30) has a solution given by Using Eq. (31) in Eq. (29) and reminding that p = ex, one obtains where N is the number of molecules per unit volume. With Eq. (32) and (27), one finds expression for the optical responses (ω) or χ(ω) .

Dielectric universality
As pointed out at the beginning of the section, there are several models for polarization, such as the Clausius-Massoti and Cole-Cole models. Interested reader should refer to [4] and [27] for deeper analysis on the subject. However, it is interesting to mention the important model due to Debye [4,7]. In this model, the electrical dipole p is supposed to suffer viscous resistance, being allowed to rotate in the medium, but it is not able to interact to other dipoles and there is not Brownian motion in such "fluid" medium. This model leads to the following movement equation [4]: in which a is a coupling constant and τ is the molecular relaxation time, i.e., the characteristic time interval the dipole takes to return to a random orientation after E ceases. It can be shown [4,6] that Debye model leads to a susceptibility of the form It is important to emphasize that Debye's model can be improved [4,27], but still fails at some point. Experimental observations in this direction can be found in [28] and [29], for instance. In fact, one might build more complete models for dielectric response, but at the expense of dealing with more variables and a increasing degree of complexity. Instead, we look now at the behavior of a large class of media at high frequencies. Interestingly, very distinct materials show some power law behavior for the dielectric response. Let us denote the susceptibility by According Jonscher [6,7], the following proportionality is observed at high frequencies: Eq. (36) shows that the ratio between imaginary and real parts of χ does not depend on the frequency: In [7,8,9] some theoretical scenarios are discussed in order to explain this universal behavior. Some of them are built on generalizations of the Debye model; for instance, models where there are many relaxation times, instead of only one. In such models, one should consider not only one relaxation time, but an ensemble average on the distribution. In addition, the concept of fractionality is also considered, i.e., according to such approach, the power-law behavior could be associated to some sort of scale invariance of the medium.

Universal dielectric response as a fractional process
Among a plurality of interpretative scenarios for universal dielectric response, we turn our focus to the formal aspects subjacent to the phenomenon. Following Tarasov in [8,9,10], we start from constitutive relations for linear media and the high-frequency universal behavior, Eq. (36), moving towards a fractional integral relation. According to the discussion carried out in Section 3.1, one has to expect the following relation between P and E: Eq. (38) assumes that the susceptibility χ does not depend on r and is itself a causal function, so that the integral limits can be set at infinity. Taking temporal Fourier transform on both sides, one has P(r, ω) = 0 χ(ω)E(r, ω).
Using Eq. (36) as an equality with the proportionality constant equals one in the previous result, and denoting −α = n − 1, we obtain P(r, ω) = 0 (ω) −α E(r, ω), with α ∈]0, 1[ and ω 1 τ . (40) Samko et al. [11] show that the Fourier transform derivative property, given by Eq. (2) remains valid in the fractional scenario. The right-hand side of Eq (40) contains a term (ω) −α , leading to an α-order integral in the time domain [10]: The result shown in Eq. (41) allows an interesting interpretation: the polarization of the medium depends on E past history in a more complex way than in the usual scenario. Temporal dispersion is related to memory effects, vide Eq. (25). But the fractional relation between stimulus and effect opens door to further interpretations; for instance, to those related to fractality and multiple time scales for the relaxation times. One last comment: the previous development is due to Tarasov [9,8,10], but he was not the first to notice fractional characteristics in the universal dielectric response. In fact, in [27], Garrappa et al. investigate several models whose original formulations made use of fractional integral kernels.

Numerical approximation of fractional integrals 4.1. Fundamentals
Once FC main ideas were exposed and one of its applications to electromagnetism was detailed, it is interesting to consider how one can discretize fractional integrals and derivatives and approximate them numerically. Obviously, such study would comprise an entire work by itself, so that we will focus on the elements which provide the minimum apparatus on how to proceed with applied problems in FC using Riemann-Liouville integrals. The interested reader can refer to [22] and [17] for much deeper analysis. In order to approximate the Riemann-Liouville (hereafter identified by RL) integral, one might explore its "convolutional" structure, shown in Eqs. (8) and (9), since one has -for the left-sided case - where * is the convolution symbol. An analog expression exists for the right-hand sided RL integral. It may be tempting to approximate the convolution shown in Eq. (42) in a straight manner, i.e, through discretization of the independent variable and then performing a discrete convolution in time domain or an usual product in the reciprocal space, through use of the discrete Fourier transform -DFT. However, one should notice that the convolution kernel is singular at the origin when (α) < 1. Therefore, care must be taken in order to perform such approximation.

Linear multi-step method
In this work we want to approximate left-sided RL integrals, Eq. (8) with a = 0, through the use of fractional linear multistep method (FLMM), first proposed by Lubich [17]. This method is a generalization of the the well-known class of linear multistep methods (LMM) used to approximate solutions of ordinary differential equations (ODEs) [16]. Below, the main concepts of usual LMM are reviewed, building the path to the exposition of Lubich generalization.
Consider first the classical case in which α = 1. One may employ LMMs to approximate the following initial value problem: Let t i = i∆t be the i-th node of the mesh generated by discretization of the interval of interest and y i the approximate solution of Eq. (43) at t i . For simplicity, we consider that the increment ∆t is constant along the process. The simplest approach to approximate the considered ODE is the explicit Euler method [16,15], given by It can be shown [16,15] that explicit Euler scheme is a first order method, denoted by O(∆t). However, the main aspect to be considered here is the fact that Eq. (44) is a one-step method, i.e., the value of y i+1 depends only on the value of y i . Another example of one-step method is the class of Runge-Kutta schemes. Alternatively, one might make use of other entries of y to approximate the desired solution. In general, one would deal with the following difference equation where f k = y k and {α k } N 0 and {β k } N 0 are two sets of constants that depend on the structure of the employed method. Eq. (45) is called a linear multi-step method. There is a plethora of methods based on it, including the Adams-Bashfort (explicit) and Adams-Moulton (implicit) methods. References [16,15] are recommended to the reader who seeks for detailed analysis on the subject.
It is clear that the coefficients α k , β k play a fundamental role on both stability and convergence of the LMM shown in Eq. (45). In fact, since one deals with a difference equation with constant coefficients, discrete Laplace transform rises as a powerful tool for analysis of such systems. It is defined by 5 and allows us to obtain the generating polynomials of the considered LMM [15,31,32,18,19] Since the polynomials ρ(z) and σ(z) completely determine the LMM, it is customary to designate the method by (ρ, σ).

Classical convolution quadratures
As pointed out at the beginning of Section 4.1, one might approximate the left-sided RL integral through the use of discrete convolutions. Despite the singularity at kernel's origin, Eq. (42), convolution is, in fact, a tool to be considered. We want to approximate an integral (possibly singular at some finite set) I(t) = t 0 f (t )g(t − t )dt by means of a discrete convolution: where ω j depends on ∆t and s ≤ n . The second sum is necessary to ensure convergence near the origin, and the coefficients µ nj are the so-called starting quadrature terms [32,20,19]. They can be obtained forcing the approximation given by Eq. (66) to be exact for polynomials of degree up to s [20]. It can be shown [32,19] that ω j is the j-th coefficient of the rational series Observing Eq. (49), one notices the need for a Laplace transform inversion in order to estimate ω j . However, as pointed by Lubich in [20] and [21], ω j is proportional to the s-order approximation of the inverse Laplace transform L −1 (s) = f (t), when t is far from any singularity. In the general case, Eq. (49) requires a priori knowledge of the analytic Laplace transform of the considered LMM. If this is not the case, numerical procedures can be used to perform an inversion.

Fractional quadratures and the Grünwald-Letniknov fractional derivative
As considered in the previous section, linear multi-step methods can be used to approximate convolution integrals, in the classical sense. Interestingly, the generalization for the fractional case is quite straightforward: Lubich shows in [17] that the fractional left-sided Riemann-Liouville integral can be approximated through the coefficients of in which ω(z) is given by Eq. (49). With this result, one has the approximation One fundamental requirement about fractional linear multistep methods (FLMM) is that they must be implicit [17]. A beautiful result rises when one uses the implicit Euler scheme, given by By direct inspection of Eq. (52) and application of Lubich fractional generalization, Eq. (51), one obtains the following polynomial: Using this result in the discrete convolution given by Eq.
(51), and discarding for the moment the starting weights, one has which corresponds precisely to a truncated approximation for the Grünwald-Letnikov derivative of α-order, vide Eq. (14). From a theoretical perspective, Eq. (54) builds the path to the proof on the equivalence between Grünwald-Letnikov (GL) and Riemann-Liouville fractional derivatives; see, for example, [5] and [33]. This connection allows us to approach theoretically a certain fractional problem using RL integrals, but numerically approximating it through GL derivatives of real negative orders [22].

Gamma and factorial functions overflow and the short-memory principle for Grünwald-Letnikov approximation
The approximation obtained by means of Eq. (54) has an important drawback: the sequence decays rapidly to zero as α or j increase, as can be seen in Figure 1. Simultaneously, the coefficients −α j diverge. In this sense, one loses any control over precision and stability of GL approximation. Another error source due the form of Eq. (55) is related to the behavior of gamma and factorial functions for large numbers in computational environment. For instance, SCILAB language [34] returns infinity for j! = Γ(j + 1); j > 170. In practical terms, this overflow sets a limit on the mesh refinement. Such limitation can be visualized in Figure 2, where it is approximated the solution of I 0.5 + [f ](t) for f (t) = 1, t ∈]0, 10], whose exact solution is given by It is shown in Figure 3 how the relative error of the GL approximation increases artificially above t ≈ 3.4.   An immediate solution for this issue consists in making the grid coarser; vide Figures 4 and 5. However, it should be noticed that Euler's scheme, even in the fractional case [17], is an order one method. Therefore, making ∆t large leads to a progressive loss of the accuracy.
The observed errors are due solely to number representation limitations on computers, and have nothing to do with the nature of the fractional operations discussed so far. On the other hand, the rapid decay of GL approximation coefficients introduces some interesting aspects. It shows that, at least for the case of fractional derivative (α > 0), the lower terminal a (taken equal to zero here) contributions are not always fundamental for the evaluation of D α [·]. In other words, there is some sort of "numerically induced locality" in a non-local operator. That means that one could use only few past entries for each n in order to approximate D α n and, therefore, I −α n , according Section 2.1. This is what Li and Zeng call short-memory principle [22]. In their work, authors digress on how one can estimate a truncation term. We, however, follow a different path, in which zero and one order approximations of the generalized Newton-Cotes are used [22,35]. Such approach provides similar results than those obtained by FLMM methods (up to order one) and does not suffer from the former approach practical complications, such as the need for Laplace transform inversions and truncation terms estimation.

Interpolation methods and the Fractional Newton-Cotes formula
Instead of using FLMM to approximate I α n , one can appeal to interpolation methods, whose formalism is built on orthogonal polynomials local expansions [14]. The main idea consists in approximating the integrand by an easily integrable function within a given partition of the original interval. As shown by Li and Zheng [22], this classical method can be extended without further difficulties to the fractional scenario. For example, let us consider a zero order approximation for f , i.e., when one wants to approximate f (t) in a given partition t ∈ [t k , t k+1 ] by f k+1/2 ≡ 0.5(f k + f k+1 ).
In the case of f (t k ) ≈ f k , the left-sided RL integral reads in which The case for f (t k ) ≈ f k+1 is completely analogous, leading to where c k is given by Eq. (57). Taking the average between Eqs (56) and (58), one deals with the fractional composite trapezoid rule [22]: In principle, one might choose higher order approximations for f and apply the same procedure above in order to obtain better estimations. This is possible since the integrand can be locally approximated by a polynomial. In the general case, one can represent f (t) by its Lagrange polynomial [15,14], given by with In the expressions above, P is the number of sub-partitions of the interval [t k , t k+1 ].
In a formal sense, there is no impediment to the derivation of higher orders schemes, provided by the Fractional Newton-Cotes formula [22]: It may be tempting to increase in order to achieve better results; however, one should notice that for equidistant points, Lagrange interpolation presents oscillatory behavior near terminals [14], leading to bad representations of the integrand. In fact, for sufficiently smooth integrands, low order methods provide satisfactory results.

Numerical examples
To illustrate the use of fractional linear multi-step and interpolation methods to approximate left-sided RL integrals, let us consider two integrands: f (t) = 1 and g(t) = e t , whose exact integrals are given respectively by where γ(·, ·) is the lower incomplete gamma function, given by whose numerical evaluation can be performed through Taylor series expansion, at least for small t [36]: It is interesting to notice that both Grünwald-Letnikov approximation and zero order Newton-Cotes method (NC-0) can be represented by a discrete causal convolution: in which the coefficients are given respectively by Therefore, both methods can be efficiently implemented by means of use of FFT algorithms, largely available [34]. Fig. 6 shows results for I[f ] α (t) using GL approximation, where it can be noticed a good agreement among exact and approximate results. In Fig. 7 are shown the errors of the GL approximation, corroborating to the good accuracy statement, since far from origin (a point of singularity), relative deviations fall below 0.1% for all tested α, in the worst case scenario, at t 6.0. NC-0 method provides similar results to those shown in Fig. 6. In fact, both methods differ little, as can be seen in Fig. 8, where absolute difference between approximations are shown. In either case, it is important to keep in mind that the integrand is a zero-order polynomial, which by construction of the methods, must lead to a high fidelity approximation.     9 shows results for I α [g](t) using GL approximation, presenting again a good agreement between exact and approximate solutions. Errors are shown in Fig. 10. NC-0 method also led to good approximations, but since neither of methods are built to match higher order polynomials, it was verified an increasing of the error magnitude for larger  values of t, vide Fig. 11. This is also related to the limited extent of the series representation for the lower incomplete gamma function, vide Eqs. (64) and (65). Figure 11: Absolute deviation among GL and NC-0 approximations for g(t). One notices an error amplification for values near superior terminal.

Conclusions
The present paper introduced Fractional Calculus main concepts, performing a connection between this theoretical background to applied classical electrodynamics. Such connection is feasible due mainly to the "convolutional" structure of both fractional operators -such Riemman-Liouville integrals and Grünwald-Letnikov derivativesand constitutive relations of linear non-local media.
In order to provide means for numerical experimentation on this theoretical development, a review on two of the most simple numerical tools was carried out, in which were pointed out practical aspects, such as gamma function overflow and implementation through FFT algorithms. With respect to the former occurrence, it is clear that order zero Newton-Cotes approximation should be chosen when fine grids are necessary.
Despite the fact that only low order schemes were used in this work, general formulations were exposed, allowing one to build more complex models, at the expense of either performing a Laplace transform inversion (FLMM), or using higher order polynomial approximations of the integrand (fractional Newton-Cotes method). Both choices have their virtues and drawbacks and it is on the operator to decide which one is better for a given problem.