An Implementation Solution for Fractional Partial Differential Equations

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.

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 ⟨𝑟²⟩∼𝑡^𝛾 with 𝛾 between 0 and 2 (for classical diffusion ⟨𝑟²⟩∼𝑡).

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 Section 2, 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) 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): ℎ(𝑡)=𝐿⁻¹𝑠^−𝛾=lim_𝜔→∞12𝑗𝜋^𝑐+𝑗∞_{𝑐−𝑗∞}𝑒^𝑠𝑡𝑠^−𝛾𝑑𝑠,(2) where 𝑐 is greater than the abscissa of the singular points of 𝐻(𝑠). Using poles definition that can be found in [6, 7], this system does not generate poles and its impulse response is thus given by ℎ(𝑡)=sin𝛾𝜋𝜋^∞₀𝑥^−𝛾𝑒^−𝑡𝑥𝑑𝑥.(3)

Response of system (1) to an input 𝑢(𝑡) is defined as the convolution product of the impulse response ℎ(𝑡) with the input 𝑢(𝑡): 𝑦(𝑡)=^𝑡₀ℎ(𝑡−𝜏)𝑢(𝜏)𝑑𝜏,(4) and thus using relation (3) and through an integral permutation 𝑦(𝑡)=^∞₀sin𝛾𝜋𝜋𝑥^−𝛾𝑡₀𝑒^{−(𝑡−𝜏)𝑥}𝑢(𝜏)𝑑𝜏𝑑𝑥.(5)

Let 𝑤(𝑡,𝑥)=^𝑡₀𝑒^{−(𝑡−𝜏)𝑥}𝑢(𝜏)𝑑𝜏(6) the following state space representation can be obtained for system (1): 𝜕𝑤(𝑡,𝑥)()(𝑡),𝑦(𝑡)=𝜕𝑡=−𝑥𝑤𝑡,𝑥+𝑢sin𝛾𝜋𝜋∫^∞₀𝑥^−𝛾𝑤(𝑡,𝑥)𝑑𝑥.𝑥∈ℝ⁺.(7)

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) and representations (7) and (8) are denoted by diffusive representation. For a fractional transfer function defined by 𝐹𝑝=𝐵𝑝𝐴𝑝(9) with 𝐵(𝑝)=∑^𝑟_𝑙=0𝑞_𝑙𝑝^𝛽𝑙 and 𝐴(𝑝)=∑^𝑚_𝑘=0𝑟_𝑘𝑝^𝛼𝑘 where 𝛽_𝑙+1≥𝛽_𝑙≥0, function 𝜇(𝑥) is defined by [9] 𝜇(𝑥)=1𝐹)2𝑖𝜋((−𝑥⁻)−𝐹(−𝑥)⁺=1𝜋∑^𝑚_𝑘=0∑^𝑞_𝑙=0𝑎_𝑘𝑞_𝑙sin𝛼_𝑘−𝛽_𝑙𝜋𝑥^{𝛼𝑘+𝛽𝑙}∑^𝑚_𝑘=0𝛼²_𝑘𝑥^2𝛼𝑘+∑_{0≤𝑘<𝑙<𝑚}2𝑎_𝑘𝑎_𝑙cos𝛼_𝑘−𝛼_𝑙𝜋𝑥^{𝛼𝑘+𝛼𝑙}.(10)

Initial conditions are defined for system (7) by 𝑤(0,𝑥)=𝜌(𝑥) and thus permits giving the exact expression of the system response with initial conditions [7] 𝑦(𝑡)=^∞₀𝜇(𝑥)𝑤(0,𝑥)𝑒^−𝑥𝑡+^𝑡₀𝑒^{−𝑥(𝑡−𝜏)}𝑢(𝜏)𝑑𝜏𝑑𝑥.(11)

Through several changes of variables described in [7], system (1), (but also a large number of fractional systems) can be described by 𝜕𝜙(𝑡,𝜁)=𝜕𝜕𝑡²𝜙(𝑡,𝜁)𝜕𝜁²(𝑡)𝛿(𝜁),𝑦(𝑡)=+𝑢_∞−∞𝑚(𝜁)𝜙(𝑡,𝜁)𝑑𝜁.(12)

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 input 𝑢(𝑡) applied at the abscissa 𝜁=0,(ii)the real distributed state 𝜙(𝑡,𝜁),(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:𝑦(𝑡)=^𝜁1_−𝜁2𝑚(𝜁)𝜙(𝑡,𝜁)𝑑𝜁.(13)

In (13) 𝑚(𝜁)=𝔉⁻¹{4𝜋²𝑥𝜇(4𝜋²𝑥²)} where 𝔉⁻¹ 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 𝜇(𝑥) function (𝜇(𝑥)=sin(𝛾𝜋)𝑥^−𝛾/𝜋 for a fractional integrator such as (1)), relation (3) becomes 𝐻(𝑠)=₀^+∞𝜇(𝑥)(𝑠+𝑥)𝑑𝑥.(14)

Using change of variable 𝑥=𝑒^−𝑧, relation (15) becomes 𝐻(𝑠)=^+∞_−∞𝜇𝑒^−𝑧𝑒^−𝑧(𝑠+𝑒^−𝑧)𝑑𝑧.(15)

Implementation of such a transfer function requires the integral truncation, namely, 𝐻(𝑠)≈^𝑍𝑓_𝑍𝑖𝜇𝑒^−𝑧(𝑠/𝑒^−𝑧+1)𝑑𝑧.(16)

Note that 𝑍_𝑖∈ℝ, 𝑍_𝑓∈ℝ are homogenous to the logarithm of a frequency. Now Let 𝑥(𝑧,𝑡) be a function of the space variable 𝑧 of finite dimension (𝑧∈⌊𝑍_𝑖⋯𝑍_𝑓⌋) and of the time variable 𝑡. This function satisfies the class of partial differential equations 𝛽(𝑧)𝜕𝑤(𝑧,𝑡)(𝑧)𝜕𝜕𝑧+𝛾²𝑤(𝑧,𝑡)(𝑡).𝜕𝑡𝜕𝑧=𝑢(17) Also, let the system output 𝑦(𝑡) be given by 𝑦(𝑡)=𝑤𝑍_𝑓,𝑡−𝑤𝑍_𝑖,𝑡=^𝑍𝑓_𝑍𝑖𝜕𝑤(𝑧,𝑡)𝜕𝑧𝑑𝑧.(18)

This partial differential equation class has been studied in [10]. Transfer function that links the system input and output is defined by 𝐺(𝑠)=𝑌(𝑠)𝑈(𝑠)=^𝑍₀𝛽⁻¹(𝑧)1+𝛾(𝑧)𝛽⁻¹(𝑧)𝑠𝑑𝑧.(19)Now if𝛽⁻¹(𝑧)=𝜇𝑒^−𝑧,𝛾⁻¹(𝑧)𝛽(𝑧)=𝑒^−𝑧(𝑧)=𝑒andthus𝛾^𝑧𝜇(𝑒^−𝑧).(20)

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: 𝜕^1−𝛾𝐶₁(𝑡,𝑥)𝜕𝑡^1−𝛾𝜕=−𝐷²𝐶₁(𝑡,𝑥)𝜕𝑥²,(21) with the following boundary and initial conditions: 𝜕𝐶₁(𝑡,0)=𝜕𝑥𝜕𝐶₁(𝑡,𝐿)𝜕𝑥=0,𝐶₁(0,𝑥)=𝑔(𝑥).(22)

It is an ADIb type anomalous diffusion equation that can also be rewritten as 𝑌()𝜕𝑡,𝑥=−𝐷²𝐶₁(𝑡,𝑥)𝜕𝑥²,(23) with 𝑌()=𝜕𝑡,𝑥^1−𝛾𝐶₁(𝑡,𝑥)𝜕𝑡^1−𝛾,(24)𝜕𝐶₁(𝑡,0)=𝜕𝑥𝜕𝐶₁(𝑡,𝐿)𝜕𝑥=0,𝐶₁(0,𝑥)=𝑓(𝑥).(25) According to the demonstration in Section 2 and in [11], (1−𝛾) fractional derivative of 𝐶₁(𝑡,𝑥) can be approximated by 𝑈(𝑡,𝑥)=𝜕𝐶₁(𝑡,𝑥),𝑈()=()𝜕𝑡(26)𝑡,𝑥𝜕𝑤𝑧,𝑥,𝑡1𝜕𝑧𝑅(𝑧)(𝑧)𝜕+𝐶²𝑤(𝑧,𝑥,𝑡),𝑌()𝑍𝜕𝑧𝜕𝑡(27)𝑡,𝑥=𝑤_𝑓,𝑥,𝑡−𝑤𝑍_𝑖,𝑥,𝑡=^𝑍𝑓_𝑍𝑖𝜕𝑤(𝑧,𝑥,𝑡)𝜕𝑧𝑑𝑧,(28) where functions 𝐶(𝑧) and 𝑅(𝑧) are, respectively, defined by 𝐶(𝑧)=𝐶(0)𝑒^−𝐴𝑧,𝑅(𝑧)=𝑅(0)𝑒^𝐵𝑧.(29)

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.

Values 𝑍_𝑖 and 𝑍_𝑓 for the 2D geometry along the 𝑧 axis are used to define the range of frequency [𝜔_𝑖,𝜔_𝑓] for which the approximation of the fractional differentiation is expected using the relations 𝜔_𝑖=𝑒^𝑍𝑖,𝜔_𝑓=𝑒^𝑍𝑓.(30) Note that frequencies 𝜔_𝑖 and 𝜔_𝑓 depend, respectively, on the simulation duration 𝑇_𝑑 and the sampling time 𝑇_𝑠. These frequencies can be defined using the following rules: 𝜔_𝑖≪2𝜋/𝑇_𝑑 and 𝜔_𝑓≫2𝜋/𝑇_𝑠.

Information produced in the two geometries are then exchanged as described in Figure 4.

4. Example

The following diffusion system with 𝑥∈[0⋯𝐿] is considered: 𝜕^1−𝛾𝐶₁(𝑡,𝑥)𝜕𝑡^1−𝛾𝜕=−𝐷²𝐶₁(𝑡,𝑥)𝜕𝑥²,(31) with the following initial and boundaries conditions: 𝐶₁(𝑡,0)=ℎ(𝑡),𝐶₁()(𝑡),𝐶𝑡,𝐿=𝑗₁(0,𝑥)=0.(32) For 𝛾=0.5, 𝐷=0.1 m/s, 𝑗(𝑡)=0 and ℎ(𝑡)=𝐾⋅𝑡−𝑡₀heaviside𝑡−𝑡₀−𝐾𝑡−𝑡₁heaviside𝑡−𝑡₁=ℎ₁(𝑡)−ℎ₂(𝑡)𝑡≥0.(33) As shown in the appendix, system (31) and (32) solution is defined by 𝐶₁(𝑥,𝑡)=0for𝑡≤𝑡₀,𝐶(34)₁()=𝑥,𝑡_∞02⋅1−𝑒^{(𝑛𝜋)4𝐷2(𝑡−𝑡0)}erfc(𝑛𝜋)²𝐷𝑡−𝑡₀(𝑛𝜋)⁵𝐷²−4⋅𝑡−𝑡₀^0.5(𝑛𝜋)³()()+()𝐷⋅Γ0.5×sin𝑛𝜋𝑥1−𝑥⋅ℎ₁(𝑡)for𝑡₀≤𝑡≤𝑡₁,𝐶(35)₁()=𝑥,𝑡_∞02⋅1−𝑒^{(𝑛𝜋)4𝐷2(𝑡−𝑡0)}erfc(𝑛𝜋)²𝐷𝑡−𝑡₀(𝑛𝜋)⁵𝐷²−4⋅𝑡−𝑡₀^0.5(𝑛𝜋)³𝐷⋅Γ(0.5)−2⋅1−𝑒^{(𝑛𝜋)4𝐷2(𝑡−𝑡1)}erfc(𝑛𝜋)²𝐷𝑡−𝑡₁(𝑛𝜋)⁵𝐷²+4⋅𝑡−𝑡₁^0.5(𝑛𝜋)³()()+()⋅ℎ𝐷⋅Γ0.5×sin𝑛𝜋𝑥1−𝑥₁(𝑡)−ℎ₂(𝑡)for𝑡≥𝑡₁.(36) For 𝑡₁=5 s and 𝑡₂=10 s, this solution is represented by Figure 5.

The analytical solution is compared with the results produced by COMSOL Multiphysics for function 𝐶₁(𝑥,𝑡) using the implementation and geometries described in Figures 2–4. For the implementation, 𝜔_𝑖=0.01 rd/s and 𝜔_𝑓=1000 rd/s. Figure 6 is a representation of the absolute error between function 𝐶₁(𝑥,𝑡) 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 

Function 𝐶1(𝑥,𝑡) produced by COMSOL Multiphysics.

Figure 6 

Absolute error between function 𝐶1(𝑥,𝑡) 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.

Appendix

This appendix demonstrates how, using material provided in [12], the analytical solution of the following fractional diffusion equation: 𝜕^1−𝛾𝐶₁(𝑡,𝑥)𝜕𝑡^1−𝛾𝜕=−𝐷²𝐶₁(𝑡,𝑥)𝜕𝑥²,𝑥∈[0⋯𝐿](A.1) with the following initial and boundaries conditions 𝐶₁(𝑡,0)=ℎ(𝑡),𝐶₁(𝑡,𝐿)=𝑗(𝑡),𝐶₁(0,𝑥)=0,(A.2) is obtained. To obtain homogenous conditions at 𝑥=0 and 𝑥=1, the following change of variable is used: 𝑉(𝑥,𝑡)=𝐶₁(𝑥,𝑡)+𝑈(𝑥,𝑡)(A.3) with 𝑈()=𝑥𝑥,𝑡1−𝐿ℎ(𝑡)+𝑥𝐿𝑗(𝑡).(A.4) System (A.1) thus becomes 𝜕^1−𝛾𝑉(𝑡,𝑥)𝜕𝑡^1−𝛾𝜕=−𝐷²𝑉(𝑡,𝑥)𝜕𝑥²(),𝑉()(𝐿)()(𝑥),()+𝑓𝑥,𝑡𝑡,0=0,𝑉=0,𝑉0,𝑥=𝑔A.5 with 𝑓()𝑥𝑥,𝑡=−1−𝐿𝜕^1−𝛾ℎ(𝑡)𝜕𝑡^1−𝛾−𝑥𝐿𝜕^1−𝛾𝑗(𝑡)𝜕𝑡^1−𝛾,𝑔(𝑥)𝑥=−1−𝐿ℎ(0)−𝑥𝐿𝑗(0).(A.7) Separation variable method leads to writing 𝑉(𝑥,𝑡) as: 𝑉(𝑥,𝑡)=^∞₁𝑇_𝑛(𝑡)sin(𝑛𝜋𝑥)(A.9) with 𝜕^1−𝛾𝑇_𝑛(𝑡)𝜕𝑡^1−𝛾+(𝑛𝜋)²𝐷𝑇_𝑛(𝑡)=𝑓_𝑛(𝑡),()𝑇A.10_𝑛(0)=2¹₀𝑔(𝜁)()()𝑓sin𝑛𝜋𝜁𝑑𝜁,A.11_𝑛(𝑡)=2¹₀𝑓(𝑥,𝑡)sin(𝑛𝜋𝑥)𝑑𝑥.(A.12) Now let   𝛾=0.5, 𝑗(𝑡)=0, and ℎ(𝑡)=0,𝑡≤𝑡₀,ℎ(𝑡)=𝐾𝑡−𝑡₀,𝑡≤𝑡≤𝑡₁,ℎ(𝑡)=𝐾𝑡₁−𝑡₀,𝑥≥𝑡₁.(A.13) Solution of system (A.1) is given, according to (A.9), (A.10), and (A.11)𝑉(𝑥,𝑡)=^∞₁𝑇_𝑛(𝑡)(),()𝜕sin𝑛𝜋𝑥A.16_1−𝛾𝑇_𝑛(𝑡)𝜕𝑡^1−𝛾+(𝑛𝜋)²𝐷𝑇_𝑛(𝑡)=𝑓_𝑛(𝑡),()𝑇A.17_𝑛(0)=0(A.18) with 𝑓_𝑛(𝑡)=2¹₀𝑓()()()𝑓()()𝜕𝑥,𝑡sin𝑛𝜋𝑥𝑑𝑥,A.19𝑥,𝑡=−1−𝑥^1−𝛾ℎ(𝑡)𝜕𝑡^1−𝛾,()𝑔(𝑥)()ℎ(0)()𝑉()A.20=−1−𝑥=0,A.21𝑥,𝑡=𝐶₁()(),()𝑈()=()ℎ(𝑡).()𝑥,𝑡+𝑈𝑥,𝑡A.22𝑥,𝑡1−𝑥A.23 Combining relations (A.19) and (A.20) leads to 𝑓_𝑛(𝑡)𝜕=−2^1−𝛾ℎ(𝑡)𝜕𝑡^1−𝛾1₀(1−𝑥)sin(𝑛𝜋𝑥)𝑑𝑥,(A.24) and thus 𝑓_𝑛(𝑡)4=−𝜕𝑛𝜋^1−𝛾ℎ(𝑡)𝜕𝑡^1−𝛾.(A.25) Laplace transform applied to (A.17) leads to 𝑝^0.5𝑇_𝑛𝑝+(𝑛𝜋)²𝐷𝑇_𝑛𝑝=4𝑝𝑛𝜋^0.5ℎ𝑝,(A.26) and thus 𝑇_𝑛𝑝=(4/𝑛𝜋)𝑝^0.5ℎ𝑝𝑝^0.5+(𝑛𝜋)²𝐷.(A.27) If ℎ(𝑡) is written as ℎ(𝑡)=𝐾𝑡−𝑡₀heaviside𝑡−𝑡₀−𝐾𝑡−𝑡₁heaviside𝑡−𝑡₁=ℎ₁(𝑡)−ℎ₂(𝑡)(),𝑡≥0,A.28 relation (A.27) becomes 𝑇_𝑛𝑝=(4/𝑛𝜋)𝑝^0.5ℎ₁𝑝𝑝^0.5+(𝑛𝜋)²𝐷−(4/𝑛𝜋)𝑝^0.5ℎ₂𝑝𝑝^0.5+(𝑛𝜋)²𝐷=𝑇_𝑛1𝑝−𝑇_𝑛2𝑝.(A.29) Within the interval 0<𝑡<𝑡₀, ℎ₁(𝑡)=0 and thus using (A.27)  𝑇_𝑛1(𝑡)=0.

If 𝑡>𝑡₀, function ℎ(𝑡) is a ramp and thus 𝑇_𝑛1𝑝=(4/𝑛𝜋)𝑝^0.5𝐾/𝑝²𝑝^0.5+(𝑛𝜋)²𝐷.(A.30) Inverse Laplace transform then permits 𝑇_𝑛1(𝑡)=2⋅1−𝑒^{(𝑛𝜋)4𝐷2(𝑡−𝑡0)}erfc(𝑛𝜋)²𝐷𝑡−𝑡₀(𝑛𝜋)⁵𝐷²−4⋅𝑡−𝑡₀^0.5(𝑛𝜋)³().()𝐷⋅Γ0.5A.31 Within the interval 0<𝑡<𝑡₁, ℎ₂(𝑡)=0, and thus 𝑇_𝑛2(𝑡)=0.

Using a similar method, 𝑇_𝑛2(𝑡)   where   𝑡>𝑡₁ is given by 𝑇_𝑛2(𝑡)=2⋅1−𝑒^{(𝑛𝜋)4𝐷2(𝑡−𝑡1)}erfc(𝑛𝜋)²𝐷𝑡−𝑡₁(𝑛𝜋)⁵𝐷²−4⋅𝑡−𝑡₁^0.5(𝑛𝜋)³().()𝐷⋅Γ0.5A.32 Finally, using (A.3), system (A.1) and (A.2) solution is defined by, using (34), (35), and (36).

Conflict of Interests

The authors declare no conflict of interests.