A stabilized multiple time step method for coupled Stokes-Darcy flows and transport model

1.   Introduction

Fractional differential equation models have been founded in a lot of fields of science and engineering, such as physics, chemistry, biology, dynamics, and control ^[1,2,3,4]. In these models, the fractional wave model plays an important role in many practical application fields including transmission and modeling propagation of electrical signals, neural conduction, weak current propagation in the animal nervous system, wave phenomena, and wave propagation. However, it is often difficult to get the analytic solutions of these complex problems. In view of the importance of this kind of model, more and more scholars have focused on solving them numerically by developing a lot of numerical methods including finite element method ^{[5,6,7,8,9,10,11,12,13,14]}, wavelet method ^[15], finite difference method ^{[16,17,18,19,20,21,22,23]}, meshless method ^[24], collocation method ^[25,26], and B-spline method ^[27,28].

In this article, we consider the initial-boundary problem of the following nonlinear fractional hyperbolic wave model

where $\Omega = (a, b)$ is the spatial domain and $J = (0, T]$ with $0 < T < \infty$ is the time interval. $u_0(x)$ and $u_1(x)$ are given initial functions, $f(x, t)$ is the given source term and the nonlinear term $g(u) \in C^2(R)$, fractional parameter $\beta = \alpha+1$, and $_0^RD_t^\gamma w(x, t)$ is the Riemann-Liouville fractional-order derivative defined by

and

The fractional hyperbolic wave model (1.1), which includes both propagation and diffusion of the wave, can be degenerated into the pseudo-hyperbolic equation for $\beta = 2$ and diffusion equation for $\beta = 1$.

In the following, for formulating our numerical method we need to introduce numerical techniques including the weighted and shifted Grünwald difference (WSGD) formula, BDF2-$\theta$, $H^1$-Galerkin MFE method, and time two-mesh (TT-M) finite element algorithm. The WSGD formula, which was proposed by Tian et al. in ^[29], is a useful approximate method for the Riemann-Liouville fractional derivative. Due to its high-order approximation characteristics, many scholars have developed efficient numerical methods based on the WSGD formula; see ^{[30,31,32,33,34,35]}. The $H^1$-Galerkin MFE method is an important numerical method, which was proposed by Pani ^[36]. Due to several advantages of this method, many scholars have begun to use it to solve evolution partial differential equation (PDE) models, such as integer PDE models ^[37,38,39], fractional PDE models ^[40], and distributed-order PDE models ^[41]. The TT-M finite element method was proposed by Liu et al. in ^[42] to quickly solve the fractional water wave model, which can also combine many other numerical methods, such as the finite difference method and the finite volume element method, to solve evolution differential equation models ^{[43,44,45,46]}.

In this article, considering the characteristics of the nonlinear fractional hyperbolic wave equation, we introduce an auxiliary function with a fractional derivative, and formulate a fast high-order fully discrete $H^1$-Galerkin MFE method, where the time direction is discretized by the BDF2-$\theta$ with the WSGD operator, the space direction is approximated by the $H^1$-Galerkin MFE method, and the fast TT-M algorithm is used to reduce calculation time. The main works and contributions of this article are as follows:

(Ⅰ) Propose a fast TT-M mixed element method with the WSGD operator to numerically solve the nonlinear pseudo-hyperbolic wave equation with two term fractional derivatives.

(Ⅱ) Introduce a special auxiliary function, transform the original high-order equation into the coupled system of equations with lower order space-time derivatives, and directly formulate a second-order fully discrete BDF2-$\theta$ $H^1$-Galerkin MFE system, which can avoid difficulties in numerical calculations and theoretical analysis by directly discretizing fractional derivatives. Further, develop the fast fully discrete TT-M MFE system, and derive optimal a priori error estimates for two functions.

(Ⅲ) Provide the detailed numerical algorithm by taking smooth and weakly regular solutions. Validate the correctness of the theoretical results and the effectiveness of the numerical algorithm, and illustrate that the TT-M MFE method has good computational efficiency by comparing the calculation results with the standard nonlinear MFE method.

The rest of the article is outlined as follows: In Section 2, the fully discrete scheme based on the combination of an MFE method and the BDF2-$\theta$ with the WSGD formula is derived. In Section 3, the optimal error estimates in both $L^2$-norm and $H^1$-norm for the fully discrete TT-M MFE scheme are derived. In Section 4, the numerical algorithm is shown. Some experiments in Section 5 are conducted to further confirm our theoretical results. Finally, in Section 6, conclusions and advancements are provided.

2.   Weak form and TT-M MFE scheme

Letting $q = _0^R\!\!D_t^\alpha u_{x}(x, t)+u_x(x, t)$, we rewrite equation (1.1) as

We multiply the first equation of (2.1) by $v_x$ and the second equation of (2.1) by $-\omega_x$, respectively, then make the inner product on the spatial domain $\bar{\Omega} = [a, b]$ to have

For the second equation of (2.2), by the integration by part and the boundary condition

we obtain

Now, we can get the following mixed weak form

For obtaining the fully discrete TT-M MFE scheme, we introduce the nodes $t_n = n\tau_c(n = 0, 1, 2, \cdots, N)$ in the time interval $[0, T]$, where $t_n$ satisfies $0 = t_0 < t_1 < t_2 < \cdots < t_N = T$ with fine mesh length $\tau = T/NM$ and coarse mesh length $\tau_c = M\tau$ for some positive integer $N$. Define $u^n = u(\cdot, t_n), q^n = q(\cdot, t_n)$ for smooth functions $u$ and $q$ on $[0, T]$. Some useful lemmas will also be introduced as follows.

Lemma 2.1. (^[35]) With $v(t) \in C^3[0, T]$, at time $t_{n-\theta}$, the following formula with second-order accuracy for approximating the first-order derivative holds

where

for any $\theta \in [0, \frac{1}{2}]$.

Lemma 2.2. At time $t_{n-\theta}$, the following important results hold for any $\theta \in [0, 1]$ and $v(t) \in C^2[0, T]$,

Lemma 2.3. (^[29,35]) At time $t_n$, the second-order approximate formula for the Riemann-Liouville fractional derivative with parameter $\gamma \in (0, 1)$ holds

with

where series $w_i^\gamma$ are defined as $w_0^\gamma = 1, \; w_l^\gamma = (-1)^l\big{(}\begin{matrix} \gamma \\ l \end{matrix}\big{)} = \frac{\Gamma(l-\gamma)}{\Gamma(-\gamma)\Gamma(l+1)}, l\geq1,$ which satisfy $w_l^\gamma < 0, w_l^\gamma = (1-\frac{\gamma+1}{l})w_{l-1}^\gamma, (l = 1, 2, \cdots), \sum\limits_{l = 1}^{\infty}w_l^\gamma = -1.$

Based on the weak form (2.3) and the numerical approximate formulas above, we can get the following equivalent weak form

Case $n = 1$:

Case $n\geq 2$:

where

We now formulate the fully discrete TT-M MFE system at time $t_{n-\theta}$ for handling the computational time-consuming problems of implicit finite element systems, and we denote $U_c^n, Q_c^n$ as solutions of the system on the time coarse mesh and $U_f^m, Q_f^m$ as solutions of the system on the time fine mesh. The TT-M MFE algorithm can be implemented as the following three steps.

$\mathbf{STEP 1}$: First, we arrive at the following nonlinear coupled system based on the time coarse mesh $\tau_c$: Find $(U_c^n, Q_c^n):[0, T]\times[0, T]\mapsto V_h\times W_h$ such that

Case $n = 1$:

Case $n\geq 2$:

$\mathbf{STEP 2}$: Second, we can get all the interpolated values $U_I^m (m = 0, 1, \cdots, M, M+1, \cdots, 2M, \cdots, NM)$ by using an interpolation formula

where $\; \lambda_m = n-\frac{m}{M} \in [0, 1) (n = \lceil\frac{m}{M}\rceil)$ and $U_I^0 = U_c^0$. Values of $Q_I^m$ can be obtained similarly.

$\mathbf{STEP 3}$: Finally, we establish the following linear system on the time fine mesh $\tau$ based on the solutions $U_I^m, Q_I^m$; that is, to find $(U_f^m, Q_f^m):[0, T]\times[0, T]\mapsto V_h\times W_h$ for any $(v, \omega)\in V_h\times W_h$ such that

Case $m = 1$:

Case $m\geq 2$:

where finite element spaces are defined as

Remark 2.4. Here, we provide two other equivalent linearized techniques besides the one mentioned in (2.15)-(2.16).

3.   Error estimates

For subsequent analysis, we introduce some useful lemmas.

Lemma 3.1. (^[30,35]) Let $\mathcal{A}_\gamma(i)$ be defined in $\textrm{(2.8)}$, then for any positive integer $L$ and real vector $(v^0, v^1, \cdots, v^L) \in R^{L+1}$, the following inequality holds

Lemma 3.2. (^[23,35]) For series ${\chi^n}$ $(n\geq 2)$, the following inequality holds

and

Lemma 3.3. (^[47,48]) For any function $v \in H_0^1(\Omega)$, we have

For considering a priori error estimates for the TT-M MFE system, the projection operator and the inequality should be introduced.

Lemma 3.4. (^[36]) Define an elliptic-projection operator $\Upsilon_{h}: H^{1}_{0}(\Omega)\rightarrow V_{h}$, for any $\phi_{h} \in V_{h}$ such that

with an estimate inequality

Lemma 3.5. (^[36]) Define a Ritz-projection operator $\Pi_{h}:H^{1}(\Omega) \rightarrow W_{h}$ by

where $\mathcal{A}(q, \phi)\doteq (q_x, \phi_x)+\lambda(q, \phi)$, $\mathcal{A}(\phi, \phi) \geq \mu_{0}\|\phi\|^{2}_{1}, \mu_{0} > 0$ is a constant. Further, the following estimate inequality holds

Theorem 3.6. Let $u(\cdot, t_n)$, $q(\cdot, t_n)$ be the solutions of system $(1.1)$ and suppose $U_c^n$, $Q_c^n$ and $U_f^m$, $Q_f^m$ are the solutions of TT-M MFE systems $(2.9)$-$(2.10)$ and $(2.15)$-$(2.16)$, respectively, then there exists a constant $C > 0$ that depends only on $u(\cdot, t_n)$, $q(\cdot, t_n)$, such that

and

Proof. For convenience, we write error as

(1) Error estimate on the time coarse mesh.

Applying the projection operators in Lemmas 3.4 and 3.5, the error equation on the time coarse mesh is as follows:

Case $n = 1$:

Case $n\geq 2$:

Set $\omega_h = \sigma_c^{n-\theta}$ in (3.12), and use Lemma 3.3, the Cauchy-Schwarz inequality, and the Young inequality to obtain

Multiply (3.13) by $4\tau_c$, replace $n$ with $l$, and sum for $l$ from $2$ to $n$ to arrive at

Setting $v_h = \xi_c^{n-\theta}$ in (3.12), summing the resulting equation from 1 to $n$, and using the Cauchy-Schwarz inequality as well as the Young inequality, we have

Applying Lemma 3.1 and the Poincaré inequality, we obtain, for $n\geq 1$,

For the term $\mathbb{H}(\sigma_c^1)$, we take $\omega_h = \sigma_c^{1-\theta}$ in (3.11) and apply the Cauchy-Schwarz inequality as well as the Young inequality to have

Omitting the nonnegative term on the left hand side of (3.17), we obtain

Substitute (3.18) into (3.14), apply (3.16), and use the Gronwall inequality to have

Notice that the inequalities (3.6) and (3.8) hold; combine (3.16) and (3.19) with the triangle inequality to finish the proof of the first result of Theorem 3.6.

(2) Error estimate on the time fine mesh.

Based on Lemmas 3.4 and 3.5, the error equation on the time fine mesh is as follows:

Case $m = 1$:

Case $m\geq 2$:

For the nonlinear term on the right hand side of (3.21), we use Taylor's formula to get

Set $\omega_h = \sigma_f^{m-\theta}$ in (3.21) and use (3.22), the Cauchy-Schwarz inequality, and the Young inequality to arrive at

Using a similar derivation to (3.14), we have

To estimate $\mathbb{H}(\sigma_f^1)$, we set $\omega_h = \sigma_f^{1-\theta}$ in (3.20) and apply Taylor's formula to deal with the nonlinear term to arrive at

Combining (3.25) with (3.24), we have

Setting $v_h = \xi_f^{m-\theta}$ in (3.12) and using a derivation similar to (3.16), we get

We now estimate the error $C\tau \sum\limits_{l = 1}^{m}\|u^{l-\theta}-U_I^{l-\theta}\|_{L^4}^4$. Denote $n = \lceil\frac{l}{M}\rceil$ as the smallest integer that is equal to or greater than $\frac{l}{M}$, then by the notations introduced in (2.11), we get

where $\bar{t}_{l-\theta} \in (t_{n-\theta-1}, t_{n-\theta})$. For $\lambda_l \in [0, \frac{1}{2}]$, follow the idea from ^[46] and use (3.9) and (3.28) to obtain the following result

Using the techniques applied to (3.29), we easily get the inequality

Making use of Lemma 3.3, (3.30), and (3.29), we can obtain

Substitute (3.27) and (3.31) into (3.26) and apply the Gronwall inequality to obtain

Combine (3.27), (3.30), (3.32) and (3.6) with (3.8) and use the triangle inequality to finish the proof of the second result of Theorem 3.6.

4.   Numerical algorithm

In this section, we provide a numerical algorithm for solving the examples with smooth solutions and weakly regular solutions. For the solution $u$ with weak regularity, referring to ^[49,50], we split it into the smooth part and the weak regular part as the following

where $c_k = c_k(x)$ are coefficient functions, parameters $\sigma_{k}$ satisfy $0\leq\sigma_1 < \cdots < \sigma_{j+1}$, $\sigma_j < 3$ and $\sigma_{j+1}\geq 3$ and $\varrho$ is sufficiently smooth with respect to $t$. Thus, we can think of $u_1$ as the nonsmooth part of the $u$, which may cause a loss of accuracy in time. For solving this problem, based on the idea presented in ^[51], we develop a corrected technique by adding correction parts. We now discretize the spatial domain $\bar{\Omega}$ as $a = x_0 < x_1 < \cdots < x_L = b$, where the nodes are $x_k = x_0+kh$ with the uniform spatial step size $h = \frac{b-a}{L}$. Next, considering mixed linear element spaces with linear basis functions $\{\phi_i(x)\}_{i = 0}^{L}$ and $\{\varphi_i(x)\}_{i = 0}^{L}$, we can write numerical solution $U_c$ and $Q_c$ as: $U_c^n = \sum\limits_{i = 0}^Lu_i^n\phi_i$, $Q_c^n = \sum\limits_{i = 0}^Lq_i^n\varphi_i$, respectively. Based on the numerical scheme (2.9)-(2.10) combined with the corrected technique, we formulate a numerical algorithm in the matrix form.

Case $n = 1$:

Case $n\geq 2$:

where

In the above algorithm, $\{\omega_{n, k}^{(\alpha)}\}_{k = 1}^j$, $\{\omega_{n, k}^{(0)}\}_{k = 1}^j$ are correction weights of $I_\alpha^{n-\theta}[U_{cx}^{n-\theta}]$, $U_{cx}^{n-\theta}$, respectively. $\{\bar{\omega}_{n, k}^{(0)}\}_{k = 1}^{j+1}$ are correction weights of $Q_{c}^{n-\theta}$. $\{\tilde{\omega}_{1, k}^{(1)}\}_{k = 1}^{j+1}$ are correction weights of $\partial_t[u^1]$. $\{\bar{\omega}_{n, k}^{(1)}\}_{k = 1}^{j+1}$ are correction weights of $\partial_t[u^{n-\theta}]$. The correction weights $\{\omega_{n, k}^{(\alpha)}\}_{k = 1}^j$ can be obtained by the following formula ^[51]

where

Similarly, we can get the correction coefficients $\{\omega_{n, k}^{(0)}\}_{k = 1}^j$, $\{\bar{\omega}_{n, k}^{(0)}\}_{k = 1}^{j+1}$, $\{\bar{\omega}_{n, k}^{(1)}\}_{k = 1}^{j+1}$, and $\{\tilde{\omega}_{1, k}^{(1)}\}_{k = 1}^{j+1}$ by using (4.5).

We divide the calculation process into two parts. First, we calculate $\textbf{u}_c^k$ and $\textbf{q}_c^k$ by (4.4)-(4.5), where $k = 1, 2, \cdots, j+1$, then we can obtain $\textbf{u}_c^m$ and $\textbf{q}_c^m, m > j+1$ by $\textbf{u}_c^k$ and $\textbf{q}_c^k$. The process of computation on the fine mesh is similar to that on the coarse mesh, so we will not introduce details here.

5.   Numerical tests

Here, for showing the feasibility and validity of our numerical method and the efficiency of the TT-M MFE system, we consider a linear element and provide the computing results by our numerical procedure.

5.1.   Example 1

In this example, we use the linearized method (2.17)(a) to finish our calculations. Considering the space domain $\bar\Omega = [0, 1]$ and the time interval $\bar{J} = [0, 1]$, we take the nonlinear term $g(u) = u^3-u$, the following given source term

and then easily validate that the exact solution is $u = t^{3+\alpha}\sin\pi x$ and the corresponding auxiliary function is $q = \frac{\Gamma(4+\alpha)} {\Gamma(4)}t^3\pi\cos\pi x+t^{3+\alpha}\pi\cos\pi x$.

In Table 1, with the fixed time step length $\tau = \frac{1}{{10000}}$, changed space step length $h = \frac{1}{200}, \frac{1}{300}, \frac{1}{400}$, fractional parameter $\alpha = 0.3, 0.8, 0.99$, and shifted parameter $\theta = 0.1, 0.3, 0.5$, we calculate the spatial convergence results of the standard nonlinear mixed element algorithm under different parameters and record the Central Processing Unit (CPU) time required for the algorithm. One can see that for the currently selected exact solution, the spatial convergence results are optimal, which is consistent with the theoretical results of the linear element ($k = r = 1$) we selected. Further, in Table 2, based on the chosen changed parameters as in Table 1 with the fixed time step length $\tau = \tau_c/M = 1/NM = 1/10000$ ($N = M = 100$) for the TT-M MFE method, we get the optimal convergence results and CPU time. Comparing the data in Tables 1–2, one can see that the fast TT-M MFE algorithm can greatly reduce the CPU time while maintaining the same convergence accuracy as the standard nonlinear mixed element method.

In Tables 3–4, by taking the space step length $h = 1/5000$, time step length $\tau = 1/144$, $1/256$, $1/400$ ($\tau = \tau_c^2$ for TT-M method), time fractional parameter $\alpha$ as $0.3, 0.8, 0.99$, and shifted parameter $\theta$ as $0.1, 0.3, 0.5$, we implement the numerical calculations by using standard nonlinear MFE method and fast TT-M MFE method, respectively. From this, one can see that these two methods have almost the same error results and time convergence rate, and that our TT-M MFE algorithm can save the CPU time.

5.2.   Example 2

In this example, we continue to use the linearization technique (2.17)(a) to verify the efficiency of the current TT-M MFE algorithm. Considering the space-time domain $\bar\Omega\times\bar{J} = [0, 1]\times[0, 1]$, we choose the nonlinear term $g(u) = \arctan u$ and the source term $f(x, t) = 100\sin^2(5\pi t)\sin^2(3\pi x)\big{[}0.15-(t-\frac{1}{2})^2-(x-\frac{1}{2})^2\big{]}^2$. Here, we just consider the TT-M MFE algorithm with $M = 4$. Because of the unknown exact solution, we choose the numerical solution with $h = \tau = \frac{1}{1200}$ as the approximating exact solution.

In Table 5, with the fixed time step length $\tau = \tau_c/M = 1/NM = 1/1200$, changed space step length $h = 1/30, 1/40, 1/50$, fractional parameter $\alpha = 0.3, 0.5, 0.99$, and shifted parameter $\theta = 0.1, 0.3, 0.5$, we can get the errors and spatial convergence results of the TT-M MFE system. In Table 6, considering the fixed space step length $h = 1/1200$, fine time step length $\tau = \tau_c/M = 1/NM = 1/80, 1/100, 1/120$ ($N = 20, 25, 30$), fractional parameter $\alpha = 0.3, 0.5, 0.99$ and shifted parameter $\theta = 0.1, 0.3, 0.5$, we calculate the error results and time convergence rate for the TT-M MFE algorithm. The computed data shows that the TT-M MFE algorithm can also maintain a good calculation effect for the selected numerical example with an unknown exact solution.

Further, in order to check the behaviors of numerical solution, we provide the comparison figures of numerical solutions between different time step length sizes. In Figure 1, we show the comparison surfaces of numerical solutions $U_f$ with the fixed space step length $h = 1/1200$, fractional parameter $\alpha = 0.3$, shifted parameter $\theta = 0.1$, and changed time step length $\tau = 1/120, 1/1200$. We also provide the comparison surfaces of numerical solutions $Q_f$ in Figure 2. The comparison results tell us the corresponding numerical solutions have similar behavior. Moreover, in Figure 3, for fixed fractional parameter $\alpha = 0.3$ and parameter $\theta = 0.1$, we depict the figures of difference in $L^2$-norm between reference solution with $h = \tau = 1/1200$ and numerical solution with $h = 1/1200$ and $\tau = 1/120$, from which one can see the performances of $\|u^n-U_f^n\|$ and $\|q^n-Q_f^n\|$. It is easy to see the changes of actual errors at different time nodes from the figures, which can reveal the overall distribution of errors.

5.3.   Example 3

For comparison and validation of feasibility, we still carry out the numerical calculation by taking Example 1. Here, we apply the linearized technique (2.16) to deal with the nonlinear term. One can see from the numerical results in Table 7 that the optimal spatial convergence data is almost consistent with the calculation results in Example 1, which uses the linearized method (2.17)(a). It indicates that the linearization technique adopted in this paper is feasible. Further, comparison of CPU time in Table 2 and Table 7 shows that computing time in this example is slightly slower, which may be caused due to the linearization for the $m-1$ layer.

5.4.   Example 4

In this example, we choose the solution $u = t^{2+\alpha}\sin(\pi x)$, which has weaker regularity with comparison to the case in Example 1. We choose the same space-time domain and the nonlinear term $g(u)$ as in Example 1. We provide the source term $f$, which we omit here, such that the equation has the current exact solution $u$. For this case with weak regularity, we continue to apply the linearized technique (2.16) to deal with the nonlinear term.

By taking $\tau = 1/16, 1/25, 1/36, 1/49$, $\alpha = 0.1, 0.3, 0.5$, $\theta = 0.1, 0.3, 0.5$, and the fixed space step $h = 1/1000$, we implement numerical tests and obtain the numerical results shown in Table 8, from which one can see that most data cannot achieve second-order approximation results in time. For solving this issue, under the same parameters, we consider the corrected scheme with correction parts, and arrive at optimal time second-order convergence results listed in Table 9, which imply that the numerical scheme by adding the correction parts can effectively solve the problem of accuracy loss and restore the optimal convergence order. Further, based on the data from Tables 8–9, we show the case of the convergence rate in Figures 4–5 for $U_f$ and $Q_f$, from which one can see intuitively, that with comparison to the case without adding the correction parts, the optimal convergence rate can be achieved by adding the correction parts.

6.   Conclusions and Advancements

In this article, we developed a fast TT-M MFE method for solving the nonlinear fractional hyperbolic wave model. We derived optimal a priori error results for the fully discrete TT-M MFE scheme. To verify the correctness of theoretical results and the computational efficiency of the algorithm, we implemented four numerical tests. For the cases with smooth solutions, one can see from the computing results that our TT-M MFE algorithm can obtain the similar convergence results as that computed by using the nonlinear MFE algorithm, while the computing time was reduced. Further, for the case with a weakly regular solution, the considered numerical scheme under certain parameters may lose computational accuracy. For handling this problem, we designed the corrected TT-M MFE method by adding the correction term to restore calculation accuracy.

In future works, we will design other TT-M MFE methods to solve some nonlinear fractional PDE models.

Acknowledgments

The authors would like to thank the editor and all the anonymous referees for their valuable comments, which greatly improved the presentation of the article. This work is supported by the National Natural Science Foundation of China (12061053, 12161063), Natural Science Foundation of Inner Mongolia (2022LHMS01004), Young innovative talents project of Grassland Talents Project, Program for Innovative Research Team in Universities of Inner Mongolia Autonomous Region (NMGIRT2413, NMGIRT2207).

Use of AI tools declaration

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

Conflict of interest

The authors declare there is no conflict of interest.