Numerical method for solving the subdiffusion differential equation with nonlocal boundary conditions

1.   Introduction

Fractional differential equations are widely used in various fields of science and engineering ^[1,2,3]. Among the popular fields, we can mention physics models, such as anomalous diffusion ^[4], viscoelastic media ^[5], propagation of acoustic waves in porous media ^[6], etc. Fractional calculus is used in control problems, to model dynamic systems in which the current state depends on the whole history ^[7], in signal processing ^[8,9], in biology for the modeling of neurons ^[10], in epidemiology ^[11], in environmental sciences ^[12], and in economic models with memory ^[13].

Differential equations with nonlocal boundary conditions arise in problems of mathematical modeling of various processes such as heat transfer, chemical diffusion, hydrology, and biochemistry.

Boundary conditions of Samarskii-Ionkin type were originally considered in works ^[14] for the heat transfer problem in a tenuous plasma:

The last condition means that the total energy of the system is constant. This problem was further investigated in other works ^[15,16]. The nonlocal condition was shown to be equivalent to

which means equal heat flows at the ends of the interval. The existence of the solution was proved.

Later, the heat equation with two-point boundary conditions of the general type was considered in ^[17]:

Restricting the class of such conditions to strongly regular boundary conditions allows one to use the well-established theory for self-adjoined operators and well-known methods ^[18,19,20]. In the case of not strongly regular boundary conditions, the system of eigenfunctions of the spectral problem do not form a Riesz basis. Thus, these methods cannot be used. In this case, new methods or modifications of existing methods must be developed.

In ^[21], the heat problem with nonlocal boundary condition

was considered. It was shown that these conditions are not strongly regular for $a\ne0$. The existence, uniqueness, and stability of the solution was established.

Work ^[22] was devoted to solving an initial boundary value problem for a heat equation:

with regular but not strongly regular boundary conditions of the general type (1.1). It was shown that for the case of the even potential function $q(x) = q(1-x)$, the considered class of problems can always be reduced to the sequential solution of two similar problems with strongly regular boundary conditions. The proof does not depend on whether the system of eigenfunctions of a spectral problem forms a basis.

Note that the problem of establishing the basis properties of a system of eigenfunctions remains open so far. For special cases of not strongly regular boundary conditions, the basis property was confirmed in ^[23].

Problems for fractional differential equations with local and nonlocal boundary conditions have been the topic of many works ^{[24,25,26,27,28]}. There are many methods for numerical solutions of such problems, such as the finite difference methods ^[29], spectral methods ^[30,31], etc. For various types of fractional derivatives, stability, convergence, and solvability of difference schemes were researched in ^[32].

One of the methods for solving the initial boundary value problem for a difference equation with nonlocal boundary conditions was developed in the works ^[22,33] and applied to the nonlocal heat equation with not strongly regular boundary conditions. It consists in reducing the nonlocal difference problem to a sequential solution of two local difference problems.

In ^[34], a parallel numerical algorithm was constructed for solving the initial boundary value problem for the subdiffusion equation with the homogeneous Dirichlet condition. The stability of the difference scheme was established. Numerical experiments were performed to study the performance of parallel implementations. These results were used in ^[35] to develop a parallel numerical algorithm for solving the inverse problem of identifying the space-dependent source term in the two-dimensional fractional diffusion equation. To solve the inverse problem, a regularized iterative conjugate gradient method is used. At each iteration of the method, a finite difference scheme is used to solve the auxiliary direct initial boundary value problem.

Research on correctness and stability of direct and inverse problems of mathematical physics was performed in numerous works ^{[36,37,38,39,40,41]}.

In this work, we construct a numerical method for solving the initial boundary value problem for the subdiffusion equation with nonlocal boundary conditions. We solve the problem by reducing it to two problems with local boundary conditions. The solution of the original problem is found as a sum of the solutions of the subproblems. These subproblems are solved by using implicit finite difference schemes. The schemes produce systems of linear algebraic equations (SLAE) that must be solved for each successive time step. To solve these SLAEs, we utilize the Thomas algorithm. We research the stability and convergence of the constructed difference schemes and obtain the stability estimates for the proposed method. To confirm the stability and convergence order of the developed method, we perform numerical experiments.

Our approach to the numerical solution of the nonlocal initial boundary value problem is used for the first time to solve the problem for a time-fractional equation. For simplicity, we consider the problem for the case of a zero initial value. The proposed numerical method and algorithm presented below may be applied for more generalized statements of the problem.

The rest of the article is structured as follows: Section 2 describes the statement of the initial boundary value problem for the subdiffusion differential equation with nonlocal boundary conditions. Section 3 is devoted to the numerical algorithm for solving the problem. It describes the construction of the finite difference scheme for solving the subproblems. It also contains the results on researching the stability and convergence of the constructed finite difference scheme. Section 4 describes the performed numerical experiments and discusses their results. Section 5 concludes the article.

2.   Problem statement

In this work, we consider numerical algorithms for solving the initial boundary value problem for the subdiffusion equation with nonlocal boundary conditions. Consider the equation

with initial conditions

and the following nonlocal boundary condition ^[21]

where $a > 0$.

Here, we consider the Caputo fractional derivative with order $\alpha$ in the form ^[42]

with $m = \lceil \alpha \rceil: \alpha \in (m-1, m)$, $m\in N, x > 0$.

By attempting to solve the initial boundary value problem (2.1)–(2.3) by using the Fourier method, we obtain the spectral problem for operator $l$, which is defined by the differential expression and boundary conditions:

Boundary conditions (2.5) are regular but not strongly regular. Therefore, the system of eigenfunctions of operator $l$ is a complete system, but it does not form a basis in $L_2(0, 1)$. The spectral problem (2.4) and (2.5) has two series of eigenvalues:

where $\beta_k$ are roots of the equation

They satisfy the inequalities $\pi k < \beta_k < \pi k + \pi/2, \ k = 0, 1, 2...$. For $\delta_k = \beta_k - \pi k$ with sufficiently large $k$, the following inequalities hold:

The eigenvalues of the problem are

The system is almost normalized, but it does not even form a normal basis in $L_2(0, 1)$. Therefore, the direct application of the Fourier method is impossible.

3.   Numerical algorithm for solving the problem

3.1.   Reducing the nonlocal problem to two local problems

For the numerical solution of problem (2.1)–(2.3), we adapt the algorithm ^[22]. It was proposed for nonlocal problems with classical derivatives. The algorithm is based on reducing the nonlocal differential problem to the successive solving of two problems with Sturm-type boundary conditions. This method allows one to avoid studying correctness and stability of difference schemes for the original nonlocal problem, and replace it with studying local difference schemes for the auxiliary subproblems.

Let us represent the solution $u(x, t)$ as a sum of two functions $C(x, t)$ and $S(x, t)$, such as

where

$C(x, t)$ is an even part of function $u(x, t)$ and $S(x, t)$ is an odd part of $u(x, t)$ on the interval $0\le x \le 1$. Thus, these functions have the following properties:

Then, for the functions $C(x, t)$ and $S(x, t)$, we solve two initial boundary value problems.

Let us obtain the boundary conditions for the subproblems for $C(x, t)$ and $S(x, t)$. From (2.3), we have

Using the properties of $C(x, t)$ and $S(x, t)$, we transform these conditions to contain only values for $x = 0$:

Now, let us transform the conditions to contain only values for $x = 1$:

Now, for $C(x, t)$, we have a problem with homogeneous boundary conditions of the third kind and a homogeneous initial condition:

where $f^C(x, t) = \Big(f(x, t)+f(1-x, t)\Big)/2$.

For $S(x, t)$, we have the following problem with the Dirichlet boundary condition and the homogeneous initial condition:

where $f^S(x, t) = \Big(f(x, t)-f(1-x, t)\Big)/2$.

By construction, we can formulate the following statement:

Statement 1. A solution of problem (2.1)–(2.3) can always be equivalently reduced to a sequential solution of two boundary value problems (3.2) and (3.3) with local boundary conditions.

Thus, to solve problem (2.1)–(2.3), we need first to solve problem (3.2) for $C(x, t)$. Then, knowing the function $C(x, t)$, we solve problem (3.3), because its boundary conditions depend on $C(0, t)$. Finally, we obtain the sought function $u(x, t)$ by formula (3.1).

3.2.   Discretization and difference schemes

Let us introduce regular grids for $x$ and $t$ with $M+1$ and $N+1$ points, respectively: $i = 0, ..., M,$ $h = 1/M, \ x_i = ih$, $n = 0, ..., N, \ \tau = 1/N, \ t_n = n\tau$. Denote the values of grid functions at grid points as $u_{i, n} = u(x_i, t_n)$.

To approximate the Caputo fractional derivative of function $u(x, t)$ (and similarly for functions $C(x, t)$ and $S(x, t)$) at the time level $n$, we use the L1 formula ^[32,43]:

Using formula (3.4) of order $O(\tau^{2-\alpha})$ and a central difference scheme of order $O(h^2)$ to construct an implicit difference scheme for Eq (3.2), we obtain

Omitting the small term $O\left (\tau^{2-\alpha}+h^2\right)$, we obtain the difference equations

Then, let us rearrange

The boundary condition at the point $x_0 = 0$ is

To approximate this condition with the order of approximation $O(h^2)$, let us obtain the derivative $C_x(0, t)$. Consider the Taylor series expansion

Substituting $C_{xx}(x_0, t) = D^{\alpha}_{t} C(x_0, t) - f^C(x_0, t)$ from (3.2), we obtain

Substitute this expression for derivative $C_x(x_0, t)$ and formula (3.4) into the boundary condition

Let us rearrange the terms:

Similarly, for the boundary condition at the point $x_M = 1$, we obtain a difference equation:

Combining all the difference equations for points $x_i, \ i = 0, ..., M$, we can form the system of $(M+1)$ linear algebraic equations in the matrix form

where

Similarly, we can obtain the difference equations for the function $S(x, t)$ and problem (3.3):

Due to the Dirichlet boundary conditions, we have a SLAE of smaller size, $M-1$:

where

Thus, the numerical algorithm for solving problem (2.1)–(2.3) is the following:

At each sequential time level $n = 1, 2, ..., N$:

$(1)$ Solve system (3.10) for $C_{i, n}, \quad i = 0, ..., M$.

$(2)$ Solve system (3.12) for $S_{i, n}, \quad i = 1, ..., M-1$.

$(3)$ Compute the values $u_{i, n} = C_{i, n}+S_{i, n}, \quad i = 0, ..., M.$

3.3.   Solving the SLAEs

To solve SLAEs (3.10) and (3.12) with tridiagonal matrices, the Thomas algorithm (also known as the sweep method or elimination method) ^[44] is used. It is a direct method for systems with special matrices. For the tridiagonal systems (3.10) and (3.12), we denote the elements of the main diagonal of the matrices as $q_i$, and the elements of the lower and upper diagonals as $p_i$ and $r_i$, respectively. The algorithm may be written as follows (in the case of system (3.10) for $C$):

● The forward elimination phase:

● The backward substitution phase:

The conditions for correctness of this algorithm are the diagonal dominance in the matrix

This property obviously holds for matrices of systems (3.10) and (3.12).

3.4.   Stability of the difference schemes

To justify the proposed algorithm, let us derive estimates of the stability of schemes (3.10) and (3.12) with respect to the initial data and the right-hand side.

Theorem 1. Suppose $\{C_{i, n}\ |\ 0\le i\le M, \ 0\le n \le N\}$ is the solution of the difference scheme (3.10). Then, the following estimate holds:

where $C_n = \left \{C_{0, n}, C_{1, n}, ..., C_{M, n}\right \}$.

Proof. Rewrite (3.7) as follows:

Similarly, we can rewrite the difference equations for boundary conditions (3.8) and (3.9):

Suppose $\left \| C_n \right \|_{\infty} = |C_{l, n}|$ for some inner point $0 < l < M$. Note that $\left (w_j^{\left(\alpha\right)}-w_{j+1}^{\left(\alpha\right)} \right) > 0$ and $\sigma_{\alpha, \tau} > 0$. Taking the absolute value of both sides and using the triangle inequality, we obtain

We can simplify it to

For the cases of $\left \| C_n \right \|_{\infty} = |C_{0, n}|$ or $\left \| C_n \right \|_{\infty} = |C_{M, n}|$, we can obtain the same inequality:

Next, let us use the mathematical induction method. Let

Note that since $w_{n}^{\left(\alpha\right)} < w_{n-1}^{\left(\alpha\right)}, \quad A_n > A_{n-1}$.

The base case. For $n = 1$, we have

The induction step. Suppose $\left \| C_{k} \right \|_\infty \le A_k$ for $1\le k \le n-1$. Consider the case of $k = n$.

Thus, by induction, inequality (3.16) is true for $k = 1, 2, ..., n$. □

Similarly, we can obtain an estimate for $S$.

Theorem 2. Suppose $\{S_{i, n}\ |\ 0\le i\le M, \ 0\le n \le N\}$ is the solution of the difference scheme (3.12). Then, the following estimate holds:

Proof. The proof is similar to Theorem 1. For the boundary cases, we utilize the inequality (3.16). □

Now, we can obtain an estimate for approximate solution $u_{n} = C_{n} + S_{n}, \ 0 < n \le N$, to the base problem (2.1)–(2.3).

Theorem 3. Suppose $\{C_{i, n}, \ S_{i, n}\ |\ 0\le i\le M, \ 0\le n \le N\}$ are the solutions of the difference schemes (3.10) and (3.12).

Then, the following estimate holds:

Proof. Let us utilize the triangle inequality. Since $u_{n} = C_{n} + S_{n}$, then $\left \| u_n \right \|_\infty \le \left \| C_n \right \|_\infty + \left \| S_n \right \|_\infty.$

Now, by inequalities (3.16) and (3.17), we have

□

Thus, the difference schemes (3.10) and (3.12) are unconditionally stable with respect to the right-hand side.

3.5.   Convergence of the difference schemes

In ^[32], it is shown that the approximation schemes used in (3.10) and (3.12) have a truncation error of order $O(\tau^{2-\alpha}+h^2)$, i.e., for the error $r_{i, n}$, there exists a constant $c > 0$, such that

Theorem 4. Suppose $\left \{ C_{i, n} | \ 0\le i\le M, \ 0\le n\le N \right \}$ is the exact solution of problem (3.2) and $\left \{ \overline{C_{i, n}} | \ 0\le i\le M, \ 0\le n\le N \right \}$ is the solution of difference scheme (3.10). Consider the error

Then, for some $c_1 > 0$, it holds that

Proof. Let us subtract Eq (3.6) from (3.5), and, similarly, subtract boundary difference Eqs (3.8) and (3.9) from the boundary conditions. We obtain the difference scheme for error ${e^C}_n$:

where

We can apply Theorem 1 to this scheme and obtain

By (3.19), for some constant $c_1 > 0$, we have

For the coefficients $w^{(\alpha)}_{j}$, the following property holds (see [32, p. 32]):

Thus,

□

Similarly, we can obtain an estimate for the error ${e^S}_{i, n} = S_{i, n} - \overline{S_{i, n}}, \ 0\le i\le M, \ 0\le n\le N.$

Theorem 5. Suppose $\left \{ S_{i, n} | \ 0\le i\le M, \ 0\le n\le N \right \}$ is the exact solution of problem (3.3) and $\left \{ \overline{S_{i, n}} | \ 0\le i\le M, \ 0\le n\le N \right \}$ is the solution of the difference scheme (3.12).

Then, for some $c_2 > 0$, it holds that

Now, consider the error between the solution of the original problem and the approximate solution obtained by the proposed numerical algorithm.

Theorem 6. Suppose $\left \{ u_{i, n} | \ 0\le i\le M, \ 0\le n\le N \right \}$ is the exact solution of problem (2.1)–(2.3) and $\left \{\overline{C_{i, n}}, \ \overline{S_{i, n}} | \ 0\le i\le M, \ 0\le n\le N \right \}$ are the solutions of the difference schemes (3.10) and (3.12). Let

where $u_{i, n} = C_{i, n} + S_{i, n}\quad$ and $\quad \overline{u_{i, n}} = \overline{C_{i, n}} + \overline{S_{i, n}}$.

Then, for some $c_3 > 0$, it holds that

Proof. By Theorems 4 and 5, we have

The proof is complete. □

In the next section, we perform numerical experiments to confirm the convergence order.

4.   Numerical experiments

The test problem is as follows:

where

The analytical solution is $u\left(x, t\right) = t^m\left[x^2\left(1-x\right)^2-ax^2+\left(2+a\right)x\right].$

For subproblems $C(x, t)$ and $S(x, t)$, the right-hand terms are

The analytical solutions for $S(x, t)$ and $C(x, t)$ are

4.1.   Experiment 1: convergence in time

Experiment 1 was performed for the parameter values $\alpha = 0.5, \ m = 2$, and $a = 1$. To estimate the order of convergence in time, we solved the test problem for grid sizes $N = 64,128,256,512$, and $1024$. To reduce the impact of the space step $h = 1/M$, we used the fixed value $M = 4096$. Thus, $\tau^{2-\alpha} > h^2$ for all values of $\tau = 1/N$.

Table 1 shows the absolute errors $\delta_1(\tau) = \left \| u(1, t) - \overline{u_\tau(1, t)} \right \|_\infty$ of the approximate solution $\overline{u_\tau(1, t)}$ obtained for step size $\tau$, as well as the mixed relative errors $\delta_2(\tau) = \frac{ \left \| u(1, t) - \overline{u_\tau(1, t)} \right \|_\infty} { 1+ \left \|u(1, t)\right \|_\infty}$ for various dimensions $N$ of the computational grid. It also contains the convergence orders $\frac{\log(\delta(\tau_1)/\delta(\tau_2))}{log(\tau_1/\tau_2)}$ for both errors.

4.2.   Experiment 2: convergence in space

Experiment 2 was performed for the parameter values $\alpha = 0.5, \ m = 2$, and $a = 1$. To estimate the order of convergence in space, we solved the test problem for grid sizes $M = 64,128,256,512$, and $1024$. To reduce the impact of the time step $\tau = 1/N$, we used the fixed value $N = 16{, }384$. Thus, $h^2 > \tau^{2-\alpha}$ for all values of $h = 1/M$.

Table 2 shows the absolute errors $\delta_1(h) = \left \| u(x, 1) - \overline{u_h(x, 1)} \right \|_\infty$ of the approximate solution $\overline{u_h(x, 1)}$ obtained for step size $h$, as well as the mixed relative errors $\delta_2(h) = \frac{ \left \| u(x, 1) - \overline{u_h(x, 1)} \right \|_\infty} { 1+ \left \|u(x, 1)\right \|_\infty}$ for various dimensions $M$ of the computational grid. It also contains the convergence orders $\frac{\log(\delta(h_1)/\delta(h_2))}{log(h_1/h_2)}$ for both errors.

Figure 1 shows the exact and approximate solutions $u\left(x, 1\right)$ at the final instant $t = T = 1$, as well as the auxiliary solutions $S(x, 1)$ and $C(x, 1)$ obtained for the grid sizes $M = 4096$ and $N = 1024$. Solid graphs are exact solutions and dots are approximate ones.

Figures 2(a)–(c) show the surface plots for $C(x, t)$, $S(x, t)$, and $u(x, t)$ obtained for the grid sizes $M = 4096$ and $N = 1024$. Solid graphs are the exact solutions and dots are the approximate ones. Figure 2(d) shows the plot of error $\delta(x, t) = u(x, 1) - \overline{u(x, 1)}$.

Let us discuss the results of Experiments 1 and 2.

$(1)$ In Experiments 1 and 2, the finer the grid is for $x$ or $t$, the lower the errors of the solutions. This confirms the stability and convergence of the difference schemes (3.10) and (3.12).

$(2)$ In Experiment 1, the computed convergence order is $1.5$ in time for $\alpha = 0.5$.

$(3)$ In Experiment 2, the convergence order is $2$ in space.

$(4)$ This coincides with the theoretical convergence order for the difference schemes based on the L1 formula for the Caputo time derivative and central difference scheme for space.

$(5)$ In the future, we plan to utilize higher-order approximations for the Caputo fractional derivative ^[45,46,47].

5.   Conclusions

In this work, we have constructed a numerical algorithm for solving the initial boundary value problem for the subdiffusion equation with nonlocal boundary conditions. The algorithm consists of reducing the nonlocal problem to the successive solving of two subproblems with local boundary conditions. The solution of the original problem has been found as a sum of the solutions of the subproblems. To solve the subproblems, the implicit difference schemes have been constructed. For solving systems of linear algebraic equations, we have applied the Thomas algorithm. We have established the unconditional stability and obtained the convergence estimates of the difference schemes. We have conducted numerical experiments to study the constructed algorithm. The numerical results show the stability and convergence of the algorithm. The experiments confirm that the order of convergence of the proposed method coincides with the theoretical one.

In the future, we plan to adapt our approach to subdiffusion equations with more generalized formulations (1.1) of not strongly regular boundary conditions and nonzero initial conditions. We plan to adapt our approach to two-dimensional problems. Our approach for solving the initial boundary value problems may also be utilized in numerical algorithms for solving the inverse problems for nonlocal subdiffusion equations.

Author contributions

Murat A. Sultanov: Methodology, investigation, writing–review and editing, supervision, project administration, funding acquisition; Vladimir E. Misilov: Methodology, investigation, writing–original draft preparation, writing–review and editing; Makhmud A. Sadybekov: Methodology, validation, investigation, writing–review and editing, supervision. All authors have read and agreed to the published version of the manuscript.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was financially supported by the Ministry of Science and Higher Education of the Republic of Kazakhstan (project No. AP19676663).

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.