Dynamics and Control of a Novel Discrete Internet Rumor Propagation Model in a Multilingual Environment

Abstract

In the Internet age, the development of intelligent software has broken the limits of multilingual communication. Recognizing that the data collected on rumor propagation are inherently discrete, this study introduces a novel SIR discrete Internet rumor propagation model with the general nonlinear propagation function in a multilingual environment. Then, the propagation threshold $R_0$ is obtained by the next-generation matrix method. Besides, the criteria determining the spread or demise of rumors are obtained by the stability theory of difference equations. Furthermore, combined with optimal control theory, prevention and refutation mechanisms are proposed to curb rumors. Finally, the validity and applicability of the model are demonstrated by numerical simulations and a real bilingual rumor case study.

1. Introduction

Rumor usually refers to statements without corresponding factual basis, which are fabricated and spread by certain means [1]. With the arrival of the Internet era, socializing through the network has become a mainstream way. Social networks have made the spread of rumors no longer limited by language, time and space. As a result, it has brought greater harm to the productive life of society. For example, in the beginning of 2020, rumors emerged in some parts of India that certain foods or herbs could prevent or cure COVID-19. Soon after, the rumor appeared on some social media platforms in many countries, which hampered the fight against the epidemic. Therefore, it is of important and realistic significance to investigate the rumor propagation mechanism in a multilingual environment.

Due to the similarity of a rumor spreading mechanism, most studies on rumors are based on infectious disease models. In 1965, Daley and Kendal proposed a classical rumor spreading model, i.e., the DK model, by applying compartment modeling ideas [2]; the model lays the foundation for the dynamic analysis of rumor spreading. Subsequently, Maki and Thomson [3] developed the MK model based on the DK model. The classic DK model and MK model pioneered the study of rumor propagation dynamics. Since then, many rumor models have been proposed, such as the SIR [4,5] model, the SEIR [6,7] model, the ISCR [8,9] model, and the SIDRW [10] model. With the further research, many realistic factors have been added to the rumor model, such as the educational mechanism [11], the refutation mechanism [12], the trust mechanism [13] and correction mechanism [14]. Some excellent results have also been achieved with regard to multilingual environments. For example, in 2021, Yu et al. [15] investigated two new 2I2SR rumor propagation models on OSN. In 2023, Ye et al. [16] investigated a fractional-order reaction diffusion rumor propagation model in a multilingual environment.

It is noted that data on rumor propagation are collected in discrete time units, but current research about multilingual rumor propagation is limited to building continuous models. Continuous modeling requires a large amount of data support, and acquire the data may be confined, resulting in higher difficulty and cost of model construction. Meanwhile, due to the simplification of the model and the deviation of the parameters, it is difficult for the rumor continuous model to accurately predict the results of rumor propagation. In this paper, we will make up for the shortcomings of the existing work with the use of differential equations to build a discrete rumor propagation model in a multilingual environment, and analyze its stability to obtain the propagation dynamics and propagation regularity.

Rumor propagation is a complex phenomenon typically involving multiple factors such as social psychology, information dissemination channels, and individual cognitive differences. These influencing factors often result in the propagation of rumors exhibiting nonlinear characteristics. For instance, the effect of increased education levels may vary across different social groups and regions, thereby manifesting nonlinear features. For this purpose, some excellent nonlinear functions have been proposed. Concretely, saturation incidence ${\displaystyle \frac{\beta SI}{1+\alpha I}}$ has been proposed to characterize psychological effects [9], where S and I represent the number of susceptible individuals and infected individuals respectively, $\beta$ is the disease transmission coefficient and $\alpha$ is the parameter measures the psychological or inhibitory effect. In addition, the nonlinear propagation function $IS(1+I)$ has been raised to characterize the multi-body interaction [17]. General nonlinear models are more applicable and have achieved good results. In 2021, Xia et al. [18] investigated the ILSR rumor propagation model with general nonlinear propagation rate, and found that considering nonlinear spreading rate in the rumor spreading model can better adapt to the realistic situation. Therefore, it is necessary to consider general nonlinear propagation rates in rumor propagation models.

With the aim of controlling the propagation of rumors, many control methods have been proposed, such as target immunization [19], artificial immunization [20], and efficient immunization [21]. However, most of the controls are continuous methods, which may have some shortcomings, mainly including waste of resources and inefficiency. Therefore, this paper will design a discrete optimization control method, including prevention mechanism and rumor refutation mechanism. Among them, the preventive mechanism is to convey the civilized concept of “not believing in rumors and not spreading rumors” by holding various open classes on online rumor refutation, so as to improve one’s own discernment ability and reduce the spread of rumors. The rumor refutation mechanism is the governmental control and establishment of the rumor debunking mechanism, which stops rumors from spreading. By applying the optimal control theory, we obtain the optimal control intensity, and the conclusions that we get will provide theoretical support for the strategy development of government departments.

At present, there is very little work to establish a discrete model to analyze the spread of rumors. As far as we know, only a discrete rumor propagation model has been established, and no complete propagation dynamics analysis is carried out in [22]. Moreover, the development of online social networks makes multilingual communication possible. The current research on rumor propagation in multilingual environment only explores the continuous model. In fact, discrete models are often better matched with the sampling frequency of data sources such as social media because the data itself is discrete. Specifically, the discrete features of data are more pronounced in multilingual environments. Therefore, we consider a novel SIR discretel rumor propagation model in a multilingual environment, which also considers general nonlinear propagation rates. The contributions of this paper are as follows:

• Different from the continuous rumor models in the multilingual environment [14,18], an SIR discrete rumor propagation model is developed;
• Contrary to the linear propagation rate in [23,24], we introduce a rumor propagation model with a general nonlinear propagation rate in a multilingual environment;
• The optimal control problem of discrete multilingual SIR rumor propagation model is analyzed, the optimal control intensities for the prevention mechanism $u^{{\mathsf{\Gamma }}_k}\left(n\right)$ and the refutation mechanism $v^{{\mathsf{\Gamma }}_k}\left(n\right)$ are provided.

The paper is organized as follows: Section 2 gives the model description and prerequisite knowledge. Section 3 discusses the dynamics of the rumor-free equilibrium and rumor-prevailing equilibrium of the SIR rumor propagation model. Section 4 analyzes the optimal control problem. In Section 5, numerical simulations are performed. In Section 6, numerical simulation examples are given to illustrate the validity of the results obtained. Finally, conclusions are given in Section 7.

2. Model Instruction and Preliminaries

Some definitions and lemmas are given below for the convenience of subsequent work.

2.1. Preliminaries

Definition 1

([25]). A matrix $A={\left(a_{ij}\right)}_{n\times n}$ is called irreducible if there does not exist a permutation matrix P such that

$$P^{T}AT=\left(\begin{array}{cc}{A}_{11}& A_{12}\\ 0& A_{22}\end{array}\right),$$

(1)

where $A_{11}$ and $A_{22}$ are square matrices of order r and $n-r$ respectively, with $1\le r<n$. If such a permutation exists, the matrix A is said to be reducible.

Lemma 1

([26]). A weighted digraph $(G,A)$ is strongly connected if and only if and ongly if the weight matrix A is irreducible. The Laplacian matrix $L=\left[l_{ij}\right]$ of $(G,A)$ is defined as

$$\begin{array}{c}{l}_{ij}=\left\{\begin{array}{c}-a_{ij}\mathrm{for}i\ne j,\hfill \\ \sum _{k\ne i}{a}_{ij}\mathrm{for}i=j,\hfill \end{array}\right.\hfill \end{array}$$

The following result gives a graph-theoretic description of the cofactors of the diagonal entries of L. Suppose $n\ge 2$ and denote $c_i$ be the cofactor of $l_{ii}$ in L. Then $c_i$ = ${\displaystyle {\sum }_{\mathcal{T}\in {\mathbb{T}}_i}\omega \left(\mathcal{T}\right)}$, $i=1,2,\cdots ,n$ where ${\mathbb{T}}_i$ is the set of all spanning trees $\mathcal{T}$ of $(G,A)$ that are rooted at vertex i, and $\omega \left(\mathcal{T}\right)$ is the weight of $\mathcal{T}$. If $(G,A)$ is strongly connected, then $c_i>0$ for $1\le i\le n$.

Suppose $c_i$ be as given in the Kirchhoff’s matrix tree theorem, and denote ${\left\{H_i\left(y\right)\right\}}_{i=1}^m$ be any family of functions with $y={(y_1,\cdots ,y_n)}^T\in {\mathbb{R}}^n$, then

$$\sum _{i,j=1}^m{c}_i{a}_{ij}{H}_i\left(y\right)=\sum _{i,j=1}^m{c}_i{a}_{ij}{H}_j\left(y\right)$$

Assumption A1.

Let $F_{ij}(S_i,I_j)$ represent the general nonlinear propagation function. Assume that the function $F_{ij}(S_i,I_j)(i,j=1,2,\cdots ,k)$ satisfies

(a) 

$F_{ij}(S_i\left(n\right),I_j\left(n\right)=f_i\left(S_i\left(n\right)\right)g_j\left(I_j\left(n\right)\right)$,

(b) 

$f_i\left(S\right)$ is monotonically non-decreasing, and $f_i\left(0\right)=0$,

(c) 

I/$g_j\left(I\right)$ is monotone non-decreasing.

(d) 

$f_k\left(S_k\left(n\right)\right)g_j\left(I_j\left(n\right)\right)\le S_k\left(n\right)I_j\left(n\right)$

2.2. Model Instruction

In this section, we consider a multilingual model that assumes that the total population is divided into different groups based on different languages. We will assume in the SIR rumor spreading model that rumors can be spread through i languages.

It is assumed that the population can be divided into m groups: Group mth (persons spreading rumors through the mth language). Each group is divided into three categories: those who were never informed of the rumor (Spreaders), those who were informed of the rumor and spread it (Ignorants), and those who already knew the truth or chose not to spread the rumor due to media reports, official statements, or forgetting mechanisms (Removers). Where $S_i\left(n\right),I_i\left(n\right),R_i\left(n\right)$ denote the density of Ignorants, Spreaders and Removers in group i at time n,respectively. We assume that the coming probability is equal to the leaving probability. Therefore, at each time n, it follows that $e_k=b_k^S{S}_k\left(n\right)+b_k^I{I}_k\left(n\right)+b_k^R{R}_k\left(n\right)$. At the same time, $S_k\left(n\right)$ + $I_k\left(n\right)$ + $R_k\left(n\right)$ = 1. The parameter symbols are shown in

Table 1. Meaning of required parameters in model (2).

Symbols	Implications
$e_k$	the coming probability of the susceptibility ($k=1,2,\cdots$, m).
$b_k^S$	the leaving probability of the $S_k$ ($k=1,2,\cdots$, m).
$b_k^I$	the leaving probability of the $I_k$ ($k=1,2,\cdots$, m).
$b_k^R$	the leaving probability of the $R_k$ ($k=1,2,\cdots$, m).
${\beta }_{kj}$	the cross-transmitted probability from the susceptible $S_k$ to the infected $I_j$
	($k\ne j,k,j=1,2,\cdots$, m).
${\gamma }_k$	the transfer probability from spreader $I_k$ to stiflers $R_k$
	due to forgetting mechanism ($k=1,2,\cdots$, m).
$〈k〉$	the average degree of the homogeneous network.

. Since some people speak multiple languages, different groups can also exchange information. Define ${\beta }_{kj}$ as the rumor propagation coefficient between $S_k$ and $I_j$, when $1\le k,j\le m$. Suppose that the state transformation rule of the rumor is shown in Figure 1, where the parameters are shown in Table 1. Based on the above analysis, we can build the following SIR discrete rumor propagation model in a multilingual environment.

$$\left\{\begin{array}{c}{S}_k(n+1)=S_k\left(n\right)+e_k-b_k^S{S}_k\left(n\right)-〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{F}_{kj}\left(S_k\left(n\right),I_j\left(n\right)\right),\hfill \\ I_k(n+1)=I_k\left(n\right)+〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{F}_{kj}\left(S_k\left(n\right),I_j\left(n\right)\right)-\left(b_k^I+{\gamma }_k\right)I_k\left(n\right),\hfill \\ R_k(n+1)=R_k\left(n\right)+{\gamma }_k{I}_k\left(n\right)-b_k^R{R}_k\left(n\right),\hfill \end{array}\right.$$

(2)

where the initial conditions is $S_k\left(0\right)>0,I_k\left(0\right)>0,R_k\left(0\right)>0$ and $S_k\left(0\right)$ + $I_k\left(0\right)$ + $R_k\left(0\right)$ = 1.

Since all the parameters are the description of the probability of some specific behavior, the values of each parameter are in the interval [0, 1].

Remark 1.

The positive invariant set $\mathsf{\Omega }$ is defined as

$$\mathsf{\Omega }=\left\{(S_1\left(n\right),I_1\left(n\right),S_2\left(n\right),I_2\left(n\right),\cdots ,S_m\left(n\right),I_m\left(n\right))\in R^{2m}\left|\begin{array}{c}{\displaystyle 0\le S_k\left(n\right)\le \frac{e_k}{b_k^S}}\hfill \\ {\displaystyle 0\le S_k\left(n\right)+I_k\left(n\right)\le \frac{e_k}{b_k^{*}}}\hfill \end{array}\right.\right\},$$

(3)

where $k=1,2,\cdots ,mand{b}_k^{*}=min\{b_k^S,b_k^I+{\gamma }_k\}$. Denote the interior of $\mathsf{\Omega }$ by ${\mathsf{\Omega }}^0$.

3. Theoretical Analysis

In this section, the positivity of the solution to model (2), the stability of the rumor-free equilibrium and the rumor-prevailing equilibrium will be shown.

3.1. The Positive and Boundedness of Solutions

Based on Assumption 1, the model (2) can also be rewritten:

$$\left\{\begin{array}{c}{S}_k(n+1)=S_k\left(n\right)+e_k-b_k^S{S}_k\left(n\right)-〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k\left(n\right)\right)g_j\left(I_j\left(n\right)\right),\hfill \\ I_k(n+1)=I_k\left(n\right)+〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k\left(n\right)\right)g_j\left(I_j\left(n\right)\right)-\left(b_k^I+{\gamma }_k\right)I_k\left(n\right),\hfill \\ R_k(n+1)=R_k\left(n\right)+{\gamma }_k{I}_k\left(n\right)-b_k^R{R}_k\left(n\right),\hfill \end{array}\right.$$

(4)

where the initial conditions is $S_k\left(0\right)>0, I_k\left(0\right)>0, R_k\left(0\right)>0$ and $S_k\left(0\right)$ + $I_k\left(0\right)$ + $R_k\left(0\right)$ = 1.

Because the state R in model (4) only appears in the third equation, we can decouple the R in the equations. Therefore, the dynamic equation of model (4) is equivalent to that of model (5).

$$\left\{\begin{array}{c}{S}_k(n+1)=S_k\left(n\right)+e_k-b_k^S{S}_k\left(n\right)-〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k\left(n\right)\right)g_j\left(I_j\left(n\right)\right),\hfill \\ I_k(n+1)=I_k\left(n\right)+〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k\left(n\right)\right)g_j\left(I_j\left(n\right)\right)-\left(b_k^I+{\gamma }_k\right)I_k\left(n\right),\hfill \end{array}\right.$$

(5)

Theorem 1.

Suppose $1-{\displaystyle \sum _{j=1}^m}{\beta }_{kj}\ge 0$. Let $(S_k\left(n\right),I_k\left(n\right))$ be the solution to model (4) with initial condition $S_k\left(0\right)>0$, $I_k\left(0\right)>0$, $R_k\left(0\right)>0$ and $S_k\left(0\right)$ + $I_k\left(0\right)$ + $R_k\left(0\right)$ = 1, then we have

(1) 

the solution is positive, namely $S_k\left(n\right)>0,I_k\left(n\right)>0$ for any $n>0$,

(2) 

$$\left\{\begin{array}{c}{\displaystyle {lim}_{n\to +\infty }}sup{S}_k\left(n\right)\le {\displaystyle \frac{e_k}{b_k^S}},\hfill \\ {\displaystyle {\displaystyle {lim}_{n\to +\infty }}sup(S_k\left(n\right)+I_k\left(n\right))\le \frac{e_k}{b_k^{*}},}\hfill \end{array}\right.$$

namely, the solutions are bounded, where $k=1,2,\cdots ,m$, $d_k^{*}=max\{b_k^S,b_k^I+{\gamma }_k\}$ and $b_k^{*}=min\{b_k^S,b_k^I+{\gamma }_k\}$.

Proof of Theorem 1.

First of all, according to the first equation of the model (4) and the initial condition $S_k\left(0\right)>0, I_k\left(0\right)>0, R_k\left(0\right)>0$, we can obtain that

$$I_k\left(1\right)=I_k\left(0\right)-\left(b_k^I+{\gamma }_k\right)I_k\left(0\right)+〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k\left(0\right)\right)g_j\left(I_j\left(0\right)\right).$$

That is

$$\begin{array}{cc}\hfill I_k\left(1\right)& =I_k\left(0\right)-\left(b_k^I+{\gamma }_k\right)I_k\left(0\right)+〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k\left(0\right)\right)g_j\left(I_j\left(0\right)\right),\hfill \\ & >I_k\left(0\right)\left[1-\left(b_k^I+{\gamma }_k\right)\right].\hfill \end{array}$$

If $0\le b_k^I+{\gamma }_k\le 1$, then $I_k\left(1\right)>I_k\left(0\right)\left[1-\left(b_k^I+{\gamma }_k\right)\right]\ge 0$. By induction, when $0\le b_k^I+{\gamma }_k\le 1$ and $I_k\left(0\right)>0$, for all $n\ge 0$, $I_k\left(n\right)>0$, namely $I_k\left(n\right)>0$.

Secondly, according to the first equation of the model (4), we can obtain that

$$S_k(n+1)=S_k\left(n\right)+e_k-b_k^S{S}_k\left(n\right)-〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k\left(n\right)\right)g_j\left(I_j\left(n\right)\right).$$

According to (d) of Assumption 1, we can obtain that

$$\begin{array}{cc}\hfill S_k(n+1)& =S_k\left(n\right)+e_k-b_k^S{S}_k\left(n\right)-〈k〉\sum _{j=1}^m{\beta }_{kj}{f}_k\left(S_k\left(n\right)\right)g_j\left(I_j\left(n\right)\right)\hfill \\ & \ge e_k+(1-b_k^S)S_k\left(n\right)-〈k〉\sum _{j=1}^m{\beta }_{kj}{S}_k\left(n\right)I_j\left(n\right),\hfill \end{array}$$

Then, we can get

$$\begin{array}{cc}\hfill S_k(n+1)& \ge e_k+(1-b_k^S)S_k\left(n\right)-〈k〉\sum _{j=1}^m{\beta }_{kj}{S}_k\left(n\right)·1\hfill \\ & =e_k+(1-b_k^S-\sum _{j=1}^m{\beta }_{kj})S_k\left(n\right).\hfill \end{array}$$

According to the assumption $1-{\displaystyle \sum _{j=1}^m}{\beta }_{kj}\ge 0$ and $e_k-b_k^S{S}_k\left(n\right)>0$, we have $S_k(n+1)>0$ if $S_k\left(n\right)>0$.

At the same time, we know the initial condition $S_k\left(0\right)>0$, then for each time n, we can deduce

$$\begin{array}{cc}{S}_k(n+1)& \ge e_k+(1-b_k^S-{\displaystyle \sum _{j=1}^m}{\beta }_{kj})S_k\left(n\right)>0.\end{array}$$

namely, $S_k$(n) > 0.

Finally, from the third equation of model (4), we can get

$$\begin{array}{cc}\hfill R_k\left(n\right)& ={(1-b_k^R)}^n{R}_k\left(0\right)+{(1-b_k^R)}^{n-1}\left[{\gamma }_k{I}_k\left(0\right)\right]+{(1-b_k^R)}^{n-2}\left[{\gamma }_k{I}_k\left(1\right)\right]\hfill \\ & +\cdots +{\gamma }_k{I}_k(n-1)>0.\hfill \end{array}$$

with nonnegative initial value for $n\ge 0$. According to Theorem 1.1 of reference [26], the positivity of our solution is proved.

Nextly, we prove the second conclusion. Let $V_k\left(n\right)=S_k\left(n\right)+I_k\left(n\right)$. Then, we have

$$\begin{array}{c}\hfill V_k(n+1)-V_k\left(n\right)=e_k-b_k^S{S}_k\left(n\right)-\left(b_k^I+{\gamma }_k\right)I_k\left(n\right)\le e_k-b_k^{*}{V}_k\left(n\right).\end{array}$$

So,

$$\underset{n\to +\infty }{lim}sup(S_k\left(n\right)+I_k\left(n\right))\le {\displaystyle \frac{e_k}{b_k^{*}}},$$

Similarly, from the first equation of the model (5), we further have

$$\underset{n\to +\infty }{lim}sup{S}_k\left(n\right)\le {\displaystyle \frac{e_k}{b_k^S}}.$$

 □

3.2. Rumor-Free Equilibrium

In this section, our goal is to explore the dynamical behavior of rumor-free equilibrium. First, it is not difficult to find that model (5) has a rumor-free equilibrium

$$E^0=(S_1^0,0,S_2^0,0,\cdots ,S_m^0,0),$$

where ${\displaystyle S_k^0=\frac{e_k}{b_k^S},}$$k=1,2,\cdots ,m$.

For the convenience of the narrative, define a function T: $R^m\to R^{m\times m}$,

$$T(S_1\left(n\right),S_2\left(n\right),\cdots ,S_m\left(n\right))={({\displaystyle \frac{〈k〉{\beta }_{kj}{f}_k\left(S_k\left(n\right)\right)g_j^{\prime }\left(0\right)}{b_k^I+{\gamma }_k}})}_{m\times m}.$$

Then, by the method of next-generation matrix, we obtain the propagation threshold $R_0$ that

$$R_0=\rho \left(T_0\right),$$

where $\rho$ is the spectral radius, ${\displaystyle T_0={(\frac{〈k〉{\beta }_{kj}{f}_k\left(S_k^0\right)g_j^{\prime }\left(0\right)}{b_k^I+{\gamma }_k})}_{m\times m}.}$

Theorem 2.

Let $B={\left({\beta }_{kj}\right)}_{m\times m}$ and Assumption 1 holds.

(1) 

If $R_0\le 1$, then the rumor-free equilibrium $E^0$ is globally asymptotically stable on Ω.

(2) 

If $R_0>1$, then the rumor-free equilibrium $E^0$ is unstable.

Proof of Theorem 2.

Since B is irreducible, it follows that $T_0$ is irreducible. Based on the Perron-frobenius theorem, there exists a positive principal eigenvector $\lambda =\left\{({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_m)\right\}$ such that ${\lambda }_k>0$ for $k=1,2,\cdots ,m$ and $\lambda \rho \left(T_0\right)=\lambda T_0$. Consider the following Lyapunov function:

$$\begin{array}{c}{M}_n={\displaystyle \sum _{k=1}^m}{\displaystyle \frac{{\lambda }_k}{b_k^I+{\gamma }_k}}\left[S_k\left(n\right)-S_k^0-{\int }_{S_k^0}^{S_k\left(n\right)}{\displaystyle \frac{f_k\left(S_k^0\right)}{f_k\left(\xi \right)}d\xi +I_k\left(n\right)}\right].\hfill \end{array}$$

It can be obtained that

$${\int }_{S_k\left(n\right)}^{S_k(n+1)}{\displaystyle \frac{f_k\left(S_k^0\right)}{f_k\left(\xi \right)}}d\xi \ge {\displaystyle \frac{(S_k(n+1)-S_k\left(n\right))f_k\left(S_k^0\right)}{f_k\left(S_k(n+1)\right)}},S_k\left(n\right)>0,S_k(n+1)>0,\mathrm{for}k=1,2,\cdots ,m.$$

Note that $g_j\left(I\right)\le g_j^{\prime }\left(0\right)I$ for all $I>0$. Then, one has

$$\begin{array}{cc}& M_{n+1}-M_n\hfill \\ \hfill =& {\displaystyle \sum _{k=1}^m}{\displaystyle \frac{{\lambda }_k}{b_k^I+{\gamma }_k}}\left[S_k(n+1)-S_k\left(n\right)-{\int }_{S_k\left(n\right)}^{S_k(n+1)}{\displaystyle \frac{f_k\left(S_k^0\right)}{f_k\left(\xi \right)}}d\xi +I_k(n+1)-I_k\left(n\right)\right],\hfill \\ \hfill \le & {\displaystyle \sum _{k=1}^m}{\displaystyle \frac{{\lambda }_k}{b_k^I+{\gamma }_k}}\left[\left(1-{\displaystyle \frac{f_k\left(S_k^0\right)}{f_k\left(S_k(n+1)\right)}}\right)(e_k-b_k^S{S}_k\left(n\right)-〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k\left(n\right)\right)g_j\left(I_j\left(n\right)\right))\right.\hfill \\ & \left.+〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k\left(n\right)\right)g_j\left(I_j\left(n\right)\right)-(b_k^I+{\gamma }_k)I_k\left(n\right)\right],\hfill \\ \hfill \le & {\displaystyle \sum _{k=1}^m}{\displaystyle \frac{{\lambda }_k{b}_k^S}{f_k\left(S_k(n+1)\right)\left(b_k^I+{\gamma }_k\right)}}\left(f_k\left(S_k(n+1)\right)-f_k\left(S_k^0\right)\right)\left(S_k^0-S_k\left(n\right)\right)\hfill \\ & +{\displaystyle \sum _{k=1}^m}{\displaystyle \frac{{\lambda }_k}{b_k^I+{\gamma }_k}}\left[{\displaystyle \frac{〈k〉{\sum }_{j=1}^m{\beta }_{kj}{f}_k\left(S_k^0\right)f_k\left(S_k\left(n\right)\right)g_j\left(I_j\left(n\right)\right)}{f_k\left(S_k(n+1)\right)}}-I_k\left(n\right)(b_k^I+{\gamma }_k)\right],\hfill \end{array}$$

$$\begin{array}{cc}\hfill \le & \sum _{k=1}^m{\displaystyle \frac{{\lambda }_k{b}_k^S}{f_k\left(S_k(n+1)\right)\left(b_k^I+{\gamma }_k\right)}}\left(f_k\left(S_k(n+1)\right)-f_k\left(S_k^0\right)\right)\left(S_k^0-S_k\left(n\right)\right)\hfill \\ & +\sum _{k=1}^m{\displaystyle \frac{{\lambda }_k}{b_k^I+{\gamma }_k}}\left[〈k〉\sum _{j=1}^m{\beta }_{kj}{f}_k\left(S_k^0\right)g_j\left(I_j\left(n\right)\right)-I_k\left(n\right)(b_k^I+{\gamma }_k)\right],\hfill \\ \hfill \le & \sum _{k=1}^m{\displaystyle \frac{{\lambda }_k{b}_k^S}{f_k\left(S_k(n+1)\right)\left(b_k^I+{\gamma }_k\right)}}\left(f_k\left(S_k(n+1)\right)-f_k\left(S_k^0\right)\right)\left(S_k^0-S_k(n+1)\right)\hfill \\ & +\sum _{k=1}^m{\lambda }_k\left({\displaystyle \frac{〈k〉{\sum }_{j=1}^m{\beta }_{kj}{f}_k\left(S_k^0\right)g_j^{\prime }\left(0\right)I_j\left(n\right)}{b_k^I+{\gamma }_k}}-I_k\left(n\right)\right),\hfill \\ \hfill =& \sum _{k=1}^m{\displaystyle \frac{{\lambda }_k{b}_k^S}{f_k\left(S_k(n+1)\right)\left(b_k^I+{\gamma }_k\right)}}\left(f_k\left(S_k(n+1)\right)-f_k\left(S_k^0\right)\right)\left(S_k^0-S_k(n+1)\right)\hfill \\ & +\left({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_m\right)\left(T_0{I}_n-I_n\right),\hfill \end{array}$$

$$\begin{array}{cc}\hfill =& \sum _{k=1}^m{\displaystyle \frac{{\lambda }_k{b}_k^S}{f_k\left(S_k(n+1)\right)\left(b_k^I+{\gamma }_k\right)}}\left(f_k\left(S_k(n+1)\right)-f_k\left(S_k^0\right)\right)\left(S_k^0-S_k(n+1)\right)\hfill \\ & +\left(\rho \left(T_0\right)-1\right)\left({\lambda }_1,{\lambda }_2\cdots ,{\lambda }_m\right)I_n,\hfill \end{array}$$

where $I_n={\left(I_1\left(n\right),I_2\left(n\right),\dots I_m\left(n\right)\right)}^T$. Thus is $M_{n+1}-M_n\le 0$. It follows that ${\left\{M_n\right\}}_{n\in N}$ is a decreasing sequence when $R_0\le 1$. Thus, there is a constant M such that $\underset{n\to \infty }{lim}{M}_n=\tilde{M}$ and there are $\underset{n\to \infty }{lim}\left(M_{n+1}-M_n\right)=0$. In addition, we have

(1)

When $R_0<1$

$$\underset{n\to \infty }{lim}\left(M_{n+1}-M_n\right)=0⟺\underset{n\to \infty }{lim}{S}_k\left(n\right)=S_k^0⟺\underset{n\to \infty }{lim}{I}_k\left(n\right)=0,$$

for $k=1,2,\cdots ,m$.

(2)

When $R_0=1$

$$\underset{n\to \infty }{lim}\left(M_{n+1}-M_n\right)=0⟺\underset{n\to \infty }{lim}{S}_k\left(n\right)=S_k^0,$$

for $k=1,2,\cdots ,m$.

We refer to $\underset{n\to \infty }{lim}{I}_k\left(n\right)=0$ for $k=1,2,\cdots ,m$. Otherwise, it must have the subsequence $\left\{n_p\right\}$ and $\tilde{I}={(\widetilde{I_1},\widetilde{I_2},\cdots ,{\tilde{I}}_m)}^T\ge 0$ which makes $\underset{n\to \infty }{lim}{I}_n=\tilde{I}\ne 0$. Choose the subsequence $\left\{n_p\right\}$ of the first equation in model (5).

$$\begin{array}{c}{S}_k\left(n_{q+1}\right)-S_k\left(n_q\right)=e_k-b_k^S{S}_k\left(n_q\right)-〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k{S}_k\left(n_q\right)\left(n_q\right)\right)g_j\left(I_j\left(n_q\right)\right),\hfill \end{array}$$

$k=1,2,\dots ,m.$,

Let $n_q$ tend to positive infinity in the previous equation

$$0=0-〈k〉\sum _{j=1}^m{\beta }_{kj}{f}_k\left(S_k^0\right)g_j\left(\widetilde{I_j}\right)\mathrm{for}k=1,2,\cdots ,m.$$

Since $B={\left({\beta }_{kj}\right)}_{m\times m}$ is non-negative and irreducible, one obtains $g_j\left({\tilde{I}}_j\right)=0$, it is a contradiction to say that ${\tilde{I}}_j=0$ when $j=1,2,\cdots ,m$. So $\underset{n\to \infty }{lim}{I}_k\left(n\right)=0$ for $k=1,2,\cdots ,$m. In summary, which is contradictory when $E^0$ is globally asymptotically stable when $R\le 1$ is present.

When $R_0>1$ and $I\ne 0$, there are $({\lambda }_1,{\lambda }_2,\dots {\lambda }_m)T_0-({\lambda }_1,{\lambda }_2,\dots {\lambda }_m)E_n=\left[\rho \left(T_0\right)-1\right]·({\lambda }_1,{\lambda }_2,\dots {\lambda }_m)>0,$ where $E_n$ is the identity matrix of order n. Due to the continuity, we can obtain that $M_{n+1}-M_n>0$ in a neighborhood $E^0$ of ${\mathsf{\Omega }}^0$, which shows that $E^0$ is unstable.

3.3. The Rumor-Prevailing Equilibrium

In this section, we explore the dynamical behavior of the rumor-prevailing equilibrium. It is easy to find that there is a rumor-prevailing equilibrium $E^{*}$ for model (4),

$$E^{*}=(S_1^{*},I_1^{*},S_2^{*},I_2^{*},\cdots ,S_m^{*},I_m^{*}),$$

satisfing

$$\left\{\begin{array}{c}{e}_k=b_k^S{S}_k^{*}+〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right),\hfill \\ 〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)=\left(b_k^I+{\gamma }_k\right)I_k^{*}.\hfill \end{array}\right.$$

 □

Theorem 3.

Suppose $B={\left({\beta }_{kj}\right)}_{m\times m}$ is irreducible and Assumption 1 holds. And if $R_0>1$, and $g_j\left(I_j\right)(1\le j\le m)$ satisfies $\left({\displaystyle \frac{g_j\left(I_j\left(n\right)\right)}{g_j\left(I_j^{*}\right)}-\frac{I_j\left(n\right)}{I_j^{*}}}\right)\left({\displaystyle 1-\frac{g_j\left(I_j^{*}\right)}{g_j\left(I_j\left(n\right)\right)}}\right)\le 0$ for any $I>0,$ then $E^{*}$ is globally asymptotically stable in ${\mathsf{\Omega }}^0$.

Proof of Theorem 3.

Defining Lyapunov function:

$$\begin{array}{cc}\hfill P_n& =\sum _{k=1}^m{v}_k\left[S_k\left(n\right)-S_k^{*}-{\int }_{S_k^{*}}^{S_k\left(n\right)}{\displaystyle \frac{f_k\left(S_k^{*}\right)}{f_k\left(t\right)}}dt+I_k\left(n\right)-I_k^{*}ln{\displaystyle \frac{I_k(n+1)}{I_k\left(n\right)}}\right.\hfill \\ & +\left.〈k〉\sum _{j=1}^m{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)\left(\mathsf{\Psi }\left({\displaystyle \frac{I_k(n+1)}{I_k^{*}}}\right)-\mathsf{\Psi }\left({\displaystyle \frac{I_k\left(n\right)}{I_k^{*}}}\right)\right)\right]\hfill \end{array}$$

where $\mathsf{\Psi }\left(x\right)=x-1-lnx.$

We have

$\begin{array}{cc}& P_{n+1}-P_n\hfill \\ \hfill =& \sum _{k=1}^m{v}_k\left[S_k(n+1)-S_k\left(n\right)+I_k(n+1)-I_k\left(n\right)-I_k^{*}ln{\displaystyle \frac{I_k(n+1)}{I_k\left(n\right)}}\right.\hfill \\ & \left.-{\int }_{S_k^{*}}^{S_k\left(n\right)}{\displaystyle \frac{f_k\left(S_k^{*}\right)}{f_k\left(t\right)}}dt+〈k〉\sum _{j=1}^m{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)\left(\mathsf{\Psi }\left({\displaystyle \frac{I_k(n+1)}{I_k^{*}}}\right)-\mathsf{\Psi }\left({\displaystyle \frac{I_k\left(n\right)}{I_k^{*}}}\right)\right)\right],\hfill \\ \hfill \le & \sum _{k=1}^m{v}_k\left[\left(1-{\displaystyle \frac{f_k\left(S_k^{*}\right)}{f_k\left(S_k(n+1)\right)}}\right)\left(S_k(n+1)-S_k\left(n\right)\right)+\left(1-{\displaystyle \frac{I_k^{*}}{I_k(n+1)}}\right)\left(I_k(n+1)\right.\right.\hfill \end{array}$

$$\begin{array}{c}\left.-I_k\left(n\right)\right)+\left.〈k〉\sum _{j=1}^m{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)\left(\mathsf{\Psi }\left({\displaystyle \frac{I_k(n+1)}{I_k^{*}}}\right)-\mathsf{\Psi }\left({\displaystyle \frac{I_k\left(n\right)}{I_k^{*}}}\right)\right)\right].\hfill \end{array}$$

Since

$$\left\{\begin{array}{c}{e}_k=b_k^S{S}_k^{*}+〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right),\hfill \\ 〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)=\left(b_k^I+{\gamma }_k\right)I_k^{*}.\hfill \end{array}\right.$$

It can be further acquired that,

$\begin{array}{cc}& P_{n+1}-P_n\hfill \\ \hfill \le & {\displaystyle \sum _{k=1}^m}{v}_k\left\{\left(1-{\displaystyle \frac{f_k\left(S_k^{*}\right)}{f_k\left(S_k(n+1)\right)}}\right)\left[b_k^S\left(S_k^{*}-S_k(n+1)\right)+〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)\right.\right.\hfill \\ & -\left.〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k(n+1)\right)g_j\left(I_j\left(n\right)\right)\right]+\left(1-{\displaystyle \frac{I_k^{*}}{I_k(n+1)}}\right)\times \hfill \end{array}$

$\begin{array}{cc}& \left[〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k(n+1)\right)g_j\left(I_j\left(n\right)\right)-\left(b_k^I+{\gamma }_k\right)I_k(n+1)\right]\hfill \\ & \left.+〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)\left(\mathsf{\Psi }\left({\displaystyle \frac{I_k(n+1)}{I_k^{*}}}\right)-\mathsf{\Psi }\left({\displaystyle \frac{I_k\left(n\right)}{I_k^{*}}}\right)\right)\right\},\hfill \\ \hfill =& {\displaystyle \sum _{k=1}^m}{v}_k\{{\displaystyle \frac{b_k^S}{f_k\left(S_k(n+1)\right)}}\left(S_k^{*}-S_k(n+1)\right)\left[f_k\left(S_k(n+1)\right)-f_k\left(S_k^{*}\right)\right]\hfill \\ & +〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)[2-{\displaystyle \frac{f_k\left(S_k^{*}\right)}{f_k\left(S_k(n+1)\right)}}+{\displaystyle \frac{g_j\left(I_j\left(n\right)\right)}{g_j\left(I_j^{*}\right)}}\hfill \\ & -{\displaystyle \frac{f_k\left(S_k(n+1)\right)g_j\left(I_j\left(n\right)\right)I_k^{*}}{f_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)I_k(n+1)}}-{\displaystyle \frac{I_k\left(n\right)}{I_k^{*}}}+ln{\displaystyle \frac{I_k^{*}}{I_k(n+1)}}+ln{\displaystyle \frac{I_k\left(n\right)}{I_k^{*}}}\left]\right\},\hfill \\ \hfill =& {\displaystyle \sum _{k=1}^m}{v}_k\left\{{\displaystyle \frac{b_k^S}{f_k\left(S_k(n+1)\right)}}\left(S_k^{*}-S_k(n+1)\right)\left[f_k\left(S_k(n+1)\right)-f_k\left(S_k^{*}\right)\right]\right.\hfill \\ & +〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)\left[-\mathsf{\Psi }\left({\displaystyle \frac{f_k\left(S_k^{*}\right)}{f_k\left(S_k(n+1)\right)}}\right)-\mathsf{\Psi }\left({\displaystyle \frac{f_k\left(S_k(n+1)\right)g_j\left(I_j\left(n\right)\right)I_k^{*}}{f_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)I_k(n+1)}}\right)\right.\hfill \\ & \left.\left.+ln{\displaystyle \frac{g_j\left(I_j^{*}\right)}{g_j\left(I_j\left(n\right)\right)}}+{\displaystyle \frac{g_j\left(I_j\left(n\right)\right)}{g_j\left(I_j^{*}\right)}}-{\displaystyle \frac{I_k\left(n\right)}{I_k^{*}}}+ln{\displaystyle \frac{I_k\left(n\right)}{I_k^{*}}}\right]\right\},\hfill \\ \hfill =& {\displaystyle \sum _{k=1}^m}{v}_k\left\{{\displaystyle \frac{b_k^S}{f_k\left(S_k(n+1)\right)}}\left(S_k^{*}-S_k(n+1)\right)\left[f_k\left(S_k(n+1)\right)-f_k\left(S_k^{*}\right)\right]\right.\hfill \\ & +〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)\left[-\mathsf{\Psi }\left({\displaystyle \frac{f_k\left(S_k^{*}\right)}{f_k\left(S_k(n+1)\right)}}\right)-\mathsf{\Psi }\left({\displaystyle \frac{f_k\left(S_k(n+1)\right)g_j\left(I_j\left(n\right)\right)I_k^{*}}{f_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)I_k(n+1)}}\right)\right.\hfill \\ & +ln{\displaystyle \frac{g_j\left(I_j^{*}\right)}{g_j\left(I_j\left(n\right)\right)}}+\left({\displaystyle \frac{g_j\left(I_j\left(n\right)\right)}{g_j\left(I_j^{*}\right)}}-{\displaystyle \frac{I_j\left(n\right)}{I_j^{*}}}\right)\left(1-{\displaystyle \frac{g_j\left(I_j^{*}\right)}{g_j\left(I_j\left(n\right)\right)}}\right)+1-{\displaystyle \frac{I_j\left(n\right)g_j\left(I_j^{*}\right)}{I_j^{*}{g}_j\left(I_j\left(n\right)\right)}}\hfill \\ & \left.\left.+{\displaystyle \frac{I_j\left(n\right)}{I_j^{*}}}-{\displaystyle \frac{I_k\left(n\right)}{I_k^{*}}}+ln{\displaystyle \frac{I_k\left(n\right)}{I_k^{*}}}\right]\right\},\hfill \end{array}$

$\begin{array}{cc}\hfill =& {\displaystyle \sum _{k=1}^m}{v}_k\left\{{\displaystyle \frac{b_k^S}{f_k\left(S_k(n+1)\right)}}\left(S_k^{*}-S_k(n+1)\right)\left[f_k\left(S_k(n+1)\right)-f_k\left(S_k^{*}\right)\right]\right.\hfill \\ & +〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)\left[-\mathsf{\Psi }\left({\displaystyle \frac{f_k\left(S_k^{*}\right)}{f_k\left(S_k(n+1)\right)}}\right)-\mathsf{\Psi }\left({\displaystyle \frac{f_k\left(S_k(n+1)\right)g_j\left(I_j\left(n\right)\right)I_k^{*}}{f_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)I_k(n+1)}}\right)\right.\hfill \\ & -\mathsf{\Psi }\left({\displaystyle \frac{I_j\left(n\right)g_j\left(I_j^{*}\right)}{I_j^{*}{g}_j\left(I_j\left(n\right)\right)}}\right)+\left({\displaystyle \frac{g_j\left(I_j\left(n\right)\right)}{g_j\left(I_j^{*}\right)}}-{\displaystyle \frac{I_j\left(n\right)}{I_j^{*}}}\right)\left(1-{\displaystyle \frac{g_j\left(I_j^{*}\right)}{g_j\left(I_j\left(n\right)\right)}}\right)+{\displaystyle \frac{I_j\left(n\right)}{I_j^{*}}}-{\displaystyle \frac{I_k\left(n\right)}{I_k^{*}}}\hfill \\ & \left.\left.+ln{\displaystyle \frac{I_j^{*}{I}_k\left(n\right)}{I_j\left(n\right)I_k^{*}}}\right]\right\},\hfill \end{array}$

Set

$$V_1={\displaystyle \sum _{k=1}^m}{\displaystyle \sum _{j=1}^m}〈k〉v_k{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(l_j^{*}\right)\left({\displaystyle \frac{I_j\left(n\right)}{I_j^{*}}}-{\displaystyle \frac{I_k\left(n\right)}{I_k^{*}}}\right),$$

$$V_2={\displaystyle \sum _{k=1}^m}{\displaystyle \sum _{j=1}^m}〈k〉v_k{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)ln{\displaystyle \frac{I_j^{*}{I}_k\left(n\right)}{I_j\left(n\right)I_k^{*}}},$$

we can get $P_{n+1}-P_n\le V_1+V_2$. Next it will be shown that for all $I_1\left(n\right),I_2\left(n\right),\cdots ,$ $I_m\left(n\right)>0,V_1=0$, it follows that ${\displaystyle \sum _{j=1}^m}{\overline{\beta }}_{kj}{v}_j={\displaystyle \sum _{j=1}^m}{\overline{\beta }}_{kj}{v}_k$, from $\overline{B}v=0$.

Because

$${\overline{\beta }}_{kj}={\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right),$$

ones has

$$〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{jk}{f}_j\left(S_j^{*}\right)g_k\left(I_k^{*}\right)v_j=〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)v_k,$$

that’s why the introduction of the

$$\begin{array}{cc}& 〈k〉\sum _{k,j=1}^m{v}_k{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right){\displaystyle \frac{I_j\left(n\right)}{I_j^{*}}}\hfill \\ \hfill =& \sum _{k=1}^m{\displaystyle \frac{I_k\left(n\right)}{I_k^{*}}}〈k〉\sum _{j=1}^m{\beta }_{jk}{f}_j\left(S_j^{*}\right)g_k\left(I_k^{*}\right)v_j,\hfill \\ \hfill =& \sum _{k=1}^m{\displaystyle \frac{I_k\left(n\right)}{I_k^{*}}}〈k〉\sum _{j=1}^m{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)v_k,\hfill \\ \hfill =& 〈k〉\sum _{k,j=1}^m{v}_k{\beta }_{kj}{f}_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right){\displaystyle \frac{I_k\left(n\right)}{I_k^{*}}}.\hfill \end{array}$$

Thus, $V_1=0$ for all $I_1\left(n\right),I_2\left(n\right),\cdots ,I_m\left(n\right)>0$, similar to the proof of Theorem 4.2 of the midpoint [27], $V_2=0$ can be deduced for all $I_1\left(n\right),I_2\left(n\right),\cdots ,I_m\left(n\right)>0$. So, $P_{n+1}-P_n\le 0$. This shows that ${\left\{P_n\right\}}_{n\in N}$ is a decreasing sequence, so there is a constant $\tilde{P}$ such that ${lim}_{n\to +\infty }{P}_n=\tilde{P}$. So we have ${lim}_{n\to +\infty }(P_{n+1}-P_n)=0$, which means

$$\underset{n\to +\infty }{lim}(S_k^{*}-S_k(n+1))(f_k\left(S_k(n+1)\right)-f_k\left(S_k^{*}\right))=0,$$

$$\underset{n\to +\infty }{lim}\left(\mathsf{\Psi }\left({\displaystyle \frac{f_k\left(S_k^{*}\right)}{f_k\left(S_k(n+1)\right)}}\right)+\mathsf{\Psi }\left({\displaystyle \frac{f_k\left(S_k(n+1)\right)g_j\left(I_j\left(n\right)\right)I_k^{*}}{f_k\left(S_k^{*}\right)g_j\left(I_j^{*}\right)I_k(n+1)}}\right)+\mathsf{\Psi }\left({\displaystyle \frac{I_j\left(n\right)g_j\left(I_j^{*}\right)}{I_j^{*}{g}_j\left(I_j\left(n\right)\right)}}\right))=0,\right.$$

$$\underset{n\to +\infty }{lim}\left({\displaystyle \frac{g_j\left(I_j\left(n\right)\right)}{g_j\left(I_j^{*}\right)}}-{\displaystyle \frac{I_j\left(n\right)}{I_j^{*}}}\right)\left(1-{\displaystyle \frac{g_j\left(I_j^{*}\right)}{g_j\left(I_j\left(n\right)\right)}}\right)=0.$$

Hence, it can be concluded that

$$\underset{n\to +\infty }{lim}{S}_k\left(n\right)=S_k^{*}\mathrm{and}\underset{n\to +\infty }{lim}{I}_j\left(n\right)=I_j^{*},\mathrm{for}j=1,2,\cdots ,m.$$

Based on the decreasing trend of $P_n$ and the LaSalle’s invariance principle, it can be concluded that the rumor-prevailing equilibrium $E^{*}$ is globally asymptotically stable in $\mathsf{\Omega }$ when ${\mathcal{R}}_0>1$. □

4. Optimal Control

How to take control when rumors are prevailing is a very important issue, for this reason, this paper considers two control strategies, including prevention mechanism and rumor refutation mechanism. Specially, the prevention mechanism refers to the dissemination of the civilized concept of “don’t believe in rumors, don’t spread rumors” by holding various online public lecture on rumor refutation, so as to improve their ability to distinguish. The rumor-refuting mechanism is government control and the establishment of a rumor-refuting mechanism to stop the spread of rumors.

We assume that there are m languages of the domain studied $\mathsf{\Gamma }$, where denotes ${\mathsf{\Gamma }}_k$ the density of speakers of the kth language. The ${\mathsf{\Gamma }}_k$ population was divided into three compartments, $S^{{\mathsf{\Gamma }}_k}\left(n\right),I^{{\mathsf{\Gamma }}_k}\left(n\right),$ and $R^{{\mathsf{\Gamma }}_k}\left(n\right)$.

4.1. Presentation of the Model

The multilingual discrete time SIR model associated with ${\mathsf{\Gamma }}_k$ can be derived from (4) as follows:

$$\begin{array}{c}\left\{\begin{array}{c}{S}^{{\mathsf{\Gamma }}_k}(n+1)=S^{{\mathsf{\Gamma }}_k}\left(n\right)+e_k-b_k^S{S}^{{\mathsf{\Gamma }}_k}\left(n\right)-〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S^{{\mathsf{\Gamma }}_k}\left(n\right)\right)g_j\left(I^{{\mathsf{\Gamma }}_j}\left(n\right)\right),\hfill \\ I^{{\mathsf{\Gamma }}_k}(n+1)=I^{{\mathsf{\Gamma }}_k}\left(n\right)+〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S^{{\mathsf{\Gamma }}_k}\left(n\right)\right)g_j\left(I^{{\mathsf{\Gamma }}_j}\left(n\right)\right)-\left(b_k^I+{\gamma }_k\right)I^{{\mathsf{\Gamma }}_k}\left(n\right),\hfill \\ R^{{\mathsf{\Gamma }}_k}(n+1)=R^{{\mathsf{\Gamma }}_k}\left(n\right)+{\gamma }_k{I}^{{\mathsf{\Gamma }}_k}\left(n\right)-b_k^R{R}^{{\mathsf{\Gamma }}_k}\left(n\right).\hfill \end{array}\right.\hfill \end{array}$$

In the model, we introduce two control variables $u^{{\mathsf{\Gamma }}_k}\left(n\right)$ and $v^{{\mathsf{\Gamma }}_k}\left(n\right)$ to characterize the prevention mechanism and the rumor refutation mechanism. Then, for a given group ${\mathsf{\Gamma }}_k$, the model can be further shown as the following equation:

$$\left\{\begin{array}{c}{S}^{{\mathsf{\Gamma }}_k}(n+1)=S^{{\mathsf{\Gamma }}_k}\left(n\right)+e_k-b_k^S{S}^{{\mathsf{\Gamma }}_k}\left(n\right)-〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S^{{\mathsf{\Gamma }}_k}\left(n\right)\right)g_j\left(I^{{\mathsf{\Gamma }}_j}\left(n\right)\right)-u^{{\mathsf{\Gamma }}_k}\left(n\right)S^{{\mathsf{\Gamma }}_k}\left(n\right),\hfill \\ I^{{\mathsf{\Gamma }}_k}(n+1)=I^{{\mathsf{\Gamma }}_k}\left(n\right)+〈k〉{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S^{{\mathsf{\Gamma }}_k}\left(n\right)\right)g_j\left(I^{{\mathsf{\Gamma }}_j}\left(n\right)\right)-\left(b_k^I+{\gamma }_k\right)I^{{\mathsf{\Gamma }}_k}\left(n\right)-v^{{\mathsf{\Gamma }}_k}\left(n\right)I^{{\mathsf{\Gamma }}_k}\left(n\right),\hfill \\ R^{{\mathsf{\Gamma }}_k}(n+1)=R^{{\mathsf{\Gamma }}_k}\left(n\right)+{\gamma }_k{I}^{{\mathsf{\Gamma }}_k}\left(n\right)-b_k^R{R}^{{\mathsf{\Gamma }}_k}\left(n\right)+u^{{\mathsf{\Gamma }}_k}\left(n\right)S^{{\mathsf{\Gamma }}_k}\left(n\right)+v^{{\mathsf{\Gamma }}_k}\left(n\right)I^{{\mathsf{\Gamma }}_k}\left(n\right).\hfill \end{array}\right.$$

(6)

For this control problem, our goal is to reduce the density of the ignorants $S^{{\mathsf{\Gamma }}_k}\left(n\right)$ and the spreaders $I^{{\mathsf{\Gamma }}_k}\left(n\right)$ with the lowest cost, and try to increase the density of removers individuals in ${\mathsf{\Gamma }}_k$. Suppose that the control variables $u^{{\mathsf{\Gamma }}_k}\left(n\right)$ and $v^{{\mathsf{\Gamma }}_k}\left(n\right)$ are bounded. Define the maximum value of $u^{{\mathsf{\Gamma }}_k}\left(n\right)$ and $v^{{\mathsf{\Gamma }}_k}\left(n\right)$ is $u^{{\mathsf{\Gamma }}_k}\left(max\right)$ and $v^{{\mathsf{\Gamma }}_k}\left(max\right)$, the minimum value of $u^{{\mathsf{\Gamma }}_k}\left(n\right)$ and $v^{{\mathsf{\Gamma }}_k}\left(n\right)$ is $u^{{\mathsf{\Gamma }}_k}\left(min\right)$ and $v^{{\mathsf{\Gamma }}_k}\left(min\right)$, and satisfies $0<u^{{\mathsf{\Gamma }}_k}\left(min\right)<u^{{\mathsf{\Gamma }}_k}\left(max\right)<1$ and $0<v^{{\mathsf{\Gamma }}_k}\left(min\right)<v^{{\mathsf{\Gamma }}_k}\left(max\right)<1$, where $k=1,2,\dots ,m$.

4.2. Optimal Control Strategy

Based on the above description, we can get the objective function of model (6) is

$$\begin{array}{cc}\hfill J_k\left(u^{{\mathsf{\Gamma }}_k},v^{{\mathsf{\Gamma }}_k}\right)& =\left({\alpha }_k^S{S}^{{\mathsf{\Gamma }}_k}\left(N\right)+{\alpha }_k^I{I}^{{\mathsf{\Gamma }}_k}\left(N\right)-{\alpha }_k^R{R}^{{\mathsf{\Gamma }}_k}\left(N\right)\right)\hfill \\ \hfill & +{\displaystyle \sum _{k=1}^m}{\displaystyle \sum _{n=1}^{N-1}}\left[{\alpha }_k^S{S}^{{\mathsf{\Gamma }}_k}\left(n\right)+{\alpha }_k^I{I}^{{\mathsf{\Gamma }}_k}\left(n\right)-{\alpha }_k^R{R}^{{\mathsf{\Gamma }}_k}\left(n\right)+{\displaystyle \frac{A_k}{2}}{\left(u^{{\mathsf{\Gamma }}_k}\left(n\right)\right)}^2+{\displaystyle \frac{B_k}{2}}{\left(v^{{\mathsf{\Gamma }}_k}\left(n\right)\right)}^2\right],\hfill \end{array}$$

(7)

where $A_k>0, B_k>0, {\alpha }_k^S>0, {\alpha }_k^I>0, {\alpha }_k^R>0$ are the weight coefficient. In other words, we are looking for an optimal control $u^{{\mathsf{\Gamma }}_k^{*}}\mathrm{and}{v}^{{\mathsf{\Gamma }}_k^{*}}$ such that,

$$\begin{array}{c}{J}_k\left(u^{{\mathsf{\Gamma }}_k^{*}},v^{{\mathsf{\Gamma }}_k^{*}}\right)=min\left\{J_k\left(u^{{\mathsf{\Gamma }}_k},v^{{\mathsf{\Gamma }}_k}\right)\left|u^{{\mathsf{\Gamma }}_k}\in U_k,v^{{\mathsf{\Gamma }}_k}\in V_k\right.\right\},\hfill \end{array}$$

where $U_k$ and $V_k$ is defined by the control set

$$\begin{array}{c}\left\{\begin{array}{c}\begin{array}{cc}\hfill U_k& =\left\{u^{{\mathsf{\Gamma }}_k}\left|u^{{\mathsf{\Gamma }}_k}\left(min\right)\le u^{{\mathsf{\Gamma }}_k}\le u^{{\mathsf{\Gamma }}_k}\left(max\right)\right.,n=0,1,\cdots ,N-1\right\},\hfill \\ \hfill V_k& =\left\{v^{{\mathsf{\Gamma }}_k}\left|v^{{\mathsf{\Gamma }}_k}\left(min\right)\le v^{{\mathsf{\Gamma }}_k}\le u^{{\mathsf{\Gamma }}_k}\left(max\right)\right.,n=0,1,\cdots ,N-1\right\},\hfill \end{array}\hfill \end{array}\right.\hfill \end{array}$$

(8)

where $0<u^{{\mathsf{\Gamma }}_k}\left(min\right)<u^{{\mathsf{\Gamma }}_k}\left(max\right)<1$ and $0<v^{{\mathsf{\Gamma }}_k}\left(min\right)<v^{{\mathsf{\Gamma }}_k}\left(max\right)<1$.

Therefore, we define the Hamiltonian function as:

$$\begin{array}{cc}\hfill H\left({\mathsf{\Gamma }}_k\right)& ={\displaystyle \sum _{k=1}^m}{\alpha }_k^S{S}^{{\mathsf{\Gamma }}_k}\left(N\right)+{\displaystyle \sum _{k=1}^m}{\alpha }_k^I{I}^{{\mathsf{\Gamma }}_k}\left(N\right)-{\alpha }_k^R{R}^{{\mathsf{\Gamma }}_k}\left(N\right)+{\displaystyle \frac{A_n}{2}}{\left(u^{{\mathsf{\Gamma }}_k}\left(n\right)\right)}^2+{\displaystyle \frac{B_n}{2}}{\left(v^{{\mathsf{\Gamma }}_k}\left(n\right)\right)}^2\hfill \\ & +{\zeta }_1^k(n+1)\times [S^{{\mathsf{\Gamma }}_k}\left(n\right)+e_k-b_k^S{S}^{{\mathsf{\Gamma }}_k}\left(n\right)-{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S^{{\mathsf{\Gamma }}_k}\left(n\right)\right)g_j\left(I^{{\mathsf{\Gamma }}_j}\left(n\right)\right)-u^{{\mathsf{\Gamma }}_k}\left(n\right)S^{{\mathsf{\Gamma }}_k}\left(n\right)]\hfill \\ & +{\zeta }_2^k(n+1)\times [I^{{\mathsf{\Gamma }}_k}\left(n\right)+{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S^{{\mathsf{\Gamma }}_k}\left(n\right)\right)g_j\left(I^{{\mathsf{\Gamma }}_j}\left(n\right)\right)-(b_k^I+{\gamma }_k)I^{{\mathsf{\Gamma }}_k}\left(n\right)-v^{{\mathsf{\Gamma }}_k}\left(n\right)I^{{\mathsf{\Gamma }}_k}\left(n\right)]\hfill \\ & +{\zeta }_3^k(n+1)[R^{{\mathsf{\Gamma }}_k}\left(n\right)+{\gamma }_k{I}^{{\mathsf{\Gamma }}_k}\left(n\right)-b_k^R{R}^{{\mathsf{\Gamma }}_k}\left(n\right)+u^{{\mathsf{\Gamma }}_k}\left(n\right)S^{{\mathsf{\Gamma }}_k}\left(n\right)+v^{{\mathsf{\Gamma }}_k}\left(n\right)I^{{\mathsf{\Gamma }}_k}\left(n\right)]\hfill \end{array}$$

Theorem 4

(sufficient condition).  For the optimal control problem given by (7) and equation of state (6), there exists an optimal control $u^{{\mathsf{\Gamma }}_k^{*}}\in U_k$ such that

$$\begin{array}{c}{J}_k\left(u^{{\mathsf{\Gamma }}_k^{*}},v^{{\mathsf{\Gamma }}_k^{*}}\right)=min\left\{J_k\left(u^{{\mathsf{\Gamma }}_k},v^{{\mathsf{\Gamma }}_k}\right)\left|u^{{\mathsf{\Gamma }}_k}\in U_k,v^{{\mathsf{\Gamma }}_k}\in V_k\right.\right\}.\hfill \end{array}$$

Proof of Theorem 4.

The specific proof steps are similar to Theorem 1 in reference [28], here, we no longer repeat. □

Theorem 5

(necessary condition). Given two optimal controls $u^{{\mathsf{\Gamma }}_k^{*}}$ $v^{{\mathsf{\Gamma }}_k^{*}}$ and solution $S^{{\mathsf{\Gamma }}_k^{*}}$, $I^{{\mathsf{\Gamma }}_k^{*}}$, and $R^{{\mathsf{\Gamma }}_k^{*}}$, there exists ${\zeta }_{j,k}$ that the accompanying variables satisfy the following equations: ($k=1,\cdots ,N;j=1,2,3.$)

$$\begin{array}{c}\left\{\begin{array}{c}\begin{array}{cc}\hfill {\zeta }_1^k\left(n\right)=& -[{\alpha }_k^S+{\zeta }_1^k(n+1)(1-b_k^S-{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{g}_j\left(I^{{\mathsf{\Gamma }}_j}\left(n\right)\right)f_k^{\prime }\left(S^{{\mathsf{\Gamma }}_k}\left(n\right)\right)-u^{{\mathsf{\Gamma }}_k}\left(n\right))\hfill \\ & +{\zeta }_2^k(n+1)({\displaystyle \sum _{j=1}^m}{\beta }_{kj}{g}_j\left(I^{{\mathsf{\Gamma }}_j}\left(n\right)\right)f_k^{\prime }\left(S^{{\mathsf{\Gamma }}_k}\left(n\right)\right))+{\zeta }_3^k(n+1)({\mu }_k+u^{{\mathsf{\Gamma }}_k}\left(n\right))],\hfill \\ \hfill {\zeta }_2^k\left(n\right)=& -[{\alpha }_k^I+{\zeta }_1^k(n+1)\times (-{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S^{{\mathsf{\Gamma }}_k}\left(n\right)\right)g_j^{\prime }\left(I^{{\mathsf{\Gamma }}_j}\left(n\right)\right))+{\zeta }_2^k(n+1)\times \hfill \\ & (1+{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S^{{\mathsf{\Gamma }}_k}\left(n\right)\right)g_j^{\prime }\left(I^{{\mathsf{\Gamma }}_j}\left(n\right)\right)-b_k^I-{\gamma }_k-v^{{\mathsf{\Gamma }}_k}\left(n\right))+{\zeta }_3^k(n+1)({\gamma }_k+v^{{\mathsf{\Gamma }}_k}\left(n\right))],\hfill \\ \hfill {\zeta }_3^k\left(n\right)=& -[-{\alpha }_k^R+{\zeta }_3^k(n+1)(1-b_k^R)],\hfill \end{array}\hfill \end{array}\right.\hfill \end{array}$$

where ${\zeta }_1^k\left(N\right)={\alpha }_k^S,{\zeta }_2^k\left(N\right)={\alpha }_k^I,{\zeta }_3^k\left(N\right)=-{\alpha }_k^R$ are transversal conditions, in addition, there are

$$\begin{array}{c}\left\{\begin{array}{c}{u}^{{\mathsf{\Gamma }}_k^{*}}\left(n\right)=min(max(u^{{\mathsf{\Gamma }}_k}\left(min\right),{\displaystyle \frac{({\zeta }_1^k(n+1)-{\zeta }_3^k(n+1))S^{{\mathsf{\Gamma }}_k^{*}}\left(n\right)}{A_k}},u^{{\mathsf{\Gamma }}_k}\left(max\right))),\hfill \\ v^{{\mathsf{\Gamma }}_k^{*}}\left(n\right)=min(max(v^{{\mathsf{\Gamma }}_k}\left(min\right),{\displaystyle \frac{({\zeta }_2^k(n+1)-{\zeta }_3^k(n+1))I^{{\mathsf{\Gamma }}_k^{*}}\left(n\right)}{B_k}},v^{{\mathsf{\Gamma }}_k}\left(max\right))),\hfill \end{array}\right.\hfill \end{array}$$

$n=0,1,\cdots ,N-1$.

Proof of Theorem 5.

Let $S^{{\mathsf{\Gamma }}_k}=S^{{\mathsf{\Gamma }}_k^{*}},I^{{\mathsf{\Gamma }}_k}=I^{{\mathsf{\Gamma }}_k^{*}},R^{{\mathsf{\Gamma }}_k}=R^{{\mathsf{\Gamma }}_k^{*}}$ and $u^{{\mathsf{\Gamma }}_k}=u^{{\mathsf{\Gamma }}_k^{*}}$, based on the Pontryagin maximum principle [29], we get the adjoint equation:

$$\begin{array}{c}\left\{\begin{array}{c}\begin{array}{cc}\hfill {\zeta }_1^k\left(n\right)=-{\displaystyle \frac{\partial H}{\partial S^{{\mathsf{\Gamma }}_k}\left(n\right)}}& =-[{\alpha }_k^S+{\zeta }_1^k(n+1)(1-b_k^S-{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{g}_j\left(I^{{\mathsf{\Gamma }}_j}\left(n\right)\right)f_k^{\prime }\left(S^{{\mathsf{\Gamma }}_k}\left(n\right)\right)-u^{{\mathsf{\Gamma }}_k}\left(n\right))\hfill \\ & +{\zeta }_2^k(n+1)({\displaystyle \sum _{j=1}^m}{\beta }_{kj}{g}_j\left(I^{{\mathsf{\Gamma }}_j}\left(n\right)\right)f_k^{\prime }\left(S^{{\mathsf{\Gamma }}_k}\left(n\right)\right))+{\zeta }_3^k(n+1)({\mu }_k+u^{{\mathsf{\Gamma }}_k}\left(n\right))],\hfill \\ \hfill {\zeta }_2^k\left(n\right)=-{\displaystyle \frac{\partial H}{\partial I^{{\mathsf{\Gamma }}_k}\left(n\right)}}& =-[{\alpha }_k^I+{\zeta }_1^k(n+1)\times (-{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S^{{\mathsf{\Gamma }}_k}\left(n\right)\right)g_j^{\prime }\left(I^{{\mathsf{\Gamma }}_j}\left(n\right)\right))\hfill \\ & +{\zeta }_2^k(n+1)(1+{\displaystyle \sum _{j=1}^m}{\beta }_{kj}{f}_k\left(S^{{\mathsf{\Gamma }}_k}\left(n\right)\right)g_j^{\prime }\left(I^{{\mathsf{\Gamma }}_j}\left(n\right)\right)-b_k^I-{\gamma }_k-v^{{\mathsf{\Gamma }}_k}\left(n\right))\hfill \\ & +{\zeta }_3^k(n+1)({\gamma }_k+v^{{\mathsf{\Gamma }}_k}\left(n\right))],\hfill \\ \hfill {\zeta }_3^k\left(n\right)=-{\displaystyle \frac{\partial H}{\partial R^{{\mathsf{\Gamma }}_k}\left(n\right)}}& =-[-{\alpha }_k^R+{\zeta }_3^k(n+1)(1-b_k^R)],\hfill \end{array}\hfill \end{array}\right.\hfill \end{array}$$

with ${\zeta }_1^k\left(N\right)={\alpha }_k^S,{\zeta }_2^k\left(N\right)={\alpha }_k^I,{\zeta }_3^k\left(N\right)=-{\alpha }_k^R$, $k=1,2,\cdots ,m$ are transversal conditions. To obtain the optimal control conditions, we take the change relative to the control $u^{{\mathsf{\Gamma }}_k},v^{{\mathsf{\Gamma }}_k}$ and set it to zero

$$\begin{array}{c}\left\{\begin{array}{c}\begin{array}{cc}\hfill {\displaystyle \frac{\partial H}{\partial u^{{\mathsf{\Gamma }}_k}\left(n\right)}}& =A_k{u}^{{\mathsf{\Gamma }}_k}\left(n\right)-{\zeta }_1^k(n+1)S^{{\mathsf{\Gamma }}_k}\left(n\right)+{\zeta }_3^k(n+1)S^{{\mathsf{\Gamma }}_k}\left(n\right)=0,\hfill \\ \hfill {\displaystyle \frac{\partial H}{\partial v^{{\mathsf{\Gamma }}_k}\left(n\right)}}& =B_k{v}^{{\mathsf{\Gamma }}_k}\left(n\right)-{\zeta }_2^k(n+1)I^{{\mathsf{\Gamma }}_k}\left(n\right)+{\zeta }_3^k(n+1)I^{{\mathsf{\Gamma }}_k}\left(n\right)=0,\hfill \end{array}\hfill \end{array}\right.\hfill \end{array}$$

And then we have

$$\begin{array}{c}\left\{\begin{array}{c}{u}^{{\mathsf{\Gamma }}_k}\left(n\right)={\displaystyle \frac{({\zeta }_1^k(n+1)-{\zeta }_3^k(n+1))S^{{\mathsf{\Gamma }}_k}\left(n\right)}{A_k}},\hfill \\ v^{{\mathsf{\Gamma }}_k}\left(n\right)={\displaystyle \frac{({\zeta }_2^k(n+1)-{\zeta }_3^k(n+1))I^{{\mathsf{\Gamma }}_k}\left(n\right)}{B_k}}.\hfill \end{array}\right.\hfill \end{array}$$

Then, we combine the Equation (7) to obtain the optimal solution as follows

$$\begin{array}{c}\left\{\begin{array}{c}{u}^{{\mathsf{\Gamma }}_k^{*}}\left(n\right)=min(max(u^{{\mathsf{\Gamma }}_k}\left(min\right),{\displaystyle \frac{({\zeta }_1^k(n+1)-{\zeta }_3^k(n+1))S^{{\mathsf{\Gamma }}_k}\left(n\right)}{A_k}},u^{{\mathsf{\Gamma }}_k}\left(max\right)))\hfill \\ v^{{\mathsf{\Gamma }}_k^{*}}\left(n\right)=min(max(v^{{\mathsf{\Gamma }}_k}\left(min\right),{\displaystyle \frac{({\zeta }_2^k(n+1)-{\zeta }_3^k(n+1))I^{{\mathsf{\Gamma }}_k}\left(n\right)}{B_k}},v^{{\mathsf{\Gamma }}_k}\left(max\right)))\hfill \end{array}\right.\hfill \end{array}$$

where k = 0, 1, …, N − 1. □

5. Numerical Simulations

In this section, we will perform numerical simulations for model (4) with several classical nonlinear spreading rates, such as the bilinear diffusion function $\beta SI$, the saturated diffusion function ${\displaystyle \frac{\beta SI}{1+\alpha I}}$ [30], and the nonmonotonic diffusion function ${\displaystyle \frac{\beta SI}{1+\alpha I^2}}$ [31], to demonstrate the reliability of the theoretical results. The initial value of model (4) is chosen as $S_i\left(t\right)=0.7$, $I_i\left(t\right)=0.2$, $R_i\left(t\right)=0.1$. Without losing generality, we choose the bilinear diffusion function $F_{kj}\left(S_k\left(n\right),I_j\left(n\right)\right)=f_k\left(S\right)g_k\left(I\right)=SI$, the saturated diffusion function $F_{kj}\left(S_k\left(n\right),I_j\left(n\right)\right)=f_k\left(S\right)g_k\left(I\right)=S{\displaystyle \frac{S}{1+I}}$ and the nonmonotonic diffusion function $F_{kj}\left(S_k\left(n\right),I_j\left(n\right)\right)=f_k\left(S\right)g_k\left(I\right)=S{\displaystyle \frac{S}{1+I^2}}$ satisfying the model (4) to conduct some numerical simulations in bi-lingual environments.

For the first case, we choose $〈k〉=1$, $e_k=0.24$, $b_k^S=0.15$, $b_k^I=0.3$, $b_k^R=0.5$, ${\beta }_{11}=0.15$, ${\beta }_{12}=0.1$, ${\beta }_{21}=0.15$, ${\beta }_{22}=0.15$, ${\gamma }_k=0.15$ in the model (4) with the bilinear diffusion function, the saturated diffusion function and the nonmonotonic diffusion function, respectively. With the parameter configuration in this case, it follows that $R_0\le 1$. The numerical simulations of the bilinear diffusion function, saturated diffusion function, and nonmonotonic diffusion function respectively correspond to Figure 2, Figure 3 and Figure 4.

Through observing Figure 2, Figure 3 and Figure 4, we find that in all three cases, regardless of the specific from of the nonlinear diffusion function, rumors will prevail in the network when $R_0\le 1$. It verifies that the result of Theorem 2 is correct.

For the second case, we choose $〈k〉=1$, $e_k=0.16$, $b_k^S=0.05$, $b_k^I=0.05$, $b_k^R=0.5$, ${\beta }_{11}=0.4$, ${\beta }_{12}=0.3$, ${\beta }_{21}=0.4$, ${\beta }_{22}=0.5$, ${\gamma }_k=0.3$ in the model (4) with the bilinear diffusion function, the saturated diffusion function and the nonmonotonic diffusion function, respectively. With the parameter configuration in this case, it follows that $R_0>1$. The numerical simulations of the bilinear diffusion function, saturated diffusion function, and nonmonotonic diffusion function respectively correspond to Figure 5, Figure 6 and Figure 7.

By observing Figure 5, Figure 6 and Figure 7, we find that in all three cases, regardless of the specific from of the nonlinear diffusion function, rumors will prevail in the network when $R_0>1$. The validity of the result of Theorem 2 is shown.

Remark 2.

From Figure 5, Figure 6 and Figure 7, we can sell that the propagation rate of the bilinear diffusion function is the fastest and reaches the peak first. After reaching equilibrium, the recovery density of the bilinear diffusion function is also the largest. The saturated diffusion function has the smallest recovery density and the smallest propagation density among the three functions.

Remark 3.

As we all know, the purpose of scholars analyzing the dynamic behavior of rumor spreading is to guide and control the spread of rumors. As can be seen from Figure 5, Figure 6 and Figure 7, the number of spreaders $I\left(t\right)$ will peak somewhere, indicating that at this time we can take measures to effectively prevent the spread of rumors.

We now give numerical simulations related to the above optimal control problem. For example, to show the significance of our work, we choose m = 40, i.e., we consider forty groups, and we try to control $\mathsf{\Omega }$ with control variables $u^{{\mathsf{\Gamma }}_k^{*}}$ and $v^{{\mathsf{\Gamma }}_k^{*}}$. Let $〈k〉=1$, $e_k=0.16$, $b_k^S=0.05$, $b_k^I=0.05$, $b_k^R=0.5$, ${\beta }_{11}=0.4$, ${\beta }_{12}=0.3$, ${\beta }_{21}=0.4$, ${\beta }_{22}=0.5$, ${\gamma }_k=0.3$ in the model (4) with the bilinear diffusion function. We obtain the change trend diagram of the states with bilinear diffusion function before control in Figure 8. Figure 9 is the change trend diagram of each state with bilinear diffusion function after control.

It can be seen from Figure 8 and Figure 9 that the spreader $I\left(t\right)$ is controlled, which better illustrates that the theoretically proposed strategy is effective and proves the applicability of the model.

6. Numerical Simulation of the Case

In this section, we give a real bilingual rumor case for numerical simulation and verify the applicability of the proposed model.

The main content of the rumor is that in the early morning of 26 August 2018, a multi-person brawl broke out in the city of Chemnitz, Germany. At the same time, there is a lot of false information on the network. For example, there are rumors that the perpetrator attacked the victim because he wanted to protect women from sexual assault that the perpetrator did not want. The incident involved people of multiple nationalities, so the rumors related to the incident were spread in many languages, mainly in English and German.

We collected rumor spreading data within 217 h after the event from the literature [32], and used 217 h of data as a test set. In order to fit the real data as accurately as possible, we use Matlab R2017a to simulate model (4) to evaluate the parameter values. On this basis, the parameter values shown in Figure 10 and Figure 11 are obtained. We use red and yellow to represent the true density curves of English and German rumor spreaders, respectively. Figure 10 shows the density evolution curve of English rumor spreaders represented by blue. It is found that the peaks are basically the same and the fitting effect is good. Blue is used to represent the density evolution curve of German rumor spreaders, as shown in Figure 11. We found that in this case, especially in the second half, the fitting effect is not good, but it also shows that the rumor is in a popular state. The reason may be due to the influence of random factors, which also inspires us to establish a random model to further explore this issue in the future. It is worth noting that since the incident occurred in Germany, rumors usually occur in the spread of Germany.

The first case is modeling the data in English. Let $〈k〉=6$, $e_k=0.4$, $b_k^S=0.001$, $b_k^I=0.8$, $b_k^R=0.22$, ${\beta }_{11}=0.2$, ${\beta }_{12}=0.1$, ${\beta }_{21}=0.1$, ${\beta }_{22}=0.2$, ${\gamma }_k=0.3$ in the model (4) with the bilinear diffusion function.

The second case is modeling the data in German. Let $〈k〉=5$, $e_k=0.6$, $b_k^S=0.001$, $b_k^I=0.2$, $b_k^R=0.2$, ${\beta }_{11}=0.2$, ${\beta }_{12}=0.2$, ${\beta }_{21}=0.02$, ${\beta }_{22}=0.12$, ${\gamma }_k=0.3$ in the model (4) with the bilinear diffusion fuction.

Remark 4.

Note that the previous multilingual rumor propagation models used single-language cases for case simulation [16,24,33], which lacks representativeness and cannot reflect the generalization ability of the model. It is worth affirming that our case is multilingual, and the data and the model are consistent in discrete forms, which enhances the reliability of the data and better demonstrates the usefulness of the model.

7. Conclusions

In this paper, a discrete multilingual SIR Rumor propagation model with general nonlinear propagation rate is established. The results show that when $R_0\le 1$, the rumor-free equilibrium point is globally asymptotically stable. When $R_0>1$, the rumor-prevailing equilibrium is globally asymptotically stable. Secondly, the optimal control problem of discrete multilingual SIR rumor propagation model is analyzed, the optimal control intensities for the prevention mechanism $u^{{\mathsf{\Gamma }}_k}\left(n\right)$ and the refutation mechanism $v^{{\mathsf{\Gamma }}_k}\left(n\right)$ are provided. Finally, through some numerical simulations, we verified the rationality of the theory and the practicality of the proposed model.