The concept of the evolutionary stable state of the water pollution particle swarm game is given. At the same time, the replication dynamics model with bounded continuous distribution time delay is deduced in detail under a special kind of game model of single water pollution particle swarm of two-strategy water pollution particle swarm.

Definition 1. PG=(X, F) is a single water pollution particle swarm game. If the social state x ¯ ∈ X satisfies the following conditions:

(1) x ¯ is the Nash equilibrium of the single water pollution particle swarm game PG.

Then, the social state x ¯ is called the evolutionary stable state (ESS) of the single-water pollution particle swarm game PG.

Note 1. It can be seen from the above definition that the evolutionary stable state is actually a refinement of the Nash equilibrium, in which the first condition is the equilibrium condition and the second condition is the stability condition. It should be noted that although the ESS is still a static concept, condition 2) already reflects the dynamic idea. If social state y is the optimal response of x ¯ , condition 2) means that when social state changes to y, social state x ¯ is better than y.

The legacy is that not all water pollution particle swarm games have evolutionary stable states.

Example 1. PG is a water pollution particle swarm game model with three strategies of single water pollution particle swarm. Therefore, the social state space is = x 1 , x 2 , x 3 ∈ R + | x 1 + x 2 + x 3 = 1 . For any social state ∈ X , the payoff function F : X → R 3 of the water pollution particle swarm game PG is defined as: F x = F 1 x F 2 x F 2 x = x 1 + 2 x 2 x 2 + 2 x 3 2 x 1 + x 3

It is easy to verify that the social state x = 1 3 , 1 3 , 1 3 is the only Nash equilibrium of PG by the definition of Nash equilibrium in the water pollution particle swarm game. In the following, we show that condition 2) in Definition 1 does not hold.

In fact, for any social state = y 1 , y 2 , y 3 ∈ X , it is easy to see that: y − x ¯ , F x ¯ = y 1 − 1 3 , y 2 − 1 3 , y 3 − 1 3 , 1,1,1 = y 1 + y 2 + y 3 − 1 = 0

If we take y = y 1 , y 2 , y 3 = 1,0,0 , then there is F(y)=(1,0,2), thus y − x ¯ , F y = y 1 − 1 3 , y 2 − 1 3 , y 3 − 1 3 , 1,0,2 = 0 ≮ 0 .

Therefore, condition 2) in Definition 1 cannot be established.

We assume that there is a water pollution particle swarm composed of a sufficiently large number of individuals, in which the individuals repeatedly interact with each other through random matching. In each strategy interaction, each individual has two strategies O 1 and O 2 for him to choose, and the corresponding interaction result is represented by the following matrix: ⋯ O 1 O 2 O 1 a b O 2 c d

Among them, c>a and b > d.

At time t, x(t) represents the individual share of choosing strategy O 1 ; F i t represents the payment corresponding to strategy O i at time t and the social state is (x(t), 1-x(t)), where i = 1,2. From a biological point of view, the above payoff function can also be understood as the fitness function of the organism, which is simply regarded as the number of its offspring. Therefore, we have: F 1 t = a x t + b 1 − x t , F 2 t = c x t + d 1 − x t .

To sum up, we induce a water pollution particle swarm game PG with two strategies of single water pollution particle swarm from the classical non-cooperative two-person symmetric game. We set α 1 = c − a , α 2 = b − d , α = α 1 + α 2 . From the definition of Nash equilibrium and evolutionary stable state in the water pollution particle swarm game, it is easy to know that there is a unique evolutionary stable state y ¯ = x ¯ , 1 − x ¯ in the above water pollution particle swarm game PG, where x ¯ = α 2 / α .For the convenience of notation, we use the individual share x ¯ of the selection strategy O 1 to represent. Because after the individual share of the selection strategy O 1 is determined, the individual share of the selection strategy O 2 is naturally determined.

Next, under the above-mentioned water pollution particle swarm game framework, considering the existence of bounded continuous distribution time delay, we deduce the replication dynamics model with bounded continuous distribution time delay in detail.

At the same time, we assume that f i t + ϵ is related to the payoff function or fitness function at each moment before time τ , that is, we consider the bounded continuous distribution delay phenomenon. Therefore, we have f i t + ϵ = 1 − ϵ f i t + ϵ f i t ∫ t − τ t e − t − s F i s d s , (2)

Among them, i = 1,2.

Note 2. When F i s ≡ F i t − τ and ∫ t − τ t e − t − s d s = 1 , formula (2) is f i t + ϵ = 1 − ϵ f i t + ϵ f i t F i t − τ ,

Among them, i = 1,2. It should be noted that from the above equation, we can obtain the fixed-delay replication kinetic model of T a o , so our model is more general. Since e − t − s is increasing with respect to s at t − τ , t , the payoff function F i s closer to time t has a greater impact on f i t + ϵ . In fact, e − t − s can also be simply regarded as a weight factor.

In the equation, it is assumed that the kernel function e − t − s is always equal to 1. From the explanation of the kernel function in the previous section, it can be known that the time delay information at each moment is equally important. Therefore, the equation has the simple form as d x t d t = − α x t 1 − x t ∫ t − τ t ( x ( s ) − x ¯ ) d s . (3)

Next, we study the stability of the evolutionary stable state x ¯ . To this end, by substituting the variable t = x t − x ¯ , Eq. 3 becomes d y t d t = − α y t + x ¯ 1 − y t − x ¯ ∫ t − τ t y s d s . (4)

Therefore, in order to discuss the stability of the constant solution x ¯ in Eq. 3, only the stability of the zero solution in Eq. 4 needs to be discussed. Linearizing Eq. 4 at y t ≡ 0 can be obtained: d z t d t = − α x ¯ 1 − x ¯ ∫ t − τ t z s d s = − δ ∫ t − τ t z s d s , (5)

Among them, = α x ¯ 1 − x ¯ .

According to the lemma, the zero solution of Eq. 5 is asymptotically stable if and only if all the eigenvalues of its corresponding eigenvalue equation have negative real parts. The corresponding feature equation can be obtained by simply changing the elements of the integral of Eq. 5. λ = − δ ∫ − τ 0 e λ s d s (6)

That is, λ 2 = − δ 1 − e − λ τ . (7)

Theorem 1. The evolutionary stable state g is asymptotically stable under replication dynamics 3) with bounded continuous distributed delays if and only if 0 < τ < π 2 2 δ .

Certification. Sufficiency. If 0 < τ < π 2 2 δ , the eigenvalues of the feature Eq. 6 all have negative real parts.

Proof by contradiction. We set λ = ρ + i ω to be a characteristic root of the feature Eq. 6, where ρ , ω ∈ R , i 2 = − 1 and ρ ≥ 0 .Obviously λ = 0 is not a solution of the feature equation. From the knowledge of linear algebra, we know that complex roots always appear in pairs, so λ = ρ − i ω is also the root of the feature Eq. 6.

Without loss of generality, we assume > 0 . When λ = ρ + i ω is substituted into the feature Eq. 9, the real part and the imaginary part are separated to obtain ρ = − δ ∫ 0 τ e − ρ s cos ⁡ ω s d s ω = δ ∫ 0 τ e − ρ s sin ⁡ ω s d s

When performing variable substitution v = ws, we get ρ = − δ ω ∫ 0 ω τ e − ρ v ω cos ⁡ v d v , ω = δ ω ∫ 0 ω τ e − ρ v ω sin ⁡ v d v (8)

Further, we set the upper limit of the integral to be = 2 k π + β , where k ≥ 0 and 0 ≤ β < 2 π , Then ω = δ ω ∫ 0 ω τ e − ρ v ω sin ⁡ v d v = δ ω ∫ 0 2 k π + β e − ρ v ω sin ⁡ v d v ≤ δ ω ∫ 0 2 k π + β sin ⁡ v d v = δ ω ∫ 0 β sin ⁡ v d v ≤ δ ω ∫ 0 π sin ⁡ v d v = 2 δ ω . (9)

Note that the first inequality holds because e − p u w decreases over the interval 0,2 k π + β with respect to v. Since the sine function is non-negative on the interval 0 , π and non-positive on the interval π , 2 π , the second inequality holds. Furthermore, there is ≤ 2 δ .

Since 0 < τ < π 2 2 δ , 2 k π + β = ω τ < 2 δ ⋅ π 2 2 δ = π (10)

The above inequality states that k = 0 and < π .

i) If 0 < ω τ ≤ π 2 , then ∫ 0 ω τ e − p u ω cos ⁡ v d v > 0 . Note that ω > 0 ; > 0 , and then i) ρ = − δ ω ∫ 0 ω τ e − ρ v ω cos ⁡ v d v < 0 .

This contradicts hypothesis ≥ 0 .

(ii)If π 2 < ω τ < π , then ∫ 0 ω τ e − ρ v ω cos ⁡ v d v = ∫ 0 π 2 e − ρ v ω cos ⁡ v d v + ∫ π 2 ω τ e − ρ v ω cos ⁡ v d v ≥ e − ρ π 2 ω ∫ 0 π 2 cos ⁡ v d v + e − ρ π 2 ω ∫ π 2 ω τ cos ⁡ v d v > 0

Note that ∫ π 2 ω τ e − ρ v ω cos ⁡ v d v ≥ e − ρ π 2 ω ∫ π 2 ω τ cos ⁡ v d v in the first inequality is because π 2 , ω τ is non-positive on the interval cos ⁡ v . Thus ρ = − δ ω ∫ 0 ω τ e − ρ v ω cos ⁡ v d v < 0 .

This contradicts hypothesis ≥ 0 .

To sum up, we always have ρ ≥ 0 contradiction, so the characteristic roots of feature Eq. 9 all have negative real parts.

Necessity. It is only necessary to prove that when ≥ π 2 2 δ , the zero solution of the feature Eq. 10 is unstable, that is, there are characteristic roots with positive real parts.

Note that when τ = π 2 2 δ and = π , it is easy to verify that λ = i ω is the root of the feature Eq. 10. The following shows that when τ continues to increase, the characteristic root of the feature Eq. 10 will enter the right half plane of the complex plane and it is impossible to cross the imaginary axis again and return to the left half plane. 2 λ d λ d τ + δ τ e − λ τ d λ d τ + δ λ e − λ τ = 0 . (11)

Note that since the feature Eq. 10 has e − λ τ = δ + λ 2 , the above Eq. 11 can be simplified as: d λ d τ = − λ 3 + δ λ τ λ 2 + τ δ + 2 λ .

Then, d λ d τ λ = i ω = i ω 3 − i δ ω τ δ − ω 2 τ + 2 i ω = i ω 3 − i δ ω τ δ − ω 2 τ − 2 i ω τ δ − ω 2 τ 2 + 4 ω 2 = 2 ω 2 ω 2 − δ − i τ ω 2 − δ 2 τ δ − ω 2 τ 2 + 4 ω 2

Furthermore, it can be known from δ = 1 2 ω 2 that:   R   e   d λ d τ λ = i ω = 2 ω 2 ω 2 − δ ω 2 τ − δ τ 2 + 4 ω 2 > 0

Similarly, we can get: R e d λ d τ λ = − i ω = 2 ω 2 ω 2 − δ ω 2 τ − δ τ 2 + 4 w 2 > 0

Therefore, when τ increases slightly from π 2 2 δ , the root of the feature Eq. 10 on the imaginary axis will enter the right half-plane.

Furthermore, when τ continues to increase, if the characteristic root returns to the imaginary axis at τ = τ ∼ , that is, there is ω ∼ > 0 such that λ ∼ = ± i ω ∼ is a pure imaginary root of the feature Eq. 10. Exactly like the above discussion, we have: R e d λ d τ λ ∼ = ± i ω ∼ = 2 ω ∼ 2 ω ∼ 2 − δ ω ∼ 2 τ − δ τ 2 + 4 ω ∼ 2 > 0 .

This shows that with the increase of, and if the root of the feature Eq. 10 returns to the imaginary axis, it is impossible to pass through the imaginary axis and enter the left half plane.

To sum up, as long as ≥ π 2 2 δ , the feature Eq. 10 always has characteristic roots with positive real parts.

In the hawk-dove game, two agents compete for resources with value V, and each agent has two strategies for him to choose, namely Hawk and Dove. If all individuals adopt the eagle strategy, the fight will end with one side injured, and the injured side will be at the cost of a decrease in fitness C.At the same time, individuals with two eagle strategies often have a 50% probability of defeating their opponents and a 50% probability of being injured. If one individual chooses the eagle strategy and the other chooses the dove strategy, the individual with the dove strategy withdraws and ends the competition. At this time, the eagle strategy individual enjoys the resource exclusively. In addition, if both agents choose the pigeon strategy, the individuals share the resource together. The eagle-dove game can be represented by the following matrix.   H D H V − C 2 V D 0 V 2

Among them, C > 0, V > 0. Further, we assume V < C, which means that the cost of injury is higher than the value of resources.

Obviously, the hawk-dove game model is a special case of the model discussed in our first section. At this point, the evolutionary stable state is x ¯ = V C .From Theorem 1, the evolutionary stable state x ¯ is asymptotically stable if and only if τ is less than π 2 2 δ . Specifically, under the Eagle-Dove Boyi model, the upper bound of τ is C π 2 V C − V .

Taking C = 7 and V = 3, we numerically simulate the replication dynamics 6) with a bounded continuous distribution delay under a constant kernel function.

It can be seen from Figure 1 that the numerical solution of Eq. 7 converges to the evolutionary stable state x ¯ when τ is smaller than the critical value. At the same time, when the value of τ is close to the critical value, the numerical solution of Eq. 7 converges to the evolutionary stable state x ¯ , and the amplitude of oscillation increases with the closeness of g to the critical value. When τ keeps increasing and exceeds the critical value, as shown in Figure 2, the solution of Eq. 6 keeps oscillating and does not converge to x ¯ , so the evolutionary stable state x ¯ is not asymptotically stable at this time.

The stability of replication dynamics with bounded continuous time delays under constant kernel function is discussed. This section examines the stability of the evolutionary steady state when the kernel function has an exponential form. Note that the conclusions in this section are only a sufficient condition for the evolutionary stable state to be asymptotically stable. In fact, for general delay differential equations, it is difficult to obtain sufficient and necessary conditions for the stability of the solutions as in the previous section. This is the so-called conservative problem in the field of differential equations.

Below, we consider a replication dynamics model with bounded continuous distributed delays of the form. d x t d t = − α x t 1 − x t ∫ t − τ t e − t − s x s − x ¯ d s (12)

Among them, = α 1 + α 2 , x ¯ = α 2 / α , α 1 = c − a , α 2 = b − d .

Note that in the above-mentioned time-delay evolution dynamics, the time-delay information at each moment is not equally important. The time lag information that is closer to the current time t has a greater impact on the current state. It can be said that the exponential kernel function is more general and more realistic than the constant kernel function. Similar to the discussion in the previous section, in order to examine the stability of the evolutionary stabilization strategy x ¯ , this paper first substitutes the variable y t = x t − x ¯ to get: d y t d t = − α y t + x ¯ 1 − y t − x ¯ ∫ t − τ t e − t − s y s d s (13)

Then the stability study of the solution x t ≡ x ¯ of Eq. 12 is transformed into the stability study of the zero solution of Eq. 13.

Further, linearizing Eq. 13 around t ≡ 0 , we get: d z t d t = − α x ¯ 1 − x ¯ ∫ t − τ t e − t − s z s d s = − δ ∫ t − τ t e − t − s z s d s (14)

Among them, δ = α x ¯ 1 − x ¯ .

Similar to the discussion in the previous section, it is easy to obtain the feature equation of the linearized Eq. 14 as follows: λ + δ ∫ − τ 0 e λ + 1 s d s = 0 , (15)

That is, λ 2 + λ + δ − δ e − τ λ + 1 = 0 . (16)

When the time delay τ is sufficiently small, the first-order Taylor expansion of the exponential term e − τ λ + 1 in the feature Eq. 16 is carried out, and the above equation is converted into the following quadratic equation about. λ 2 + 1 + δ τ λ + δ τ = 0 (17)

Theorem 2. For sufficiently small delays τ, the evolutionary stable state 𝑥¯ is asymptotically stable under replica dynamics (12) with bounded continuously distributed delays.

The two roots λ 1 and λ 2 of Eq. 17 satisfy: λ 1 + λ 2 = − 1 + δ τ < 0 , λ 1 ⋅ λ 2 = δ τ > 0

Therefore, the evolutionary stable state x ¯ is asymptotically stable under the replication dynamics (12) when the time lag is sufficiently small.

Next, for a general time-delay τ, we obtain a sufficient condition for the evolutionary stable state 𝑥¯ to be asymptotically stable under replication dynamics (12). We look forward to finding a larger upper bound for delay τ in future work, which involves the conservation of equations and is also one of the key issues in the study of delay differential equations.

Theorem 3. Under replication dynamics (12) with bounded continuously distributed delays, the evolutionary stable state x ¯ is asymptotically stable if < π 2 δ .

Certification: λ = u + i v is the root of the feature Eq. 15, where , v ∈ R . Then, according to the feature Eq. 15, we have λ + δ ∫ − τ 0 e λ + 1 s d s = u + i v + δ ∫ − τ 0 e u + i v + 1 s d s = u + i v + δ ∫ − τ 0 e u + 1 s ⋅ e i v s d s = u + i v + δ ∫ − τ 0 e u + 1 s cos ⁡ v s + i sin ⁡ v s d s = 0

Separating the real and imaginary parts in the above formula, we get u = − δ ∫ − τ 0 e u + 1 σ cos ⁡ v σ d σ v = − δ ∫ − τ 0 e u + 1 σ sin ⁡ v σ d σ

If the real part of the characteristic root λ = u + i v is assumed to be ≥ 0 , then v = δ ∫ − τ 0 e u + 1 σ sin ⁡ v σ d σ ≤ δ ∫ − τ 0 e u + 1 σ d σ ≤ δ < π 2 τ (18)

Further, from the integral median theorem, we know that: 0 ≥ − 1 δ ⋅ u = ∫ − τ 0 e u + 1 σ cos ⁡ v σ d σ = τ e u + 1 σ * cos ⁡ v σ *

Among them, σ * ∈ − τ , 0 . Thereby 2 k π + π 2 ≤ v σ * ≤ 2 k π + 3 π 2 ,

Among them, ∈ Z .

Note that by combining Eq. 18, we get: v σ * ≤ v τ < π 2 .

This is a contradiction. Therefore, the real part is u < 0. This shows that when τ < π 2 δ , the eigenvalues of the characteristic Eq. 15 all have negative real parts, so that the evolutionary stable state x ¯ is asymptotically stable under the replication dynamics (12) with bounded continuous distributed delays.

At the end of the previous section, we validated our results by numerical simulations of replication dynamics with a bounded continuous distributed delay under a constant kernel function against the background of the hawk-dove game model. Next, we give a simple example of replication dynamics (12) with bounded continuous time-delay under the exponential kernel function under the eagle-dove game model with the same parameters. We briefly restate the eagle-dove game model, and its game matrix is as follows.   H D H V − C 2 V D 0 V 2

Among them, C > 0, V > 0, V < C. Continuing from the previous section, when C = 7, V = 3, it is easy to calculate the evolutionary stable state x ¯ = 0.428 . It can be seen from Figure 3A that when the time delay τ is small, the solution of Eq. 12 converges to the evolutionary stable policy x ¯ . It is worth noting that the convergence rate is relatively slow at this time. It can be seen from Theorem two that the evolutionary stable state x ¯ is asymptotically stable when < π 2 δ , and Figure 3B verifies our conclusion. However, since our conclusion is only a sufficient condition for asymptotic stability. Therefore, in the future work, how to theoretically enlarge the upper bound of the delay, that is, to reduce the conservatism of the integro-differential Eq. 12 will be a very meaningful work in terms of both the application of game theory and the theory of time-delay differential equations.

In this paper, a method combining leaching surface and flux concentration is proposed, and a numerical simulation method is applied to simulate and analyze transient leakage monitoring of similar point and line source pollution in a planar two-dimensional heterogeneous aquifer. On the stratigraphic characteristic scale, the stochastic simulation method is used to generate the profile two-dimensional heterogeneous anisotropic aquifer. Under the conditions of convective dispersion and adsorption, the migration and evolution of polluted plumes in this special heterogeneous formation are analyzed by numerical simulation method, which provides a reference for the monitoring and management of such polluted aquifers. The prediction of groundwater pollution diffusion path based on the total multi-source data fusion in this paper is shown in Figure 4.

Taking the identification of sudden pollution incidents as the ultimate goal, the four factors of environmental impact factors, emergency response factors, social impact factors, and groundwater system factors are used as the criterion layers. The identification index system of sudden pollution incidents is established by refining the supported indicators respectively, as shown in Figure 5. The selection of indicators in the indicator system is based on the principles of scientificity, systematicness, pertinence and operability.

In order to further verify the reliability of the established mathematical model and model parameters, the experimental analysis is carried out, and the study is carried out in a certain area, as shown in Figure 6. From the concentration fitting diagram in Figure 7, a total of 23 observation points is taken in the model identification stage. It can be seen that the isolines of the measured and calculated concentrations also achieve a good fit as a whole. It shows that the established hydrogeological conceptual model and mathematical model are reasonable.

It is divided into six stress periods, and each 360 days is a stress period. From the inspection results, select three representative observation holes. It can be seen that the model identification and verification results prove that the established mathematical model, boundary conditions, hydrogeological parameters and source-sink terms are all reasonable, and the model can be used for the prediction of groundwater solute system (Figure 8).

The above research verifies that the multivariate data fusion proposed in this paper can play an important role in the prediction of groundwater pollution diffusion paths.