1 Introduction

1.1 State of the art and motivation

The games with a continuum of agents have been widely studied during the last decade thanks to the mean field games (MFG) theory, introduced by Lasry–Lions [34,35,36] and Huang–Caines–Malham [26,27,28]. An application of this framework is the modeling of the vaccination decision in a population where an individual has the choice between only two pure strategies: to get vaccinated or not.

This decision is a part of the behavioral epidemiology, which takes into account the interplay between human behavior and spread of infectious diseases (see [38] for a detailed presentation and [30, 40, 47] for other examples). In particular, risk perception and spread of information (on vaccination risk or probability of infection) can influence at a large scale the population behavior and thus the spread of the disease. As detailed in [38], a well-known example is Measles–Mumps–Rubella (MMR) vaccination coverage decrease between 1998 and 2003. This caused a decline in herd immunity and a measles resurgence after the publication of an hypothesis linking MMR vaccination and autism. Behavioral responses to the threat of a pandemic event can influence the spread of the disease in several ways: people can choose preventive protection (by vaccination as presented below) or limit the disease spread by isolation and quarantine (e.g., school closures). Since then, many epidemiological models have been developed to include individual response (see [47] for a review); on the contrary see [1, 6, 39, 42] for the treatment of compulsory (as opposed to voluntary) vaccination.

We use the family of compartmental models in order to describe the spread of the disease in the population. This supposes that all individuals are the same and the population is divided into several exclusive classes; the transmission between two classes is modeled through a time-dependent rate. As each compartment represents the status of an individual, it means that the structure of model with immunity is not the same as that of a disease without immunity. \({ SIR}\) terminology usually represents a disease with immunity against re-infection, indicating that the passage of an individual is from the susceptible class S to the infectious class I and then to the recovered class R. If the re-infection is possible, then the terminology used can be \({ SIS}\), which means that an individual in the infected class I goes back to the susceptible class S. If temporary immunity is given by the disease, the compartmental structure can be \({ SIRS}\), meaning that after infection, individual stays in R class before going back to the susceptible class. Other well-known compartmental models can include an exposed class (E) that represents the period between being infected and becoming contagious, like in the \({ SEIR}\) or \({ SEIS}\) models. The evolution of the size of the compartments is given by ordinary differential equations (for more details see for instance [10, 15, 25]).

In order to model individual vaccination decision in the spread of the disease, vaccination coverage can depend on available information (see [11, 18, 19]) or be determined by a vaccination game. The evolution of the vaccination coverage could be determined by an imitation dynamics depending (or not) of the available information (see [3, 16, 16, 17, 32, 41]), or be the solution of a Nash equilibrium determined at the beginning of the game, as described below. The Nash–MFG equilibrium corresponds in this case to a strategy where nobody has interest to change his own vaccination strategy. Therefore the existence and the determination of an equilibrium is central.

Before the development of the Nash–MFG theory, some work have been proposed to couple game and epidemiology. Three well-known methods that introduced specific models are presented in the next paragraph with their advantages and limitations. The compartmental model used in [5] includes five classes (susceptible, exposed, infected, recovered, vaccinated), the dynamics is presented in Eq. (1) where the population is divided into the following classes: S (susceptible), I (infected), R (recovered) and V (vaccinated). In this case preventive vaccination for the smallpox is available: the model is used to determine the probability of an individual to choose preventive vaccination. To do so, explicit individual payoffs for the two pure strategies (vaccination before or after the outbreak) are defined. Then the existence and uniqueness of the Nash equilibrium is obtained using properties of such payoffs.

$$\begin{aligned} \begin{aligned} \frac{{ dS}}{{ dt}}&= - \beta S I - f(S,t), \\ \frac{dE}{dt}&= \beta S I - \rho E, \\ \frac{{ dI}}{{ dt}}&= \rho E - \gamma I, \\ \frac{{ dR}}{{ dt}}&= \gamma I, \\ \frac{{ dV}}{{ dt}}&= f(t,S), \\ \end{aligned} \end{aligned}$$
(1)

Another method introduced in [5] uses the same methodology as previously but on the stationary state of the system (2) (or equivalently system (2bis) that distinguishes immunity given by the disease and immunity given by the vaccine) because of the possible non extinction of the disease. This allows to determine the probability of vaccination when disease is already present in the population with a compartmental system where vital dynamics is considered. As in previous case, the probability p to vaccinate is not time dependent and properties of payoff are used to prove existence and uniqueness of the equilibrium. Finally, solutions for other incidence rate (term \(\beta I\) is replaced by a function \(\lambda (I)\)) are given in [43] by using Markov decision process. Analytical solutions for the vaccination probability at the Nash equilibrium can be found with these approaches but it supposes constant probability of vaccination.

$$\begin{aligned} \begin{aligned} \frac{{ dS}}{{ dt}}&= \mu (1-p) - \beta S I - \mu S \\ \frac{{ dI}}{{ dt}}&= \beta S I - \gamma I - \mu I \\ \frac{{ dR}}{{ dt}}&= \mu p + \gamma I- \mu R \end{aligned} \end{aligned}$$
(2)
figure a

The third and last model being detailed here describes numerical approach when two \({ SEIR}\) models (see for example (1)) are used in order to separate two age classes (see [24]) to obtain the vaccination probability of each age class. The evolution of the probability of vaccination can be modeled by introducing a rule of thumb (see [3, 4]) to study the epidemiological stationary state and the role of parameters like vaccination and infection costs. But reaching the stationary state of the epidemiological system can take a long time as explained in [29]. To avoid this problem, [20] considers a \({ SIRS}\) model where the transmission rate \(\beta \) is a periodic function, causing the same property for the incidence. Economic approach (see [21]) based on the epidemiological model (3) (where notations are the same as above and \(\mathcal {V}(t)\) is the instantaneous vaccination rate) allows to determine the vaccination and no vaccination region (depending on the state of the disease) by equalization of the two payoffs but no proofs of optimality (existence or uniqueness) are provided.

$$\begin{aligned} \begin{aligned} \frac{{ dS}(t)}{{ dt}}&= - \beta S(t) I(t) - \mathcal {V}(t) \\ \frac{{ dI}(t)}{{ dt}}&= \beta S(t) I(t) - \gamma I(t) \\ \frac{{ dR}(t)}{{ dt}}&= \gamma I(t) dt \\ \frac{{ dV}(t)}{{ dt}}&= \mathcal {V}(t) \end{aligned} \end{aligned}$$
(3)

Within the MFG framework, article [33] proposes analytic solutions in some particular cases. Furthermore, finding an equilibrium becomes more complicated if the model has specific behaviour with network interaction or epidemic spread for example [14, 23]. Thus finding the equilibrium can quickly become time consuming. One commonly used method is based on finding a fixed point of some function, which can be formulated as the attractor of an evolution system.

The existence of an equilibrium in the general case of MFG is of utmost interest (see for instance [7, 9] for an entry point to this literature). Several approaches were developed, inspired by the general framework of gradient flows (see [31]) where the procedure is iterative.

1.2 Mathematical framework

Let us recall that in a gradient flow framework (which is not necessarily related to a MFG model !) the equation to be solved is written symbolically (see [2] for an introduction):

$$\begin{aligned} \frac{d}{d \tau } \xi + \nabla F(\xi ) = 0 \end{aligned}$$
(4)

where the function F is given but \(\nabla F\) cannot be computed. The numerical counterpart of (4) is the famous JKO numerical scheme (introduced by Jordan et al. in [31]) which, for \(\xi _k\) an approximation of \(\xi _{\tau _{k}}\), can be written as:

$$\begin{aligned} \xi _{k+1} \in \underset{\eta }{\text {argmin }} \left\{ \frac{d(\xi _k, \eta )^2}{2 {\varDelta }} \tau + F(\eta ) \right\} . \end{aligned}$$
(5)

This scheme can be seen as the variational implicit Euler scheme and has been used for the heat equation as a gradient flow of the entropy in the Wasserstein space. Second order in time methods have been proposed very recently (see [37]).

On the other hand, a MFG equilibrium has no associated function ‘F’ and it is likely that no such function exists. For each \(\tau , \xi _\tau \) is a probability law. Rather, the search for an equilibrium takes the form of a fixed point:

$$\begin{aligned} \xi = J(\xi ). \end{aligned}$$
(6)

In practice, inspired by the fictitious game (see [12]) and Best Reply (see [8, 9]) procedures, a new algorithm has been introduced in [44, Eq 3.2] and [46] with the following form:

$$\begin{aligned} \xi _{k+1} \in \underset{ \eta }{\text {argmin }} \left\{ \frac{ d(\eta , \xi _k)^2}{2 {\varDelta }\tau } + \mathcal {P}(\eta , \xi _k) \right\} . \end{aligned}$$
(7)

In the limit \(\tau \rightarrow 0\) the Eq. (7) describes a curve which under some assumptions on \(\mathcal {P}\) will be called a solution of the evolution equation:

$$\begin{aligned} \frac{d}{d \tau } \xi _{\tau } + \nabla _1 \mathcal {P}(\xi _{\tau }, \xi _{\tau }) = 0. \end{aligned}$$
(8)

On the other hand, the MFG framework provides a natural function \(\mathcal {P}(\cdot , \cdot )\) even when \(\xi \) does not belong to a Hilbert space but only to a metric space. Thus we are led to consider evolution equations in metric spaces with semi-explicit numerical schemes (see [9, 46]).

1.3 Scope and structure of the paper

The purpose of this work is to introduce variational Runge–Kutta explicit methods of high order in a metric space using a generalization of linearity (presented in [13]).

These approaches are applied to the control of epidemic spread with voluntary vaccination as in [44].

The paper is structured as follows: Sect. 2 is dedicated to the Runge–Kutta method recalling its general form (Sect. 2.1) and then the introduction of variational Runge–Kutta methods (Sect. 2.2). The equivalence of the two approaches in a Hilbert space is proved in Sect. 2.3. Then Sect. 3 gives the epidemiological model used in the numerical application given in Sect. 4. An economical application is proposed in Sect. 5. Some considerations are discussed in Sect. 6.

2 Runge–Kutta methods

2.1 In a Hilbert space

Let \(\mathcal {H}\) be a Hilbert space, \(f: [0, T] \times \mathcal {H} \rightarrow \mathcal {H}\) a regular function and \(y: [0,T] \rightarrow \mathcal {H}\) satisfying (in some sense to be specified):

$$\begin{aligned} \frac{d}{d \tau }y(\tau ) = f(\tau , y(\tau )) \end{aligned}$$
(9)

The Runge–Kutta approach is an iterative method using temporal discretization in order to obtain a numerical solution of an ordinary differential equation of type (9). The time horizon T is supposed to be finite and can be discretized in \((N_T+1)\) time instants \( \tau _0 = 0, \tau _1 = {\varDelta }\tau , \tau _2 = 2 {\varDelta }\tau , \ldots , \tau _{N_T} = T\) where \({\varDelta }\tau \) is the time step. Let \(y_{k} \in \mathcal {H}\) be an approximation of \(y(\tau _k)\).

Runge–Kutta method is based on the evaluation at intermediate points in time. Equation (10) presents the general form of the Runge–Kutta method for (9).

Fig. 1
figure 1

Runge–Kutta method. Left Butcher tableau. Right equation

A method is defined by the values of the coefficients \(b_i, c_i\) and \(a_{i,j}\), often presented in the Butcher tableau (see table in Fig. 1). Consistency of such a method is ensured if and only if \(\sum _{i=1}^s b_i = 1\).

We will present below some particular schemes that will later be formulated in a metric space.

2.1.1 Explicit Euler scheme

The explicit Euler (denoting EE) scheme is a Runge–Kutta method with only one step. It is defined in Fig. 2 by the Butcher tableau or equivalently by Eq. (11). In this case, this scheme is of order one in time.

Fig. 2
figure 2

Explicit Euler method. Left Butcher tableau. Right equation

2.1.2 Heun scheme

Heun scheme uses an approximate value \(p_1\), in order to compute \(y_{k+1}\), see the Fig. 3. This scheme is of order two in time.

Fig. 3
figure 3

Heun method. Left Butcher tableau. Right equation

2.1.3 Explicit Runge–Kutta 3 scheme

In the same way, we introduce the Runge–Kutta 3 (denoted RK3) method. It calculates two intermediate values (\(p_2, p_3\)) to obtain \(y_{k+1}\). See Fig. 4 for the value of coefficients in the Butcher table and Eq. (13) for the explicit equations. In this case, the method is of order three in time.

Fig. 4
figure 4

Runge–Kutta 3 method. Left Butcher tableau. Right equation

2.1.4 Explicit Runge–Kutta 4 scheme

The last method is the celebrated Runge–Kutta 4 (denoted RK4) where three intermediate values (\(p_2, p_3, p_4\)) are computed to determine \(y_{k+1}\). The value of coefficients is presented in Fig. 5 and its application to Eq. (9) is given in Eq. (14). This method is of order four in time.

Fig. 5
figure 5

Runge–Kutta 4 method. Left Butcher tableau. Right equation

2.2 Variational approach

Unfortunately, none of the schemes presented in Sect. 2.1 can be used for finding a MFG equilibrium. Firstly because the problem is not really presented that way, and secondly because the space of the unknowns is, in general, not a vector space but a metric space. Therefore we need first to generalize the Runge–Kutta schemes to metric spaces and then to apply them to the MFG setting.

Let \((\mathcal {A}, d)\) be a metric space. The Definition 1 below gives the properties of a geodesic space (see [37] for additional details).

Definition 1

A curve \(\gamma : [0,1] \rightarrow \mathcal {A}\) is called a (constant speed) geodesic provided that \(d(\gamma (s), \gamma (\tau )) = |s- \tau | \cdot d(\gamma (0), \gamma (1))\) and the space is called geodesic if, for any couple of points \((X,Y) \in \mathcal {A}^2\), there exists at least a geodesic \(\gamma \) connecting them, that is, \(\gamma (0) = X, \gamma (1) = Y\).

From now on we assume that \((\mathcal {A}, d)\) is a geodesic space. Consider \(\mathcal {P}: \mathcal {A} \times \mathcal {A} \rightarrow \mathbb {R}\) a function with the following properties:

(H1) :

\(\forall n \le 4, \forall r>0, \forall \xi , Y_1, \dots , Y_n \in \mathcal A\), the set of vectors \(\{ (\mathcal {P}(z, Y_1), \dots ,\mathcal {P}(z, Y_n)); z \in {\mathcal A}, d(z,\xi ) \le r \} \) is compact as a subset of \(\mathbb {R}^n\),

(H2) :

for any point \(Y \in \mathcal {A}\) and any \(X_0, X_1 \in \mathcal {A}\) a constant speed geodesic \(\gamma \) exits such that \(\gamma (0) = X_0, \gamma (1) = X_1\) and the function \(\tau \mapsto \mathcal {P}(\gamma (\tau ), Y)\) from [0, 1] to \(\mathbb {R}\) is linear.

Examples:

  • If \(\mathcal {H}\) is a Hilbert space, \(f: \mathcal {H} \rightarrow \mathcal {H}\) a function then \((x,y) \rightarrow \left\langle x,f(y)\right\rangle \) satisfies hypothesis (H1) and (H2), see proof of Proposition 2.2.

  • Let \(\mathcal {A}\) be the set of probability measures on [0, T]; \(\mathcal {F}\) be the set of Lipschitz functions with constant less than 1 and the metric of weak convergence d is defined, for \(a_1, a_2 \in \mathcal {A}\) by (see [45]): \(d(a_1, a_2) = \sup _{f \in \mathcal {F}} | \int _0^T f(t) a_1(dt) - \int _0^T f(t) a_2(dt) |\). For \(Y \in \mathcal {A}, t \in [0,T]\) we define \(g_Y(t) = \left[ \sin (t) \int _0^T \cos (s) Y(ds) \right] ^2\) and for \(z, Y \in \mathcal {A}, \mathcal {Q}(z,Y) = \int _0^T g_Y(t) z(dt)\).

    Let \((z_k)_{k \in \mathbb {N}}\) be a sequence of elements in \(\mathcal {A}\) converging weakly to \(z \in \mathcal {A}\) then, for \(f \in \mathcal {F}, \int _0^T f(t) z_k(dt)\) converges to \(\int _0^T f(t) z(dt)\). As \((\mathcal {A},d)\) is compact (see [45]) and the application \(z \in \mathcal {A} \rightarrow (\mathcal {Q}(z,Y_1), \dots ,\mathcal {Q}(z,Y_n)) \in \mathbb {R}^n\) continuous, \(\mathcal {Q}\) satisfies (H1).

    Consider \(a_1\) and \(a_2 \in \mathcal {A}\) and for \(s \in [0,1], \gamma _s = sa_1 + (1-s) a_2\). Then \(d(\gamma _{s}, \gamma _{\tau }) = \sup _{f \in \mathcal {F}} | \int _0^T f(u) \gamma _{s}(du) - \int _0^T f(u) \gamma _{\tau }(du) | = \sup _{f \in \mathcal {F}} | \int _0^T f(u) (\tau -s)a_1(du) - \int _0^T f(u) (\tau -s)a_2(du)| = |\tau -s| \sup _{f \in \mathcal {F}} | \int _0^T f(u) a_1(du) - \int _0^T f(u) a_2(du)| = |\tau -s| \cdot d(\gamma (0),\gamma (1))\). It shows that \(\gamma \) is a (constant speed) geodesic. Then, the function \(s \rightarrow \mathcal {Q}(\gamma _s, Y)\) of the form \(s \rightarrow s \int _0^T g_Y(u) a_1(du) + (1-s) \int _0^T g_Y(u) a_1(du) \) is linear.

  • Note that previous example also works with the total variation distance.

The evolution equation that we want to solve has the form of Eq. (8) (but note that this is only a formal expression as \(\frac{d}{d \tau } \xi _{\tau }\) does not have a well-defined meaning in a general metric space, see [46]).

This section develops a variational Runge–Kutta method in a metric space. We propose a method that does not assume any vector calculus but uses the hypothesis (H2) which defines the notion of linear application in a metric space (see also [13, Section 8, Definition 8.1, p. 480] for more details). As \({\varDelta }\tau \) tends to zero Eq. (7) describes a curve. The purpose of the scheme is to describe this evolution faster and more precisely that is to say to obtain high order schemes.

2.2.1 Variational explicit Euler scheme

The first scheme of that form was introduced in [44, Eq 3.2]. Formally, we define the variational explicit Euler scheme (VEE) by:

$$\begin{aligned} \xi _{k+1} \in \underset{\eta \in \mathcal {A}}{\text {argmin }} \left\{ \frac{ d(\eta , \xi _k)^2}{2 {\varDelta }\tau } + \mathcal {P}(\eta , \xi _k) \right\} . \end{aligned}$$
(15)

2.2.2 Variational Heun scheme

The variational Heun scheme (VH) introduces, as in vector space, an intermediate value \(\widetilde{\xi }_{k+1}\) (Eq. 16a) and uses it to compute \(\xi _{k+1}\) (Eq. 16b).

$$\begin{aligned} \widetilde{\xi }_{k+1}&\in \underset{\eta \in \mathcal {A}}{\text {argmin }} \left\{ \frac{d(\eta , \xi _k)^2}{2 {\varDelta }\tau } + \mathcal {P}(\eta , \xi _k)\right\} , \end{aligned}$$
(16a)
$$\begin{aligned} \xi _{k+1}&\in \underset{\eta \in \mathcal {A}}{\text {argmin }} \left\{ \frac{d(\eta , \xi _k)^2}{2 {\varDelta }\tau } + \frac{1}{2}\mathcal {P}(\eta , \xi _k)+\frac{1}{2}\mathcal {P}(\eta , \widetilde{\xi }_{k+1}) \right\} . \end{aligned}$$
(16b)

Note that two minimizations are required in order to obtain \(\xi _{k+1}\).

2.2.3 Variational Runge–Kutta 3 scheme

The variational Runge–Kutta 3 scheme (VRK3) is defined by Eq. (17). Two intermediate values are computed in this case [\(\widetilde{\xi }^1_{k+1}\) with Eq. (17a) and \(\widetilde{\xi }^2_{k+1}\) with Eq. (17b)] in order to obtain \(\xi _{k+1}\) using Eq. (17c).

$$\begin{aligned} \widetilde{\xi }^1_{k+1}&\in \underset{\eta \in \mathcal {A}}{\text {argmin }} \left\{ \frac{d(\eta , \xi _k)^2}{2 {\varDelta }\tau } + \frac{1}{2}\mathcal {P}(\eta , \xi _k) \right\} , \end{aligned}$$
(17a)
$$\begin{aligned} \widetilde{\xi }^2_{k+1}&\in \underset{\eta \in \mathcal {A}}{\text {argmin }} \left\{ \frac{d(\eta , \xi _k)^2}{2 {\varDelta }\tau } - \mathcal {P}(\eta , \widetilde{\xi }^1_{k+1})+2\mathcal {P}(\eta , \widetilde{\xi }^2_{k+1}) \right\} , \end{aligned}$$
(17b)
$$\begin{aligned} \xi _{k+1}&\in \underset{\eta \in \mathcal {A}}{\text {argmin }} \left\{ \frac{d(\eta , \xi _k)^2}{2 {\varDelta }\tau } + \frac{1}{6}\mathcal {P}(\eta , \xi _k) + \frac{2}{3} \mathcal {P}(\eta , \widetilde{\xi }^1_{k+1}) \right. \nonumber \nonumber \\&\quad \left. + \,\frac{1}{3} \mathcal {P}(\eta , \widetilde{\xi }^2_{k+1}) \right\} . \end{aligned}$$
(17c)

2.2.4 Variational Runge–Kutta 4 scheme

The variational Runge–Kutta 4 scheme is presented in Eq. (18). Three intermediate values [\(\widetilde{\xi }^1_{k+1}\) with Eq. (18a), \(\widetilde{\xi }^2_{k+1}\) with Eq. (18b) and \(\widetilde{\xi }^3_{k+1}\) defined by Eq. (18c)] are computed in order to obtain \(\xi _{k+1}\) by using Eq. (18d).

$$\begin{aligned} \widetilde{\xi }^1_{k+1}&\in \underset{\eta \in \mathcal {A}}{\text {argmin }} \left\{ \frac{d(\eta , \xi _k)^2}{2 {\varDelta }\tau } + \frac{1}{2}\mathcal {P}(\eta , \xi _k) \right\} , \end{aligned}$$
(18a)
$$\begin{aligned} \widetilde{\xi }^2_{k+1}&\in \underset{\eta \in \mathcal {A}}{\text {argmin }} \left\{ \frac{d(\eta , \xi _k)^2}{2 {\varDelta }\tau } + \frac{1}{2}\mathcal {P}(\eta , \widetilde{\xi }^1_{k+1}) \right\} , \end{aligned}$$
(18b)
$$\begin{aligned} \widetilde{\xi }^3_{k+1}&\in \underset{\eta \in \mathcal {A}}{\text {argmin }} \left\{ \frac{d(\eta , \xi _k)^2}{2 {\varDelta }\tau } + \frac{1}{2}\mathcal {P}(\eta , \widetilde{\xi }^2_{k+1}) \right\} , \end{aligned}$$
(18c)
$$\begin{aligned} \xi _{k+1}&\in \underset{\eta \in \mathcal {A}}{\text {argmin }} \left\{ \frac{d(\eta , \xi _k)^2}{2 {\varDelta }\tau } + \frac{1}{6}\mathcal {P}(\eta , \xi _k)+\frac{1}{3}\mathcal {P}(\eta , \widetilde{\xi }^1_{k+1}) \right. \nonumber \\&\quad \left. +\,\frac{1}{3}\mathcal {P}(\eta , \widetilde{\xi }^2_{k+1})+\frac{1}{6}\mathcal {P}(\eta , \widetilde{\xi }^3_{k+1}) \right\} . \end{aligned}$$
(18d)

Note that four minimizations are needed in order to obtain \(\xi _{k+1}\).

See Sect. 4.2 for a numerical application and the study of the order in time of these schemes.

2.3 Property of variational scheme

In this subsection we prove that using hypotheses (H1) and (H2), the variational schemes introduced above are well-defined (Proposition 2.1) and correspond, when the ambient space is Hilbert, to the Runge–Kutta method presented in Sect. 2.1 (Proposition 2.2).

Proposition 2.1

Under the hypothesis (H1) and (H2), the Eqs. (15), (), () and () defining respectively the VEE, VH, VRK3 and VRK4 schemes admit a solution, i.e., the schemes are well defined.

Proof

As the argumentation is very similar for all schemes, we present the general case where \(\xi _{k+1} \in \text {argmin}_{\eta \in \mathcal {A}} \{ \mathscr {F}(\eta ) \}\) with \(\mathscr {F}: \mathcal {A}\rightarrow \mathbb {R} \), for fixed \(\xi \in \mathcal {A}, \mathscr {F}(\eta ) = \frac{d(\eta ,\xi )^2}{2 {\varDelta }\tau } + l(\eta )\) where l satisfies hypothesis (H2) and the following hypothesis (H1bis):

(H1bis) :

\(\forall r>0, \forall \xi \in \mathcal A\), the set \(\{ l(z); z \in {\mathcal A}, d(z,\xi ) \le r \} \) is compact as a subset of \(\mathbb {R}\).

We will prove that the application \(\mathscr {F}\) attains its minimum in \(\mathcal {A}\).

Let \((z_k)_{k \in \mathbb {N}} \in \mathcal {A}\) be a sequence such that \(\lim _{k\rightarrow \infty } \mathscr {F} (z_k) = m \) with \( m = \inf _{\eta \in \mathcal {A}} \mathscr {F}(\eta )\). We show that \(d(z_k, \xi )\) is bounded then we find a minimizer.

Suppose \(d( z_k, \xi ) = \lambda > 1\) and let \(\gamma \) be the geodesic connecting \(\xi \) and \(z_k\) (\(\gamma (0) = \xi \) and \(\gamma (1) = z_k\)) satisfying (H2). Consider an element \(z = \gamma (t)\) on the geodesic \(\gamma \) such that \(d(z, \xi ) = 1\). By definition of a geodesic, \(d(z, \xi ) = d(\gamma (t), \gamma (0)) = |t-0|\cdot d(\gamma (0), \gamma (1)) = t d(z_k, \xi ) = \lambda \), implies necessarily \(t = 1/\lambda \) and using hypothesis (H2), we obtain:

$$\begin{aligned} \begin{aligned} l(z)&= l(\gamma (t)) = t l(\gamma (1)) + (1-t) l(\gamma (0)), \\&= \left( 1- \frac{1}{\lambda }\right) l(\xi ) + \frac{1}{\lambda } l(z_k). \end{aligned} \end{aligned}$$
(19)

For k large enough, \(z_k\) is such that \(\mathscr {F}(z_k) < 2m\) that is \(\frac{\lambda ^2}{2 {\varDelta }\tau } + l(z_k) < 2 m\). Rearranging Eq. (19) gives \(l(z_k) = \lambda l(z) - (\lambda -1) l(\xi )\), note that \(l(\xi )\) is constant and as \(d(z, \xi ) = 1\), hypothesis (H1bis) insures that l(z) is bounded, so \(l(z_k)\) is also bounded (by \(M \in \mathbb {R}\)). Consequently, there exists a constant C, depending on \(m, \tau \) and M such that \(\lambda < C\).

As \(d(z_k, \xi )\) is bounded, there exists a converging subsequence, \((x_k)_{k \in \mathbb {N}} \in \mathcal {A}\) such that \(d(x_k, \xi ) \rightarrow d_1\), and \(l(x_k)\) converges to \(m - \frac{d_1^2}{2 {\varDelta }\tau }\). Starting from some index \(k, x_k \in B(\xi , d_1)\) (the ball with center \(\xi \) and radius \(d_1\)) and if it is not the case, we replace \(x_k\) by the element at the intersection of the ball and the geodesic \(\gamma \) linking \(x_k\) and \(\xi \). This operation does not change the limit of \(l(z_k)\): noting temporarily \(\tilde{x}_k\) the obtained sequence, as l satisfies hypothesis (H2), we have \(l(\tilde{x}_k) = (1- \frac{d_1}{d(x_k, \xi )}) l(\xi ) + \frac{d_1}{d(x_k, \xi )}l(x_k)\) with \(d(x_k, \xi ) \rightarrow d_1\). By hypothesis (H1bis), the set \({l(x_k)}\) is closed, so there exists an element \(Z \in B(\xi , d_1)\) such that \(l(Z) = m - \frac{d_1^2}{2 {\varDelta }\tau }\) with, of course, \(d(Z, \xi ) \le d_1\). In this case, \(\mathscr {F}(Z) = \frac{d(Z, \xi )^2}{2 {\varDelta }\tau } + m - \frac{d_1^2}{2 {\varDelta }\tau } \le m\), so Z is a minimizer.

Let \(n \le 4, Y_1, \dots , Y_n \in \mathcal {A}, \lambda _1, \dots , \lambda _n \in \mathbb {R}\) and \(\mathcal {P}\) satisfying (H1) and (H2). We define the linear combination l of \(\mathcal {P}\) by: \(l(\eta ) = \sum _{i=1}^n \lambda _i \mathcal {P}(\eta , Y_i)\). It is immediate that l satisfies (H2); we prove below that if \(\mathcal {P}\) satisfies (H1) then l satisfies (H1bis). As all elements of the sum are in a compact set, the sum is bounded. Let \((z_k)_{k \in \mathbb {N}} \in \mathcal {A}\) be a sequence such that \(\lim _{k\rightarrow \infty } l(z_k) = \alpha \). Using (H1), there exists a subsequence of \(z_k\) (also noted \(z_k\)) such that the element \((\mathcal {P}(z_k, Y_1), \dots , \mathcal {P}(z_k, Y_n)) \in \mathbb {R}^n\) converges to \((\alpha _1, \dots \alpha _n)\) so \(l(z_k)\) converges to l(z) with \(l(z) = \sum _{i=1}^n \lambda _i \alpha _i = \alpha \). \(\square \)

Proposition 2.2

Let \(\mathcal {H}\) be a Hilbert space, \(f: \mathcal {H} \rightarrow \mathcal {H}\) be a function and \(\mathcal {P}(x, y) = \left\langle x,-f(y)\right\rangle \). Then the schemes VEE (Eq. 15), VH (Eq. 16), VRK3 (Eq. 17), and VRK4 (Eq. 18) correspond to Explicit Euler (Eq. 11), Heun (Eq. 12), Runge–Kutta 3 (Eq. 13) and Runge–Kutta 4 (Eq. 14) schemes respectively for the evolution Eq. (9).

Proof

We prove hypothesis (H1) in the case \(n=1\), the proof for the general case is left to the reader. We show that, for \(\xi \in \mathcal {H}, \mathcal {D} = \{ \left\langle z,-f(\xi )\right\rangle ; z \in {\mathcal H}, d(z,\xi ) \le r \} \) is bounded and closed. Denoting \(||\cdot ||\) the norm of \(\mathcal {H}\), by the triangular inequality when \(x \in B(\xi , r), ||x|| \le r + ||\xi ||\). For \(x \in \mathcal {D}, ||\left\langle x,f(\xi )\right\rangle || \le ||x||\cdot ||f(\xi )|| \le (r + ||\xi ||) \cdot ||f(\xi )|| < \infty \), so \(\mathcal {D}\) is bounded. Let z be an element of \(B(\xi ,r)\) and \(x \in \mathcal {H}\) such that \(z = \xi + x\) with \(||x|| \le r\), by linearity \(\left\langle z,-f(\xi )\right\rangle = \left\langle \xi ,-f(\xi )\right\rangle + \left\langle x,-f(\xi )\right\rangle \). Furthermore, \(\left\langle x,-f(\xi )\right\rangle \le ||x||\cdot ||f(\xi )|| \le r ||f(\xi )||\) so \(\{\left\langle z,-f(\xi )\right\rangle , z \in B(\xi , r) \}\) is included in the closed interval \([\left\langle \xi ,-f(\xi )\right\rangle - r ||f(\xi )||, \left\langle \xi ,-f(\xi )\right\rangle + r ||f(\xi )||]\). If we take \(x = - \lambda \frac{f(\xi )}{||f(\xi )||}\) with \(\lambda \in [-r,r]\) then \(\left\langle x,-f(\xi )\right\rangle = \lambda ||f(\xi )||\) which implies \(\left\langle z,-f(\xi )\right\rangle = \left\langle \xi ,-f(\xi )\right\rangle + \lambda ||f(\xi )||\). Therefore \(\{\left\langle z,f(\xi )\right\rangle , z \in B(\xi , r) \}\) is exactly \([\left\langle \xi ,f(\xi )\right\rangle - r ||f(\xi )||, \left\langle \xi ,f(\xi )\right\rangle + r ||f(\xi )||]\) and is closed as a closed interval of \(\mathbb {R}\).

To prove hypothesis (H2), consider \(x_0, x_1 \in \mathcal {H}\) and for \(t \in [0,1], \gamma (t) = t x_1 + (1-t) x_0\) the segment between \(x_0\) and \(x_1\). For \(t, s \in [0,1]\), with straightforward computation, we have \(d(\gamma (t), \gamma (s)) = |t-s| \cdot d(\gamma (0), \gamma (1))\), showing that \(\gamma \) is a (constant speed) geodesic. In this case, we can also prove the uniqueness of the geodesic. In fact, assume that there exists another geodesic \(\gamma _2\) linking \(x_0\) and \(x_1\). Let \(z = \gamma _2(t)\) be an element only on \(\gamma _2\); using the triangle inequality and that \(\gamma \) is the segment linking \(x_0\) and \(x_1\) we have \(d(x_0, z) + d(z,x_1) > d(x_0, x_1)\). Furthermore as \(\gamma _2\) is a geodesic, left member of previous equation is \(t d(x_0, x_1) + (1-t) d(x_0, x_1) = d(x_0, x_1)\). We obtain a contradiction thus there does not exist another geodesic. The function \(t \rightarrow \mathcal {P}(\gamma (t), Y)\) has the form \(t \rightarrow \left\langle t x_1 + (1-t) x_0,-f(Y)\right\rangle \), so is linear.

We now prove the equivalence of schemes. The scheme VEE is defined by the following equation:

$$\begin{aligned} \xi _{k+1} \in \underset{\eta \in \mathcal {H}}{\text {argmin }} \left\{ \frac{ ||\eta - \xi _k||^2}{2 {\varDelta }\tau } + \left\langle \eta ,f(\xi _k)\right\rangle \right\} . \end{aligned}$$

The application \(\mathscr {F}: \mathcal {H} \rightarrow \mathbb {R}, \mathscr {F} ( \eta ) = \frac{ ||\eta - \xi _k||^2}{2 {\varDelta }\tau } + \left\langle \eta ,f(\xi _k)\right\rangle \) is differentiable and its derivative \(\mathscr {F}': \mathcal {H} \rightarrow \mathcal {H} \) is \(\mathscr {F}' ( \eta ) = \frac{ \eta - \xi _k}{{\varDelta }\tau } + f(\xi _k)\). The minimizer of \(\mathscr {F}\) is a critical point, the minimum noted \(\xi _{k+1}\), satisfies \(\mathscr {F}'(\xi _{k+1})\) = 0. After trivial rearrangement, we obtain \(\xi _{k+1} = \xi _k - {\varDelta }\tau f(\xi _k)\), that is the Explicit Euler scheme presented in (11) for Eq. (9). Arguments are similar for other schemes. \(\square \)

3 Epidemiological model

From now on, all variables can have a double dependence: on the real time t which will be subscript and on the iteration time \(\tau \) (previously k in Eq. (7)) which will be a superscript l.

This section reproduces a short presentation of the epidemiological model and of the cost computing as introduced in [44, Sect. 2]. Please refer to the original article for more details and proof.

3.1 Spread of the disease

The model simulates the dynamics of an epidemic in a population. The final time period T is supposed finite, and the time horizon can be discretized in \((N+1)\) time instants, noted \(t_0 = 0, t_1 = {\varDelta }T, \ldots , t_N = T\).

The model is compartmental, and the population is divided into several classes: susceptible individuals \((S_n)\) is the proportion of individuals susceptible to catch the disease at time instant \(t_n\); infected individuals regroups individuals who are infected, more precisely, \(I_n^{\omega }\) is the proportion of individuals who have been infected at time instant \(t_{n-\omega }\) (with \(\omega \in \{0, 1, \ldots , {\varOmega }\}\)) and \(I_n\) is the sum of all \(I_n^{\omega }\); recovered individuals \((R_n)\) is the proportion of individuals who have exited from one infected class; vaccinated individuals (\(V_n^{\theta }\) where \(\theta \) counts the time between vaccination time instant and current time instant with \(\theta \in \{ 0, 1, \ldots , {\varTheta }\}\)) is the proportion of vaccinated individuals at time instant \(t_n\) and not infected since; failed vaccinated individuals (\(F_n\)) is the proportion of individuals who were vaccinated at some time \(t \le t_n\) but whose vaccination failed and have not been infected since. As individuals can only take vaccination once, \(V^{{\varOmega }}\) represents individuals who have lost immunity. Lastly \(U_n\) represents the proportion of individuals vaccinated between \(t_n\) and \(t_{n+1}\), a proportion f of them will never develop any immunity and go to the failed vaccination class.

The specificities of disease are represented by the vector \(\gamma = (\gamma ^0, \ldots , \gamma ^{{\varOmega }}) \in (\mathbb {R}_+)^{{\varOmega }+1}\) (which reflects the recovery of an infected individual); the function \(\beta (t)\) (which characterizes the contact between an infected individual and a susceptible one at time t); and the vector \(\alpha \in [0,1]\) which quantifies the protection provided by the vaccine in terms of the probability of infection.

The model is defined by the following system of equations:

$$\begin{aligned} S_{n+1}&= \left( 1- \beta _{{\varDelta }T}^nI_n\right) (S_n- U_n) \end{aligned}$$
(20a)
$$\begin{aligned} I_{n+1}^0&= \beta _{{\varDelta }T}^n\left[ F_n+ S_n+\sum _{\theta =0}^{N-1} \alpha _\theta V^\theta _n \right] I_n \end{aligned}$$
(20b)
$$\begin{aligned} I_{n+1}^{\omega +1}&= (1-\gamma _{{\varDelta }T}^\omega )I_{n}^{\omega } \qquad \omega =0, \dots , {\varOmega }-1 \end{aligned}$$
(20c)
$$\begin{aligned} V_{n+1}^0&= (1-f)\cdot \left( 1- \beta _{{\varDelta }T}^nI_n\right) U_n \end{aligned}$$
(20d)
$$\begin{aligned} V_{n+1}^{\theta +1}&= \left( 1- \beta _{{\varDelta }T}^n\alpha _{\theta }I_n \right) V_{n}^{\theta }, \qquad \theta =0, \dots , {\varTheta }-2 \end{aligned}$$
(20e)
$$\begin{aligned} V_{n+1}^{{\varTheta }}&= \left( 1- \beta _{{\varDelta }T}^n\alpha _{{\varTheta }-1}I_n \right) V_{n}^{{\varTheta }-1} + \left( 1- \beta _{{\varDelta }T}^nI_n \right) V_{n}^{{\varTheta }} \end{aligned}$$
(20f)
$$\begin{aligned} F_{n+1}&= f\cdot \left( 1- \beta _{{\varDelta }T}^nI_n\right) U_n + F_{n} \left( 1- \beta _{{\varDelta }T}^nI_n \right) \end{aligned}$$
(20g)

with initial conditions:

$$\begin{aligned} S_0&=S_{0^-},\quad I_0^\omega = I_{0^-}^\omega , \quad V^\theta _0=0, \quad \forall \theta \ge 0. \end{aligned}$$
(20h)

The continuous version of the previous discrete model is given by the following equations:

$$\begin{aligned}&S(t) = - \beta S(t) I(t) - U'(t) \end{aligned}$$
(21a)
$$\begin{aligned}&I(t) = \beta \left[ S(t) + F(t) + \int _0^{+ \infty } A(\theta ) V(t, \theta )d\theta \right] I(t) - \gamma I(t) \end{aligned}$$
(21b)
$$\begin{aligned}&\partial _t V(t, \theta ) + \partial _\theta V(t, \theta ) = - \beta A(\theta )V(t,\theta )I(t) \end{aligned}$$
(21c)
$$\begin{aligned}&F'(t) = fU'(t) - \beta F(t)I(t) \end{aligned}$$
(21d)

with initial and boundary conditions:

$$\begin{aligned} \begin{aligned}&S(0)=S_{0^-},\quad I(0) = I_{0^-}, \quad \forall t \ge 0 : S(t) \ge 0, \quad F(0^-) = 0, \\&\forall \theta \ge 0 : V(0^-, \theta ) = 0, \quad V(t, 0) = (1-f)U'(t). \end{aligned} \end{aligned}$$
(21e)

In this continuous time model, we use the following notation (refer to [44, Appendix] for details):

  • \(U : [0, \infty ] \rightarrow [0,1]\) represents the vaccination, U(t) is the fraction of the population that has been vaccinated by time t;

  • constants \(\beta , \gamma , f\) and function A, have the same meaning as in the discrete model;

  • \(V(t, \theta )\) is the fraction of individuals who have been vaccinated at time \(t-\theta \) and have not been infected since vaccination.

3.2 Individual vaccination: as minimization of individual cost

The system (20) presents the spread of the disease at the population level. For an individual, we consider the Markov chain \(M_n\) which describes, in terms of transition probabilities (see [44, Eq 2.10]), the state of the individual (\(M_n \in \{ S, R, F, I^j, V^l \}\) with \(j\in \{0, \dots , {\varOmega }\}\) and \(l \in \{ 0, \dots , {\varTheta }\}\)).

In this subsection, we look for the global decision of vaccination \(U_n\) which corresponds to the sum of individual decisions, as shown by Eq. (24). In order to reflect the decision of an individual, we introduce the probability law \(\mu \) defined on \(\{ t_0, \ldots , t_{N-1} \} \cup \infty \). In practice, before the epidemic starts, the individual chooses his probability of vaccination at each time step \(t_k\), assuming his non-infected status.

The collection of conditional rates \(\lambda = (\lambda _n)_{n=1}^N\) is given by the law \(\mu \):

$$\begin{aligned} \forall n\le N-1: \ \lambda _n = {\left\{ \begin{array}{ll} \displaystyle \frac{\mu _n}{\mu _n + \dots + \mu _\infty }, &{}\text { if } \mu _n + \dots + \mu _\infty > 0 \\ 0, &{}\text { otherwise} \end{array}\right. } \end{aligned}$$
(22)

The individual in the susceptible class at time \(t_0\) (\(M_0 = S\)) has the following cost: \(J_{indi}(\mu ; U) = \left\langle \mu ,\mathscr {C}_U\right\rangle \) where \(\mathscr {C}_U\) (the vector representing the cost of all pure strategies “vaccination sure at time t”) is defined by:

$$\begin{aligned} \mathscr {C}_U (t_n) = {\left\{ \begin{array}{ll} r_I\varphi ^I_n + (1- \varphi ^I_n)(r_V+ (1-f)r_I\varphi ^{V,I}_n) + r_If (\varphi ^I_\infty -\varphi ^I_n) &{} \text {for } n < N, \\ r_I\varphi ^I_\infty &{} \text {for } n = N. \end{array}\right. } \end{aligned}$$
(23)

Here \(\varphi ^{V,I}_n = 1- \prod _{k=n}^{{\varTheta }} \left( 1- \beta _{{\varDelta }T}^k\alpha _{k-n-1} I_k \right) ,\) with the convention \(\alpha _{-1}=1\) and \(\varphi ^I_n= 1-\prod _{k=0}^{n} \left( 1- \beta _{{\varDelta }T}^kI_k \right) (\text {for } n=0,\ldots ,N-1), \varphi ^I_\infty = 1-\prod _{k=0}^{N} \left( 1- \beta _{{\varDelta }T}^kI_k \right) \).

The equilibrium between individual dynamics and global dynamics (20) is attained when:

$$\begin{aligned} U_n = \lambda _n S_n. \end{aligned}$$
(24)

The purpose of an individual is to minimize \(J_{indi}(\mu ; U) \) under the constraint \(\mu \in {\varSigma }_{N+1}\) with \({\varSigma }_{N+1} = \left\{ x \in \mathbb {R}^{N+1} \left| \sum _{k=0}^N \right. x_k =1 \text { and } x_k \ge 0, 0 \le k \le N \right\} \).

The numerical schemes give a sequence of probability laws \(\mu ^l\); here the index “l” means that \(\mu ^l\) is a numerical approximation of the solution \(\xi \) of (8) at the time \(\tau _l = l\cdot {\varDelta }\tau \). For instance, if we consider the discrete model (20), each \(\mu ^l = (\mu ^l_0, \mu ^l_1, \dots , \mu ^l_{N-1}, \mu ^l_\infty ) \in {\varSigma }_{N+1}\) is a probability law on \(\{t_0; \dots , t_{N-1}\} \cup \{\infty \}\), while if we consider the continuous time model (21), each \(\mu ^l\) is a probability law on \( [0,T] \cup \{\infty \}\). Note that here a step of the numerical schemes increments the “l” index in order to converge to a Nash–MFG equilibrium.

3.3 Numerical values

This section describes the numerical values used in the two cases presented in this paper: Short persistence, large efficacy and Long persistence, 100% efficacy. These values are sensibly similar to the ones used in [44, Subsect. 4.2 and 4.3].

For the epidemic parameters, we consider a total simulation time at 1 year (\(T = 1\)); an initial proportion of susceptibles \(S_0 = 0.94\) and infected \(I_0 = 2.0 ^{-6}\); three time instants by day (\(N = 365*3\)); a recovery rate \(\gamma ^{\omega } = 365/3.2\) (\({\varOmega }= 20\)); the reproduction number \(R_0 = 1.35\) thus \(\beta = \gamma R_0\); \(t_2^{\beta } = 1/2\) such that \(\beta (t) = \beta \) for \(t \le t_2^{\beta }\) and then \(\beta (t) = \beta _{{ min}}\) for \(t > t_2^{\beta }\) where \(\beta _{\min } = \gamma /S_0\); the relative cost of the epidemic is \(r_I= 1\).

Finally, we introduce \(t_1^{\alpha }, t_2^{\alpha }\) and set \(\alpha _{\theta } = 1 - \mathbf {1}_{[t_1^{\alpha }, t_2^{\alpha }]}\). For the case Short persistence, large efficacy, \(t_1^{\alpha } = 5/365\) and \(t_2^{\alpha }=1/12\) (\({\varOmega }= 93\)). For the case Long persistence, 100% efficacy, \(t_1^{\alpha } = 0\) and \(t_2^{\alpha }=1/2\) (\({\varOmega }= 549\)). We suppose a failure rate \(f=0\), and the relative cost of the vaccination \(r_V= 0.005\).

4 Numerical illustration

4.1 Framework

Salvarani and Turinici introduced an epidemiological model with possibility of vaccine which has imperfect efficiency and limited persistence recalled in the previous Sect. 3. To model the spread of the disease in the population, they use a compartmental model and a probability distribution to reflect the individual decisions.

To find stable individual decision they define the problem as a Nash equilibrium. They find the probability distribution \(\mu \) making an individual indifferent to change his vaccination decision if all individuals have the same \(\mu \) (see [44, Theorem 2.1] for the proof of equilibrium existence).

Let \(\mathscr {C}_{\mu }\) be the cost of pure strategies “vaccination happens at time t” when all individuals choose as vaccination strategy \(\mu \). In that case \(\left\langle \eta ,\mathscr {C}_{\mu }\right\rangle \) represents the cost of an individual with strategy \(\eta \) when others use strategy \(\mu \).

The definition of the mapping \(E(\mu )\) as introduced in [44] is the maximum gain obtained by an individual if he changes unilaterally his strategy and everybody else remains with the strategy \(\mu \); with mathematical notation:

$$\begin{aligned} E(\mu ) = \left\langle \mu ,\mathscr {C}_{\mu }\right\rangle - \min _{\eta \in {\varSigma }_{N+1}} \left\langle \eta ,\mathscr {C}_{\mu }\right\rangle . \end{aligned}$$

A minimum of the mapping \(\mu \rightarrow E(\mu )\) is a Nash equilibrium. To find it they introduce an iterative method depending on a step \({\varDelta }\tau \) (algorithm 1). In this algorithm, Eq. (25) is used. The following intuitive interpretation is also provided: an individual in a population with strategy \(\mu ^l\) will, if necessary, adjust his strategy to minimize his cost function \(\eta \rightarrow \left\langle \eta ,\mathscr {C}_{\mu ^l}\right\rangle \) while at the same time keeping small the distance between the previous strategy \(\mu ^l\) and the new \(\mu ^{l+1}\). But Algorithm 1 can also be seen as numerical resolution of an evolution e.g., describing a curve in the metric space of the admissible strategies. That is why we apply variational methods in the metric space \({\varSigma }_{N+1}\) with the standard Euclidean distance to obtain faster convergence to the Nash equilibrium. Here \(\mathcal {P}(x, y) = \left\langle x,\mathscr {C}_y\right\rangle \) and for the same reasons as in the proof of Proposition 2.2, \(\mathcal {P}\) satisfies hypotheses (H1) and (H2). Note that with high order schemes VH, VRK3 and VRK4, we lose the intuitive idea but we increase the Nash equilibrium computing.

figure b

4.2 Results

In this section we test the variational schemes proposed in Sect. 2.2.

The reader is invited to refer to [44, Sect. 4.6] or Sect. 3 for the numerical values of the two studied cases Short persistence, large efficacy (with vaccine persistence at 1 month with a delay of action at 5 days) and Long persistence, 100% efficacy (with persistence of the vaccine at 6 months with no delay of action).

In order to appreciate the convergence scheme for the case Short persistence, large efficacy the evolution of the mapping \(E(\cdot )\) is presented (see Fig. 6 for VEE with \({\varDelta }\tau = 0.1\)), where the evolution for the four schemes with lower \({\varDelta }\tau \) value is added to stabilize result of mapping \(E(\cdot )\).

Fig. 6
figure 6

The evolution of the mapping \(E(\cdot )\) for the four schemes introduced in Sect. 2.2 in the case Short persistence, large efficacy for the model (20a)–(20h). In [44] VEE with \({\varDelta }\tau = 0.1\) is used

For the Short persistence, large efficacy case, Fig. 7 presents the numerical estimation of the scheme error and order: only the VEE is of order one, the others are of order two in \({\varDelta }\tau \). Furthermore, Fig. 8 introduces the obtained errors for each schemes depending of the number of minimizations. Recall that a high order scheme needs a very regular functional to provide high order convergence, but in our case the regularity of the function is completely unknown. However, this case remains very interesting because it shows that even if the regularity of the function is not enough to obtain order three or four, the VRK3 and VRK4 still improve the regularity of the function and have a better order than the VEE scheme.

Fig. 7
figure 7

Numerical estimated error and order obtained for the four schemes for the case Short persistence, large efficacy of the model (20a)–(20h). Reference solution (\(\mu ^{{\varDelta }\tau _{{ ref}}}\)) is given by a VRK4 scheme with \({\varDelta }\tau _{{ ref}} = 0.001\). Other steps used are \({\varDelta }\tau = 0.04, 0.06, 0.08, 0.10\) and 0.12. Denote \(\mu _n^{{\varDelta }\tau }\) the probability of vaccination at time n / N given by the last iteration (\(10/{\varDelta }\tau \)) of the numerical scheme (VEEVHVRK3 or VRK4) with parameter \({\varDelta }\tau \). We plot: (left) the error \( \sqrt{\sum _{n = 1}^{N} \left( \mu _n^{{\varDelta }\tau } - \mu _n^{{\varDelta }\tau _{{ ref}}}\right) ^2 + \left( \mu _{\infty }^{{\varDelta }\tau } - \mu _{\infty }^{{\varDelta }\tau _{{ ref}}}\right) ^2 }\). Right: order of convergence

Fig. 8
figure 8

Numerical estimated error obtained (see Fig. 7 for the definition of the error) for the Short persistence large efficacy case of the model (20a)–(20h). This figure shows for each of the four detailed schemes the error depending on number of minimisations. Reference solution is given by a VRK4 scheme with \({\varDelta }\tau _{{ ref}} = 0.001\). Other steps used are 0.04, 0.06, 0.08, 0.10 and 0.12

Fig. 9
figure 9

Mapping \(E( \cdot )\) for the case Short persistence, large efficacy and Long persistence, 100% efficacy (of the model (20a)–(20h)), generated with VEE scheme and \({\varDelta }\tau = 0.1\)

Fig. 10
figure 10

Numerical estimated error and order obtained for the four schemes and the case Long persistence, 100% efficacy of the model (20a)–(20h). Reference solution is given by a VRK4 scheme with \({\varDelta }\tau _{{ ref}} = 0.005\). Other steps used are \({\varDelta }\tau = 0.04, 0.06, 0.08, 0.10\) and 0.12. Denote \(\mu _n^{{\varDelta }\tau , l}\) the probability of vaccination at time n / N given by the lth iteration of the numerical scheme (VEEVHVRK3 or VRK4) with parameter \({\varDelta }\tau \). We plot: (left) the error \(\max _{l\le 10/ {\varDelta }\tau } \sqrt{\sum _{n = 1}^{N} \left( \mu _n^{{\varDelta }\tau , l} - \mu _n^{{\varDelta }\tau _{{ ref}}, l}\right) ^2 + \left( \mu _{\infty }^{{\varDelta }\tau , l} - \mu _{\infty }^{{\varDelta }\tau _{{ ref}}, l}\right) ^2 }\). Right: order of convergence

For the other example Long persistence, 100% efficacy, the convergence is faster: Fig. 9 compares the evolution of mapping \(E(\cdot )\) for the two cases. In the second case, the graph in Fig. 10 shows that the scheme has high order as it was expected previously. VH is indeed of order 2 and VRK3 and VRK4 are respectively order 3 and 4. Furthermore, Fig. 11 introduces the error for each scheme depending on the number of minimisations: this shows the accuracy of the schemes.

Fig. 11
figure 11

Numerical estimated error obtained (see Fig. 10 for the definition of the error) for the Long persistence, 100% efficacy case of the model (20a)–(20h). This figure shows for each of the four detailed schemes the error depending on number of minimisations. Reference solution is given by a VRK4 scheme with \({\varDelta }\tau _{{ ref}} = 0.005\). Other steps used are 0.04, 0.06, 0.08, 0.10 and 0.12

5 Second application: an economical approach

5.1 Epidemiological and economic model

This numerical application is based on [22]. In order to model the spread of the disease, three classes (as explained in Sect. 1.1) are introduced: Susceptible (S), Infected (I), Recovered (R) with the following time evolution:

$$\begin{aligned} \begin{aligned} \dot{S}(t)&= - \beta \frac{S(t) I(t)}{N} - r(t) \\ \dot{I}(t)&= \beta \frac{S(t) I(t)}{N} - \gamma I(t) \\ \dot{R}(t)&= r(t) + \gamma I(t) \end{aligned} \end{aligned}$$
(26)

where S(t), I(t), R(t) represent respectively the number of susceptible, infected and recovered individuals at time tN the total number of individuals and r(t) the rate of vaccination at time t. We also introduce the utility \(u_H\) of a healthy individual and the utility \(u_S\) when sick.

The individual cost introduced is represented by the expected loss if individual is falling ill \(\left( EL = \frac{u_H-u_S}{\gamma } \left( 1 - \frac{S(\infty )}{S(0)} \right) \right) \) and the cost of vaccination \(\theta \).

The computation that follows, taken from [22], is valid when the dynamics (26) traverses at most two regions of vaccination/non vaccination. Under this assumption, the expected payoff of the two pure strategies (vaccination or not) is equal on the frontier of the vaccination/non vaccination region, that is: \(\theta = \frac{u_H-u_S}{\gamma } \left( 1 - \frac{S(\infty )}{S(0)} \right) \).

Set \(w = \theta \gamma /(u_H-u_S)\). Rewrite previous equation under the following form: \(w = 1 - S(\infty )/S(0)\), that gives \(S(\infty ) = S(0) (1-w)\). Using that, with \(\rho = N\gamma /\beta , \frac{d}{dt} S(t)e^{R(t)\rho } = 0\) for \(r(t) = 0 \ \forall t \ge 0, S(0)e^{R(0)\rho } = S(\infty )e^{R(\infty )\rho }\) (see [1]), and after replacing of \(S(\infty )\) by \(S(0) (1-w)\), we obtain: \(R(0) = \rho \ln (1-w) + R(\infty )\). As the total number of the population is constant equal to N, we obtain Eq. (27). This last equation can be seen as a “switching curve” in the plane (SI):

$$\begin{aligned} {\varGamma }= \left\{ (S, I) \left| I = - N \frac{\gamma }{\beta } \ln \left( 1-\frac{\theta \gamma }{u_H - u_S}\right) - \frac{\theta \gamma }{u_H - u_S} S \right. \right\} , \end{aligned}$$
(27)

where the individual is indifferent between being vaccinating and doing nothing.

With this economical approach, there is no insurance that the curve \({\varGamma }\) generates a Nash–MFG equilibrium. In fact, [22] does not use Nash–MFG framework but assumes equality between two payoffs, false in general cases (see for example [24]).

Fig. 12
figure 12

Numerical error and order for the four schemes and the model (26). Reference solution (\(\mu ^{{\varDelta }\tau _{{ ref}}}\)) is given by a VRK4 scheme with \({\varDelta }\tau _{{ ref}} = 0.005\). Other steps used are \({\varDelta }\tau = 0.04, 0.06, 0.08, 0.10\) and 0.12. Denote \(\mu _n^{{\varDelta }\tau , l}\) the probability of vaccination at time n / N given by the lth iteration of the numerical scheme (VEEVHVRK3 or VRK4) with parameter \({\varDelta }\tau \). We plot: (left) the error \(\max _{l\le 10/{\varDelta }\tau } \sqrt{\sum _{n = 1}^{N} \left( \mu _n^{{\varDelta }\tau , l} - \mu _n^{{\varDelta }\tau _{{ ref}}, l}\right) ^2 + \left( \mu _{\infty }^{{\varDelta }\tau , l} - \mu _{\infty }^{{\varDelta }\tau _{{ ref}}, l}\right) ^2 }\). Right: order of convergence

Nevertheless we seek to replicate these results in order to benchmark the performance of our numerical schemes. For each initial point, the algorithm reports if vaccination is optimal or not, which can be compared to the region given by the switching curve \({\varGamma }\).

We use the same approach as in Sect. 3.2 to determine the global vaccination response r(t). The proportion of individuals vaccinated up to time t is the cumulative function V(t) defined by:

$$\begin{aligned} V(t) = \int _0^t r(u) du, \end{aligned}$$
(28)

where r(u) is the instantaneous vaccination rate at time u with \(r(u) \in [0, r_{{ max}}]\) and \(r_{{ max}} < \infty \).

5.2 Numerical values

For the simulation, we use the following numerical values: an initial proportion of susceptibles \(S_0 = 0.75\) and infected \(I_0 = 0.1\); a cost of the vaccine \(\theta = 0.5\); an expected utility loss associated if illness \(\frac{\gamma }{u_H - u_S} = 1\); a recovery rate for an infected individual \(\gamma = 365/10\) and a transmission rate of the disease \(\beta = 73\). As in the Sect. 4, total simulation time is 1 year (\(T=1\)) with \(L = 365 \times 3\) total time steps.

5.3 Results

We obtain coherent results with the analytic solution by using the same approach as in Sect. 4.1. The equilibrium individual strategy is a mixed strategy with \(\xi _0 = 66 \%\) and \(\xi _{\infty } = 33 \%\) and its cost is 0.51. This results are coherent with the analytic result that provides a mixed strategy with \(\xi _0 = 34\%, \xi _\infty = 66\%\) and a cost at 0.5. Figure 12 presents the numerical estimation errors and order for the four schemes and Fig. 13 the errors for each schemes depending of the number of minimizations.

Fig. 13
figure 13

Numerical error obtained (see Fig. 12 for the definition of the error) for the four schemes depending on the number of minimizations for the model (26). Reference solution is given by a VRK4 scheme with \({\varDelta }\tau _{{ ref}} = 0.005\). Other steps used are 0.04, 0.06, 0.08, 0.10 and 0.12

6 Perspectives

This work introduces three high order schemes to find a Nash equilibrium and illustrates their use in an epidemiological application.

We show numerically that the schemes VH, VRK3 and VRK4 exhibit better order than VEE and we can obtain, depending on the regularity of the function, up to order four convergence. The approach based on the Runge–Kutta method can be applied for other numerical schemes, for instance Mid Point, Leap–Frog or Adams–Bashford methods. The presented algorithms are not optimized, as minimizations can take a long time and need a high numerical precision, so the time execution can be significant. A perspective of this work is to provide an extension to bi-dimensional or tri-dimensional problems.