Positivity-Preserving Numerical Method and Relaxed Control for Stochastic Susceptible–Infected–Vaccinated Epidemic Model with Markov Switching

Abstract

The stochastic susceptible–infected–vaccinated (SIV) epidemic model includes a nonlinear term, making it difficult to obtain analytical solutions. Thus, numerical approximation schemes are an important tool for predicting the dynamics of infectious diseases and establishing optimal control strategies. However, the convergence rate of the existing numerical methods [e.g., Euler-Maruyama (EM) and truncated EM scheme] is only 1/2 order of the time step $Δ t$ . This article describes the construction of a logarithmic truncated EM scheme that achieves order-1 convergence and ensures positive numerical solutions of the stochastic SIV epidemic model. The existence of an invariant measure is proved for the stochastic SIV epidemic model with Markov switching. In addition, relaxed controls for the stochastic SIV epidemic model are investigated by using the Markov chain approximation method. It is demonstrated that the approximation schemes converge to the optimal strategy as the mesh size goes to zero. Finally, the results of numerical examples are presented to illustrate the theoretical results derived in this article.

1. INTRODUCTION

Medical research has shown that some diseases, particularly tuberculosis, hepatitis B, hand, foot, and mouth diseases, and COVID-19, do not produce permanent immunity, and so subjects can become reinfected after vaccination (Anderson et al., 1985; Kaufmann and Mcmichael, 2005; Keeling, 2003; Tildesley et al., 2006; Wei and Chen, 2014). Such diseases can be simulated by susceptible–infected–vaccinated (SIV) epidemic models (Arino et al., 2003; Safan and Rihan, 2014). As SIV epidemic models predict the dynamics of an epidemic, they enable doctors to provide better treatment (Deng and Liu, 2020) and have thus attracted considerable attention (Arino et al., 2003; Lin et al., 2014; Liu et al., 2018; Magpantay et al., 2014; Zhao and Jiang, 2015).

To date, there have been several studies on the deterministic SIV epidemic model (Arino et al., 2003; Magpantay et al., 2014). Considering that the spread of an epidemic is always affected by random noise, deterministic models cannot capture all the inherent fluctuations (Liu et al., 2018). Thus, stochastic SIV epidemic models incorporating white noise have been established. For example, Zhao and Jiang (2015) studied the ergodicity of a stochastic SIV epidemic model with vaccination, whereas Lin et al. (2014) proved the stationary distribution of a stochastic SIV epidemic model.

As well as being affected by temperature, humidity, climate, and other factors, epidemic models are also subject to telegraph noise, which may lead to the parameters of the system switching from one environmental state to another (Li et al., 2017; Mao and Yuan, 2006; Zhang et al., 2016). In most cases, the switching between environmental regimes is memoryless, and the waiting time for the next switching event follows an exponential distribution (Li et al., 2017). Many studies (Mao and Yuan, 2006; Li et al., 2017; Zhang et al., 2016) have shown that Markov switching can describe this phenomenon. For instance, Zhang et al. (2016) analyzed the following system:

Let $S (t)$ , $I (t)$ , and $V (t)$ denote the density of SIV individuals at time t, respectively. Let $p (Λ_{t})$ represents the proportion of the population that is vaccinated immediately after birth, and let $u_{1} (Λ_{t})$ denotes the vaccination rate; $β (Λ_{t})$ , $μ (Λ_{t})$ , and $b (Λ_{t})$ denote the per capita transmission coefficient, per capita natural death rate, and per capita natural birth rate, respectively. $α (Λ_{t})$ and $σ (Λ_{t})$ denote the per capita recovery rate and the white noise intensity, respectively.

For model (1), considerable attention has focused on the long-time behavior; see Li et al. (2017) and Zhang et al. (2016) and the references therein. Note that, to better predict the number of infected individuals and design controls for infectious diseases, it is necessary to obtain exact solutions of model (1). As this model includes a nonlinear drift term, it is fairly difficult to obtain an exact solution. Therefore, the construction of accurate and efficient numerical methods for approximating models of this form is an urgent topic. As we know, the explicit Euler-Maruyama (EM) numerical method has a simple algebraic structure and is computationally inexpensive (Berres and Ruiz-Baier, 2011; Mao and Yuan, 2006).

For most stochastic differential equations, the EM scheme is a valuable tool (Lamba, 2007; Mao and Yuan, 2006; Yuan and Mao, 2004). However, the positivity of the numerical solution cannot be guaranteed using the classical EM method for model (1), which satisfies the local Lipschitz condition and has a complex structure. The convergence rate is half the order of the time step for most existing numerical methods (Lamba, 2007; Mao and Yuan, 2006). Therefore, the study of positive numerical solutions with high convergence rates for the stochastic SIV epidemic model with Markov switching remains incomplete.

The novel coronavirus epidemic outbreak has had a great impact on people's lives and on the global economy. The control of infectious diseases has aroused widespread concern among scholars, such as isolation, inoculation, treatment, and social distancing (Telenti et al., 2021; Worthmann et al., 2021). Therefore, studying the optimal control of model (1) is of great significance. Most previous studies have applied Pontryagin's maximum principle and dynamic programming to solve optimal control problems. For example, Buonomo et al. (2014) studied a stochastic SIV epidemic model with regimen switching in which the aim was to minimize the average costs, whereas another study investigated the near-optimal control of the SIV epidemic model with Markov switching (Kar and Batabyal, 2011).

The above references (Buonomo et al., 2014; Kar and Batabyal, 2011) investigated the optimal control of an epidemic model under an admissible set of convex solutions (Buonomo et al., 2014; Laarabi et al., 2011; Kar and Batabyal, 2011). In the real world, convexity is not always satisfied, and the optimal control may not be determined by Pontryagin's maximum principle and dynamic programming. To better control the number of infected individuals using model (1), it is crucial to establish optimal control theory and present more challenging analysis.

This article presents a logarithmic truncated EM (log-TEM), scheme combined with a log transform (Lamba, 2007) and truncated EM numerical methods (Mao and Yuan, 2006; Yang and Huang, 2021). Our key idea is to transform model (1) into another low-dimensional system using the Itô formula and log transformation. We can preserve the positivity of numerical solutions for model (1) by transforming back to the original space. Furthermore, we propose relaxed controls to study the existence of optimal control with nonconvex assumptions. The main innovations of our article are as follows:

The existing literature on stochastic SIV epidemic models primarily focuses on dynamic behavior; however, numerical solutions are rarely studied. The convergence accuracy of the positivity of numerical solutions to stochastic epidemic models (Berres and Ruiz-Baier, 2011; Li et al., 2017) is only 1/2 order of the time step. In this article, we present the log-TEM scheme, which provides approximate solutions that achieve order-1 of the time step and ensure the positivity of the numerical solution for model (1).

For the study of near-optimal control theory of stochastic differential equations, existing literatures (Zong and Zhang, 2023; Zong et al., 2019) assume that the admissible control set U is convex. This article is based on a non-convex U, and a relaxed control strategy is provided for the stochastic SIV epidemic model with Markov switching. At the same time, a numerical approximation scheme, which is converging to the optimal control strategy, has been constructed. Our results can be considered as extensions of existing literatures on optimal control of the stochastic SIV epidemic model.

The remainder of this article is organized as follows. Some necessary preliminary knowledge and a novel SIV epidemic model are introduced in Section 2. The invariant measure of the numerical solution to our SIV epidemic model is given in Section 3, before the convergence of the log-TEM numerical method is proved in Section 4. In Section 5, the numerical approximation scheme for relaxed control by the Markov chain approximation method is established. Numerical examples are given to illustrate our theoretical results in Section 6. Finally, some concluding remarks are presented in Section 7.

2. PRELIMINARIES

We assume the Markov chain $Λ (\cdot)$ and the scalar standard Brownian motion $B (\cdot)$ are independent and defined on a complete filtered probability space $(Ω, ℱ, {ℱ_{t}}_{t \geq 0}, ℙ)$ with a filtration ${ℱ_{t}}_{t \geq 0}$ satisfying the usual conditions (i.e., it is increasing and right continuous while $ℱ_{0}$ contains all $ℙ$ -null sets). If $0 \leq t < s < \infty$ , we have $ℱ_{t} \subset ℱ_{s} \subset ℱ$ . Let $ℰ$ denotes the expectation corresponding to $ℙ$ and $| \cdot |$ denotes the Euclidean norm in R. Let $Λ_{t},$ $t > 0$ , be a right-continuous Markov chain on the probability space taking values in a finite state $S = {1, 2, \dots, N}$ for some positive integer $N < \infty$ . The generator of ${Λ_{t}}_{t > 0}$ is specified by such that for a sufficiently small $Δ$ , $ℙ (Λ_{t + Δ} = j | Λ_{t} = i) = (\begin{matrix} q_{i j} Δ + o (Δ), & i \neq j, \\ 1 + q_{i i} Δ + o (Δ), & i = j, \end{matrix}$ (2)

where $Δ > 0$ , $o (Δ)$ satisfies ${lim}_{Δ \to 0} \frac{o (Δ)}{Δ} = 0$ , here $q_{i j}$ is the transition rate from i to j satisfying $q_{i i} = - \sum_{i \neq j} q_{i j}$ . The discrete Markov chain ${Λ_{k}, k = 0, 1, \dots,}$ can be depicted as follows. Given a stepsize $Δ > 0$ and let $t_{k} = k Δ, Λ_{k} = Λ (t_{k})$ for $k \geq 0$ , and one-step transition probability matrix (Mao and Yuan, 2006).

Based on the idea of the epidemic control policy, the SIV epidemic model that incorporates the treatment of the infected population as a control measure is extended. As in Anderson et al. (1985), Wei and Chen (2014), and Keeling (2003), we assume that the treatment control can containment of disease outbreaks. That is, the treatment control $u_{2} (\cdot)$ is incorporated into model (1), and $u_{2} (\cdot)$ takes a value in the control set $U$ . Then, the following dynamic model is obtained:

where $h (I (t))$ denotes the treatment function. By standardizing the equation (3), it is easy to get the positive invariant set $D = {(S (t), I (t), V (t), Λ_{t}) | S (t) \geq 0, I (t) \geq 0, V (t) \geq 0, S (t) + I (t) + V (t) = 1} .$ (4)

To proceed further, our work bases on the following assumption.

Assumption. For any $ϕ, φ \in C ([0, T]; R)$ , the functional $h (\cdot) : C ([0, T]; R) \to R$ satisfies the following conditions $\begin{matrix} (i) | h (ϕ) | \leq c (1 + | ϕ |), \\ (i i) | h (ϕ) - h (φ) | \leq K | ϕ - φ |, \\ (i i i) h (ϕ) > 0, f o r a n y ϕ > 0 a . e . o n [0, T], \end{matrix}$

where c and K denote different positive constants.

Remark 1. Conditions (i)–(ii) guarantee the existence and uniqueness of the solution for system (3). The biological background significance of condition (iii) means that treatment control always has a positive effect on the disease. Conditions (i)–(iii) can be verified by various types of treatment functions, for example, the saturation treatment functional (Wang, 2012) $h (I (t)) : = \frac{m I (t)}{1 + η I^{q} (t)}, m, η, q \in R_{+},$

and linear treatment function $h (I (t)) = I (t)$ (Zaman et al., 2007) can also satisfied the Assumption.

3. INVARIANT MEASURE

From application and biological perspectives, the long-time behavior of an epidemic is one of its most important properties. A fundamental question in studying long-time behavior is whether an invariant measure exists. In this section, the existence and uniqueness of an invariant measure for model (3) is discussed. Similar to the study by Gray et al. (2012), one can define the basic reproduction number of the stochastic SIV epidemic model with Markov switching (3) as $R_{0}^{s} = \frac{1}{μ (Λ_{t}) + α (Λ_{t})} [\frac{β (Λ_{t}) [(1 - p (Λ_{t})) b (Λ_{t})]}{μ (Λ_{t}) + ψ (Λ_{t})} - \frac{σ^{2} (Λ_{t})}{2} {(\frac{(1 - p (Λ_{t})) b (Λ_{t})}{μ (Λ_{t}) + ψ (Λ_{t})})}^{2}] .$

Theorem 1. For the SIV epidemic model of (3) (i.e., when $σ (Λ_{t}) = 0$ ), if $R_{0}^{s} \leq 1,$ (5)

the disease-free equilibrium E₀ is globally stable (Beretta and Takeuchi, 1995) $E_{0} = (\frac{(1 - p (Λ_{t})) b (Λ_{t})}{μ (Λ_{t}) + u_{1} (Λ_{t})}, 0, \frac{u_{1} (Λ_{t}) + p (Λ_{t}) b (Λ_{t})}{μ (Λ_{t}) + u_{1} (Λ_{t})}) .$

Moreover, if $R_{0}^{s} > 1$ the endemic equilibrium of SIV epidemic model $E^{*} = (S^{*}, I^{*}, V^{*})$ is globally stable (Enatsu et al., 2012; Zhen et al., 2006), where $\begin{matrix} S^{*} = \frac{(μ (Λ_{t}) + α (Λ_{t}) + u_{2} (Λ_{t}))}{β (Λ_{t})}, \\ I^{*} = \frac{β (Λ_{t}) (1 - p (Λ_{t})) b (Λ_{t}) - (μ (Λ_{t}) + u_{1} (Λ_{t})) (μ (Λ_{t}) + α (Λ_{t}) + u_{2} (Λ_{t}))}{(μ (Λ_{t}) + u_{2} (Λ_{t})) β (Λ_{t})}, \\ V^{*} = \frac{β (Λ_{t}) b (Λ_{t}) (p (Λ_{t}) μ (Λ_{t}) + u_{2} (Λ_{t})) + (μ (Λ_{t}) + α (Λ_{t}) + u_{2} (Λ_{t})) (u_{1} (Λ_{t}) - u_{2} (Λ_{t})) μ (Λ_{t})}{β (Λ_{t}) μ (Λ_{t}) (μ (Λ_{t}) + u_{2} (Λ_{t}))} . \end{matrix}$ (6)

Proposition 1. There exists a unique positive solution $X^{x} (t) = (S (t), I (t), V (t))$ to model (3) with initial condition $X_{0} = x : = (S_{0}, I_{0}, V_{0}) \in ℱ_{0}$ . Then, the function $P_{t} (x, Λ_{t}) : = P^{x} (X_{t}^{x} \in Λ_{t}), t \geq 0, (x, Λ_{t}) \in R^{3} \times S,$

is Markovian and Feller continuous.

Definition 1. Let ${P_{t}}_{t \geq 0}$ be a Markovian semigroup. A probability measure $μ$ is called invariant measure of ${P_{t}}_{t \geq 0}$ if $\int_{E} P_{t} (y, Λ_{t}) μ (d y) = μ (Λ_{t}), \forall t \geq 0 .$

Define a set of probability measures ${Q_{T}^{x}}_{T \geq 0}$ by $Q_{T}^{x} (Λ_{t}) : = \frac{1}{T} \int_{0}^{T} P_{t} (x, Λ_{t}) d t .$

Theorem 2. If $R_{0}^{s} > 1, M > 0$ and $σ (Λ_{t})$ is small enough, then model (3) has a unique invariant measure and it is ergodic.

The proof of Theorem 2 is given in Section 8.1.

Remark 2. There exists a bounded open domain $U \in R^{3}$ with smooth boundary $\partial U$ . For $σ (Λ_{t})$ sufficiently small, we define

the distance $ρ (D_{σ (Λ_{t})}, \partial R_{+}^{3}) > 0$ . Then, one can take U as any neighborhood of the region $D_{σ (Λ_{t})}$ such that $Ū \subset R_{+}^{3}$ , where $Ū$ is the closure of U. Hence, for some $κ > 0, L (V_{1} + V_{2} + V_{3}) < - κ$ for any $(S, I, V) \in R_{+}^{3} ∖ U$ . By the proposition 2.1 (Fu, 2019), one can derive that the model (3) has a unique invariant measure and it is ergodic.

4. LOGARITHMIC TRUNCATED EM SCHEME

Give a stepsize $Δ > 0$ and $t_{k} = k Δ, Λ_{k} = Λ (t_{k})$ for $k \geq 0$ , and transition probability matrix In general, the data of vaccinated and infected individuals are available directly from the Centers for Disease Control, whereas the data regarding susceptible individuals are difficult to be obtained. It is obvious that $S (t) = 1 - I (t) - V (t)$ , the transformed model (3) is described by

then $\begin{matrix} d I (t) & = [M (Λ_{t}) I (t) - β (Λ_{t}) I^{2} (t) - β (Λ_{t}) V (t) I (t) - h (I (t)) u_{2} (Λ_{t})] d t \\ + σ (Λ_{t}) (1 - I (t) - V (t)) I (t) d B (t) \\ =Δ f (I (t), V (t)) d t + σ (Λ_{t}) (1 - I (t)) I (t) d B (t) - σ (Λ_{t}) V (t) I (t) d B (t), \end{matrix}$ (8)

where $M (Λ_{t}) = β (Λ_{t}) - (μ (Λ_{t}) + α (Λ_{t}))$ . One of the most essential difficulties in the study of positivity and the improvement of calculation accuracy is to select the appropriate logarithmic transform. Motivated by Lamperti transform, we give the following logarithmic transformation: $y (t) = log (I (t)) - log (1 - I (t)) = log \frac{I (t)}{1 - I (t)}, z (t) = log (V (t)),$

or equivalently $I (t) = \frac{e^{y (t)}}{1 + e^{y (t)}}, \frac{I (t)}{1 - I (t)} = e^{y (t)}, V (t) = e^{z (t)} .$

By the Itô formula, one sees that

Applying the Itô formula, we get $\begin{matrix} d z (t) & = \frac{1}{V (t)} (p (Λ_{t}) b (Λ_{t}) + u_{1} (Λ_{t}) (1 - I (t) - V (t)) - μ (Λ_{t}) V (t) + h (I (t)) u_{2} (Λ_{t})) d t \\ = e^{- z} (p (Λ_{t}) b (Λ_{t}) + u_{1} (Λ_{t})) + u_{1} (Λ_{t}) - μ (Λ_{t}) - e^{- z} u_{1} (Λ_{t}) \frac{e^{y}}{1 + e^{y}} + e^{- z} h (\frac{e^{y}}{1 + e^{y}}) u_{2} (Λ_{t}) \\ =Δ G (y (t), z (t)) d t, \end{matrix}$ (10)

with $z (0) = log (V_{0})$ . Then, by EM scheme of (9) and (10), we deduce that

where $Δ B_{k} = B (t_{k + 1}) - B (t_{k})$ . Transforming that, that is $I_{k}^{E, Δ} = \frac{e^{Y_{k}^{E}}}{1 + e^{Y_{k}^{E}}}, \frac{I_{k}^{E, Δ}}{1 - I_{k}^{E, Δ}} = e^{Y_{k}^{E}}, V_{k}^{E, Δ} = e^{Z_{k}^{E}},$ (12)

gives a strictly positive approximation of the system (7).

4.1. Positivity of numerical solutions

To define appropriate numerical solutions, an explicit scheme for approximating the exact solution of the stochastic SIV epidemic model (3) is proposed. For any given $Δ \in (0, 1]$ , each $Λ_{k} \in S$ , defining a truncation mapping $π_{Δ} : R \to R$ by $π_{Δ} (x) = (| x | \land log (K Δ^{- θ})) \frac{x}{| x |}, θ \in (0, \frac{1}{2}] .$

Then, we define a truncated EM scheme $\{\begin{matrix} Ȳ_{k + 1} = Y_{k} + F (Y_{k}, Z_{k}) Δ + σ (Y_{k}, Z_{k}) Δ B_{k}, {\bar{Z}}_{k + 1} = Z_{k} + G (Y_{k}, Z_{k}) Δ, \\ Y_{k + 1} = π_{Δ} (Ȳ_{k + 1}), Z_{k + 1} = π_{Δ} ({\bar{Z}}_{k + 1}), \end{matrix}$ (13)

where $K > (1 + e^{Y_{0}}) \lor (1 + e^{Z_{0}})$ . One has the convention $\frac{x}{| x |} = 0$ when $x = 0 \in R$ . Transforming back, that is, $Ī_{k}^{Δ} = \frac{e^{Ȳ_{k}}}{1 + e^{Ȳ_{k}}}, I_{k}^{Δ} = \frac{e^{Y_{k}}}{1 + e^{Y_{k}}}, {\bar{V}}_{k}^{E, Δ} = e^{{\bar{Z}}_{k}^{E}}, V_{k}^{E, Δ} = e^{Z_{k}^{E}},$ (14)

where $Ī_{k}^{Δ}$ and ${\bar{V}}_{k}^{Δ}$ are the strictly positive approximation for Equations (8) and (10), respectively.

Remark 3. In this article, we are devoted to studying the existence and uniqueness of invariant measure for analytic solutions of stochastic SIV epidemic systems, and the numerical solution invariant measure tends to the analytical solution invariant measure as $Δ \to 0$ in the Wasserstein metric.

Using the above proposition, we will obtain the existence and uniqueness of positive solutions to system (8). In particular, the positivity of the solution ensures the practical constraints of epidemic models.

Theorem 3. For $F \in C^{2} ([0, N]; R)$ , the model (8) has a unique strong solution $(I (t), V (t))$ , which satisfies $P (I (t) \lor V (t) \in (0, N), t \in [0, T]) = 1 .$

Proof. First, picking $m_{0} > 0$ sufficiently large such that $\frac{1}{m_{0}} < I (0) + V (0) < N - \frac{1}{m_{0}}$ . For each integer $m \geq m_{0}$ , we define a stopping time: $τ_{m} = inf {t \in [0, T] : I (t) + V (t) \notin (\frac{1}{m}, N - \frac{1}{m})} .$ (15)

By similar method with Gray et al. (2012), we can proof the theorem.

Lemma 1. For any $p \geq 0$ , we have

where $C_{p} = (I_{0}^{- p} \lor V_{0}^{- p}) exp (p K_{0} T + \frac{p (p + 1)}{2} σ^{2} (Λ_{t}) N^{2} T),$ (16)

and $K_{0} = | β (Λ_{t}) - (μ (Λ_{t}) + α (Λ_{t})) | + c u_{2} (Λ_{t}) + 3 β (Λ_{t}) N .$

Proof. For any $y \in (0, N)$ , we compute

Let $τ_{m}$ be defined by (15). By the Itô formula, we then derive that $ℰ [Y_{t \land τ_{m}}^{- p}] = Y_{0}^{- p} + ℰ \int_{0}^{t \land τ_{m}} L (Y_{s}^{- p}) d s \leq Y_{0}^{- p} + (p K_{0} + \frac{p (p + 1)}{2} σ^{2} (Λ_{t}) N^{2}) \int_{0}^{t} ℰ [Y_{s \land τ_{m}}^{- p}] d s .$

Letting $m \to \infty$ , the Gronwall inequality implies ${sup}_{t \in [0, T]} ℰ [Y_{t}^{- p}] \leq C_{p} .$

For any $y \in (0, N)$ , one can see that

By similar method as above, can be obtained.

The proof is completed.

Remark 4. Theorem 3 proves that model (8) has a unique positive strong solution. The following lemmas prove the boundedness of its inverse moment for the SIV epidemic model (7) with Markov switching. The inverse moment plays an important role in analyzing the convergence rate of the numerical scheme.

Lemma 2. For any $p \in R$ , the transformed (9) has the following property ${sup}_{t \in [0, T]} ℰ [e^{p y (t)} + e^{p z (t)}] \leq N^{| p |} C_{| p |},$ (18)

where C_p is given in Equation (16).

Proof. Note that $y (t) = ln \frac{Y_{t}}{N - Y_{t}}$ , $\forall p \geq 0$ , we can show that

Similarly, for any $p < 0$ , we obtain ${sup}_{t \in [0, T]} ℰ [e^{p y (t)}] \leq N^{| p |} C_{| p |} .$

Noticing that $z (t) = log (V (t))$ , by Theorem 3, one has ${sup}_{t \in [0, T]} ℰ [e^{p z (t)}] = {sup}_{t \in [0, T]} ℰ [V^{p} (t)] \leq N^{p} C_{p} .$

Therefore, ${sup}_{t \in [0, T]} ℰ [e^{p y (t)} + e^{p z (t)}] \leq N^{| p |} C_{| p |} .$

The proof is completed.

Remark 5. Lemma 2 concerns the exponential integrability property of $y (t), t \in [0, T]$ , which will be used to analyze the convergence of the numerical scheme.

4.2. Convergence

This section derives the rate of convergence for the log-TEM explicit scheme. To prove the main results, two lemmas are first presented. For stochastic differential equations with Markov switching (Mao and Yuan, 2006), it is obvious that the EM scheme is an approximation of the model (9).

Lemma 3. For any $q > 0$ , there exists a constant $C > 0$ such that

where $Δ \in (0, 1]$ , $T > 0$ , and $⌊ T ∕ Δ ⌋$ represent the integer part of $T ∕ Δ$ .

Remark 6. According to the above lemma, the EM scheme of (9) is an approximation of the original model (9) with order-1 convergence.

Lemma 4. For any $p > 0$ , $T > 0$ , the truncated EM scheme (13) has the property that $\begin{matrix} {sup}_{Δ \in (0, 1)} {sup}_{0 \leq k \leq ⌊ T ∕ Δ ⌋} ℰ [e^{p Y_{k}}] \leq {sup}_{Δ \in (0, 1)} {sup}_{0 \leq k \leq ⌊ T ∕ Δ ⌋} ℰ [e^{p Ȳ_{k}}] \leq C, \end{matrix}$ (20)

and $\begin{matrix} {sup}_{Δ \in (0, 1)} {sup}_{0 \leq k \leq ⌊ T ∕ Δ ⌋} ℰ [e^{p Z_{k}}] \leq {sup}_{Δ \in (0, 1)} {sup}_{0 \leq k \leq ⌊ T ∕ Δ ⌋} ℰ [e^{p {\bar{Z}}_{k}}] \leq C, \end{matrix}$ (21)

where C is a constant dependents on p.

Proof. Letting $\hat{M} = {min}_{Λ_{k} \in S} M (Λ_{k}), \overset{ˇ}{M} = {max}_{Λ_{k} \in S} M (Λ_{k}),$

and $\hat{σ} = {min}_{Λ_{k} \in S} σ (Λ_{k}), \overset{ˇ}{σ} = {max}_{Λ_{k} \in S} σ (Λ_{k}),$

the following inequalities can be derived

Using similar arguments developed above, it is easy to prove ${sup}_{Δ \in (0, 1)} {sup}_{0 \leq k \leq ⌊ T ∕ Δ ⌋} ℰ [e^{p Z_{k}}] \leq {sup}_{Δ \in (0, 1)} {sup}_{0 \leq k \leq ⌊ T ∕ Δ ⌋} ℰ [e^{p {\bar{Z}}_{k}}] \leq C .$

The proof is completed.

Lemma 4 proves that the truncated EM scheme (13) is bounded. To investigate the rate of convergence of the log-TEM numerical solutions, the boundedness of scheme (13) is first proved.

Lemma 5. For any $L \geq 0$ and integer $m \geq 0$ , we have

where $Γ (\cdot)$ is the Gamma function.

Proof. For any integer $m \geq 0$ , and $Δ B_{k} \sim N (0, Δ)$ , we have

with

The proof is completed.

The following lemma proves the boundedness of the inverse moment of the truncated EM scheme (13).

Lemma 6. For any $p > 0$ , the logarithmic truncated EM scheme of (13) has $\begin{matrix} {sup}_{Δ \in (0, 1)} {sup}_{0 \leq k \leq ⌊ T ∕ Δ ⌋} ℰ [e^{- p Y_{k}}] \leq C a n d {sup}_{Δ \in (0, 1)} {sup}_{0 \leq k \leq ⌊ T ∕ Δ ⌋} ℰ [e^{- p Z_{k}}] \leq C, \end{matrix}$

where $Δ \in (0, 1]$ and $T > 0$ .

The proof of Lemma 6 is given in Section 8.2.

Let $τ_{L} = inf {k Δ \geq 0 : y (t_{k}) \geq e^{L}}, ρ_{1, Δ} = inf {k Δ \geq 0 : Ȳ_{Δ} (t_{k}) \geq log (K Δ^{- θ})}$ , by Equation (18) yields that $\begin{matrix} e^{p L} P {τ_{L} \leq T} \leq ℰ [e^{p y (T \land τ_{L})} 1_{{τ_{L} \leq T}}] + ℰ [e^{p y (T \land τ_{L})} 1_{{τ_{L} > T}}] = ℰ [e^{p y (t \land τ_{L})}] \leq C . \end{matrix}$

It is not difficult to get $\begin{matrix} P {τ_{L} \leq T} \leq C e^{- p L}, \forall p > 0, \end{matrix}$ (22)

where C is a positive constant. Moreover, by Lemma 5, we conclude that

Using similar method, let $τ_{M} = inf {k Δ \geq 0 : z (t_{k}) \geq Z (t_{k})}, ρ_{2, Δ} = inf {k Δ \geq 0 : {\bar{Z}}_{Δ} (t_{k}) \geq Z (t_{k})}$ , we can get

and implies $\begin{matrix} P {ρ_{1, Δ} \leq T} \leq C K^{- p} Δ^{p θ}, \forall p > 0, \end{matrix}$ (23)

where C is a positive constant independent of $Δ$ .

To show that numerical scheme (13) accurately reproduces the dynamical behaviors of exact solutions, the rate of convergence of the numerical solutions needs to be studied.

Theorem 4. For any $p > 0$ , $Δ \in (0, 1)$ , and $T > 0$ , the truncated EM scheme of Equation (13) has

Proof. Define ${\bar{θ}}_{Δ} = τ_{log (K Δ^{- θ})} \land ρ_{Δ}, Ω_{1} : = {ω : {\bar{θ}}_{Δ} > T}, ū_{k} = Y_{k} - y (t_{k}),$ for any $k Δ \in [0, T]$ , where $τ_{L}$ and $ρ_{Δ}$ are defined by Equations (22) and (23). For any $q > p$ , we have $\begin{matrix} ℰ | ū_{k} |^{p} = ℰ (| ū_{k} |^{p} 1_{Ω_{1}}) + ℰ (| ū_{k} |^{p} 1_{Ω_{1}^{c}}) \leq ℰ (| ū_{k} |^{p} 1_{Ω_{1}}) + \frac{p Δ^{p}}{q} ℰ (| ū_{k} |^{q}) + \frac{q - p}{q Δ^{\frac{p^{2}}{q - p}}} P (Ω_{1}^{c}) . \end{matrix}$

It follows Lemma 2, Equations (13) and (20), we yield that $\frac{p Δ^{p}}{q} ℰ (| ū_{k} |^{q}) \leq C Δ^{p} .$

It follows from Equations (22) and (23) that $\frac{q - p}{q Δ^{\frac{p^{2}}{q - p}}} P (Ω_{1}^{c}) \leq \frac{q - p}{q Δ^{\frac{p^{2}}{q - p}}} (P {τ_{K Δ^{- θ} \leq T}} + P {ρ_{Δ} \leq T}) \leq \frac{2 (q - p)}{q Δ^{\frac{p^{2}}{q - p}}} \frac{C}{K^{δ} Δ^{- δ θ}} \leq C Δ^{δ - \frac{p^{2}}{q - p}} \leq C Δ^{p},$

where $δ \geq \frac{q}{θ (q - 2)}$ . Then, we have $\begin{matrix} ℰ | ū_{k} |^{p} \leq ℰ (| ū_{k} | 1_{Ω_{1}}) + C Δ^{p} \leq {sup}_{0 \leq k \leq ⌊ \frac{t \land ρ_{Δ}}{Δ} ⌋} ℰ [| Y_{k} - y (t_{k}) |^{p}] + C Δ^{p} \\ = {sup}_{0 \leq k \leq ⌊ \frac{t \land ρ_{Δ}}{Δ} ⌋} ℰ [| Y_{k}^{E} - y (t_{k}) |^{p}] + C Δ^{p} \leq C Δ^{p} . \end{matrix}$

By a similar analysis, we also have $| Z (t_{k + 1}) - z (t_{k + 1}) | \leq | Z (t_{k}) - z (t_{k}) | + | G (Y (t_{k}), Z (t_{k})) - G (y (t_{k}), z (t_{k})) | Δ .$

Obviously, ${sup}_{0 \leq k \leq ⌊ \frac{t \land ρ_{Δ}}{Δ} ⌋} ℰ [| Z_{k}^{E} - z (t_{k}) |^{p}] + C Δ^{p} \leq C Δ^{p} .$

Hence, the desired conclusion holds.

Remark 7. By virtue of Theorem 4, the truncated EM scheme defined by Equation (13) has an order-1 convergence rate for the approximate stochastic SIV epidemic model. As a consequence, the same convergence order can be applied to the approximate stochastic SIV epidemic model with Markov switching (8).

Theorem 5. For any $0 < p < 1$ , $Δ \in (0, 1)$ , and $T > 0$ if the Assumption is satisfied, then there exists a constant $C > 0$ such that $\begin{matrix} {sup}_{0, \dots, ⌊ T ∕ Δ ⌋} ℰ [| I_{t_{k}} - I_{k}^{Δ} |^{p}] \leq C Δ^{p}, \end{matrix}$

and $\begin{matrix} {sup}_{0, \dots, ⌊ T ∕ Δ ⌋} ℰ [| V_{t_{k}} - V_{k}^{Δ} |^{p}] \leq C Δ^{p} . \end{matrix}$

Proof. By the log-transformation, one has $\begin{matrix} {sup}_{0, \dots, ⌊ T ∕ Δ ⌋} ℰ [| I_{t_{k}} - I_{k}^{Δ} |^{p}] & = N {sup}_{0, \dots, ⌊ T ∕ Δ ⌋} ℰ [| \frac{e^{y (t_{k})}}{1 + e^{y (t_{k})}} - \frac{e^{Y_{k}}}{1 + e^{Y_{k}}} |^{p}] \\ \leq N {sup}_{0, \dots, ⌊ T ∕ Δ ⌋} ℰ [| y_{t_{k}} - Y_{k}^{Δ} |^{p}] . \end{matrix}$

By Theorem 4, for any $Δ \in (0, 1]$ , we get $\begin{matrix} ℰ [{sup}_{0, \dots, ⌊ T ∕ Δ ⌋} | I_{t_{k}} - I_{k}^{Δ} |^{p}] \leq N {sup}_{0, \dots, ⌊ T ∕ Δ ⌋} ℰ [| y_{t_{k}} - Y_{k} |^{p}] \leq C Δ^{p} . \end{matrix}$

The proof is completed.

Remark 8. As a consequence, the order-1 convergence rate for the approximation of the SIV epidemic model (3) is obtained by $I_{k}^{E, Δ} = \frac{e^{Y_{k}^{E}}}{1 + e^{Y_{k}^{E}}}$ , where $Y_{k}^{E}$ is defined by Equation (11). Again, for the numerical solutions ${Z_{k}^{E}}_{k \geq 0}$ , it follows from Theorem 4 that $\begin{matrix} ℰ [{sup}_{0, \dots, ⌊ T ∕ Δ ⌋} | I_{t_{k}} - I_{k}^{E, Δ} |^{p}] \leq C Δ^{p}, \end{matrix}$

for any $Δ \in (0, 1]$ and $T > 0$ . A similar computation yields $\begin{matrix} ℰ [{sup}_{0, \dots, ⌊ T ∕ Δ ⌋} | V_{t_{k}} - V_{k}^{E, Δ} |^{p}] \leq C Δ^{p} . \end{matrix}$

Compared with previous studies (Mao and Yuan, 2006; Yuan and Mao, 2004), the convergence rate of the numerical method for solving the stochastic SIV epidemic model with Markov switching has been improved. Our results show that an order-1 convergence rate can be achieved using the log-TEM numerical method.

5. RELAXED CONTROLS

To give the process of relaxed control, let the process $(S (t), I (t), V (t), Λ_{t})$ starts with initial condition $(s, i, v, l)$ . Define $\begin{matrix} a_{1} (s, i, l, u) = (1 - p_{l}) b_{l} + α_{l} i - (μ_{l} + u_{1 l}) s - β_{l} s i, b_{1} (s, i, l, u) = σ_{l}^{2} s^{2} i^{2}, \\ a_{2} (s, i, v, l, u) = β_{l} s i - (μ_{l} + α_{l}) i - h (i) u_{2 l}, b_{2} (s, i, v, l, u) = σ_{l}^{2} s^{2} i^{2}, \\ a_{3} (s, i, v, l, u) = p_{l} b_{l} + u_{1 l} s - μ_{l} v + h (i) u_{2 l}, b_{3} (i, v, l, u) = 0, \end{matrix}$

where $u = (u_{1}, u_{2})$ . Let $ξ$ be a positive constant such that $ξ < 1$ , then $τ = inf {t \geq 0 : I (t) \leq ξ}$ and $I (0) < ξ$ .

A strategy $u (\cdot)$ lies in the collection of all admissible control $A_{i, l}$ , if $u (t)$ is $ℱ (t)$ -adapted and $u (t)$ takes values in a nonempty compact set U. The cost function $ℱ : [0, 1] \times [0, 1] \times S \times [0, \infty) \times [0, \infty) \to [0, \infty)$ , we define $J (s, i, v, l, u) : = ℰ_{s, i, v, l} \int_{0}^{τ} e^{- δ t} ℱ (S (t), I (t), l, u_{1}, u_{2}) d t,$

where $δ > 0$ is the discounting factor and $ℰ_{s, i, v, l}$ denotes the expectation with respect to the probability law when the process $(S (t), I (t), V (t), Λ_{t})$ starts with the initial condition $(s, i, v, l)$ . The goal is to minimize the cost functional and find an optimal control strategy $u^{*} (\cdot)$ such that $J (s, i, v, l, u^{*}) = V (s, i, v, l) : = {inf}_{u \in U} J (s, i, v, l, u) .$ (24)

The cost function is as follows: $\begin{matrix} ℱ (s, i, l, u) : = ℱ (s, i, l, u_{1}, u_{2}) = a_{0} + a_{1} i + a_{2} s u_{1}^{2} + a_{3} i u_{2}^{2} . \end{matrix}$

The practically meaning of the cost function is detailed in the study by Tran and Yin (2021). The Hamilton–Jacobi–Bellman (HJB) equation of the stochastic SIV epidemic model (3) is

for all $(s, i, v, l) \in [ξ, 1] \times [ξ, 1] \times [ξ, 1] \times S$ , with the boundary condition $V (ξ, ξ, ξ, l) = 0$ for any $l \in S$ .

Since the convexity is not satisfied, the optimal control may not be determined by Pontryagin's maximum principle and dynamic programming. To overcome the non-convexity of the control set U, one has to develop the Markov chain approximation method to study relaxed controls. The definition of relaxed control is given as follows.

Definition 2. (Relaxed control) (Tran and Yin, 2021) Let $ℬ (U \times [0, \infty))$ be the $σ$ -algebra of Borel subsets of $U \times [0, \infty)$ . An admissible relaxed control or simply a relaxed control $m (\cdot)$ is a measure on $ℬ (U \times [0, \infty))$ such that $\begin{matrix} m (U \times [0, t]) = t, \forall t \geq 0 . \end{matrix}$

Given a relaxed control $m (\cdot)$ , there is a probability measure $m_{t} (\cdot)$ defined on the $σ$ -algebra $ℬ (U)$ such that $m (d u d t) = m_{t} (d u) d t$ .

Remark 9. For a sequence of controls $u^{h} = {u_{n}^{h} : n \in ℤ_{+}},$ define a sequence of equivalent relaxed controls as follows. First, set $m_{t_{n}^{h}}^{h} (d u) = δ_{u_{n}^{h}} (d u)$ , where $δ_{u_{n}^{h}} (\cdot)$ is the probability measure concentrated at $u_{n}^{h}$ . Then, $m^{h} (\cdot)$ is defined by $m^{h} (d u d t) = m_{t} (d u) d t$ . That is, $m^{h} (B \times [0, t]) = \int_{0}^{t} (\int_{B} δ_{u^{h} (s)} (d u)) d s, B \in ℬ (U), t \geq 0 .$

Using the Markov chain approximation method, the controlled Markov chain in discrete time can be constructed to approximate the controlled switching diffusions. One of the advantages of this is that little regularity or prior information of the stochastic SIV epidemic system is needed.

5.1. Approximation algorithm

Let $h > 0$ be a discretization parameter for the continuous state variable. Define $T_{h} : = {s = k_{1} h, i = k_{2} h, v = k_{3} h : k_{1}, k_{2}, k_{3} \in ℤ_{+}, k_{1} + k_{2} + k_{3} \leq 1 ∕ h} \times S .$

The sequence u^h is said to be admissible if it satisfies the following conditions:

(a) u^h is $σ {S_{0}^{h}, \dots, S_{n}^{h}, I_{0}^{h}, \dots, I_{n}^{h}, V_{0}^{h}, \dots, V_{n}^{h}, α_{0}^{h}, \dots, α_{n}^{h}, u_{0}^{h}, \dots, C_{n - 1}^{h}} - a d a p t e d .$

(b) For any $(s, i, v, l) \in T_{h}$ ,

Let $A_{s, i, v, l}^{h}$ denotes the set of admissible control sequences. For each $(s, i, v, l, u_{1}, u_{2}) \in T_{h} \times U$ , one define $Δ t^{h} (s, i, v, l, u_{1}, u_{2})$ is the interpolation intervals. Define

For $(s, i, v, l) \in T_{h}$ and $u^{h} \in A_{i, l}^{h}$ , we define the cost functional as

with $η_{h} = inf {n \geq 0 : I_{n}^{h} \leq ξ} .$ The value function is $V^{h} (s, i, v, l) = {inf}_{u^{h} \in A_{i, l}^{h}} J^{h} (s, i, v, l, u^{h}) .$

The discrete HJB equation for the stochastic SIV epidemic model (3)

for any $(s, i, v, l) \in T_{h}$ . Define $Δ S_{n}^{h} = S_{n + 1}^{h} - S_{n}^{h}, Δ I_{n}^{h} = I_{n + 1}^{h} - I_{n}^{h}, Δ V_{n}^{h} = V_{n + 1}^{h} - V_{n}^{h},$

with ${S_{k}^{h}, I_{k}^{h}, V_{k}^{h}, α_{k}^{h}, u_{k}^{h}, k \leq n, S_{n}^{h} = s, I_{n}^{h} = i, V_{n}^{h} = v, α_{n}^{h} = l, u_{n}^{h} = u}$ . And $\begin{matrix} Δ t_{k}^{h} (s, i, v, l, u) = \frac{h^{2}}{Q_{h} (s, i, v, l, u)}, \\ Q_{h} (s, i, v, l, u) = a_{1} (s, i, v, l, u) + a_{2} (s, i, v, l, u) + a_{3} (s, i, v, l, u) + h | b_{1} (s, i, v, l, u) | \\ + h | b_{2} (s, i, v, l, u) | + h | b_{3} (s, i, v, l, u) | - h^{2} Λ_{l l} + h, \\ p^{h} ((s, i, v, l), (s, i, v, l') | u) = \frac{h^{2} Λ_{l l'}}{Q_{h} (s, i, v, l, u)}, l \neq l', \\ p^{h} ((s, i, v, l), (s, i, v, l) | u) = \frac{h}{Q_{h} (s, i, v, l, u)}, \\ p^{h} ((s, i, v, l), (s + h, i, v, l) | u) = \frac{a_{1} (s, i, v, l, u) ∕ 2 + {(b_{1} (s, i, v, l, u))}^{+} h}{Q_{h} (s, i, v, l, u)}, \\ p^{h} ((s, i, v, l), (s - h, i, v, l) | u) = \frac{a_{1} (s, i, v, l, u) ∕ 2 + {(b_{1} (s, i, v, l, u))}^{-} h}{Q_{h} (s, i, v, l, u)}, \\ p^{h} ((s, i, v, l), (s, i + h, v, l) | u) = \frac{a_{2} (s, i, v, l, u) ∕ 2 + {(b_{2} (s, i, v, l, u))}^{+} h}{Q_{h} (s, i, v, l, u)}, \\ p^{h} ((s, i, v, l), (s, i - h, v, l) | u) = \frac{a_{2} (s, i, v, l, u) ∕ 2 + {(b_{2} (s, i, v, l, u))}^{-} h}{Q_{h} (s, i, v, l, u)}, \\ p^{h} ((s, i, v, l), (s, i, v + h, l) | u) = \frac{a_{3} (s, i, v, l, u) ∕ 2 + {(b_{3} (s, i, v, l, u))}^{+} h}{Q_{h} (s, i, v, l, u)}, \\ p^{h} ((s, i, v, l), (s, i, v - h, l) | u) = \frac{a_{3} (s, i, v, l, u) ∕ 2 + {(b_{3} (s, i, v, l, u))}^{-} h}{Q_{h} (s, i, v, l, u)} . \end{matrix}$ (27)

Constructing a continuous-time interpolation of the approximating chain. Define $n^{h} (t) = max {n : t_{n}^{h} \leq t}, t \geq 0$ . Let

The piecewise constant interpolation processes is denoted by $S^{h} (\cdot), I^{h} (\cdot), V^{h} (\cdot), α^{h} (\cdot), B_{i}^{h} (\cdot), M_{i}^{h} (\cdot), u^{h} (\cdot),$ where $\begin{matrix} S^{h} (t) = S_{n^{h} (t)}^{h}, I^{h} (t) = I_{n^{h} (t)}^{h}, V^{h} (t) = V_{n^{h} (t)}^{h}, α^{h} (t) = α_{n^{h} (t)}^{h}, \\ u^{h} (t) = u_{n^{h} (t)}^{h}, B_{i}^{h} (t) = B_{i n^{h} (t)}^{h}, M_{i}^{h} (t) = M_{i n^{h} (t)}^{h}, i = 1, 2, 3, t \geq 0 . \end{matrix}$ (29)

Define $ℱ^{h} (t) = σ {S^{h} (s), I^{h} (s), V^{h} (s), α^{h} (s), u^{h} (s) : s \leq t}$ , we have $\begin{matrix} S^{h} (t) = s + B_{1}^{h} (t) + M_{1}^{h} (t), I^{h} (t) = i + B_{2}^{h} (t) + M_{2}^{h} (t), V^{h} (t) = v + B_{3}^{h} (t) + M_{3}^{h} (t) . \end{matrix}$

Recall that $Δ t_{k}^{h} = \frac{h^{2}}{Q_{h} (S_{k}^{h}, I_{k}^{h}, V_{k}^{h}, α_{k}^{h}, u_{k}^{h})}$ , it follows that

with ${ε_{i}^{h} (\cdot)}, i = 1, 2, 3$ being an $ℱ^{h} (t)$ -adapted process satisfying ${lim}_{h \to 0} {sup}_{t \in [0, T_{0}]} E | ε_{i}^{h} (t) | = 0$ for $T_{0} \in (0, \infty)$ . Define $τ_{h} = τ_{η_{h}}^{h}$ , we can rewrite Equation (25) as $J^{h} (s, i, v, l, u^{h}) = E \int_{0}^{τ_{h}} e^{- δ t} ℱ (I^{h} (t), α^{h} (t), u^{h} (t)) d t .$ (31)

Based on the definition of relaxed control, Equation (31) can also be expressed as: $J^{h} (s, i, v, l, u^{h}) = J^{h} (s, i, v, l, m^{h} (\cdot)) = E \int_{0}^{τ_{h}} e^{- δ t} ℱ (I^{h} (t), α^{h} (t), u) m_{t}^{h} (d u) d t .$ (32)

Note also that the value function defined in Equation (24) can be rewritten as $\begin{matrix} V (s, i, v, l) = inf {J (s, i, v, l, m (\cdot)) : m (\cdot) i s a n a d m i s s i b l e o p t i m a l c o n t r o l}, \end{matrix}$

where $J (s, i, v, l, m (\cdot)) : = E_{s, i, l} \int_{0}^{τ} e^{- δ t} ℱ (S (t), I (t), Λ_{t}, u) m_{t} (d u) d t .$

5.2. Convergence analysis for the relaxed controls

The goal of this section is to establish the convergence of the numerical algorithms for optimal control. First, a sequence for the continuous-time interpolation defined in Equations (28)–(29) is introduced. Second, the convergence of the value function is proved.

Theorem 6. Let the chain ${(S^{h} (\cdot), I^{h} (\cdot), V^{h} (\cdot), α_{n}^{h})}$ be constructed with the transition probabilities defined in Equation (27), $(S^{h} (\cdot), I^{h} (\cdot), V^{h} (\cdot), α^{h} (\cdot), B_{i}^{h} (\cdot), M_{i}^{h} (\cdot), τ_{h} (\cdot))$ be the continuous-time interpolation defined in Equations (28)–(29), ${u_{n}^{h}}$ be an admissible control and $m^{h} (\cdot)$ be the relaxed control representation of ${u_{n}^{h}}$ . Then, the following assertions hold.

(i) The sequence $\begin{matrix} H^{h} (\cdot) : = (S^{h} (\cdot), I^{h} (\cdot), V^{h} (\cdot), α^{h} (\cdot), B_{i}^{h} (\cdot), M_{i}^{h} (\cdot), τ_{h} (\cdot)) \end{matrix}$

is tight. As a result, has a weakly convergent subsequence with limit $\begin{matrix} H (\cdot) : = (S (\cdot), I (\cdot), V (\cdot), α (\cdot), B_{i} (\cdot), M_{i} (\cdot), τ (\cdot)) . \end{matrix}$

Moreover, $(S (\cdot), I (\cdot), V (\cdot), B_{i} (\cdot), M_{i} (\cdot))$ have continuous paths with probability one.

(ii) For any $t \geq 0$ , one has $B_{i} (t) = \int_{0}^{t} b_{i} (S (u), I (u), V (u), α (u), u) m_{u} (d c) d u .$

(iii) $M_{i} (\cdot)$ is a continuous $ℱ (t)$ -martingale with quadratic variation $\int_{0}^{t} a_{i} (S (u), I (u), V (u), α (u)) d u$ . Thus, there is an $ℱ (t)$ -standard Brownian motion $B (t)$ , in which we might have to augment the probability space such that $M_{i} (t) = \int_{0}^{t} σ (S (u), I (u), V (u), α (u)) d B (u) .$

(iv) The limit process satisfy $S (t) = s + B_{1} (t) + M_{1} (t), I (t) = i + B_{2} (t) + M_{2} (t), V (t) = v + B_{3} (t) + M_{3} (t), t \geq 0 .$

Proof. (i) By [(Song et al., 2006), lemma 3.4], we have ${α^{h} (\cdot)}$ is tightness. That is, for any $T_{0}, ρ \in (0, \infty)$ , it holds true that $\begin{matrix} ℰ_{t}^{h} | ζ^{h} (t + s) - ζ^{h} (t) |^{2} \leq ℰ_{t}^{h} γ (h, ρ), \forall s \in [0, ρ], t \leq T_{0} . \\ {lim}_{ρ \to 0} {lim sup}_{h \to 0} ℰ γ (h, ρ) = 0 . \end{matrix}$

By Song et al. (2006), one infer that $H^{h} (\cdot)$ is tight. Thus, $lim H^{h} (\cdot) \to^{w} H (\cdot) = (S (\cdot), I (\cdot), V (\cdot), α (\cdot), B_{i} (\cdot), M_{i} (\cdot), m (\cdot), τ) .$ Because of the jumps of $S^{h} (\cdot), I^{h} (\cdot), V^{h} (\cdot), B_{i}^{h} (\cdot)$ and $M_{i}^{h} (\cdot) \to 0$ as $h \to 0$ , then $S (\cdot), I (\cdot), V (\cdot), B_{i} (\cdot), M_{i} (\cdot)$ is continuous.

(ii) Follows from (i) and Equation (30), (ii) is obtained.

(iii) By the definition of $M_{i}^{h} (\cdot)$ , it is an ${ℱ^{h} (\cdot)}$ -martingale with quadratic variation process $\int_{0}^{t} a (S^{h} (u), I^{h} (u), V^{h} (u), α^{h} (u)) d u + ε_{3}^{h} (t),$

where ${ε_{3}^{h} (\cdot)}$ is an ${ℱ^{h} (t)}$ -adapted process satisfying ${lim}_{h \to 0} {sup}_{t \in [0, T_{0}]} ℰ | ε_{3}^{h} (t) | = 0$ for $T_{0} \in (0, \infty)$ . By the Burkholder–Davis–Gundy inequality, we have

where K is a positive constant and $h > 0$ . For any $ρ > 0$ , we get $\begin{matrix} ℰ_{t}^{h} (M_{i}^{h} (t + ρ) - M_{i}^{h} (t)) = 0, \\ ℰ_{t}^{h} (M_{i}^{h} (t + ρ) - M_{i}^{h} (t)) (M_{i}^{h} (t + ρ) - M_{i}^{h} (t))' = ℰ_{t}^{h} \int_{t}^{t + ρ} a (S^{h} (u), I^{h} (u), V^{h} (u), α^{h} (u)) d u + ε_{4}^{h} (ρ), \end{matrix}$ (33)

where $ℰ | ε_{4}^{h} (ρ) | \to 0$ as $h \to 0$ . Then in view of Equation (33), we get $ℰ ψ (H^{h} (t_{k}), k \leq r) [M_{i}^{h} (t + ρ) - M_{i}^{h} (t)] = 0,$

and

where $ℰ | ε_{5}^{h} (ρ) | \to 0$ as $h \to 0$ . Letting $h \to 0$ , we have $\begin{matrix} ℰ ψ (H (t_{k}), k \leq r) (M_{i} (t + ρ) - M_{i} (t)) = 0 . \end{matrix}$ (35)

Since $M_{i} (\cdot)$ has continuous paths with probability one, Equation (35) implies that $M_{i} (\cdot)$ is a continuous $ℱ (\cdot)$ -martingale. Moreover, Equation (34) implies that $\begin{matrix} ℰ ψ (H^{h} (t_{k}), k \leq r) ([M_{i} (t + ρ) - M_{i} (t)] (M_{i} (t + ρ) - M_{i} (t))' \\ - \int_{t}^{t + ρ} a (S (u), I (u), V (u), α (u)) d u) = 0 . \end{matrix}$ (36)

The proof is completed.

Consider the Markov chain with the transition probability in Equation (27). Using the relaxed control representation, its interpolation process can be represented by Equation (32). Following Song et al. (2006), let $V^{h} (s, i, v, l)$ and $V (s, i, v, l)$ be the value functions defined in Equations (26) and (24), respectively. Then, $V^{h} (s, i, v, l) \to V (s, i, v, l)$ as $h \to 0$ . $V^{h} (s, i, v, l)$ is the approximation of the value function in Equation (26) obtained by dynamic programming. The convergence of the value function $V^{h} (s, i, v, l)$ to $V (s, i, v, l)$ is proved in Equation (24).

6. NUMERICAL EXAMPLES

To illustrate the theoretical results obtained in previous sections, a number of examples and simulations are now considered. Let the initial values $S (0)$ , $I (0)$ , and $V (0)$ be 30,000, 100, and 500, respectively. Next, we present numerical simulation results for the SIV epidemic model assuming that the Markov chain $Λ_{t}$ is in the state space $S = {1, 2}$ with the generator $Γ = (\begin{matrix} - 3 & 3 \\ 7 & - 7 \end{matrix}) .$

The coefficients in each state are presented in Table 1. According to the studies by Lin et al. (2014) and Zong and Zhang (2023), the parameters values are chosen as follows:

Table 1.
Value of the Coefficients

p(i) $μ (i)$ $β (i)$ $η (i)$ $m (i)$ $α (i)$ $σ (i)$

i = 1 0.5/month 0.04/month 0.02/month 1.03/month 0.01/month 0.001/month 0.003/month

i = 2 0.6/month 0.05/month 0.03/month 1.04/month 0.02/month 0.002/month 0.004/month

	p(i)	$μ (i)$	$β (i)$	$η (i)$	$m (i)$	$α (i)$	$σ (i)$
i = 1	0.5/month	0.04/month	0.02/month	1.03/month	0.01/month	0.001/month	0.003/month
i = 2	0.6/month	0.05/month	0.03/month	1.04/month	0.02/month	0.002/month	0.004/month

6.1. The invariant measure

In this section, the existence of the unique ergodic invariant measure for the stochastic SIV epidemic model (3) with Markov switching is confirmed by several numerical simulations. Figure 1 shows the density kernels of solutions to (3) for the three groups $(S (t), I (t), V (t))$ . In the simulations, the density kernels are based on 10,000 sample paths. Comparing these density kernels, one can see that the density plots of each group at t = 12 months, 24 months, and 36 months remain almost the same. Therefore, the simulations strongly indicate the existence of a unique ergodic invariant measure for model (3). Figure 2 shows the density kernels of solutions to (3) for the three groups $(S (t), I (t), V (t))$ with and without control. In the simulations, the density kernels are based on 10,000 sample paths. Our initial values are presented in Table 1. Comparing these density kernels, one can see that the density plots of each group with and without the control variable remain almost the same. Therefore, the existence of a unique ergodic invariant measure for model (3) is proved.

FIG. 1.

Density plots (a–c) based on 10,000 stochastic simulations for group susceptible, infected, and vaccinated at time t = 12 months, 24 months, and 36 months.

FIG. 2.

Density plots (a–c) based on 10,000 stochastic simulations for group susceptible $S (t)$ , infected $I (t)$ , and vaccinated $V (t)$ with and without control.

6.2. Positivity of numerical solutions

Figures 3 and 4 show that the sample paths of solutions to the stochastic SIV epidemic model with Markov switching are the same as those of the classical EM scheme from 0 to 50 months. Figures 5 and 6 show the sample paths of the solutions given by the log-TEM numerical method. To compare these with the simulation results given by the classical EM method, all parameters are listed in Table 1. We find that the positivity of the simulated numerical solution of model (3) can be achieved by the log-TEM method. From these figures, one can see that the EM method is not sufficient to ensure the positivity of the numerical solutions.

FIG. 3.

The path of classical EM solution of $I (t)$ and corresponding state $Λ_{t}$ . EM, Euler-Maruyama.

FIG. 4.

The path of classical EM solution of $V (t)$ and corresponding state $Λ_{t}$ .

FIG. 5.

The path of log-TEM solution of $I (t)$ and corresponding state $Λ_{t}$ . TEM, truncated Euler-Maruyama.

FIG. 6.

The path of log-TEM solution of $V (t)$ and corresponding state $Λ_{t}$ .

6.3. Convergence rate

We regard the approximation with $Δ = e^{- 19}$ as the exact solution $I (t)$ , and $I^{(j)} (t_{k}), I_{k}^{(j)}$ represent the jth trajectories of the exact solution $I (t)$ and numerical solution I_k, respectively.

Next, we plotting the log approximation error with 500 independent trajectories. Figure 7a shows the error of the truncated logarithmic EM scheme for $I (t)$ . By similar methods above, let the exact solution $V (t)$ with regard to $Δ = e^{- 19}$ . $V^{(j)} (t_{k}), V_{k}^{(j)}$ represent the jth trajectories of the exact solution $V (t)$ and numerical solution V_k, respectively.

FIG. 7.

The solid-line trajectory depicts the approximation error of the exact solution and the log-TEM numerical solution of model (11) as the functions of stepsize $Δ \in {2^{- 5}, 2^{- 6}, 2^{- 7}, \dots, 2^{- 10}}$ .

The approximation error with 500 independent trajectories. Figure 7b shows the error of the truncated logarithmic EM scheme for $V (t)$ .

6.4. Relaxed control

In this section, we assume that the discounting factor of objective function is $δ = 0.05$ . We use the same parameters as shown in Table 1 and consider the cost function $ℱ (s, i, l, u_{1}, u_{2}) = 1 + i + 2 s u_{1}^{2} + 2 i u_{2}^{2} .$ The initial control is given by $u_{1} (0) \equiv max U_{1}, u_{2} (0) \equiv max U_{2}$ . $U_{1} = {k ∕ 5 : 0 \leq k \leq 4}$ , $U_{2} = {k ∕ 10 : 0 \leq k \leq 10}$ denote the set of controls. Thus, the vaccination control u₁ has five control levels and the treatment control u₂ has 11 control levels. The algorithm for relaxed control is given in Table 2. Figure 8a provides the optimal vaccination control and Fig. 8b provides the treatment control plan. The optimal vaccination plan and optimal treatment control appear to depend on the regimen switching.

FIG. 8.

The trajectory depicts the optimal control u₁ and u₂ of model (11).

Table 2.

Algorithm

Step 1: for

k_{1}, k_{2} \in Z_{+}, s = k_{1} h, i = k_{2} h, k_{1} + k_{2} \leq 1 ∕ h

do:

C_{0} (s, i, l) \equiv max U_{1},

D_{0} (s, i, l) \equiv max U_{2}

end for

for

k_{1} = 1, 2, \dots, 1 ∕ h

do:

V_{0}^{h} (s, i, l) \equiv \int_{0}^{\infty} e^{- δ t} ℱ (s, i, l, C_{0}, D_{0}) d t

Choose an initial

C_{0} (s, i, l), D_{0} (s, i, l)

, an initial step size s₀ and stopping tolerances

τ_{h}

We compute objective function of

ℱ

, i.e.,

ℱ (s, i, l, u_{1}, u_{2}) = 1 + i + 2 s u_{1}^{2} + 2 i u_{2}^{2},

end for

Step 2: for

k_{1} = 0, 1, \dots, n - 1

do: Equation (27)

for

j = 0, 1, 2, \dots, n

do:

\begin{matrix} V_{j + 1}^{h} (s, i, l | c) = e^{- δ Δ t^{h} (s, i, l, c)} \sum_{(s', i', l') \in S_{h}} V_{j}^{h} (s', i', l') p^{h} ((s, i, l), (s', i', l') | c) \\ + ℱ (s, i, l, c) Δ t^{h} (s, i, l, c) . \end{matrix}

end for

C_{j + 1}^{h} = a r c m i n_{c \in U} V_{j + 1}^{h} (s, i, l | c)

;

D_{j + 1}^{h} = a r c m i n_{d \in U} V_{j + 1}^{h} (s, i, l | d)

V_{j + 1}^{h} (s, i, l) = V_{j + 1}^{h} (s, i, l | C_{j + 1}^{h} (s, i, l), D_{j + 1}^{h} (s, i, l))

;

end for

Step 3: for

j = 0, 1, 2, \dots, n

do:

control, i.e.,

u_{1} = C_{j + 1}^{h}, u_{2} = D_{j + 1}^{h}

end for

7. CONCLUSIONS

In this article, the positivity of numerical solutions for a stochastic SIV epidemic model with Markov switching has been considered. A new numerical method that preserves positivity has been established under certain conditions. The proposed log-TEM scheme provides numerical solutions that are guaranteed to be positive. The existence of an invariant measure for a stochastic SIV epidemic model with Markov switching has been studied. Additionally, relaxed controls for the stochastic SIV epidemic model were investigated using the Markov chain method to approximate the continuous-time dynamics. These conditions lead to meaningful numerical approximations. Based on model (3), we have obtained the following results:

(a) The log-TEM scheme guarantees positive numerical solutions with strong order-1 convergence.

(b) Compared with previous schemes (Kar and Batabyal, 2011; Li et al., 2017; Magpantay et al., 2014), the assumptions for control problems have been relaxed, and the convexity requirements for the state function and cost function have been removed.

Footnotes

ACKNOWLEDGMENTS

The authors would like to express the sincerest gratitude to all the anonymous reviewers for their comments and opinions on the article, which were of great help to the article.

AUTHORs' CONTRIBUTIONS

All authors agree that they have read and approved the article.

AUTHOR DISCLOSURE STATEMENT

The authors declare they have no conflicting financial interests.

FUNDING INFORMATION

This work is supported by the National Natural Science Foundation of China (Grant No. 12161068).

APPENDIX

References

Anderson

, Medley

, May

, et al. Vaccination and herd immunity to infectious diseases. Nature, 1985; 318(6044):323–329.

Arino

, Mccluskey

, Van

. Global results for an epidemic model with vaccination that exhibits backward bifurcation. SIAM J Appl Math, 2003; 64(1):260–276.

Beretta

, Takeuchi

. Global stability of an sir epidemic model with time delays. J Math Biol, 1995; 33(3):250–260.

Berres

, Ruiz-Baier

. A fully adaptive numerical approximation for a two-dimensional epidemic model with nonlinear cross diffusion. Nonlinear Anal Real World Appl, 2011; 12(5):2888–2903.

Buonomo

, Lacitignola

, Vargas-De-Leon

. Qualitative analysis and optimal control of an epidemic model with vaccination and treatment. Math Comput Simul, 2014; 100(6):88–102.

Deng

, Liu

. Analysis of a stochastic tumor-immune model with regime switching and impulsive perturbations. Appl Math Modell, 2020; 78:482–504.

Enatsu

, Nakata

, Muroya

. Lyapunov functional techniques for the global stability analysis of a delayed sirs epidemic model. Nonlinear Analysis Real World Appl, 2012; 13(5):2120–2133.

On invariant measures and the asymptotic behavior of a stochastic delayed sirs epidemic model. Physica A Stat Mech Appl, 2019; 523:1008–1023.

Gray

, Greenhalgh

, Mao

, et al. The sis epidemic model with markovian switching. J Math Anal Appl, 2012; 394(2):496–516.

10.

Kar

, Batabyal

. Stability analysis and optimal control of an sir epidemic model with vaccination. Bio Syst, 2011; 104(2–3):127–135.

11.

Kaufmann

, Mcmichael

. Annulling a dangerous liaison: Vaccination strategies against aids and tuberculosis. Nat Med, 2005; 11(4):SS3.

12.

Keeling

Modelling vaccination strategies against foot-and-mouth disease. Nature, 2003; 421(6919):136–142.

13.

Laarabi

, Labriji

, Rachik

, et al. Optimal control of an epidemic model with a saturated incidence rate. Nonlinear Anal Modell Control, 2011; 4(4):448–459.

14.

Lamba

An adaptive euler-maruyama scheme for sdes: Convergence and stabilit. IMA J Numer Anal, 2007; 27(3):479–506.

15.

, Liu

, Cui

. Threshold dynamics and ergodicity of an sirs epidemic model with markovian switching. J Differ Equ, 2017; 263(12):8873–8915.

16.

Lin

, Jiang

, Wang

. Stationary distribution of a stochastic sis epidemic model with vaccination. Physica A Stat Mech Appl, 2014; 394:187–197.

17.

Liu

, Jiang

, Shi

. The threshold of a stochastic sis epidemic model with imperfect vaccination. Math Comput Simul, 2018; 144(2):78–90.

18.

Magpantay

, Riolo

, King

, Rohani

. Epidemiological consequences of imperfect vaccines for immunizing infections. SIAM J Appl Math, 2014; 74(6):1810–1830.

19.

Mao

, Yuan

Stochastic Differential Equations with Markovian Switching. Imperial College Press: London; 2006.

20.

Safan

, Rihan

. Mathematical analysis of an sis model with imperfect vaccination and backward bifurcation. Math Comput Simul, 2014; 96(2):195–206.

21.

Song

, Yin

, Zhang

. Numerical methods for controlled regime-switching diffusions and regime-switching jump diffusions. Automatica, 2006; 7(42):1147–1157.

22.

Telenti

, Arvin

, Corey

, et al. After the pandemic: Perspectives on the future trajectory of covid-19. Nature, 2021; 596(7873):1–14.

23.

Tildesley

, Savill

, Shaw

. Optimal reactive vaccination strategies for a foot-and-mouth outbreak in the UK. Nature, 2006; 440(7080):83–86.

24.

Tran

, Yin

. Optimal control and numerical methods for hybrid stochastic sis models. Nonlinear Anal Hybrid Syst, 2021; 41:101051.

25.

Wang

, Liu

, Zheng

, Takeuchi

. Qualitative and bifurcation analysis using an sir model with a saturated treatment function. Math Comput Modell, 2012; 55(3–4):710–722.

26.

Wei

, Chen

A delayed epidemic model with pulse vaccination. Discrete Dynam Nat Soc, 2014;2008(1).

27.

Worthmann

, Sauerteig

, Ritschel

, et al. How to coordinate vaccination and social distancing to mitigate sars-cov-2 outbreaks. SIAM J Appl Dynam Syst, 2021; 20(2):1135–1157.

28.

Yang

, Huang

. First order strong convergence of positivity preserving logarithmic euler-maruyama method for the stochastic sis epidemic model. Appl Math Lett, 2021; 121(4):107451.

29.

Yuan

, Mao

. Convergence of the euler-maruyama method for stochastic differential equations with markovian switching. Math Comput Simul, 2004; 64(2):223–235.

30.

Zaman

, Kang

, Jung

. Optimal vaccination and treatment in the sir epidemic model. Math Comput Simul, 2007; 3(3):240–249.

31.

Zhang

, Jiang

, Alsaedi

, et al. Stationary distribution of stochastic sis epidemic model with vaccination under regime switching. Appl Math Lett, 2016; 59:87–93.

32.

Zhao

, Jiang

. The asymptotic behavior and ergodicity of stochastically perturbed svir epidemic model. Int J Biomath, 2015; 9(3):177–190.

33.

Zhen

, Ma

, Han

. Global stability of an sirs epidemic model with delays. Acta Math Scient, 2006; 26(2):291–306.

34.

Zong

, Zhang

. Near-optimal control of a stochastic partial differential equation seir epidemic model under economic constraints. Eur J Control, 2023; 69:100752.

35.

Zong

, Zhang

, Li

. Markovian switching for near-optimal control of a stochastic SIV epidemic model. Math Biosci Eng, 2016; 16(3):1348–1375.

Positivity-Preserving Numerical Method and Relaxed Control for Stochastic Susceptible–Infected–Vaccinated Epidemic Model with Markov Switching

Abstract

1. INTRODUCTION

2. PRELIMINARIES

Table 1. Value of the Coefficients p(i) μ ( i ) β ( i ) η ( i ) m ( i ) α ( i ) σ ( i ) i = 1 0.5/month 0.04/month 0.02/month 1.03/month 0.01/month 0.001/month 0.003/month i = 2 0.6/month 0.05/month 0.03/month 1.04/month 0.02/month 0.002/month 0.004/month

Footnotes

ACKNOWLEDGMENTS

AUTHORs' CONTRIBUTIONS

AUTHOR DISCLOSURE STATEMENT

FUNDING INFORMATION

APPENDIX

References

Table 1.
Value of the Coefficients

p(i) $μ (i)$ $β (i)$ $η (i)$ $m (i)$ $α (i)$ $σ (i)$

i = 1 0.5/month 0.04/month 0.02/month 1.03/month 0.01/month 0.001/month 0.003/month

i = 2 0.6/month 0.05/month 0.03/month 1.04/month 0.02/month 0.002/month 0.004/month