Applications of mixed finite element method based on Bernstein polynomials in numerical solution of Stokes equations

1.   Introduction

Many natural phenomena can be modeled mathematically to obtain approximate models^[1]. Compared to the classical diffusion equation, the fractional diffusion equation may be more suitable for modelling anomalously slow transport processes with memory and inheritance. In recent years, fractional calculus has found widespread applications in many fields including turbulence, wave propagation, signal processing, porous media, and anomalous diffusion^[2,3,4].

Considering that the single-term time fractional derivative cannot adequately describe many complex physical or biological processes, recently, a multi-term time and time distributed order fractional equations have been developed. The time distribution order equation is also a generalisation of the multi-term time equation. Therefore, the study of multi-term time fractional differential equations becomes very important and useful in various applications^[5,6]. For example, multi-term fractional diffusion equation has simplified the modelling of phenomena such as diffusion processes, viscoelastic damping materials, oxygen delivery through capillaries and anomalous relaxation of magnetic resonance imaging signal magnitude^[7,8,9]. Because the numerical solution is the most important in practice, a great deal of research has been done in the study of numerical solutions of multi-term fractional diffusion equations^[10,11,12]. Qiu has analyzed numerical solutions for the Volterra integrodifferential equations with tempered multi-term kernels^[13]. Hu et al. have formulated a backward Euler difference scheme for the integro-differential equations with the multi-term kernels^[14]. Guo et al. have proposed the alternating direction implicit numerical approaches for computing the solution of multi-dimensional distributed order fractional integrodifferential problems^[15]. {Guo et al. have developed an efficient finite difference/generalized Hermite spectral method for the distributed-order time-fractional reaction-diffusion equation on one-, two-, and three-dimensional unbounded domains^[16].

In general, exact solutions of the fractional diffusion equation are rarely obtained in practical applications. Therefore, it is necessary to develop some effective numerical methods to solve the multi-term time-space fractional equations. A large number of numerical methods have been developed for two-dimensional time-space fractional order diffusion equations with a single time fractional order derivative as a special case of multi-term time-space fractional equations. Abd-Elhameed et al. have introduced a new set of orthogonal polynomials to effectively obtain numerical solutions of the nonlinear time-fractional generalized Kawahara equation by using the collocation algorithm^[17]. Moustafa et al. have created Chebyshev polynomials for the time fractional fourth-order Euler-Bernoulli pinned-pinned beam based on the Petrov-Galerkin^[18]. Peng et al. have developed a novel fourth-order compact difference scheme for the mixed-type time-fractional Burgers equation, using the $L1$ discretization formula and a nonlinear compact difference operator^[19]. Marasi and Derakhshan have proposed a hybrid numerical method based on the weighted finite difference and the quintic Hermite collocation methods for the for solving the variable-order time fractional mobile-immobile advection-dispersion model^[20].

It is well known that meshless methods are a type of point set based numerical method that considers a set of scattered and uniform data points. Due to this property, the meshless method can be applied to high-dimensional models with irregular and complex domains^[21,22,23]. However, due to the singularity of the spatial fractional operators, we only deal with problems on convex domains. A meshless method with Hermite splines of order five is used to discretize the Riesz fractional operator in the spatial direction, which gives higher accuracy with fewer points.

From the last few decades, there are many methods to solve the single time fractional diffusion equations, for instance the finite difference method, interpolation, implicit stepping methods, etc. The Laplace transform is one of the powerful tools for solving differential equations in engineering and other scientific disciplines. However, solving differential equations with the Laplace transform sometimes results in solutions in the Laplace domain that are not easily invertible to the real domain by analytical methods. Therefore, we use numerical inversion methods to transform the obtained solutions from the Laplace domain to the real domain^[24,25]. The Laplace transform overcomes the memory effect arising from the convolution integral expressions of the time fractional derivative term, and better results can be obtained in the case of more general smoothness.

The remaining sections of our paper are organised as follows. Some important preliminary, definitions and lemmas are given in Section 2. The introduction of the model and the time discretization based on the Laplace transform are given in Section 3. The simplification and approximation theory of the equation is given in Section 4. The basis construction, the meshless method for solving the simplified space fractional equation and convergence analysis of the quintic Hermite spline are presented in Section 5. Meanwhile, numerical examples are given in Section 6. Section 7 explains the analysis and the results of the research. Finally, a brief conclusion is given in Section 8.

2.   Preliminary concepts and properties

In this section we will introduce some concepts and properties. Let $\Omega$ satisfy segment conditions of the form ^[26], let $\Upsilon$ be a rectangular domain containing $\Omega$, let the symbol $\cdot |_{ \Omega}$ stand for restriction to $\Omega$, and let $\Upsilon = [a, b]\times [c, d]\supseteq \Omega$.

Definition 2.1. The left and right Caputo fractional derivational of order $\alpha$ on $[a, b]$ is defined by

where $n-1 < \alpha \leq n$, and $n = \left \lceil \alpha \right \rceil$.

Definition 2.2. The left and right Riemann-Liouvile fractional derivatives operator with respect to order $\alpha$ on $[a, b]$ is defined by

where $n-1 < \alpha \leq n$, and $n = \left \lceil \alpha \right \rceil$.

Definition 2.3. Let $\alpha > 0, \; m = \lceil \alpha \rceil$, the connection between Riemann-Liouville derivatives and Caputo derivatives is

Definition 2.4. The Riesz fractional derivative with order $\alpha > 0$ on a finite interval $[a, b]$ is defined by

where

$\alpha \neq 2k+1, k = 0, 1, \cdots,$ and for $n-1 \le \alpha \le n, n \in \mathbb{N}$.

Definition 2.5. For given $v$: $[0, \infty]\rightarrow \mathbb{R}$, the definition of Laplace transform is

Lemma 2.1. ([27,Lemma 1.2.]) Suppose that $v(t)\in C^p[0, \infty)$, the Laplace transform of Caputo fractional derivative about $v(t)$ is

We introduce some definitions of spaces. Let $\Omega$ be a domain in $\mathbb{R}^n$,

where ${\rm P}_{5}$ is the set of polynomial functions with order not greater than $5$ over $k_{j}$.

Definition 2.6. For any nonnegative integer $m$ let $C^m(\Omega)$ denote the vector space consisting of all functions $f$ which, together with all their partial derivatives $D^{\alpha}f$ of orders $\alpha\le m$, are continuous on $\Omega$. We abbreviate $C^0(\Omega)\equiv\Omega$. Let

The subspaces $C_0(\Omega)$ and $C^{\infty}_0(\Omega)$ consist of all those functions in $C(\Omega)$ and $C^{\infty}(\Omega)$, respectively, that have compact support in $\Omega$.

Definition 2.7. Give a positive integer $\tau$ and a real number $r(1 \le r < \infty)$. The Sobolev space $\mathcal{W}^{\tau, r}$ is defined by

the weak derivative $D^{\theta} u \in L^r(\Omega)$ for $0\le |\theta| \le \tau \}$, with norm

where $\theta = (\theta_1, \theta_2)$, $|\theta| = \theta_1+\theta_2$ and $\theta_1$, $\theta_2$ are non-negative integers.

Definition 2.8. (^[26]) For $\forall y \in \partial \Omega$, there exists a nonzero vector $o_y$ and a neighborhood $U_y$ such that if $z \in \bar{\Omega}\bigcap U_y$, then $z+t o_y \in \Omega$ for $0 < t < 1$, and call that $\Omega$ satisfies the segment condition.

Lemma 2.2. (^[26]) If $\Omega$ satisfies the condition of Definition 2.8, then the set of restrictions to $\Omega$ of functions in $C^{\infty}_{0}(R^2)$ is dense in $\mathcal{W}^{\tau, r}(\Omega)$.

3.   Time discretization for multi-term time-space fractional diffusion equations

3.1.   Multi-term time-space fractional diffusion equations

We consider the multi-term time-space fractional diffusion equations of the following form

subject to

where $0 < \alpha_q < \cdots < \alpha_1 < \alpha_0 < 1$ is the time fractional order, ${^{C}_{0}\!D^{\alpha_q} _{t}}$ is the Caputo fractional derivative with $\alpha_q$ order given by Definition 2.1,

are the coefficients, and $\Omega\subset \mathbb{R}^2$ is bounded convex domain.

The spatial operator $\Delta^{\beta, \gamma}$, $1 < \beta, \gamma < 2$ is a Riesz fractional order operator given by Definition 2.4,

where the constants $K_x > 0, K_y > 0$ are diffusion coefficients. And the left side and right side Riemann-Liouville derivatives on $x, y$ direction, respectively, are defined by Definition 2.2,

3.2.   Time discretization based on Laplace transform technique

Suppose that $u(x, y, t)\in C^1(\Omega)$, using the Laplace transform on Eq (3.1) and owing to the property of Lemma 2.1, we have

So this equation could be

Equation (3.2) becomes with the boundary conditions

Then the methods for the inverse Laplace transform methods are based on numerical integration of the Bromwich complex contour integral. From ^[27,29], using the strategy of Talbot, the Bromwich line can be transformed into a contour that starts and ends in the left half plane,

where $\sigma_0$ is the convergence abscissa. Two simpler types of contours have mainly been proposed mainly:

● Parabolic path: $s = \mu(1+{\rm i}z)^2$, $z = \gamma+{\rm i} c$, where $c > 0, -\infty < \gamma < \infty$, then,

● Hyperbolic path:

for $-\infty < \gamma < \infty.$

On either of the above contours, the Bromwich integral becomes

4.   Simplify equation and approximation theory

4.1.   Simplify equation

Let operator $\Delta ^{\beta, \gamma}u(x, y, t)$ on $\Upsilon = [a, b]\times[c, d]$, using the Laplace transform, it becomes as follows:

To avoid the singularity of the operator $\Delta ^{\beta, \gamma}$, let

for $s_l^{\alpha_q} \in \mathbb{C}$, the Eq (3.3) becomes

then expand the $s_l$ and denote

where ${\rm Re}\left(\mathcal{U}(x, y, s_l)\right)$ stands the real part of $\mathcal{U}$, ${\rm Im}\left(\mathcal{U}(x, y, s_l)\right)$ stands the imaginary part of $\mathcal{U}$. Then the Eq (3.3) becomes

Meanwhile, the Eq (3.4) becomes

4.2.   Approximation theory

Let $\Omega_r$ be rectangular domains containing $\Omega$. Denote by $\mathcal{S}$ the set of 2-dimension polynomial functions.

Lemma 4.1. ([28,Theorem 2.2]) $C^{\infty}_{0}(\mathbb{R}^2)|_{\Omega}$ is dense in $C^2(\Omega)$.

Lemma 4.2. ([6,Lemma 2.1]) $\mathcal{S}$ is dense in $C^{\infty}(\Omega_r)$ with the norm $\|\cdot\|_{C^2(\Omega_r)}$.

If $\Omega$ is bounded and closed, then $\Omega$ contains the segment condition, so according the Lemmas 2.2, 4.1 and 4.2, we can obtain the theorem.

Lemma 4.3. Assume that the closed domain $\Omega$ is bounded, then $\mathcal{S}$ is dense in $C^2(\Omega)$.

According to Lemma 4.3, we can obtain the polynomial dense theory.

Remark 4.1. Let $\Omega\subset\Upsilon$ be an arbitrary domain. Then the set of restrictions to $\Omega$ of functions in $S_{5, 2}(\Pi_1)\times S_{5, 2}(\Pi_2)$ is dense in $C^2(\Omega)$, which leads to the set of restrictions to $\Omega$ of functions in $\left(S_{5, 2}(\Pi_1)\times S_{5, 2}(\Pi_2) \right)^2$ is dense in $\left(C^2(\Omega)\right)^2$.

5.   Meshless method based on quintic Hermite spline functions

5.1.   Quintic Hermite spline functions

Let $\Pi_1 = [a, b]$, then the division is

$h$ is the max length of the division. $S_i(x), V_i(x)$ and $W_i(x)$ denote the Hermite splines of

For $S_i(x), V_i(x)$ and $W_i(x)$ from above, we have the following properties:

Remark 5.1. Hermite bases are derived from segmented Hermite interpolating basis functions by the division $\Pi_1$. On the kth divisions $[x_{k}, x_{k+1}], \; k = 0, \cdots, N-1$, it satisfies

where $P(x)$ is the interpolation polynomial and $f(x)$ is the interpolated function, then the number of interpolating basis functions $6$ is obtained. Thus, the total number of basis functions $6N$ on $\Pi_1$ is obtained. The number of Hermite spline functions $3(N+1)$ is obtained from the definition of $S_i(x), V_i(x), W_i(x)$ above, and at the interior points $S_i(x), V_i(x)$ and $W_i(x)$ is a function with two segments and at the endpoints is only a function with one segment. The following theorem will prove that $S_i(x), V_i(x)$ and $W_i(x)$ are the bases.

Theorem 5.1.

is linearly independent and is the base of $S_{5, 2}(\Pi_1)$.

Proof. First, we will show that $S_i(x), V_i(x)$ and $W_i(x)$ is linearly independent. Assume that,

Due to the properties of the Hermite splines, when $x = x_k$, $c_k = 0, \ k = 0, 1, \cdots, N$, then take the derivative of the above

when $x = x_k$, $d_k = 0, \ k = 0, 1, \cdots, N$, then take the derivative of the above

when $x = x_k$, $e_k = 0, \ k = 0, 1, \cdots, N$, so $S_k(x_i), V_k(x_i)$ and $W_k(x_i)$ are linearly independent.

Next, we will verify that it is a base of $S_{5, 2}(\Pi_1)$. Due to the definition of the $S_{5, 2}(\Pi_1)$, so $S_i(x)$, $V_i(x)$, $W_i(x) \in C^2 [a, b]$. On the other hand, $S_i(x)$, $V_i(x)$ and $W_i(x)$ are a piecewise quintic polynomial. Thus, $S_i(x)$, $V_i(x)$, $W_i(x) \in S_{5, 2}(\Pi_1)$.

Since

and

so, $\{S_i(x), V_i(x), W_i(x)\}$ is a base of $S_{5, 2}(\Pi_1)$. □

Then according to Theorem 5.1 and Remark 4.1, it yields a new base

on $\Upsilon$ is dense on $\Omega$. So,

then, using the inverse Laplace transform based the Talbots strategy from Eq (3.5), we could obtain the numerical solution $u_N(x, y, t)$.

5.2.   Meshless method

Definition 5.1. For any $\varepsilon > 0,$ if

then, $\mathcal{U}(x, y, s_l)$ is an $\varepsilon$-approximate solution for Eq (4.1).

We will provide the method of obtaining the $\varepsilon$-approximate solution. First, the minimum bounding rectangle

containing $\Omega$ is given. Subsequently, we will calculate residuals of two parts:

$(1)$ The residual inside $\Omega$ is defined as

$(2)$ The residual on the boundary $\partial \Omega$ is defined as

For any $\varepsilon > 0$, if there exists $\mathcal{U}_{N}(x, y, s_l)$ such that

so, $\mathcal{U}_{N}(x, y, s_l)$ is residual approximate solution of Eq (4.1) on $\Omega$. If

then $\mathcal{U}^*_{N}(x, y, s_l)$ is called the best $\varepsilon$-approximate solution.

Lemma 5.1. $\mathbb{G}$: $(C_{2}(\Upsilon))^2\to (C(\Upsilon))^2$ is a bounded operator.

Proof. For $s_l = (\kappa_l+i\omega_l)$, and denoted that $\mathcal{U}(x, y, s_l)\triangleq\mathcal{U}_l(x, y)$,

where $\theta_1$ is constants, and it could be similarly obtained that

and

Since $A(x, y)$, $(x-a)^{\beta-1}, (b-x)^{\beta-1}, (y-c)^{\gamma-1}, (d-y)^{\gamma-1}$ is continuous, it has

and

similarly,

hence,

so $\mathbb{G}$ is bounded. □

Theorem 5.2. Let $\mathcal{U}(x, y, s_l)$ be the exact solution of Eq (4.1) on $\Omega$, $\mathcal{U}^*_{N}(x, y, s_l)$ be the $\varepsilon$-approximate solution. For every $\varepsilon > 0$, there exists $N_1$, when $N\ge N_1$, coefficients $d^*_{ijl}$ of $\mathcal{U}^*_{N_1}(x, y, s_l)$ from Eq (5.1) satisfy Eq (5.2).

Proof. $\mathcal{U}(x, y, s_l)$ could be approximated by $\mathcal{U}^*_{N_1}(x, y, s_l)$ on $\Upsilon\cap\Omega$. For each fixed $\varepsilon > 0$, there exists $N_1$ such that the residual $\mathbb{L}(\mathcal{U}^*_{N_1}(x, y, s_l))$ satisfies Eq (5.2).

Let $\mathcal{U}_{N_1}(x, y, s_l)$ be residual approximate solutions, taking min $\{ \frac{\varepsilon}{4 { \| \mathbb{G} \| } }, \frac{\varepsilon}{4} \}$, in which $\|G\|$ is defined by

there exists $N_1$ such that the following two parts hold. Inside $\Omega$, we suppose that

when $(x, y) \in \Omega$,

On the $\partial \Omega$, from the boundary condition $\mathcal{U}(x, y, s_l) = 0$, we suppose that

hence, when $(x, y)$ on the $\partial \Omega$,

so,

and

so the theorem holds. □

Theorem 5.3. If Eq (4.1) is well-posed, then $\mathcal{U}^*_{N_1}(x, y, s_l)$ obtained from Theorem 5.2 is the approximate solution of Eq (4.1) on $\Omega$.

Proof. Since $\mathcal{U}^*_{N_1}(x, y, s_l)$ is the $\varepsilon$-approximate solution, for every $\varepsilon > 0$, it yields,

where $\mathbb{G}$ is bounded. It implies that $\mathcal{U}^*_{N_1}(x, y, s_l)$ is the approximate solution of Eq (4.1) on $\Omega$. □

5.3.   Error analysis of meshless method

Let $S_{5, 2}(\Pi_1)$ and $S_{5, 2}(\Pi_2)$ be two quintic spline space with partition

and

The quintic spline bases have the following properties.

Theorem 5.4. Let $u(x)\in C^m[a, b]$, $1\le m \le 5$, then there exists $z(x) \in S_{5, 2}(\Pi_1)$, such that

which $h$ is the partition of the spline space, and $K$ is the constant.

Proof. The division of $[a, b]$ is

$h$ is the max length of the division, and set subinterval

$l_{jk}(x), \; \bar{l}_{jk}(x), \; \bar{\bar{l}}_{jk}(x)$ be Hermite interpolation polynomials, $j = 0, 1, \cdots, N, k = 0, 1,$ and satisfy

First, on $[x_j, x_{j+1}]$, we prove

$1\le p \le d\le 5$, $x_{j0} = x_j, x_{j1} = x_{j+1}$, when $p = 1$, it has

Consider that, $u(y) = (y-x)^p$, it could be interpolated as follows:

Setting $y = x$, we obtain that

Next, due to property of the Hermite interpolation polynomial,

Similarly,

For any $u\in C^m[a, b]$, there has $z(x)\in S_{5, 2}(\pi)$, suppose that

so,

Then, on $\pi_j$ set

let

so

then let

so

such that $g(\xi) = 0$, and

Due to $z'(x)\in P_4$, so $z'(x)$ is the polynomial interpolation of $u'(x)$ at the point $(\xi, u'(\xi)), (x_{jk}, u'(x_{jk})), (x_{jk}, u''(x_{jk})), k = 0, 1$, so

Then, let

so

$\exists \eta_1 \in (x_{j}, \xi), \eta_2 \in (\xi, x_{j+1})$, such that

due to $z''\in P_3,$ so $z''(x)$ is the cubic polynomial interpolation of $z''(x)$, so

Finally,

□

According to ^[30], the following lemma is given.

Lemma 5.2. Let $u(x, y)\in C^m(\Omega)$, $2\le m \le 6$, then there exists

such that

which $h$ is the partition of the space, and $\lambda$ is the constant.

According to Theorem 5.4 and Lemma 5.2, we can infer that:

Remark 5.2. Let $u(x, y)\in C^4(\Omega)$, then there exists

such that

which $h$ is the partition of the space, and $\lambda$ is the constant.

Theorem 5.5. The numerical solution $\tilde{\mathcal{U}}_{N}(x, y, s_l)$ obtained from the proposed meshless method converges to the exact solution $\mathcal{U}(x, y, s_l)$.

Proof. Owing to Theorem 5.1, $\mathbb{H}_i(x)$ is the base of $S_{5, 2}(\pi)$, so numerical solution $\mathcal{U}_{N, l}(x, y)$ obtained from Eq (5.1) belongs to $S_{5, 2}(\pi) \times S_{5, 2}(\pi)$. From the Remark 5.2 and Theorem 5.3, we have

Assume that

so

Then $\exists N$, such that

where $\lambda_1, \lambda_2, \lambda_3$ are constants.

□

6.   Numerical examples

In this section we give two examples to demonstrate the effectiveness of our theoretical analysis. The examples will discuss a single time fractional term and a multiple time fractional term on different domains, respectively. Calculate the

and

where $u(x, y, t)$ is the exact solution, $u_N(x, y, t)$ is the approximate solution by our method. If $t = 1, L_{\infty} = L_{\infty}(1)$. Meanwhile, let the $L = 10$ of Bromwich be integrated by the inverse Laplace transform. The node

from Eq (5.1).

Example 6.1. Consider the single term form Eq (3.1), where $r = 1, K_x = K_y = 1$,

with

Let $\alpha = 2/3$, $\beta = 3/2$, $\gamma = 5/4$. $L_{\infty}$ of Example 6.1 on rectangular domain and circular domain are shown in Table 1.

It can be concluded that our method is valid in a verifiable way and that it gives better results in the general case of smoother time solutions.

(1) When $(x, y)$ on rectangular domains, $\Omega = [0, 1]\times [0, 1],$ the true solution is

The error figure is shown in the Figure 1a at $n = 3\times 3$. And the error $L_{\infty}$ are shown in Table 1.

(2) When $(x, y)$ on circular domains $\Omega$,

the true solution is

The figure of error $u(x, y, t)-u_N(x, y, t)$ when $t = 1$ at $n = 3\times 3$ is shown in the Figure 1b. And the error $L_{\infty}$ are shown in Table 1.

With the above two numerical examples we find that our method gets high accuracy on different regions, showing that our method can handle arbitrary convex regions. Our error convergence is second-order, and since the solution of $u$ with respect to the space is $x^2(1-x)^2y^2(1-y)^2$ or $(x-1/2)^2(y-1/2)^2$, we have fewer points to get a high accuracy error, which is in accordance with the theory. At the same time, the solution of $u$ with respect to the time is $t^{\frac{4}{5}}$, Laplace transform can be used to deal with lower order smooth solutions.

Example 6.2. ^[31] Consider the multi-term from Eq (3.1), {where $r = 4, K_x = K_y = 1$

with

where

Then the exact solution is

(1) When $(x, y)$ on rectangular domains,

When

We calculate the error $E(T)$ and compare it with ^[31] in Table 2 at $T = 1$.

(2) When $(x, y)$ on circular domains, the

When

where $\beta = 1.02$, $\gamma = 1.02.$

The numerical solution and the absolute errors when $t = 1$ at $n = 3\times 3$ are shown in Figure 2.

We calculated the $L_2$ error $E(T)$ on the rectangular domain and compared it with ^[31]. It can be seen that we obtain higher accuracy with fewer points, which proves the high efficiency of our method. We also carry out experiments with different parameters $\alpha, {\bf a}$ and $\beta, \gamma$ on the circular domain and calculate the absolute errors at the $n = 3$.

From Figure 2, it can be seen that our method also achieves high error accuracy, indicating the applicability of our method. The high error accuracies obtained by our method in different regions and also with different parameters show the stability and efficiency of our method. Because of the high smoothness of $u$ with respect to $x, y$, we get high error accuracy with fewer points, which is consistent with our theoretical analysis.

7.   Conclusions and discussion

In this paper, we proposed a meshless method of solving the minimum residual approximate solution for Eq (3.1). Different from previous methods, we use the Laplace transform method to deal with the multi-term time fractional operator, we transform the time into complex frequency domain by Laplace transform, Eq (3.1) is transformed into complex equation Eq (3.3). Then, on the spatial direction, we proposed a quintic Hermite meshless method to deal with space fractional operators on arbitrary convex region based on the theory of polynomial functions dense theorem. The approximate accuracies become higher by increasing number of Quintic Hermite spline functions. The minimum residual approximate solution of Eq (4.1) is obtained by Theorems 5.3 under the condition of well-posed equations. Meanwhile, using Theorem 5.4 and Lemma 5.2, it infers Remark 5.2, which is the convergence of the biquintic spline function. Then by using Remark 5.2 and Theorem 5.3, we can obtain Theorem 5.5 to show the convergence of the method in the spatial direction. We use numerical inversion methods to transform the obtained the minimum residual approximate solution from the Laplace domain to the real domain by using the strategy of Talbot through parabolic path.

In Numerical examples, we fix the $L = 10$ in Eq (3.5) by parabolic path to get the numerical solution. First, we handle the single term time-space fractional diffusion equations, we can deduce that the method can deal with time fractions that are not sufficiently smooth, and we can get higher precision with fewer nodes in arbitrary convex region from Table 1 and Figure 1. This also proves that Laplace transform is effective for dealing with insufficiently smooth time-fractional operators. Then, we solve the multi-term time-space fractional diffusion equations with $4$ terms. These results are compared with ^[31], and it is found that our method achieves better accuracy with fewer points. At the same time, we found that the accuracy of the single term is better than that of the multi-term. In addition, the accuracy is higher on rectangular areas than on circular areas. These experimental results are consistent with theoretical expectations and demonstrate the effectiveness and efficiency of our method.

In this paper, the use of the extension theorem allows the meshless method to be applied to arbitrary convex regions in two dimensions, and the use of the Laplace transform allows to deal with multi-term low-order time solutions. In the future, through the study of spatial Riesz operators, we will investigate meshless methods for solving equations in arbitrary regions of higher dimensions. In addition, this method can also be used to study equations of time-distributed order.

Use of AI tools declaration

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

Acknowledgments

The authors would like to thank the referees for their helpful comments and suggestions, which lead to a much improved version of this paper.

Conflict of interest

The authors declare that they have no conflicts of interest.