CC BY
2APPROXIMATE SOLUTION OF BOUNDARY VALUE PROBLEMS FOR A FRACTIONAL ORDER PARABOLIC-HYPERBOLIC EQUATION WITH THE CAPUTO OPERATOR
1Khairullaev U.B., 2S.B. Temirova
Novosibirsk State University, Novosibirsk city (Russian Federation) National University of Uzbekistan 2National University of Uzbekistan, Tashkent https://doi.org/10.5281/zenodo.13636954
Abstract. In the given article an initial boundary value problem posed for the fractional order diffusion equation, and the posedproblem is approximately reduced to a system of algebraic equations using the finite difference method. Using the Python software package there was found an approximate solution to the problem.
Keyworks: fractional differential operator, finite difference method, system of linear equations, approximate solution, nodes.
Introduction. In applications of partial differential equations, it is often necessary to find an approximate solution to specific problems of mathematical physics. In this paper, methods for constructing approximate solutions of a fractional differential equation, which are finite difference methods or grid methods, are described. It should be noted that the updated power of technologies and software packages makes it possible to determine more approximate solutions to the problems under study.
In the article, the initial boundary value problem for the equation is investigated using the finite difference method cD'^tu( x, t) — uxx (x, t) = 0 (1) in a finite zone Q,
where Q = {( x, t): (x — r )2 +12 < r2, t > 0|, r, a, b,a = const,0 <a< 1, f (x, t) - given
function, qD^u - fractional Caputo differential operator, which has the following form:
1 t
C Da0tu = —-- f (t — z)—a uz (x, z)dz, (2)
r(1 — a)0
r(a) = f e~tata—1dt.
0
Statement of the problem. Problem A.
It is required to find an approximate value of the solution of equation (1), at fixed points of the zone, satisfying the conditions: u(x, t) |(x_r)2+i2_r2 = ( (x, t), 0 < t; (3)
u(x,0) = t(x), 0 < x < 2r; (4) where (p() eC[0,f ] (i = 1,2), r(x) eC[0,/] .
Research method.
Assuming Cartesian orthogonal coordinates, we cover the plane of variables with a square network x = m • h, x = n • h, m,n = 0,1,2,... (5)
where h - a given positive number. The vertices of each square of this network are called nodes and the h number is called a step.
Below, using the software package, we demonstrate the coverage of the plane of variables by a x, t square network (5).
## Importing the required libraries ##
# Import the NumPy library and alias it 'np' for brevity and readability import numpy as np
# Import the math library for basic math functions and constants. import math
# We import the Matplotlib library to create graphs and charts, assigning it the alias 'plt' for convenience.
10| import matplotlib.pyplot as plt 11|
12| ## Importing Sympy modules for symbolic math operations ##
13||
14| # Import all functions and symbols from the Sympy library into the current namespace for full access.
15| from sympy import * 16|
17| # We import the 'parse_expr' function from the Sympy parsing module to parse mathematical expressions from strings.
18| from sympy.parsing.sympy_parser import parse_expr 19|
20| # Import the 'mathml' function from the Sympy print module to convert Sympy objects into MathML strings suitable for display in LaTEX. 21| from sympy.printing.mathml import mathml
22|
23| # Let us import the function 'display' from the library IPython for the reflection of objects in notes IPython.
24| from IPython.display import display Let m = 3, n = 2 u h = 0.5, then
1 | # Let's determine the number of nodes along the x and y axes
2 | N = 2
3 | M = 4
4 | # Sampling step
5 | h = 0.5
6 | # order of fractional derivative
7 | \alpha = 0.5
8| # Creating arrays of node coordinates 9| x_ = np.arange(0, 5, h) 10| y_ = np.arange(0, 5.5, h)
Now, we need to determine the nodal points at which approximate equations corresponding to equation (1) are written.
1 | # To do this, we create a set of points that bound the zone in which the solutions of the 2 | # equations lie..
3
4
5
6
7
8
9
10 11 12
13
14
15
16 17
def get_circle_coordinates(center_x, center_y, radius):
# Step between the points step = 9
# Coordinates of points
x_p = [] y_p = []
for angle in range(0, 189, step):
x = center_x + radius * math.cos(math.radians(angle)) y = center_y + radius * math.sin(math.radians(angle)) x_pappend(round(x,2)) y_pappend(round(y,2)) return x_p,y_p
x_l,f_=get_circle_coordinates(2,0,2)
Then, we define the internal and boundary nodal points and mark them with different lights.
1 | # Let's check if a given point (a, b) is within the rectangular area defined by the grid.
2 | def check(x,y,h):
3 | if (x - 2)**2 + (y )**2 > (2-h)**2:
4 | return False
5 | else :
6 | return True
7 | # Let us determine the closest point of the boundary to a given point (a, b) based on the Euclidean distance.
9 | def change_n(a,b): 10| r=[]
11| for i in range(len(f_)):
12| r.append(np.sqrt((x_l[i]-a)**2+(f_[i]-b)**2)) 13| return x_l[np.argmin(r)],f_[np.argmin(r)]
The purpose of this part of the code is to visualize the nodes and edge points of the solution of the equation in a rectangular zone. 1 2
3
4
5
6
7
8
import matplotlib.pyplot as plt
fig, ax = plt.subplots() ax.plot(x_l,f_,'--',color="y")
ax.grid(which='both', color='grey', linewidth=1, linestyle='-', alpha=0.2)
a_l=[]
9 | b_l=[] 10|
11| f_c=[]
12| fig.set_figwidth(15)
13| fig.set_figheight(15)
14| n=0
15| m=0
16| for a in x_:
17| for b in y_:
18| al=[]
19| bl=[]
20| if check(a,b,0) :
21| if check(a,b,0) and check(a,b+h,0) and check(a+h,b+h,0) and check(a+2*h,b+h,0):
22| n=n+1
23| ax.plot(a,b,MoM,color=MgM)
24| ax.plot(a,b,M.M,color=Mtab:orangeM)
25|
26| a_l.append([round(a,1),round(a,1),round(a,1),round(a+h,1),round(a+2*h,1)])
27|
28| b_l.append([0.0,round(b+h,1),round(b,1),round(b,1),round(b,1)])
29|
30| else :
31| ax.plot(a,b, M.M,color=Mtab:orangeM)
32|
33| if not check_in(a,b,h) :
34| a1,b1=change_n(a,b)
35| if a1==0:
36| f_c.append([round(a,1),round(b,1),round(b**2,2)])
37| if b1==0:
38| f_c.append([round(a,1),round(b,1),round(a,1)])
39| if a1>1:
40| f_c.append([round(a,1),round(b,1),round(b+2,1)])
41| ax.plot(a1,b1,M.M,color=McM)
42|
43| m=m+1
44| else :
45| ax.plot(a,b,M.M,color=MrM)
Based on the definition of partial derivatives, at each node, provided that all four points ( x, t), ( x + h, t) , ( x + 2h, t), ( x, t + h), , , belong to the Q domain of equation (1), we can assume that
1 f (t - z)~au2(x,z)dz «-1-Wx,z + h) - u(x.z) dz,
r(1 -a) f( ) - ( , ) h .r(1 -a) f (t - z)a •
u( x + h, t) - u( x, t)
ux (x, t)
Uxx (^t):
h
u(x + 2h,t) - 2u(x + h,t) + u(x, t)
h2
i.e. we have the right to approximately replace equation (1) at a = 1 u b = 0 and in each node indicated above with a linear algebraic equation in the form
h \u(x, z + h) - u(x, z)
ru(x,z + h) - u(x,z) 7 . x „ . , x .
J - a dz ~ u(x + 2h, t) - 2u(x + h, t) + u(x, t) (6)
r(1 - a) 0 (t - z)
When (x, t) runs through the nodes belonging to the zone , taking into account the monotonicity (in certain intervals) and the class of the fundamental solution of equation (1), as equation (6) we will have a system of algebraic equations with respect to the values at u( x, t) with the specified nodes. Some of these quantities are either directly determined independently of equation (6), based on the initial and boundary conditions, or the latter generate additional equations to (6), which in place of the system constitute an approximate grid replacement of the entire original problem.
1 | # We construct a system of equations for points belonging to the solution domain of the equation
2 | def integ(a,b,n):
3 | x = Symbol('x')
4 | z = Symbol('z')
5 | t = Symbol('t')
6 | h = Symbol('h')
7 | alpha = Symbol('a')
8 | f = Function('f)
9 | U = Function('u') 10| G = Function(T')
11| eq =Eq(((h**2)*((n)*h)**(1-alpha))/(2*G(2-alpha))*(U(a[1],b[1])-U(a[0],b[0])),
U(a[4],b[2])-2*U(a[3],b [2])+U(a[2],b[2])) 12 | return eq
This section of code is designed to automatically generate the equations needed to solve a boundary value problem in a given domain.
1 2
3
4
5
6
7
8
9
10 11 12
13
14
15
16
17
18
19
20 21
u=0
nc=0
mc=0
eq_l=[]
alpha = 0.5
for a,b in zip(a_l,b_l):
equation = str(integ(a,b,nc+1))
print(a[0],b[0],nc,mc)
for fc in f_c:
for i,j in zip(a,b):
if fc[0]==i and fc[1]==j: u=fc[2]
equation=equation.replace(f"u({i}, {j})",f'{u}") print(equation)
break else : continue
equation=equation.replace(fT({2-alpha})",f'{math.gamma(2-alpha)}")
eq_l. append(equation)
nc=nc+1
1 | # This code is intended to display the equations that were obtained in the process of solving a boundary value problem.
2 | for e in eq_l:
3 | display(parse_expr(e).subs('h',h).subs('a',alpha))
The obtained result, example.
Using the Python programming language, the following results were obtained: The plane of variables is covered by a x, t square network x = m • h, x = n • h,
m, n = 0,1,2,... ;
The number of nodal points included in the final zone is determined Q ; The number of approximate linear equations containing internal and boundary nodal points of the finite domain is determined Q ;
d) The values of the desired function are replaced by the values of the given functions at the boundary nodal points and at nodal points closer to the boundaries, controlling the number of unknowns in the system of equations.
e) The only solution to the resulting system of algebraic equations is found, which indicates an approximate solution to problem A at each internal node.
For demonstrating the obtained result, we find an approximate solution to problem A: Approximate system of algebraic equations: 0.061836053462222 - 0.0997355701003581 u{\left(0.5,0.0 \right)} = u{\left(0.5,0.0 \right)} - 2 u{\left(1.0,0.0 \right)} + u{\left(1.5,0.0 \right)}
0.166435927146588 - 0.141047395886939 u{\left(0.5,0.0 \right)} = - 2 u{\left(1.0,0.5 \right)} + u{\left(1.5,0.5 \right)} + 0.62
- 0.172747074735668 u{\left(1.0,0.0 \right)} + 0.172747074735668 u{\left(1.0,0.5 \right)} = u{\left(1.0,0.0 \right)} - 2 u{\left(1.5,0.0 \right)} + u{\left(2.0,0.0 \right)}
- 0.199471140200716 u{\left(1.0,0.0 \right)} + 0.199471140200716 u{\left(1.0,1.0 \right)} = u{\left(1.0,0.5 \right)} - 2 u{\left(1.5,0.5 \right)} + u{\left(2.0,0.5 \right)} 0.361285133520937 - 0.223015514519097 u{\left(1.0,0.0 \right)} = u{\left(1.0,1.0 \right)} - 2 u{\left(1.5,1.0 \right)} + u{\left(2.0,1.0 \right)}
- 0.24430125595146 u{\left(1.5,0.0 \right)} + 0.24430125595146 u{\left(1.5,0.5 \right)} = u{\left(1.5,0.0 \right)} - 2 u{\left(2.0,0.0 \right)} + u{\left(2.5,0.0 \right)}
- 0.263875515352797 u{\left(1.5,0.0 \right)} + 0.263875515352797 u{\left(1.5,1.0 \right)} = u{\left(1.5,0.5 \right)} - 2 u{\left(2.0,0.5 \right)} + u{\left(2.5,0.5 \right)} 0.535980104370368 - 0.282094791773878 u{\left(1.5,0.0 \right)} = u{\left(1.5,1.0 \right)} - 2 u{\left(2.0,1.0 \right)} + u{\left(2.5,1.0 \right)}
- 0.299206710301074 u{\left(2.0,0.0 \right)} + 0.299206710301074 u{\left(2.0,0.5 \right)} = u{\left(2.0,0.0 \right)} - 2 u{\left(2.5,0.0 \right)} + u{\left(3.0,0.0 \right)}
- 0.31539156525252 u{\left(2.0,0.0 \right)} + 0.31539156525252 u{\left(2.0,1.0 \right)} = u{\left(2.0,0.5 \right)} - 2 u{\left(2.5,0.5 \right)} + u{\left(3.0,0.5 \right)}
- 0.330785464275076 u{\left(2.0,0.0 \right)} + 0.330785464275076 u{\left(2.0,1.5 \right)} = u{\left(2.0,1.0 \right)} - 2 u{\left(2.5,1.0 \right)} + 1.4
- 0.345494149471335 u{\left(2.5,0.0 \right)} + 0.345494149471335 u{\left(2.5,0.5 \right)} = u{\left(2.5,0.0 \right)} - 2 u{\left(3.0,0.0 \right)} + u{\left(3.5,0.0 \right)}
- 0.359601711984474 u{\left(2.5,0.0 \right)} + 0.359601711984474 u{\left(2.5,1.0 \right)} = u{\left(2.5,0.5 \right)} - 2 u{\left(3.0,0.5 \right)} + 0.62
Solution of the system: u(0.5, 0.0)= 3.14023214625562, u(1.0, 0.0)=2.72218657189261, u(1.0, 0.5)=1.65037752844739, u(1.0, 1.0)=-1.42107398023909, u(1.5, 0.0)=2.05278420763754, u(1.5, 0.5)=2.40426941733157, u(1.5, 1.0)= -0.169146306794998, u(2.0, 0.0)= 1.19822996645208, u(2.0, 0.5)= 2.33170039973338, u(2.0, 1.0)= 0.836976661222425, u(2.0, 1.5)= -2.92233603775452, u(2.5, 0.0)=0.429544003443219, u(2.5, 0.5)=1.67281832256120
Conclusion. It should be summarized that in practice, exact or approximate values of unknowns are often used or required to be found. This is due to the intensive development of supercomputers and software packages. These developments and capabilities insist on the use of approximate methods in the study of more modern boundary value problems. In this work, the possibilities of using an approximate method (i.e., the grid method) and the Python software package were shown, when solving boundary value problems for a fractional parabolic equation (i.e., for the diffusion equation).
REFERENCES
1. R. I. Bagley. A theoretical basis for the application of fractional calculus to viscoelasticity. \\ Journal of Rheology, 27.3. Pp. 201—210, 1983.
2. R. Magin. Fractional calculus in bioengineering. \\Crit. Rev. Biom. Eng.321. P.104, 2004.
3. M. Ortigueira. Special issue on fractional signal processing and applications.\\ Signal Processing., 83.11,Pp.2285—2480, 2003.
4. K. B. Oldham.: Fractional differential equations in electrochemistry. \\ Advances in Engineering Software. 2009.
5. A. A. Kilbas, H. M. Srivastava, J. J. Trujillo.Theory and Applications of Fractional Differential Equations.\\ North-Holland Mathematics Studies. P.204. Elsevier Science B. V., Amsterdam. 2006.
6. K. S. Miller, B. Ross.An Introduction to the Fractional Calculus and Differential Equations.\\ John Wiley., New York.1993.
7. I. Podlubny. Fractional Differential Equations.\\ Academic Press. New York. 1999.
8. S. G. Samko, A. A. Kilbas, O. I. Marichev. Fractional Integral and Derivatives, Theory and Applications.\\Gordon and Breach, Longhorne. PA. 1993.
