# An Implementation Solution for Fractional Partial Differential Equations

Authors: Nicolas Bertrand; Jocelyn Sabatier; Olivier Briat; Jean-Michel Vinassa
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

The link between fractional differentiation and diffusion equation is used in this paper to propose a solution for the implementation of fractional diffusion equations. These equations permit us to take into account species anomalous diffusion at electrochemical interfaces, thus permitting an accurate modeling of batteries, ultracapacitors, and fuel cells. However, fractional diffusion equations are not addressed in most commercial software dedicated to partial differential equations simulation. The proposed solution is evaluated in an example.


---

## Body

## 1. Introduction

For an electrochemical system, species diffusion is often modeled by Fick’s laws [1]. However, in almost all systems, the transport mechanism is different from the classical diffusion characterized by Fick’s laws. This kind of diffusion is denoted by anomalous diffusion [2]. Anomalous diffusion is characterized by a mean squared displacement of the diffusing particles that has a power law dependence on time 〈r2〉~tγ with γ between 0 and 2 (for classical diffusion 〈r2〉~t).The theoretical approach of this type of diffusion is strongly related to fractional calculus [3]. It was indeed demonstrated that stochastic processes of random walks can be represented by fractional diffusion equations [4].Among all existing anomalous diffusion equations, three diffusion modes, respectively, called “anomalous diffusion Ia” (ADIA), “anomalous diffusion Ib” (ADIB), and “anomalous diffusion II” (ADII) are characterized by a Fick’s equations adaptation for phenomena considered [5].In this paper, the authors have only considered ADIB type diffusion equations but the proposed contribution can be extended to others classes of equation. In the sequel, the link between fractional differentiation and diffusion is used to propose a solution for the implementation of a fractional diffusion equation in software such as COMSOL Multiphysics. These software applications are now powerful tools for engineers to simulate complex systems combining several physical domains such as electrochemistry and thermal. However they are not adapted to take into account anomalous diffusion and thus to model diffusion interfaces as in batteries, ultracapacitors, or fuel cells.The link between fractional differentiation and diffusion equation is reminded in the second section of the paper. This link should be used to implement fractional differentiation in software dedicated to numerical solving of partial differential equation such as COMSOL Multiphysics software. However, as shown in Section2, the diffusion equation form of a fractional system requires the computation of an inverse Fourier transform that is in most cases impossible to get analytically. This is why this paper proposes alternative partial differential equations approximation that exhibits a fractional behavior in a given frequency band. These differential equations can be easily implemented to simulate a fractional differentiator and thus a fractional diffusion equation.


## 2. Link between Fractional Systems and Partial Differential Equations

### 2.1. Partial Differential Equation Representation and Approximation of a Fractional System

For presentation simplicity, the following fractional system (fractional integrator) is considered(1)H(s)=s-γ,
with 0<γ<1. Its link with diffusion equation can be demonstrated using the system impulse response [6] defined by the Mellin-Fourier integral of (1):
(2)h(t)=L-1{s-γ}=limω→∞12jπ∫c-j∞c+j∞ests-γds,
where c is greater than the abscissa of the singular points of H(s). Using poles definition that can be found in [6, 7], this system does not generate poles and its impulse response is thus given by
(3)h(t)=sin(γπ)π∫0∞x-γe-txdx.Response of system (1) to an input u(t) is defined as the convolution product of the impulse response h(t) with the input u(t):
(4)y(t)=∫0th(t-τ)u(τ)dτ,
and thus using relation (3) and through an integral permutation
(5)y(t)=∫0∞sin(γπ)πx-γ(∫0te-(t-τ)xu(τ)dτ)dx.Let(6)w(t,x)=∫0te-(t-τ)xu(τ)dτ
the following state space representation can be obtained for system (1):
(7)∂w(t,x)∂t=-xw(t,x)+u(t),y(t)=sin(γπ)π∫0∞x-γw(t,x)dx.x∈ℝ+.Such a representation can be generalised to a large class of fractional systems as demonstrated in [8, 9]. In these works, second relation in (7) is rewritten as
(8)y(t)=∫0∞μ(x)w(t,x)dx,
and representations (7) and (8) are denoted by diffusive representation. For a fractional transfer function defined by
(9)F(p)=B(p)A(p)
with B(p)=∑l=0rqlpβl and A(p)=∑k=0mrkpαk where βl+1≥βl≥0, function μ(x) is defined by [9]
(10)μ(x)=12iπ[F((-x)-)-F((-x)+)]=1π∑k=0m∑l=0qakqlsin((αk-βl)π)xαk+βl∑k=0mαk2x2αk+∑0≤k<l<m2akalcos((αk-αl)π)xαk+αl.Initial conditions are defined for system (7) by w(0,x)=ρ(x) and thus permits giving the exact expression of the system response with initial conditions [7]
(11)y(t)=∫0∞μ(x)(w(0,x)e-xt+∫0te-x(t-τ)u(τ)dτ)dx.Through several changes of variables described in [7], system (1), (but also a large number of fractional systems) can be described by
(12)∂ϕ(t,ζ)∂t=∂2ϕ(t,ζ)∂ζ2+u(t)δ(ζ),y(t)=∫-∞∞m(ζ)ϕ(t,ζ)dζ.Relation (10) shows that a fractional integrator can thus be seen as an infinite dimensional system described by a diffusion equation. This interpretation is represented by Figure 1 [7] in which(i)
the inputu(t) applied at the abscissa ζ=0,(ii)
the real distributed stateϕ(t,ζ),(iii)
the output (weighed sum on the state) appears. This remark can be generalized to a large number of fractional systems and thus demonstrates their link with diffusion equations.Figure 1
Representation of system (12) and thus system (1).Implementation of relation (10) requires the integral truncation that can be done as follows:(13)y(t)=∫-ζ2ζ1m(ζ)ϕ(t,ζ)dζ.In (13) m(ζ)=𝔉-1{4π2xμ(4π2x2)} where 𝔉-1 denotes the inverse Fourier transform. This relation is in practice impossible to compute analytically in most cases. To solve this problem, another partial differential equation is now proposed.


### 2.2. Another Partial Differential Equation Approximation

Using Laplace transform and introducingμ(x) function (μ(x)=sin(γπ)x-γ/π for a fractional integrator such as (1)), relation (3) becomes
(14)H(s)=∫0+∞μ(x)(s+x)dx.Using change of variablex=e-z, relation (15) becomes
(15)H(s)=∫-∞+∞μ(e-z)e-z(s+e-z)dz.Implementation of such a transfer function requires the integral truncation, namely,(16)H(s)≈∫ZiZfμ(e-z)(s/e-z+1)dz.Note thatZi∈ℝ, Zf∈ℝ are homogenous to the logarithm of a frequency. Now Let x(z,t) be a function of the space variable z of finite dimension (z∈⌊Zi⋯Zf⌋) and of the time variable t. This function satisfies the class of partial differential equations
(17)β(z)∂w(z,t)∂z+γ(z)∂2w(z,t)∂t∂z=u(t).
Also, let the system output y(t) be given by
(18)y(t)=w(Zf,t)-w(Zi,t)=∫ZiZf∂w(z,t)∂zdz.This partial differential equation class has been studied in [10]. Transfer function that links the system input and output is defined by
(19)G(s)=Y(s)U(s)=∫0Zβ-1(z)1+γ(z)β-1(z)sdz.Now if(20)β-1(z)=μ(e-z),γ-1(z)β(z)=e-zandthusγ(z)=ezμ(e-z).It turns out that the solution of the partial differential equation defined by (15) and (16) is an approximation of the response of the system whose transfer is defined by (14) (and thus by (1) for a fractional integrator). Such a representation can thus be used to approximate a fractional system or a fractional integrator. It is now used to implement a fractional diffusion equation.


## 2.1. Partial Differential Equation Representation and Approximation of a Fractional System

For presentation simplicity, the following fractional system (fractional integrator) is considered(1)H(s)=s-γ,
with 0<γ<1. Its link with diffusion equation can be demonstrated using the system impulse response [6] defined by the Mellin-Fourier integral of (1):
(2)h(t)=L-1{s-γ}=limω→∞12jπ∫c-j∞c+j∞ests-γds,
where c is greater than the abscissa of the singular points of H(s). Using poles definition that can be found in [6, 7], this system does not generate poles and its impulse response is thus given by
(3)h(t)=sin(γπ)π∫0∞x-γe-txdx.Response of system (1) to an input u(t) is defined as the convolution product of the impulse response h(t) with the input u(t):
(4)y(t)=∫0th(t-τ)u(τ)dτ,
and thus using relation (3) and through an integral permutation
(5)y(t)=∫0∞sin(γπ)πx-γ(∫0te-(t-τ)xu(τ)dτ)dx.Let(6)w(t,x)=∫0te-(t-τ)xu(τ)dτ
the following state space representation can be obtained for system (1):
(7)∂w(t,x)∂t=-xw(t,x)+u(t),y(t)=sin(γπ)π∫0∞x-γw(t,x)dx.x∈ℝ+.Such a representation can be generalised to a large class of fractional systems as demonstrated in [8, 9]. In these works, second relation in (7) is rewritten as
(8)y(t)=∫0∞μ(x)w(t,x)dx,
and representations (7) and (8) are denoted by diffusive representation. For a fractional transfer function defined by
(9)F(p)=B(p)A(p)
with B(p)=∑l=0rqlpβl and A(p)=∑k=0mrkpαk where βl+1≥βl≥0, function μ(x) is defined by [9]
(10)μ(x)=12iπ[F((-x)-)-F((-x)+)]=1π∑k=0m∑l=0qakqlsin((αk-βl)π)xαk+βl∑k=0mαk2x2αk+∑0≤k<l<m2akalcos((αk-αl)π)xαk+αl.Initial conditions are defined for system (7) by w(0,x)=ρ(x) and thus permits giving the exact expression of the system response with initial conditions [7]
(11)y(t)=∫0∞μ(x)(w(0,x)e-xt+∫0te-x(t-τ)u(τ)dτ)dx.Through several changes of variables described in [7], system (1), (but also a large number of fractional systems) can be described by
(12)∂ϕ(t,ζ)∂t=∂2ϕ(t,ζ)∂ζ2+u(t)δ(ζ),y(t)=∫-∞∞m(ζ)ϕ(t,ζ)dζ.Relation (10) shows that a fractional integrator can thus be seen as an infinite dimensional system described by a diffusion equation. This interpretation is represented by Figure 1 [7] in which(i)
the inputu(t) applied at the abscissa ζ=0,(ii)
the real distributed stateϕ(t,ζ),(iii)
the output (weighed sum on the state) appears. This remark can be generalized to a large number of fractional systems and thus demonstrates their link with diffusion equations.Figure 1
Representation of system (12) and thus system (1).Implementation of relation (10) requires the integral truncation that can be done as follows:(13)y(t)=∫-ζ2ζ1m(ζ)ϕ(t,ζ)dζ.In (13) m(ζ)=𝔉-1{4π2xμ(4π2x2)} where 𝔉-1 denotes the inverse Fourier transform. This relation is in practice impossible to compute analytically in most cases. To solve this problem, another partial differential equation is now proposed.


## 2.2. Another Partial Differential Equation Approximation

Using Laplace transform and introducingμ(x) function (μ(x)=sin(γπ)x-γ/π for a fractional integrator such as (1)), relation (3) becomes
(14)H(s)=∫0+∞μ(x)(s+x)dx.Using change of variablex=e-z, relation (15) becomes
(15)H(s)=∫-∞+∞μ(e-z)e-z(s+e-z)dz.Implementation of such a transfer function requires the integral truncation, namely,(16)H(s)≈∫ZiZfμ(e-z)(s/e-z+1)dz.Note thatZi∈ℝ, Zf∈ℝ are homogenous to the logarithm of a frequency. Now Let x(z,t) be a function of the space variable z of finite dimension (z∈⌊Zi⋯Zf⌋) and of the time variable t. This function satisfies the class of partial differential equations
(17)β(z)∂w(z,t)∂z+γ(z)∂2w(z,t)∂t∂z=u(t).
Also, let the system output y(t) be given by
(18)y(t)=w(Zf,t)-w(Zi,t)=∫ZiZf∂w(z,t)∂zdz.This partial differential equation class has been studied in [10]. Transfer function that links the system input and output is defined by
(19)G(s)=Y(s)U(s)=∫0Zβ-1(z)1+γ(z)β-1(z)sdz.Now if(20)β-1(z)=μ(e-z),γ-1(z)β(z)=e-zandthusγ(z)=ezμ(e-z).It turns out that the solution of the partial differential equation defined by (15) and (16) is an approximation of the response of the system whose transfer is defined by (14) (and thus by (1) for a fractional integrator). Such a representation can thus be used to approximate a fractional system or a fractional integrator. It is now used to implement a fractional diffusion equation.


## 3. Application to Fractional Partial Differential Equations Implementation

Implementation problem of a fractional partial differential equation using simulation software such as Comsol Multiphysics is now addressed and the following class of equation is considered:(21)∂1-γC1(t,x)∂t1-γ=-D∂2C1(t,x)∂x2,
with the following boundary and initial conditions:
(22)∂C1(t,0)∂x=∂C1(t,L)∂x=0,C1(0,x)=g(x).It is anADIb type anomalous diffusion equation that can also be rewritten as
(23)Y(t,x)=-D∂2C1(t,x)∂x2,
with
(24)Y(t,x)=∂1-γC1(t,x)∂t1-γ,(25)∂C1(t,0)∂x=∂C1(t,L)∂x=0,C1(0,x)=f(x).
According to the demonstration in Section 2 and in [11],  (1-γ) fractional derivative of C1(t,x) can be approximated by
(26)U(t,x)=∂C1(t,x)∂t,(27)U(t,x)=∂w(z,x,t)∂z1R(z)+C(z)∂2w(z,x,t)∂z∂t,(28)Y(t,x)=w(Zf,x,t)-w(Zi,x,t)=∫ZiZf∂w(z,x,t)∂zdz,
where functions C(z) and R(z) are, respectively, defined by
(29)C(z)=C(0)e-Az,R(z)=R(0)eBz.Implementation of (24) using approximations from (26) to (29) thus requires 2 geometries. In a first 1D-type geometry, (23) is implemented. To take into account relation (24), a second 2D-type geometry is created.ValuesZi and Zf for the 2D geometry along the z axis are used to define the range of frequency [ωi,ωf] for which the approximation of the fractional differentiation is expected using the relations
(30)ωi=eZi,ωf=eZf.
Note that frequencies ωi and ωf depend, respectively, on the simulation duration Td and the sampling time Ts. These frequencies can be defined using the following rules: ωi≪2π/Td and ωf≫2π/Ts.Information produced in the two geometries are then exchanged as described in Figure4.


## 4. Example

The following diffusion system withx∈[0⋯L] is considered:
(31)∂1-γC1(t,x)∂t1-γ=-D∂2C1(t,x)∂x2,
with the following initial and boundaries conditions:
(32)C1(t,0)=h(t),C1(t,L)=j(t),C1(0,x)=0.
For γ=0.5, D=0.1 m/s, j(t)=0 and
(33)h(t)=K·(t-t0)heaviside(t-t0)-K(t-t1)heaviside(t-t1)=h1(t)-h2(t)t≥0.
As shown in the appendix, system (31) and (32) solution is defined by
(34)C1(x,t)=0fort≤t0,(35)C1(x,t)=∑0∞(2·(1-e(nπ)4D2(t-t0)erfc((nπ)2D(t-t0)))(nπ)5D2-4·(t-t0)0.5(nπ)3D·Γ(0.5))×sin(nπx)+(1-x)·h1(t)fort0≤t≤t1,(36)C1(x,t)=∑0∞(2·(1-e(nπ)4D2(t-t0)erfc((nπ)2D(t-t0)))(nπ)5D2-4·(t-t0)0.5(nπ)3D·Γ(0.5)-2·(1-e(nπ)4D2(t-t1)erfc((nπ)2D(t-t1)))(nπ)5D2+4·(t-t1)0.5(nπ)3D·Γ(0.5))×sin(nπx)+(1-x)·(h1(t)-h2(t))fort≥t1.
For t1=5 s and t2=10 s, this solution is represented by Figure 5.The analytical solution is compared with the results produced by COMSOL Multiphysics for functionC1(x,t) using the implementation and geometries described in Figures 2–4. For the implementation, ωi=0.01 rd/s and ωf=1000 rd/s. Figure 6 is a representation of the absolute error between function C1(x,t) and the solution computed with COMSOL Multiphysics. This really small error permits us to validate the method we used for the implementation of a fractional diffusion equation using a partial differential equation for the approximation of a fractional derivative.Figure 2
1D geometry representation.Figure 3
2D geometry representation.Figure 4
Information exchanged between the two geometries.Figure 5
FunctionC1(x,t) produced by COMSOL Multiphysics.Figure 6
Absolute error between functionC1(x,t) and the solution computed with COMSOL Multiphysics.


## 5. Conclusion

This paper proposes a method for the implementation of a fractional diffusion equation into simulation softwares such as COMSOL Multiphysics. These software applications are now powerful tools for engineers to simulate complex systems combining several physical domains such as electrochemical and thermal. However they are not adapted to take into account anomalous diffusion and thus to model diffusion interfaces as in batteries, ultracapacitors, or fuel cells. To permit the implementation, the link between fractional systems and diffusion equation is used. The fractional diffusion equation considered is splitted into two parts and the remaining fractional equation is approximated by a partial differential equation. For the implementation of this partial differential equation, an additional geometry is created (a 1D system is transformed into a 2D system). The efficiency of the proposed method is evaluated in an example. The results obtained showed the efficiency of the proposed method.



---

*Source: 795651-2013-03-19.xml*