Numerical analysis of fractional order integral dynamical models with piecewise continuous kernels

Volterra integral equations (VIEs) [1] were introduced by Vito Volterra and then studied by Traian Lalescu in his 1908 thesis. In weakly regular case, VIEs were introduced in [2], and the theory of such equations is generalized to the case of systems of equations and to abstract operator equations in the monograph [3]. Integral equations are essential tools in various fields including power and electrical engineering, mathematical physics, control theory and mechanics.

Fractional integration and fractional differentiation are generalisations of notions of integer-order integration and differentiation, and include n-th derivatives and n-folded integrals (n denotes an integer number) as particular cases. The first mention of derivatives of non-integer order is presented in the correspondence between G. Leibniz and J. Bernoulli. In a letter written by G. Leibniz to G. Lopital, there is also an interesting mention about the paradox and possible useful practical application of differentials of order . In recent years, interest in fractional calculus grows [4–6].

Fractional derivatives have numerous applications; they are in the core of various mathematical models for viscoelastic bodies, in chemical physics, the theory of gravity, viscoplasticity, etc. [6, 7]. Fractional derivatives have practical applications in modelling the behavior of viscoplastic materials. In particular, the equations of state in the theory of viscoplasticity contain fractional derivatives and can be reduced to weakly singular Volterra integral equations of the second kind [8, 9].

This article discusses weakly regular Volterra equations of the first kind

t

I K(t, s)x(s)ds = f (t), 0 <  s <  t <  T, f (0) = 0,                  (1)

where the kernel is defined as follows:

{ K i (t, s), t, s G m i ,     m i = { t, s | a i - i (t) < s < a i (t) } ,

... ......... (2)

Kn(t, s), t,s G mn,    ao(t) = 0, an(t) = t, i = 1, n, ai(t), f (t) G С^тj, Ki(t, s) have continuous derivatives with respect to t for G mi, Kn(t,t) = 0, ai(0) = 0, 0 < ai(t) < a2(t) < ... < an-i(t) < t, ai(t),..., an-i(t) grow in a small neighborhood 0 < t < t, 0 < a‘(0) < ... < a^-i (0) < 1.

The works [2, 10] are the first to study such equations and to outline the non-uniqness of the solution. The paper [11] derives the existence of a continuous solution depending on free parameters as well as sufficient conditions for the existence of a unique continuous solution to the system of VIE of the first kind with discontinuous kernels. The class of Volterra operator equations of the first kind with piecewise continuous kernels is studied in [12]. Various numerical methods for linear and nonlinear VIE with piecewise continuous kernels and their systems with applications for power systems operation are proposed in [13–15]. In [16], the Volterra model is employed for a load leveling problem in modelling the hybrid AC/DC power systems with renewable energy sources and storage system.

For systematic studies of VIE of the first kind with piecewise continuous kernels, see the book [3] and the part 1 in [17]. For numerical solution of Volterra integral equation of Abel type, see [18]. The numerical solution of the Volterra fractional integral equations of the second kind using the Simpson 3/8 rule method is proposed in [4].

In the case of a fractional order of integration for these equations, consider the leftsided Riemann-Liouville fractional integral of the order в >  0

t

I ^е f (t) = Г(в) J (Т-Тв^ds-                       ⁽³⁾

For equation (1) in the case of fractional order of integration (3), we solve the integral equation

t

^{( t}, s e ^x(s)ds = ^{f (t)}, ⁰ <  ^s <  t <  T, ^{f (0)} = ⁰, ^{Г (e}) J (t - si

with piecewise continuous kernel (2).

1.    Polynomial Collocation Method

Let us write equation (4) in the following expanded form:

_n ^a i (t)

1 ^ Г K ( i,s i x e s ) ds = f(t), t G [0,T], 0 <в< 1.           (5)

Г(в) “Г J (t - s)e i_i ai-1 (t)

Apply the following representation (as power series segment) of the solution to equation (5):

with the desired a j , j = 0,1,..., m, .

a i ( t )

Г(в) 5 j i=1 ai—i (t)

.

m

Xm (t) = ^a3 tj j=0

. Let us substitute (6) into (5), then we get

K i (t, s) (t - s) ^e

m

52 a j j ds = f^(t)

Let us change the order of integration and summation

Гв)

a ( t )

nm aj i 1 j 0    ai-i(t)

K iMj ^—^ ds = ^{f (t)}

and transform the last equality as follows:

Гв)

ai (t)

mn aj j 0 i 1 ai-i(t)

K i ^(t, s)s ds = (t — s) e ^ds    f^(t).

We require the fulfillment of the last equality at the points of a uniform grid kT tk =---—, k = 1, 2,...,m + 1.

m + 1

Therefore, we arrive at the following system of equations for the unknown coefficients aj, j = 0, i,.

.

. , m:

m

52 Y kj a j = F k , k = 1, 2,...,m + 1, j 0

where

i

^{Y kj} = г(в)

n      ^a i ^{( t} k )

5 /

a i- 1 (t k )

K i (t k , s)s j

^(t k ^— s)

ds,   F k = f (t k ).

System of linear algebraic equations (8) is solved for unknown expansion coefficients (6) by the Jordan–Gauss method, and in order to calculate the elements of the matrix of this system, below we construct special quadrature formulas that take into account the power singularity of the integrands.

• 1.1.    Approximation of Weakly Singular Integrals

In the previous section, when constructing the matrix of system of linear algebraic equations, in order to determine the coefficients of the approximating polynomial, it is necessary to calculate the values of the functions ai (tk )                  j

I ^(i,k,j )= J

,i^ k, ds, i = 1,n, k = 1,m + 1, j = 0,m, (10) (tk — s)e ai-1(tk)

containing either a convergent improper integral (for i = n), or significant in magnitude integrands (for i < n). In this section, we propose a method to calculate integrals (10).

For each triple i, j, k , an internal mesh of nodes is constructed as follows:

T = a i - i (t k ) + r (a i (t k ) - a i — i (t k )),

l = 0,1,..., r,

a i (t k ) - a i — i (t k ) h

+ 1, (11)

were [t] stands for integer part of a real number t, and h is a sufficiently small parameter, which we conventionally call the grid step. Then we approximate integral (10) as follows:

τ l

I (i,k,j ) = g/

¹ T l- 1

~        ^ Ki (tk,    +  Tl ) (+  Tl )  ((tk — Т1-1)1-в — (tk — T )1-в) .

Therefore, for relatively small values T , we propose a sufficiently effective method for solving equations of the form (5), the error of which can be estimated by the inequality тm

Bx(t) - Xm(t)Bc[o,T]< + ,,, , L = max        (t)|

(m ~+ 1 )•                t t [o,T ]

for x(t) G C ^{m +1} [0, T].

With an increase in the length of the planning interval [0, T ], to maintain the order of accuracy, it is necessary to increase the degrees of the approximating polynomials, which leads to significant computational errors. In order to overcome this restriction, the method can be generalized using a polynomial spline approximation of the solution constructed in each section in a similar way.

It is known that exact solutions to weakly singular Volterra integral equations can have unbounded derivatives as t ^ 0 [19]. In this case, the exact solutions belong to the class of functions C ^r,e (0,T] defined as follows.

Definition 11. Let fi = (0,T], 0 С в < 1. We say that a function x(t) belongs to C ^r,e (0,T ], if for t G (0, T ] the function x(t) has continuous derivatives up to the order r estimated as follows:

A

| x ^{( k ) (t) | С} _+k i+,3 ^for all t ^{G (0,T],k} = ^0, 1,...^,г. t k 1+p

In this case, an approximate solution to equation (5) can be found as

m

X m (t) = a o + ^ a j t ^{j - e} . j =1

The method to determine coefficients expansion (14) is employed below in the similar way.

				Table 1
	Numerical results for Example 1
		Errors
m	h = 10 - 2	h = 10 - 3	h = 10 - 4	h = 10 - 5
3	8,3730 • 10 - 4	1,0826 • 10 - 3	1,1168 • 10 - 3	1,0538 10 - 3
4	1,1624 • 10 - 3	2,0991 • 10 - 4	1,7066 • 10 - 4	1,6319 10 - 4
5	2,0042 • 10 - 3	5,3862 • 10 - 5	6,9064 • 10 - 6	5,3063 10 - 6
6	1,8925 • 10 - 3	5,3393 • 10 - 5	1,1764 • 10 - 6	5,0464 10 - 7
7	1,5110 • 10 - 3	4,6339 • 10 - 5		4,3806 10 - 8
8	1,8389 • 10 - 3
9	1,8421 • 10 - 3
10	7,6494 • 10 - 4

2.    Numerical Results

In order to demonstrate how the numerical method works, we show two examples.

Example 1. As the first example, we solve the following equation:

t/2                       2t/3

(s + 2) ³ x(s) (t - s) ^e

= f (t),

Г(в)

Г e ^{( t + s )} x(s) , f tsx(s) . Г

• •        +(t-^

0                        t/22

where the right hand side of the equation is chosen such that the exact solution is x(t) = sin(t). Calculations are performed for t E [0,T], the order of the fractional integration is в = 1 . Table 1 shows the errors e = max | x(t i ) — X m (t i ) | for different values of m and h.

• ² 0 ≤ i ≤ N ^m

Example 2. As the second example, we solve the equation sin 2t                              sin 23t

t

/ sin ² ₃ ^t

² + s + 3)ф)    1

(t — s) e         I     f^W’

1 f Vt + s x(s) д ,   [ (t + s)x(s)

Г(в) ¹ 1 "¥—^^{ds +} 1 t^{ds +}

0                              sin 2t where, as well as in the previous example, the right hand side of the equation is chosen such that the exact solution is x(t) = ln^+lt , t E [0, T] and в = 2•

Table 2 shows the errors e = max | x(t i ) — X m (t i ) | for different values of m and h.

0 ≤ i ≤ N

The examples show the following interesting feature: as the step decreases from h = 10 - 4 to h = 10 ^{- 5} , the error is significantly decreased at Example 1, and remains at the same level at Example 2. Note that the error in the approximate solution to the second model problem is somewhat worse than for the first one. This is due to less smooth components K i (t, s) and a more complex structure of the exact solution. In the future, we intend to develop more accurate quadratures for approximation integrals arising when calculating the coefficients of the system.

Table 2

Errors

m h = 10 - 2      h = 10 - 3      h = 10 - 4

h = 10 - 5

3 6,2792 · 10 - 3 6,3200 · 10 - 3 6,3211 · 10 - 3 4 1,6741 · 10 - 3 1,7364 · 10 - 3 1,7370 · 10 - 3 5 4,9248 · 10 - 4 5,3815 · 10 - 4 5,3845 · 10 - 4 6 3,1440 · 10 - 4 1,9342 · 10 - 4 1,9291 · 10 - 4 7 3,1069 · 10 - 4 7,9240 · 10 - 5 8,0305 · 10 - 5 8 3,8965 · 10 - 5 3,8303 · 10 - 5 9 2,0569 · 10 - 5 2,0370 · 10 - 5 10 1,0447 · 10 - 5 1,1771 · 10 - 5 11 7,4383 · 10 - 6 12 4,5928 · 10 - 6 13 3,2945 · 10 - 6 14 2,7023 · 10 - 6 15 1,2863 · 10 - 6

6,3212 10 - 3 1,7371 10 - 3 5,3846 10 - 4 1,9292 10 - 4 8,0298 10 - 5 3,8296 10 - 5 2,0406 10 - 5 1,1840 10 - 5 7,3277 10 - 6 4,7671 10 - 6 3,2260 10 - 6 2,2571 10 - 6 1,6224 10 - 6

Numerical results for Example 2

Conclusion

In this article, we present the numerical collocation method that can be used to solve linear Volterra fractional order integral equations of the first kind with piecewise continuous kernel. The proposed quadrature formula effectively employ the restriction of the considered weakly regular equations with fractional order integrals and piecewise continuous kernels.

Acknowledgment. Authors owe a special debt of gratitude to Professor Nikolai A. Sidorov for his generous, thoughtful support and for sharing his vast know ledge.

This work was partially supported by the base part of the Government Assignment for Scientific Research from the Ministry of Science and Higher Education of Russia (project code: FZZS-2020-0039).