In this paper, we analyze the e ect of environment noise on the transmission dynamics of a stochastic hepatitis B virus (HBV) infection model with intervention strategies. By using the Markov semigroups theory, we define the stochastic basic reproduction number and find it can be used to govern disease extinction or persistence. When it is less than one, under a mild extra condition, the stochastic system has a disease-free equilibrium and the disease is predicted to die out with probability one. When it is greater than one, under mild extra conditions, the model admits a stationary distribution which means the persistence of the disease. Thus, we observe that larger intensity of noise (resulting in a smaller stochastic basic reproduction number) can suppress the emergence of hepatitis B outbreak. Numerical simulations are also carried out to investigate the influence of information intervention strategies that may change individual behavior and protect the susceptible from infection. Our analysis shows that the environmental noise can greatly a ect the long-term behavior of the system, highlighting the importance of the role of intervention strategies in the control of hepatitis B.
1.
Introduction
In this article, we consider the following time-fractional generalized Rosenau-RLW-Burgers equation:
with boundary conditions
and initial condition
where Ω=(a,b) is the spatial domain, J=(0,T] is the time interval with T∈(0,∞), and g(x,t) is a known source term function. The nonlinear term f(u) satisfies the assumption condition |f(u)|≤cf(u)|u|, where cf(u) is a positive constant on u. C0Dαtu and C0Dβtu are both Caputo fractional derivatives with 0<α,β<1. Since C0Dγtu=∂γ(u−u0)∂tγ, all of the above Caputo fractional derivatives can be converted into the Riemann-Liouville fractional derivative, note that
Specifically, when α=1, β=1, (1.1) degenerates into the generalized Rosenau-RLW-Burgers equation which can be seen as the combined system between the generalized Rosenau-RLW equation and the generalized Rosenau-Burgers equation.
The RLW equation, the Rosenau equation, and their combined systems with other equations are significant mathematical and physical equations that effectively describe nonlinear wave behaviors. These equations have become interesting topics in the study of nonlinear dispersion dynamics. Since obtaining analytical solutions for these equations is challenging, studying their numerical methods is paramount. Over the years, there has been extensive research on numerical methods for solving this type of equation. In [1], Atouani and Omrani discussed the numerical solution of the Rosenau-RLW (RRLW) equation based on the Galerkin finite element method. In [2], He and Pan developed a three-level, linearly implicit finite difference method for solving the generalized Rosenau-Kawahara-RLW equation. In [3], Wongsaijai and Poochinapan developed a pseudo-compact finite difference scheme for solving the generalized Rosenau-RLW-Burgers equation. In [4], Mouktonglang et al. analyzed a generalized Rosenau-RLW-Burgers equation with periodic initial-boundary value. For more papers on related equations, please refer to [5,6,7,8]. It is worth noting that the literature on the fractional generalized Rosenau-RLW-Burgers equation is relatively scarce, and its analytical solution is difficult to obtain. Therefore, we have to consider effective numerical methods such as finite element methods [9,10,11,12], finite difference methods [13,14,15], finite volume methods [16], spectral methods [17,18,19], and mixed finite element methods [20,21,22]. In addition, the existence of time-fractional derivatives increases the difficulty of studying numerical methods. Therefore, it is crucial to choose an appropriate high-order approximation formula for the fractional derivative to establish a stable numerical scheme for (1.1).
In 1986, Lubich [23] proposed the convolution quadrature (CQ) formula for Riemann Liouville fractional operators using the discrete convolution. In [24], Chen et al. developed an alternating direction implicit fractional trapezoidal rule type to solve a two-dimensional fractional evolution equation. In [25], Jin et al. proposed a corrected approximation formula for high-order BDFs through appropriate initial modifications to discretize fractional evolution equations. Based on the CQ formula, in [26], Liu et al. developed the shifted convolution quadrature (SCQ) theory, which extended the CQ formula at xn−θ and discussed the constraints of parameter θ. In [27], Yin et al. studied the generalized BDF2-θ with the finite element method for solving the fractional mobile/immobile transport model, and also developed a correction scheme by adding the starting part to restore convergence order. For more related papers, please refer to [28,29,30,31,32,33].
In this article, we develop the generalized BDF2-θ in time combined with the mixed finite element method in space to solve (1.1). The focuses of this article are as follows:
● It is noted that the time-fractional generalized Rosenau-RLW-Burgers equation containing two time-fractional operators is studied.
● The stability of the time-fractional generalized Rosenau-RLW-Burgers equation (1.1) based on the mixed finite element method is given.
● Based on a comprehensive analysis of some numerical examples, the numerical method's feasibility and effectiveness have been extensively validated. Specifically, the issue of decreasing the convergence rate of nonsmooth solutions is solved by adding correction terms.
The structure of this article is as follows: In Section 2, the generalized BDF2-θ is introduced, and the fully discrete mixed finite element scheme is provided. In Section 3, the existence and uniqueness theorem for the fully discrete mixed finite element scheme is given. In Section 4, the stability of the scheme is proved. In Section 5, some numerical examples with smooth and nonsmooth solutions based on the discrete scheme are presented. In Section 6, some conclusions are given.
2.
Numerical scheme
In this section, we present the fully discrete mixed finite element scheme for (1.1) in space, which combines the generalized BDF2-θ in time. The generalized BDF2-θ with the starting part is introduced in [27]. Further, we divide the time interval [0,T] into 0=t0<t1<⋯<tN−1<tN=T, and let tn=nτ(n=1,2,⋯,N), where τ is time step length size and N is a positive integer.
For the convenience of research, set ˆu:=u−u0, and assume that ˆu has the following form:
where cj=c(x), 1<σ1<σ2<⋯<σκ<σκ+1 and ϕ(x,t) is sufficiently differentiable with respect to t.
Using ˆu:=u−u0, we can write (1.1)–(1.3) as
with boundary conditions
and initial condition
where ˆg(x,t)=g(x,t)+(u0)x−(u0)xx.
Now, we introduce an auxiliary variable q=ˆuxx to obtain the following coupled system:
and
Multiplying (2.5) and (2.6) by v∈H10 and w∈H10, respectively, integrating the result equations, and using integration by parts, we obtain the following weak form:
and
To provide the fully discrete numerical scheme, we first introduce the relevant formulas and lemmas for the generalized BDF2-θ.
For smooth functions ˆu and q in [0,T], we let ˆun=ˆu(⋅,tn), qn=(⋅,tn). The approximation formula for the Riemann-Liouville fractional derivative at time tn−θ with the generalized BDF2-θ is
where |Rn−θγ|≤Cτ2.
The discrete convolution part is denoted as
and the starting part is
The convolution weights {ω(γ)j}∞j=0 in (2.10) are generated by the following generating function:
Lemma 2.1. [27] We give the convolution weights {ω(γ)j}∞j=0 of the generalized BDF2-θ as follows:
Lemma 2.2. [27] The starting weights {ω(γ)n,j}κj=1 of the generalized BDF2-θ are given as the following:
Lemma 2.3. [12,15] For ˆu∈C4[0,π], the following two approximate formulas at tn−θ hold:
where gn−θ:=(1−θ)gn+θgn−1 and f(ˆun−θ):=(2−θ)f(ˆun−1)−(1−θ)f(ˆun−2).
Next, we have the following approximate formula:
Without considering the starting part, we can obtain the weak form of (2.5) and (2.6) at tn−θ:
and
where Rn−θ1=O(τ2) and Rn−θ2=O(τ2).
To establish the fully discrete mixed finite element scheme, we introduce the following finite element space:
where Th is a subdivision of ˉΩ=[a,b] into M subintervals Ii=[xi−1,xi], with hi=xi−xi−1, h=max1≤i≤Mhi, and Pk(Ii) represent the polynomials with a degree less than or equal to k in Ii.
Next, we provide linear basis functions {φi}Mi=1 of finite element space Vh as follows:
Based on the above finite element space, we find {Un−θ,Qn−θ}∈Vh×Vh satisfying
and
3.
Existence and uniqueness of numerical solution
Theorem 3.1. The solution of the fully discrete mixed finite element scheme (2.21) and (2.22) is uniquely solvable.
Proof. Taking basis functions {φi}Mi=1 of finite element space Vh, we have
Taking V=φj and W=φj from (2.21) and (2.22), we have
and
where
Obviously, A and B are symmetric and positive definite. Further, processing the boundary and simplifying the right-hand term, we have
and
where
Multiplying (3.4) by τ˜A−1, we have
Further, rewrite (3.5) as
where Hn3=(1−θ)−1˜A−1Hn2−˜A−1˜BUn.
Substitute (3.7) into (3.6) to obtain
where
It is easy to see that (3.7) and (3.8) are equivalent to (3.4) and (3.5). Due to τ being small enough and E being an identity matrix, the matrix K is invertible. Additionally, since Uk(k=0,1,⋯,n−1) is known, after multiple iterations, (3.7) and (3.8) have a unique solution. □
Remark 3.1. Since we introduce the auxiliary variable q=ˆuxx to transform (2.2) into a first-order system (2.5) and (2.6), according to [34,35], the mixed finite element scheme (2.21) and (2.22) do not need to satisfy the LBB condition. In [36], the LBB condition is a condition for the problem to be well posed. From this perspective, typically satisfying the LBB condition is to obtain the existence and uniqueness of a solution. Although the mixed finite element scheme in this article does not need to satisfy the LBB condition, it still satisfies the existence and uniqueness of a solution.
4.
Stability
Lemma 4.1. [12,14] For ∀Um∈Vh, satisfying Um=0(m<0), we have
where
and
Lemma 4.2. [27] } {For any vector (v^0, v^1, \cdots, v^{n-1})\in \mathbb{R}^n , defining \{\omega_k^{(\gamma)}\}_{k = 0}^{\infty}\; (0 < \gamma < 1) be a sequence of coefficients of the generating function \omega^{(\gamma)}(\xi) in (2.12) and 0\le\theta\le\min\{\gamma, \frac{1}{2}\} , we have
Theorem 4.1. Let u^n_h = U^n+\bar{u}_h^0 , where \bar{u}_h^0 is an approximation of u_0 , the following stability of the fully discrete scheme (2.21) and (2.22) holds:
where C is a positive constant independent of h and \tau .
Proof. Taking V = U^{n-\theta} , W = \Psi_{\tau}^{\beta, n}Q , (2.21) and (2.22) can be written as
and
Adding (4.2) and (4.3), we have
Using Lemma 4.1, we obtain
Multiply (4.5) by 4\tau and sum it with respect to n from 1 to L to get
By the Hölder inequality and Young inequality, the three terms on the right-hand side of (4.6) can be expanded to
where we use the bounded condition \|c_f(U^{n-\theta})\|_{\infty}\leq C .
Substituting (4.7)–(4.9) into (4.6), we arrive at
In what follows, using Lemmas 4.1 and 4.2 and the Gronwall inequality, we have
Since U^0 = 0 , we obtain
Noting that U^L = u^L_h-\bar{u}_h^0 and using the triangle inequality, the conclusion of this theorem is derived. □
5.
Numerical results
In this section, we present numerical simulation results for both smooth and nonsmooth solutions to verify the effectiveness of the numerical scheme. Next, we set the nonlinear term f(u) = u^2 , the spatial domain \Omega = (0, 1) , and the time interval J = (0, 1] .
Example 5.1 The exact solution is u(x, t) = t^{2}\sin(2\pi x) satisfying u(x, 0) = 0 , and the known source function g(x, t) is given by
In Table 1, fixing \tau = 1/1000 and choosing h = 1/10, 1/20, 1/40, 1/80 , we provide the L^2 -errors and the spatial convergence rates for u and q with different parameters \alpha , \beta , and \theta , where \theta\le\min\{\alpha, \beta, \frac{1}{2}\} . Similarly, in Table 2, taking h = 1/1000 , we calculate the L^2 -errors and the time convergence rates with \tau = 1/10, 1/20, 1/40, 1/80 . From Tables 1 and 2, one can see that the convergence rates in both space and time are close to 2 when the exact solution is smooth. In Table 3, if \theta > \min\{\alpha, \beta, \frac{1}{2}\} , the convergence accuracy will be unstable, which verifies the range of \theta values from a numerical perspective. To observe the effect of numerical simulation more clearly, we provide the comparison images between numerical solutions and exact solutions. In Figure 1, we show distinct comparison images of the numerical solutions of u_h and q_h and the exact solutions of u and q with \tau = 1/1000 , h = 1/80 , \alpha = 0.2 , \beta = 0.8 , and \theta = 0.2 .
Example 5.2 In this example, we consider the case where the nonsmooth solution is taken as u = (t^{\alpha+\beta}+t^3) \sin(2\pi x) , and the known source term g(x, t) is
Table 4 presents the L^2 -errors and the spatial convergence rates of u and q before and after adding the starting parts with h = 1/10, 1/20, 1/40, 1/80 , \tau = 1/2000 , where Erroro denotes the error before adding the starting parts and Errorc denotes the error after adding the starting parts.
The spatial convergence rate is almost unaffected before and after correction, based on a comparison of the data in Table 5. In Tables 6 and 7, we present the L^2 -errors and the time convergence rates of u and q before and after adding the starting parts. Without the addition of the starting parts, the time convergence rates are unstable and cannot reach the second-order convergence results computed by the generalized BDF2- \theta . After adding the starting parts, the time convergence rates keep around 2 , indicating that the starting part plays a major role in correcting the time convergence rates.
In Figure 2, we obtain the comparison images between the numerical solution and the exact solution with \tau = 1/1000 , h = 1/80 , \alpha = 0.9 , \beta = 0.2 , and \theta = 0.2 . In Figures 3 and 4, we present the space and time convergence rate images of u_h and q_h under different parameters \alpha , \beta , and \theta . From Figure 4, one can see that the corrected scheme with the starting parts can effectively restore the second-order convergence rate for the nonsmooth problem.
Example 5.3. To better investigate the effect of changes of two fractional parameters \alpha and \beta on the convergence rates, we introduce the numerical example with two nonsmooth terms. Here, we take the nonsmooth solution u with
and the known source term
In Table 8, we provide the errors of \|u_h-u\| and \|q_h-q\| and the spatial convergence rates under different parameters, which indicate that the corrected term hardly affects the spatial convergence rate.
In Tables 9–11, fixing \tau = 1/4000 , choosing h = 1/10, 1/20, 1/40, 1/80 , and changing parameters \alpha , \beta , and \theta , we provide the L^2 -errors and the time convergence rates for u and q based on the corrected scheme and uncorrected scheme. The impact of different fractional parameters on the time convergence rates of nonsmooth problems is evident from Tables 9–11. Furthermore, one can see that the corrected scheme with the starting part can effectively restore the second-order convergence rate.
To further validate the performance of the parameter \theta in numerical simulations with nonsmooth solutions, we provide the computing data in Table 12, from which one can see that the parameter \theta still needs to satisfy \theta\le \min\{\alpha, \beta, \frac{1}{2}\} , whether before or after correction. Notably, when \theta is negative, as long as it is not much less than 0, we can still obtain second-order convergence accuracy.
The time convergence rates of u and q are compared before and after correction with different parameters \alpha , \beta , and \theta in Figure 5, where the slope of the line segment indicates the convergence rate. The slope of each line segment in the corrected images is the same regardless of the parameters chosen, indicating that the introduction of the starting part has a significant effect on the time convergence rates for the case with nonsmooth solutions.
6.
Conclusions
In this article, the spatial mixed finite element method with the generalized BDF2- \theta for solving the time-fractional generalized Rosenau-RLW-Burgers equation was presented. Detailed proofs of stability were shown. The numerical scheme's effectiveness and feasibility were verified by conducting numerical examples that included both smooth and nonsmooth solutions. The numerical examples with good regularity indicated that our algorithm with changed parameters \alpha , \beta , and \theta can maintain second-order convergence in time. Especially, the nonsmooth examples demonstrated that adding the correction term could effectively solve the problem of reduced order caused by weak singularity.
Author contributions
N. Yang: Writing–original draft, Formal analysis, Software; Y. Liu: Methodology, Validation, Formal analysis, Funding acquisition, Supervision, Writing–review & editing. 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 supported by the National Natural Science Foundation of China (12061053), Young Innovative Talents Project of Grassland Talents Project and Program for Innovative Research Team in Universities of Inner Mongolia Autonomous Region (NMGIRT2413).
Conflict of interest
The authors declare that they have no conflicts of interest.
DownLoad:
