Grünwald-Letnikov Fractional Derivative Applied to First-Order Ordinary Differential Equations

Scotton, J. W.;

DOI 10.5433/1679-0375.2024.v45.51533

Citation Semin., Ciênc. Exatas Tecnol. 2024, v. 45: e51533

Received: 1 October 2024 Received in revised for: 12 November 2024 Accepted: 28 November 2024 Available online: 14 December 2024

Abstract:

First-order ordinary differential equations (ODEs) are widely used in various fields of science and engineering to model natural phenomena. This study proposes an extension of these equations using fractional derivatives, specifically the Grünwald-Letnikov definition, to explore their impact on the behavior of solution curves. The fractional ODE considered is discretized using the finite difference method and solved numerically for different values of the derivative order (alpha). Tests were conducted to verify mesh independence and the quality of the computational implementation of the method, through which the accuracy and the absence of implementation errors were confirmed. The behavior of the solution curves for different values of alpha was analyzed, revealing a sharp decrease near the initial point (t = 0) and an almost linear growth at higher values of $t$, within the considered domain. Additionally, in solving a specific initial value problem with a known analytical solution, it was discovered that the accuracy of the numerical solutions for higher values of alpha was more dependent on the mesh refinement than the solutions for lower values.

Keywords: differential equations, fractional calculus, Grünwald-Letnikov

Introduction

First-order Ordinary Differential Equations (ODEs) are important mathematical tools, with applications in various fields of science and engineering. In physics, for instance, they are often used to describe the motion of particles; in biology, in population growth models; in engineering, in processes involving heat and mass transfer; among many other different applications.

These equations relate an unknown function to its first-order derivative and are generally represented in the form of equation (1)

$$\frac{dy}{dt}=f(t,y),$$

where y is a function of t. Such equations are said to be linear when f depends linearly on the dependent variable y and are often written in their standard form given by equation (2)

$$\frac{dy}{dt}+p(t)y=g(t),$$

where p and g are given functions of the independent variable t. Any differentiable function $y=\psi(t)$ that satisfies equation (2) for all t in some interval is called a solution .

In this work, however, instead of fixing a first-order derivative in the ODE, a non-integer order derivative $\alpha \in (0,1)$ is considered, such that equation (2) takes the form of equation (3)

$$\frac{d^\alpha y}{dt^\alpha}+p(t)y=g(t).$$

It happens that non-integer order derivatives, or fractional order derivatives as they are generally called, have different definitions, such as those of Riemann-Liouville, Caputo, Riesz, and Grünwald-Letnikov, among others, which have been applied to various problems in recent years.

Mark M. Meerschaert & Charles Tadjeran, 2004 developed a numerical method for solving the transient one-dimensional advection-diffusion equation, incorporating a fractional-order derivative in the diffusive term. Santanu, 2009 studied the problem of transient one-dimensional and two-dimensional diffusion with Riemann-Liouville space-fractional derivative.

Giulio Cottone & Mario Di Paola, 2011, for example, proposed a model with the Riesz fractional derivative to simulate wind velocity fields. A. Simmons et al., 2017 presented a new discretization, based on the Finite Volume Method, of the transient one-dimensional diffusion equation, with a fractional-order derivative in space. A. G. O. Goulart et al., 2017 applied the Caputo derivative to studies of pollutant dispersion in the atmosphere.

Davidson Moreira et al., 2019 presented a numerical solution for the three-dimensional advection-diffusion equation, employing Caputo fractional derivatives, applied to the dispersion of pollutants in the atmosphere. J. W. Scotton et al., 2020 investigated the problem of transient one-dimensional diffusion using Fractional Calculus. P. Zhang et al., 2021 employed Riemann-Liouville fractional derivatives in studies of diffusion and reaction processes in batch reactors.

M. Luft et al., 2021 applied Riemann-Liouville fractional derivatives in studies involving the dynamics of pneumatic systems. R. Banerjee & R. K. Biswas, 2022 proposed fractional models, relating Caputo and Grünwald-Letnikov derivatives, to describe the dynamics of COVID-19. C. Y. Kee et al., 2022 apply fractional derivatives to an urban growth model. M. Qayyum et al., 2023 present the fractional modeling of the flow of a non-Newtonian fluid between parallel plates. K. S. Nisar et al., 2023 reviewed recent studies based on fractional modeling of infectious and non-infectious diseases using different fractional operators.

All the cited works found that fractional models provided excellent results, generally better than models with integer-order derivatives, something that has led more and more researchers to become interested in the subject and to conduct studies applying non-integer order derivatives to classical models to describe different phenomena. Thus, studies that provide information on the impact of derivative orders on solution curves are important so that, once the limitations of classical models in representing certain phenomena are recognized, other more assertive options in the field of Fractional Calculus can be sought.

Among the aforementioned definitions, one with great numerical importance is that of Grünwald-Letnikov. The right and left fractional derivatives of order $\alpha$ in an interval $[a,b]$, according to the Grünwald-Letnikov definition, are given respectively by equation (4)

$$f_{a+}^{(\alpha)}(x)= \lim_{\Delta x \rightarrow{0^+}}{\frac{1}{\Delta x^\alpha} \sum_{k=0}^{\frac{x-a}{\Delta x}}{(-1)^k \dbinom{\alpha}{k}f(x-k\Delta x) }},$$

where $\alpha \in \mathbb{R}_+^*, x>a$ and equation (5)

$$f_{b-}^{(\alpha)}(x)= \lim_{\Delta x \rightarrow{0^+}}{\frac{1}{\Delta x^\alpha} \sum_{k=0}^{\frac{b-x}{\Delta x}}{(-1)^k \dbinom{\alpha}{k}f(x+k\Delta x) }},$$

where \(\alpha \in \mathbb{R}_+^*, x.

In this work, the Grünwald-Letnikov definition is applied to the first-order linear ordinary differential equation and solve the resulting equation numerically using the finite difference method. The objective is to verify the behavior of the solution curves for different values of the derivative order, in order to support new applications of these differential equations to situations or problems where modeling based on the classical equation with integer-order derivatives does not provide results with the required accuracy.

Materials and methods

Discretization of the equation with integer-order derivative   xxxxxxxxxxxxxxx

Consider the first-order linear ordinary differential equation given by equation (2), where $y=y(t)$, for $t \in [0,t_f]$ with $y(0)=y_0$. We can discretize the domain such that $t_i=ih, i=0, 1, ..., N$ with $h=t_f/N$. Thus, approximating the derivative by backward finite difference, the equation (2) can be rewritten as $$\frac{y(t_i)-y(t_{i-1})}{h}+p(t_i)y(t_i)=g(t_i).$$

Multiplying the entire equation (6) by $h$, we obtain $$y(t_i)-y(t_{i-1})+hp(t_i)y(t_i)=hg(t_i).$$

Finally, by rearranging the terms in equation (7), we arrive at the linear system given by equation (8) $$-y(t_{i-1})+[1+hp(t_i)]y(t_i)=hg(t_i),$$ for $i=1,2,...,N$ where $y(t_0)=y_0$.

Discretization of the equation with fractional derivative  xxxxxxxxxxxxxxxxxxx

We can approximate the order $\alpha$ derivative in equation (3) based on the Grünwald-Letnikov definition, on the left, such that $$\frac{d^\alpha y}{dt^\alpha} \approx \frac{1}{h^\alpha} \sum_{j=0}^{t/h}{w_j^\alpha y(t-jh)},$$ where the coefficients of equation (9) are given by the expressions presented in equation (10) $$w_0^\alpha=1, \quad w_j^\alpha=\left(1-\frac{\alpha+1}{j}\right)w_{j-1}^\alpha, \quad j=1,2,... .$$

Since the discretized domain, as in the case of the ODE with integer-order derivative, is $t_i=ih$, $i~=~0, 1, ..., N$ with $h=t_f/N$, the equation (9) becomes $$\frac{1}{h^\alpha} \left[w_0^\alpha y(t_i) + w_1^\alpha y(t_{i-1}) + \cdots + w_{i-1} y(t_1) + w_i^\alpha y(t_0) \right] + p(t_i) y(t_i) = g(t_i).$$

Rewriting the equation (11), we obtain $$(w_i^\alpha) y(t_0) + (w_{i-1}) y(t_1) + \cdots + (w_1^\alpha) y(t_{i-1}) + (w_0^\alpha + h^\alpha p(t_i)) y(t_i) = h^\alpha g(t_i),$$ for $i=1,2,\cdots,N$, where $y(t_0)=y_0$.

Note that for $\alpha=1$, the coefficients $w_j^\alpha$ in equation (12) becomes $$w_0^\alpha = 1, \quad w_1^\alpha = -1, \quad w_j^\alpha = 0, \quad j=2,3,\ldots .$$

Thus, from equation (12) and equation (13), we obtain the linear system given by equation (14) $$-y(t_{i-1})+[1+hp(t_i)]y(t_i)=hg(t_i),$$ for $i=1,2,...,N$ where $y(t_0)=y_0$, which is exactly the linear system resulting from the discretization of the classical ODE with integer-order derivative.

Existence and uniqueness of solution of the linear system   xxxxxxxxxxxxxxxxxx

The linear system given by equation (12) can be rewritten as

$$\begin{cases} y(t_0) = y_0 \\ (w_1^\alpha)y(t_0) + [w_0^\alpha + h^\alpha p(t_1)]y(t_1) = h^\alpha g(t_1) \\[1ex] (w_2^\alpha)y(t_0) + (w_1^\alpha)y(t_1) + [w_0^\alpha + h^\alpha p(t_2)]y(t_2) = h^\alpha g(t_2) \\[2ex] \vdots \\[2ex] (w_N^\alpha)y(t_0) + \cdots + (w_1^\alpha)y(t_{N-1}) + [w_0^\alpha + h^\alpha p(t_N)]y(t_N) = h^\alpha g(t_N) \end{cases},$$

which is a lower triangular system, as can be more clearly observed in matrix form

$$\begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\[1ex] w_1^\alpha & w_0^\alpha + h^\alpha p(t_1) & 0 & \cdots & 0 \\[1ex] w_2^\alpha & w_1^\alpha & w_0^\alpha + h^\alpha p(t_2) & \cdots & 0 \\[1ex] w_3^\alpha & w_2^\alpha & w_1^\alpha & \cdots & 0 \\[1ex] \vdots & \vdots & \ddots & \vdots \\[1ex] w_N^\alpha & w_{N-1}^\alpha & w_{N-2}^\alpha & \cdots & w_0^\alpha + h^\alpha p(t_N) \end{bmatrix} \cdot \begin{bmatrix} y(t_0) \\[1ex] y(t_1) \\[1ex] y(t_2) \\[1ex] y(t_3) \\[1ex] \vdots \\[1ex] y(t_N) \end{bmatrix} = \begin{bmatrix} y_0 \\[1ex] h^\alpha g(t_1) \\[1ex] h^\alpha g(t_2) \\[1ex] h^\alpha g(t_3) \\[1ex] \vdots \\[1ex] h^\alpha g(t_N) \end{bmatrix}.$$

Note that the determinant of the coefficient matrix is equal to the product of the elements on its main diagonal. Thus, this value will be zero if, and only if, $w_0^\alpha + h^\alpha p(t_i) = 0$ for some $i = 1, 2, \ldots, N$. Therefore, since $w_0^\alpha = 1$, choosing $h$ such that $h^\alpha \neq -1/p(t_i)$ for $i = 1, 2, \ldots, N$, ensures that the linear system has a unique solution. Furthermore, since $h$ must always be positive, if $p(t_i) > 0$ for $i = 1, 2, \ldots, N$, the existence and uniqueness of the solution are guaranteed, regardless of the choice of $h$.

Results and discussion

It is known that the number of points used in the discretization processes of differential equations impacts the accuracy of the numerical solutions obtained. For example, the approximation of the first derivative using the backward finite difference scheme has an error of order $O(h^2)$. In this context, to mitigate the influence of the mesh size on the solutions, a mesh independence test was performed using different numbers of points ($N$=5, 10, 20, 40, 80, 160). To facilitate understanding, we call $M_1$ the mesh containing 5 points, $M_2$ the mesh containing 10 points, $M_3$ the mesh containing 20 points, $M_4$ the mesh containing 40 points, $M_5$ the mesh containing 80 points, and $M_6$ the mesh containing 160 points.

For this preliminary analysis, the equation (3) is considered with $\alpha=0.999$, $t\in[0,1]$, $y(0)=1$, $p(t)=1/2$, $g(t)=(1/2)e^{t/3}$, and the results are evaluated for $y(1)$. Accordingly, in Figure 1, we can observe the results obtained for the considered problem, where we see that the variations between the obtained results decreased as the number of points increased.

The results shown in Figure 1 are even more evident in Table 1, which presents the relative errors between the meshes, calculated by equation (15)

$$(ER)_{i+1}=\frac{|M_{i+1}-M_{i}|}{|M_{i+1}|},$$

where $i=1,2,3,4,5$ and each $(ER)_{i+1}$ expresses the relative variation between the values of $y(1)$ obtained using mesh $i+1$ and mesh $i$. The criterion established to consider the mesh $M_i$ independent was $(ER)_{i+1}~<~10^{-3}$. In this context, the chosen mesh was $M_5$, with 80 points.

Therefore, to verify the possible existence of implementation errors, using the mesh $M_5$, the described problem was numerically solved, considering the classical equation ($\alpha=1$) and the fractional equation ($\alpha=0.9999$). In this process, as we can observe in Figure 2, the numerical results were very close to those provided by the analytical solution, indicating the absence of operational errors. Indeed, with the implementation correctly carried out, it was expected that for values of $\alpha$ close to 1, the obtained solutions would be close to those provided by the classical equation, which was indeed verified.

After verifying the mesh and the computational implementation of the discretized equation, the next step was to perform simulations to examine the impact of different values of $\alpha$ on the behavior of the solution curves. For this analysis, the following values of $\alpha$ were used: $0.9999, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65$, and $0.60$.

In Figure 3, the results obtained can be observed. It shows that as the value of $\alpha$ decreases, the solution curves exhibit an intense decrease for values of $t$ close to the initial point $t=0$ and an almost linear increase for higher values of $t$.

The linear behavior of the solution curves from certain values of $t$ can be observed more clearly in Figure 4, where the values of $y(t)$ for $0.5 \leq t \leq 10$ are presented.

To further highlight the linear behavior of the solutions, we can examine the coefficients $a$ and $b$ obtained through linear fittings of the type $g(t)=at+b$, performed using the least squares method. These coefficients are presented in Table 2, along with the respective values of the coefficient of determination $\text{R}^2$.

In Table 2, we can see that the $\text{R}^2$ values are all above 0.99, indicating that the data are extremely well-fitted by the linear function.

In Figure 5, the solutions obtained for different values of $\alpha$ at $t=0.25$, $t=0.5$, $t=0.75$ and $t=1$ can be observed. We see that, for all the analyzed values of $t$, the values of $y$ tend to increase as the value of $\alpha$ increases. However, it is also evident that the growth rate is less pronounced for lower values of $\alpha$ and more prominent for higher values of the order of the derivative.

(a)	(b)
(c)	(d)

Another case considered in this work was the initial value problem given by equation (3) and $y(0)=0$ where $$p(t)=1, \quad g(t)=\frac{\Gamma(2\alpha + 1)}{\Gamma(\alpha + 1)}t^{\alpha}+t^{2\alpha}.$$

In this case, by using the functions $p$ and $g$ from equation (16), the analytical solution of the fractional equation is $y(t)=t^{2\alpha}$. To numerically simulate this problem, two meshes were considered, with $N=10$ and $N=20$ points in the discretization, and the following values of $\alpha$ were used: 0.999, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, and 0.3.

Table 3 presents the maximum relative errors calculated between the analytical solutions and the numerical solutions obtained for different values of $\alpha$, calculated by equation (17)

$$(ER)_{max}=\max {\frac{|y_{AS}-y_{NS}|}{|y_{AS}|}},$$

By inspecting Table 3, we can see that the maximum relative errors decrease as the values of $\alpha$ also decrease, indicating that the results obtained for $\alpha$ close to 1 are more affected by the number of mesh points than the results obtained for lower values of $\alpha$. For $\alpha=0.3$, for example, even using the mesh with 10 points, we observed smaller errors compared to the results for $\alpha \geq 0.4$, using the mesh with 20 points, a pattern that repeats for other values of $\alpha$. Figures 6 and 7, in turn, display the graphs of the results for the meshes with $N=10$ and $N=20$, respectively.

(a)	(b)
(c)	(d)
(e)	(f)
(g)	(h)

(a)	(b)
(c)	(d)
(e)	(f)
(g)	(h)

The fact that mesh refinement has a greater impact on the results for larger values of $\alpha$ in terms of precision may initially seem surprising. However, if we consider that the non-local character of fractional derivatives intensifies for smaller values of $\alpha$, this makes more sense. This is because the coefficients $w_j^\alpha$, which can be considered as weights given to the points in the calculation of the derivatives, increase in absolute value for higher values of $j$ as $\alpha$ decreases. Therefore, we can infer that the errors tend to be smaller since the more distant points of the mesh are used with increasingly larger weights in the calculation of the derivatives for these values of $\alpha$.

Conclusions

In this work, the Grünwald-Letnikov definition was applied, generalizing the order $\alpha$ of the derivative present in the first-order ordinary differential equation from $\alpha=1$ to $\alpha \in (0, 1)$, and the resulting equation was numerically solved to verify the impact of this generalization on the solution curves.

Initially, tests were conducted to verify mesh independence and the quality of the computational implementation of the method. These tests confirmed the accuracy and the absence of implementation errors. Regarding the main results obtained, it was observed that as $\alpha$ decreased, the solution curves exhibited a decreasing behavior for values of $t$ close to the initial point $t=0$ and a linear growth behavior for higher values of $t$. This was confirmed by linear fittings, where the determination coefficients ($\text{R}^2$) calculated were greater than 0.99.

In the study, a specific initial value problem was also considered, for which the analytical solution of the fractional equation is known. Through the numerical solution of this problem, it was verified that the maximum relative errors between the analytical and numerical solutions decreased as the value of $\alpha$ decreased, indicating that the accuracy of the numerical solutions was more affected by the number of mesh points for larger values of $\alpha$. This reduction in errors was attributed to the fact that the non-local character of the fractional derivatives intensifies for smaller values of $\alpha$ since, for smaller values, greater weights are given to the more distant points of the mesh in the derivative calculations.

In summary, this work presented a numerical methodology to address ordinary differential equations with Grünwald-Letnikov derivatives and highlighted the enormous potential of fractional derivatives regarding the behaviors of the solution curves, providing support for the development of more accurate models in various areas where the limitations of classical models can be overcome by using fractional derivatives.

Author contributions

J.W. Scotton participated in conceptualization, data curation, formal analysis, investigation, methodology, visualization, validation and writing – original draft, review and editing.

Conflicts of interest

The author certify that there is not a commercial or associative interest that represents conflict of interest in relation to the manuscript.

Acknowledgments

I would like to thank the State University of Santa Catarina for providing the infrastructure necessary to conduct this research.

References

R. Banerjee & R. K. Biswas (2022). Fractional optimal control of compartmental SIR model of COVID-19: Showing the impact of effective vaccination. 55(1), 616--622. https://doi.org/10.1016/j.ifacol.2022.04.101

W. E. Boyce & R. C. DiPrima (2009). Elementary Differential Equations and Boundary Value Problems. John Wiley \& Sons.

Giulio Cottone & Mario Di Paola (2011). Fractional spectral moments for digital simulation of multivariate wind velocity fields. 741--747. https://doi.org/10.1016/j.jweia.2011.03.006

A. G. O. Goulart, M. J. Lazo, J. M. S. Suarez & D. M. Moreira (2017). Fractional derivative models for atmospheric dispersion of pollutants. 9--19. https://doi.org/10.1016/j.physa.2017.02.022

Mark M. Meerschaert & Charles Tadjeran (2004). Finite difference approximations for fractional advection-dispersion flow equations. 65--77. https://doi.org/10.1016/j.cam.2004.01.033

Santanu Saha Ray (2009). Analytical solution for the space fractional diffusion equation by two-step adomian decomposition method. 1295--1306. https://doi.org/10.1016/j.cnsns.2008.01.010

A. Simmons, Q. Yang & T. Moroney (2017). A finite volume method for two-sided fractional diffusion equations on non-uniform meshes. 747--759. https://doi.org/10.1016/j.jcp.2017.01.061

Davidson Moreira, Paulo Xavier, Anderson Palmeira & Erick Nascimento (2019). New approach to solving the atmospheric pollutant dispersion equation using fractional derivatives. 118667--118675. https://doi.org/10.1016/j.ijheatmasstransfer.2019.118667

K. S. Nisar, M. Farman, M. Abdel-Aty & J. Cao (2023). A review on epidemic models in sight of fractional calculus. 81--113. https://doi.org/10.1016/j.aej.2023.05.071

M. Qayyum, S. Afzal & E. Ahmad (2023). Fractional Modeling of Non-Newtonian Casson Fluid between Two Parallel Plates. 1--12. https://doi.org/10.1155/2023/5517617

C. Y. Kee, C. Chua, M. Zubair & L. K. Ang (2022). Fractional modeling of urban growth with memory effects. 1--14. https://doi.org/10.1063/5.0085933

J. W. Scotton, J. M. S. Suarez & A.G.O. Goulart (2020). Numerical study of transient one-dimensional diffusion employing the Fractional Calculus. 95--113. https://doi.org/10.5335/rbca.v12i1.10067

M. Luft, K. Krzysztoszek, D. Pietruszczak & A. Nowocień (2021). Analysis of Dynamic Characteristics of Selected Pneumatic Systems with Fractional Calculus. {Simulation. 15(4), 835--843. https://doi.org/10.12716/1001.15.04.16

I. Podlubny (1999). Fractional Differential Equations. Academic Press.

S. G. Samko, A. A. Kilbas & O. I. Marichev (1993). Fractional Integrals and Derivatives: Theory and Applications. Gordon and Breach Science Publishers.

P. Zhang, P. Li, G. Xiu & A. R. Rodrigues (2021). Modeling Riemann–Liouville fractional differential equations for diffusion and reaction in fractal porous media. 459--475. https://doi.org/10.1007/s10910-020-01209-z