1 Introduction

Human interactions supply many phenomena, which can be modelled and analysed with applied mathematics. One of the most famous models, developed in this context, concerns the role of the military strategy (and the linked decision problems) during a conflict between two or more armies. It is strange to think that, despite the numerous wars from ancient time up to nowadays, a rigorous mathematical modeling of warfare strategies is only fairly recent. In fact, the construction of a consistent theory reconciling the intuitive idea with the mathematical formalism only came close to World War I. At that time, Frederick Lanchester (1914) published a first system of differential equations to model the clash between two homogeneous armed forces (i.e., armed forces each of which has only one type of weapon system at their disposal). It is historically relevant to notice that, in the same years Lanchester, Chase, see Fisk (1905) and Fisk and Chase (1916), and Mikhail Osipov, see Osipov (1915) and Helmbold (1993), had produced similar differential equation models to describe the evolution of a battle, independently of each other. For further historical details, see Jaiswal (1997).

As we shall explain extensively in the next section, Lanchester provided a useful characterization of the dynamics of the battle by using two simple differential equations. Starting from these equations, he proved an interesting result, the so called Lanchester’s square law, which provides important information concerning the numerical evolution of the two forces. In the following decades, Lanchester’s law was studied by many authors as Bracken (1995), Davis (1995), Fox (2010), Taylor (1983). Over the years, various generalization were proposed, starting from the action of mixed forces due to Lepingwell (1987). Most of these attempts are summarized by MacKay in MacKay (2006). Moreover, the applications of Lanchester’s type models goes beyond the original purpose, as its core ideas were used to describe human, by Johnson and MacKay (2015), and nonhuman, by Clifton (2020), evolution. Recently, the interest of such a model made a comeback in the operational research literature as series of works proposing alternative stimulating models. In particular, a game theoretical approach is given in Kress et al. (2017, 2018). This approach exploits a matrix formalism similar to the one used in this work, as well as a network-theoretical approach used in Kalloniatis et al. (2021).

The asymptotic dynamics were thoroughly analysed for two heterogeneous (mixed) forces by MacKay (2009), Lin and MacKay (2014), in order to determine the optimal fire distribution policy, and for models that capture irregular warfare, such as insurgencies studied by Kress (2020). The authors Liu et al. (2012) moved some criticism to MacKay (2009), to which MacKay replied stressing the strengths of his approach. Lastly, we mention another sophisticated attempt which involves the use of partial differential equations with time variables developed by Spradlin and Spradlin (2007).

Given the significance of such a model, we believe it is worth providing an accurate formal analysis and a precise model in terms of dynamical systems research. Thus, the aim of this paper is to provide a mathematically consistent generalization of the classical model that sorts out the critical issues arising from the original construction. This construction allows us to describe general higher dimensional scenarios.

The paper is organized as follows. In Sect. 2, we summarize in a mathematical framework the classical Lanchester’s model and its related square and linear laws. Section 3 contains our main contribution to the modelling effort. It is devoted to the proposal of Lanchester’s type models with more than two homogeneous armed forces. The attempt here is to provide a more mathematically coherent description of the dynamics of the battles, which is consistent with a realistic interpretation. As we shall see, there are various possible combinations, which lead to different interpretations of the real-world scenarios we aim to describe. Our analysis moves between these different scenarios, supplying various results which are validate with numerical simulations. Finally, in Sect. 4, we take stock of the work, arguing possible improvements of our models and discussing leads for future developments.

2 Classical Lanchester’s model

We begin by reviewing the well-known Lanchester equations, first formulated by Lanchester (1914), and then extensively developed by Taylor (1983), MacKay (2006, 2009), Kress et al. (2017, 2018). The aim of these equations was to model the change in combat power of two enemy factions facing each other. Throughout the paper, we indicate with \(F'\) the derivative of the function F with respect to the time variable t. We do not normalize the size of the armies, as it brings no advantage to our analysis, and we pick their values at time \(t=0\) in order to produce clearly interpretable figures.

Lanchester’s model can be summarized by two systems of differential equations, one for the aimed fire:

$$\begin{aligned} {\left\{ \begin{array}{ll} \displaystyle {R'=-bB},\\ \displaystyle {B'=-rR}, \end{array}\right. } \end{aligned}$$
(1)

and one for the unaimed one:

$$\begin{aligned} {\left\{ \begin{array}{ll} \displaystyle {R'=-bRB},\\ \displaystyle {B'=-rRB}. \end{array}\right. } \end{aligned}$$
(2)

here \(B=B(t)\) and \(R=R(t)\) represent the size of the two factions (Blue army and Red army) at any time \(t>0\); whereas \(b,r>0\) are coefficients which indicate the fighting ability of each unit in the corresponding formation, supposed to be kept both constant throughout the battle. Due to their interpretation, it is clear that we can only accept \(B(t)\ge 0\) and \(R(t)\ge 0\) for all \(t\ge 0\). The analysis of these systems provide the celebrated Lanchester’s square and linear law respectively.

2.1 Lanchester’s square law

We start by recalling classical results for system (1). The interested reader can find the proofs of these results and other insights in Taylor (1974, 1983).

Theorem 2.1

(Lanchester’s Square Law). The solutions to (1) satisfy the state equation

$$\begin{aligned} bB^2(t)-rR^2(t) = K, \end{aligned}$$
(3)

where K is a constant that depends only on the fighting abilities of the two armies, b and r, and on their initial size B(0) and R(0).

The claim follows by differentiating (3) with respect to the time variable t along solutions of (1). Equation (3) allows an a priori estimate of the winning side: if \(K>0\), then only R(t) can ever be zero and thus faction B will win.

It is important to notice the presence of squares in (1). It indicates that the outcome of the battle is more sensitive to the changes in the number of units rather than to changes in their effectiveness.

While Eq. (3) is interesting for practical application, from a mathematical point of view we may be interested in solving system (1) explicitly.

Lemma 2.2

The solutions of (1) are:

$$\begin{aligned} {\left\{ \begin{array}{ll} B(t)= \displaystyle {B(0) \cdot \cosh \left( \sqrt{br} \cdot t\right) -R(0)\sqrt{\frac{r}{b}} \sinh \left( \sqrt{br} \cdot t\right) },\\ \\ R(t)= \displaystyle {R(0) \cdot \cosh \left( \sqrt{br} \cdot t\right) -B(0) \sqrt{\frac{b}{r}} \sinh \left( \sqrt{br} \cdot t\right) }. \end{array}\right. } \end{aligned}$$
(4)

In particular, in the case \(rR^2(0) = bB^2(0)\), the solutions to (1) are reduced to:

$$\begin{aligned} {\left\{ \begin{array}{ll} B(t)=B(0)\exp (-\sqrt{br}t), \\ R(t)=R(0)\exp (-\sqrt{br}t). \end{array}\right. } \end{aligned}$$

Notice that, since \(\cosh\) and \(\sinh\) are asymptotically equivalent, then the behaviour of system (4) when \(t\rightarrow +\infty\) is completely determined by the initial constants B(0), b, R(0), and r. In particular whenever \(bB(0)^2 \ne rR(0)^2\) one of the two will diverge to \(+\infty\), while the other will diverge to \(-\infty\), see Fig. 1a. This behaviour could be avoided by artificially introducing a constraint on the non-negativity of B and R in the corresponding ODEs.

2.2 Lanchester’s linear law

For the case of the linear law described in system (2), analogous results to Theorems 2.1 and 2.2 hold. An in-depth study of these results with the proofs and other applications can be found in Taylor (1973, 1983).

Theorem 2.3

(Lanchester’s Linear Law). The solutions to (2) satisfy the state equation

$$\begin{aligned} bB(t)-rR(t) = K, \end{aligned}$$
(5)

where K is a constant that depends only on the fighting abilities of the two armies, b and r, and on their initial size B(0) and R(0).

As for Theorem 2.1, the proof of this theorem can be obtained by direct derivation with respect to the time variables along solutions of (2). The interpretation of Eq. (5) is fundamentally different from the one of Eq. (3). The linear law does not gives a predominant role to the size of the army with respect to his combat power. In the case of the square law, a lack of combat expertise could be easily balanced by rising the number of soldiers by a small amount (i.e., increasing B(0) or R(0), respectively), whereas in the case of the linear law this is no longer true.

As for the case of the square law, we may be interested in the mathematical solution of (2).

Lemma 2.4

Let us suppose \(rR(0) \ne bB(0)\), then the solutions to (2) are:

$$\begin{aligned} {\left\{ \begin{array}{ll} B(t)=\displaystyle {\frac{(rB(0)R(0)-bB(0)^2)\exp {((bB(0)-rR(0))t)}}{rR(0) - bB(0)\exp {((bB(0)-rR(0))t)}}}, \\ \\ R(t)=\displaystyle {\frac{rR(0)^2-bB(0)R(0)}{rR(0) - bB(0)\exp {((bB(0)-rR(0))t)}}}. \end{array}\right. } \end{aligned}$$
(6)

Moreover, in the case \(rR(0) = bB(0)\), the solutions to (2) are reduced to:

$$\begin{aligned} {\left\{ \begin{array}{ll} B(t)=\displaystyle {\frac{B(0)}{1+bB(0)t}}, \\ \\ R(t)=\displaystyle {\frac{R(0)}{1+rR(0)t}}. \end{array}\right. } \end{aligned}$$

Notice that, in contrast with the solutions of (1), one of the solutions of (2) tends to a positive limit, while the other one tends to 0 as \(t\rightarrow \infty\), see also Fig. 1b.

Fig. 1
figure 1

Evolution in time of a solutions of (1), plotted together with (3), with B(0) and R(0) chosen randomly in (0, 10) and \(r=0.2\), \(b=0.9\); and b solutions of (2), as well as the constant of motion (5), with B(0) and R(0) as in Fig. 1a and \(r=0.2\), \(b=0.9\). Notice that in a R(t), the size of the losing army, quickly becomes negative and keeps decreasing. In turn, the negative sign of R makes B increase and asymptotically diverge to \(+\infty\). Instead, in b R(t), the size of the losing army, quickly approaches 0 as the system asymptotically converges to its equilibrium \(B_{\infty }=B(0)-rR(0)/b\), \(R_{\infty }=0\)

3 Multi-battle model

In this section, we provide multi-dimensional generalizations of systems (1) and (2), together with a detailed analysis of the latter. Such a generalization can model a multilateral warfare between n armies, if we assume all the attacks are carried out simultaneously. For ease of notation, in the following we shall use \(X_1, X_2, \dots , X_n\) to represent the size of the armies involved. A useful guide that provide all the necessary basics for the following discussion is Horn and Johnson (2013).

Consider a column vector \(X=(X_1, X_2, \dots , X_n)^T\), representing all the n armies involved; assume each army has a fighting ability \(x_i\), \(i=1,2,\dots ,n\), that we collect in the vector \(x=(x_1,x_2,\dots ,x_n)\).

Definition 3.1

We define the conflict matrix \({\tilde{C}}\) as the matrix whose entries are

$$\begin{aligned} {\tilde{c}}_{i,j}={\left\{ \begin{array}{ll} 1 &{} \text {if army }X_i\text { and }X_j\text { are in conflict},\\ 0 &{} \text {otherwise}. \end{array}\right. } \end{aligned}$$

The associated weighted conflict matrix C is obtained by re-scaling \({\tilde{C}}\), normalizing the column sum to 1. Its entries are

$$\begin{aligned} c_{i,j}=\dfrac{{\tilde{c}}_{i,j}}{\sum _{k=1}^n {\tilde{c}}_{k,j}}, \end{aligned}$$

where the re-scaling takes into account the homogeneous division of forces of each army into the various conflicts it is fighting. We discuss other re-scaling options in Sect. 4.

Notice that the matrix \({\tilde{C}}\) is always symmetricFootnote 1 and is telling us who is fighting against whom. The matrix C is generally not symmetric, and it represents the ongoing conflicts and the distributions of the various armies. A generalization of this construction taking into account an heterogeneous distribution of forces, or unilateral conflicts, i.e. a more general matrix C, could be used to describe more general scenarios. We remark that the matrix \({\tilde{C}}\) can be interpreted as the adjacency matrix of a graph with n nodes (the armies), whose edges represent existing conflicts. The irreducibility assumption on \({\tilde{C}}\) impedes a splitting of the ongoing conflicts into two or more sub-conflicts, as it coincides with assuming that the corresponding graph is connected.

Moreover, we introduce the diagonal matrix D, whose diagonal entries are the parameters \(d_{i,i}=x_i\). The n-dimensional generalization of system (1) is given by

$$\begin{aligned} X'=-CD X. \end{aligned}$$
(7)

We only consider solutions for which \(X(t)\ge 0\) for all \(t\ge 0\), meaning all the entries in the vector X must remain non-negative. With the same notation, system (2) generalizes to

$$\begin{aligned} X_i'=(-CD X)_i X_i, \quad i=1,2,\dots ,n. \end{aligned}$$
(8)

Likewise, we only consider solutions for which \(X_i(t)\ge 0\) for all \(t\ge 0\). Clearly, for \(n=2\), we get

$$\begin{aligned} C= {\tilde{C}} = \left( \begin{matrix} 0 &{} \quad 1 \\ 1 &{} \quad 0 \end{matrix} \right) , \end{aligned}$$

thus, system (7) reduces to system (1), and system (8) reduces to system (2). However, in general, neither (7) nor (8) admit constants of motion, which were used in the original models to predict which side would win the conflict, based solely on X(0) and x.

We remark that for systems (1) and (7), \(X_i(t)=0\) does not stop the evolution in time for the corresponding army. Indeed, if we let (1) evolve past the moment in which e.g. \(R=0\) and \(B>0\), we would see \(R\rightarrow -\infty\) and \(B\rightarrow +\infty\), as in Fig. 1a. This clearly represents a issue of such a simple model, which is inherited by its n-dimensional generalization (7).

On the contrary, both in system (2) and (8), \(X_i({\bar{t}})=0\) implies that, for any \(t>{\bar{t}}\), \(X_i(t)=0\). Of course, the losing army in (2) only converges to 0 as \(t\rightarrow +\infty\). However, if we set the initial condition of one of the two to 0, and the other to any positive constant, the system is at equilibrium, as we would expect from the interpretation of a conflict with only one army involved. The same holds true in the n-dimensional case: an army starting with 0 soldiers will not be affected by the conflict between the remaining armies, as is biologically reasonable to happen.

In the following, we distinguish three cases: the all-fight-all case, in which the matrix \({\tilde{C}}\) has all 1 entries except for the zeros on the diagonal; the army-of-one case, in which the matrix \({\tilde{C}}\) has 1 entries only in the first row and column, except the entry (1, 1), and every other case. We will show that, in the all-fight-all scenario, only one winner will emerge from the conflict; while in the army-of-one case there can be up to \(n-1\) winners (interpreted as asymptotically surviving armies). In all other situations, the outcome is highly case-dependent, as we will show for the case when \(n=4\).

We remark that, although we were able to recover a constant of motion for some particular cases, it is our opinion that for higher dimensions, i.e. high number of armies involved, the existence of such a constant is not guaranteed. Indeed, denoting for ease of notation \(A:=CD\), system (7) has \(X^TQX\) as a constant of motion if and only if \(A^TQ + QA = 0\). Similarly, (8), which can be written as \(X_i'=-(\sum _j A_{ij}X_j)X_i\) in this new notation, has \(\sum _i q_iX_i\) as a conserved quantity if and only if \(q_i A_{ij} + q_j A_{ji}=0\) for all ij. For \(n=2\), the well known constants of motion for the square and linear laws can be derived from these more general formulas.

3.1 All-fight-all scenario

We now look at Eqs. (7) and (8) for the case \(n=3\). Since we are in the all-fight-all scenario, each army is facing all the other armies, and we have that

$$\begin{aligned} C =\dfrac{1}{2}\begin{pmatrix} 0 &{} \quad 1 &{} \quad 1\\ 1 &{} \quad 0 &{} \quad 1\\ 1 &{} \quad 1 &{} \quad 0 \end{pmatrix}. \end{aligned}$$
(9)

Hence system (7) becomes

$$\begin{aligned} {\left\{ \begin{array}{ll} X_1'&{}= \frac{1}{2} \left( -x_2 X_2 -x_3 X_3 \right) ,\\ X_2'&{}= \frac{1}{2} \left( -x_1 X_1 -x_3 X_3 \right) ,\\ X_3'&{}= \frac{1}{2} \left( -x_1 X_1 -x_2 X_2 \right) , \end{array}\right. } \end{aligned}$$
(10)

whose only equilibrium is (0, 0, 0). System (10) inherits the same critical issue of system (1): one (see Fig. 2a) or two armies (see Fig. 2b) can become negative, while the remaining diverge to \(+\infty\).

Fig. 2
figure 2

Evolution in time of the solutions of (11), exhibiting a one negative and two positive diverging solutions, and b two negative and one positive diverging solutions. This fundamentally unrealistic behaviour can not be avoided with a simple linear system. The values relative to the Square Law (3) are displayed as the title of the figure, although they do not give useful information on the asymptotic behaviour of the middle solution

In the remainder of this article, in order to overcome the critical issue of solutions diverging to \(\pm \infty\), we will only focus on properties of generalized versions of (2). In the following, we completely characterize the asymptotic behaviour of its 3D generalization. Using (9), system (8) becomes

$$\begin{aligned} {\left\{ \begin{array}{ll} X_1'&{}= \frac{1}{2} \left( -x_2 X_2 -x_3 X_3 \right) X_1,\\ X_2'&{}= \frac{1}{2} \left( -x_1 X_1 -x_3 X_3 \right) X_2,\\ X_3'&{}= \frac{1}{2} \left( -x_1 X_1 -x_2 X_2 \right) X_3. \end{array}\right. } \end{aligned}$$
(11)

Any point of the form \((X_1^*,0,0)\), \((0,X_2^*,0)\) or \((0,0,X_3^*)\) is an equilibrium of (11). However, analysis of local stability at any of them provides one zero eigenvalue and two equal non-positive ones, namely \(-x_iX_i^*/2\); hence, local stability analysis does not provide useful information. We will now directly proceed with a global stability analysis.

Since (11) is non-increasing in all its variables, and bounded (assuming non-negative initial conditions), each \(X_i(t)\) will admit a finite non-negative limit as\(t \rightarrow +\infty\).

We are interested in predicting which army will win, given only the initial conditions \((X_1(0),X_2(0),X_3(0)) \in {\mathbb {R}}^3_{\ge 0}\) and the fighting abilities \((x_1,x_2,x_3) \in {\mathbb {R}}^3_{> 0}\). We will focus on the case \((X_1(0),X_2(0),X_3(0)) \in {\mathbb {R}}^3_{> 0}\). Indeed if the initial condition is \(X_i(0)=0\) for one i, system (11) reduces to (2); if \(X_i(0)=X_j(0)=0\) for \(i\ne j\), the system is at equilibrium, and the same in the case \(X_1(0)=X_2(0)=X_3(0)=0\), as one would expect.

Proposition 3.2

Assume that in system (11) initial conditions and parameters are such that \(x_iX_i(0)>x_jX_j(0)\ge x_kX_k(0)\), for some permutation of \(i,j,k\in \{1,2,3\}\). Then

$$\begin{aligned} \lim _{t\rightarrow +\infty }X_i(t)>0 \quad \text {and} \quad \lim _{t\rightarrow +\infty }X_j(t)=\lim _{t\rightarrow +\infty }X_k(t)=0. \end{aligned}$$

Proof

Let us suppose, without loss of generality, that \(x_1X_1(0)>x_2X_2(0)\ge x_3X_3(0)\). The proof for the remaining five cases is identical.

First of all, we remark that

$$\begin{aligned} X_1'=\frac{1}{2}(-x_2 X_2 -x_3 X_3)X_1\ge \frac{1}{2}(-x_2 X_2(0) -x_3 X_3(0))X_1, \end{aligned}$$

since all the solutions are non-increasing, meaning

$$\begin{aligned} X_1(t)\ge X_1(0) \text {exp}\bigg ((-x_2 X_2(0) -x_3 X_3(0))\frac{t}{2}\bigg ). \end{aligned}$$

The same reasoning applied to \(X_2(t)\) and \(X_3(t)\) ensure the solutions of (11) will never become negative.

Moreover, consider the first two equations of (11); they can be bounded as follows:

$$\begin{aligned} {\left\{ \begin{array}{ll} X_1'&{}= \frac{1}{2}(-x_2 X_2 -x_3 X_3)X_1\le -\frac{1}{2}x_2 X_1X_2,\\ X_2'&{}= \frac{1}{2}(-x_1 X_1 -x_3 X_3)X_2\le -\frac{1}{2}x_1 X_1X_2. \end{array}\right. } \end{aligned}$$

This means that the trajectories of \(X_1(t)\) and \(X_2(t)\) are bounded from above by the respective solution of system

$$\begin{aligned} {\left\{ \begin{array}{ll} Y_1'&{}= -\frac{1}{2}x_2 Y_1Y_2,\\ Y_2'&{}= -\frac{1}{2}x_1 Y_1Y_2, \end{array}\right. } \end{aligned}$$

which is the classical 2D Lanchester model (2) with parameters \(x_1/2\), \(x_2/2\). From our assumptions, we know that \(Y_1(t)\rightarrow Y_1(0)-x_2 Y_2(0)/x_1>0\), while \(Y_2(t) \rightarrow 0\) as \(t\rightarrow +\infty\), implying that \(X_2(t) \rightarrow 0\), as well.

The same reasoning applied to the couple \(X_1\) and \(X_3\) shows that \(X_3(t) \rightarrow 0\) as \(t\rightarrow +\infty\).

We conclude by remarking that

$$\begin{aligned} X_1(t)=X_1(0)\text {exp}\bigg (-\dfrac{x_2}{2}\int _0^t X_2(s)\text {d}s -\dfrac{x_3}{2}\int _0^t X_3(s)\text {d}s\bigg ). \end{aligned}$$

Recall (6): since both integrands in the exponential are bounded from above by exponentially vanishing functions, then

$$\begin{aligned} \lim _{t\rightarrow +\infty } X_1(t)=X_1(0)\text {exp}\bigg (-\dfrac{x_2}{2}\int _0^{+\infty } X_2(s)\text {d}s -\dfrac{x_3}{2}\int _0^{+\infty } X_3(s)\text {d}s\bigg )>0, \end{aligned}$$
(12)

which concludes the proof. \(\square\)

Remark 3.1

For our purpose, it is sufficient to know that the limit (12) is strictly positive. However, one might obtain a lower bound for such a limit by bounding \(X_2\) and \(X_3\) from above with the explicit solution of the 2D Lanchester given in (6).

The considerations above, and the proof of Proposition 3.2, allow us to formulate the following, more general result.

Theorem 3.3

Assume that the matrix \({\tilde{C}}\) is given as in the all-fight-all scenario, i.e. every entry is 1 except the entries in the diagonal, which are 0. Then, the winner is the \(i_W\)-th army, where \(i_W=\text {argmax } x_i X_i(0)\).

3.2 Army-of-one scenario

The army-of-one scenario is, in a way, opposite to the all-fight-all: \(n-1\) armies attack the remaining one, which we will assume, without loss of generality, to always be \(X_1\). It is natural to interpret \(X_2\), \(X_3\), \(\dots\), \(X_n\) as allies, since they only focus on a common enemy. A similar scenario is also studied by Lin and MacKay (2014). Hence, there are two possible outcomes: either \(X_1\) survives, killing all the other armies, or \(X_1\) eventually perishes, and the remaining armies are victorious.

The \(n\times n\) matrices \({\tilde{C}}\) and C in this case are respectively

$$\begin{aligned} {\tilde{C}}=\left( \begin{array}{c|ccc} 0 &{} 1 &{} \dots &{} 1 \\ \hline 1 &{}\\ \vdots &{} &{} {{\textbf{0}}}&{} \\ 1 &{} &{} &{} \\ \end{array} \right) \ \ \text {and}\ \ C=\left( \begin{array}{c|ccc} 0 &{} 1 &{} \dots &{} 1 \\ \hline \frac{1}{n-1} &{}\\ \vdots &{} &{}{\textbf{0}} &{} \\ \frac{1}{n-1} &{} &{} &{} \\ \end{array} \right) . \end{aligned}$$

For \(n=3\), the system of ODEs is the following:

$$\begin{aligned} {\left\{ \begin{array}{ll} X_1'&{}= (-x_2 X_2 -x_3 X_3)X_1,\\ X_2'&{}= \frac{1}{2}(-x_1 X_1)X_2,\\ X_3'&{}= \frac{1}{2}(-x_1 X_1)X_3. \end{array}\right. } \end{aligned}$$
(13)
Fig. 3
figure 3

Evolution in time of the solutions of (13), exhibiting a \(X_1\) surviving and \(X_2\), \(X_3\) vanishing, as \(K>0\), and b \(X_1\) vanishing and \(X_2\), \(X_3\) surviving, as \(K>0\). Notice how \(X_3\) in figure b does not vanish, even though it starts at very small value, due to its belonging to the winning side

In this case, it is possible to recover a constant of motion which is closely related to the Lanchester Linear Law. In fact, by direct derivation with respect to the time variable t, it is easy to see that the following holds, along solutions of (13):

$$\begin{aligned} x_1X_1(t)-2\left( x_2X_2(t)+ x_3X_3 (t)\right) = K, \end{aligned}$$
(14)

with K a constant depending on \(x_i\) and \(X_i(0)\), \(i=1,2,3\). In particular, as we will show in a more general setting in Theorem 3.5, (14) implies that, assuming the constant is positive at time \(t=0\), \(X_1\) will be the only survivor as \(t\rightarrow +\infty\), see Fig. 3a; whereas if the constant is negative, \(X_1\) will perish and all the remaining forces will survive, see Fig. 3b.

By using an analogous argument, one can easily prove a similar result for the general n-dimensional system:

$$\begin{aligned} {\left\{ \begin{array}{ll} X_1'&{}= (-x_2 X_2 -x_3X_3 -\dots -x_n X_n)X_1,\\ X_2'&{}= \frac{1}{n-1}(-x_1 X_1)X_2,\\ X_3'&{}= \frac{1}{n-1}(-x_1 X_1)X_3,\\ \ \vdots \\ X_n'&{}= \frac{1}{n-1}(-x_1 X_1)X_n. \end{array}\right. } \end{aligned}$$
(15)

The following statement summarizes the foregoing.

Proposition 3.4

(Lanchester’s n-dimensional Linear Law). The solutions to (15) satisfy the state equation

$$\begin{aligned} x_1X_1(t)-(n-1)(x_2X_2(t)+\dots + x_nX_n(t)) = K, \quad n\ge 2, \end{aligned}$$
(16)

where K is a constant that depends only on the fighting abilities of the armies, \(x_1, \dots x_n\), and on their initial size \(X_1(0),\dots , X_n(0)\).

The claim once again follows by deriving (16) with respect to the time variable t along the solutions of system (15).

The sign of the constant of motion (16) allows us to predict the outcome of the conflict, as stated in the following theorem.

Theorem 3.5

Recall state Eq. (16) and system (15); if at \(t=0\) the sign of the constant K is positive (resp. negative), \(X_1\) will win the war (resp. \(X_1\) will perish).

Proof

From the arguments used in the proof of Proposition 3.2, we know that all \(X_i\), for \(i=1,2,\dots ,n\) will remain non-negative.

We introduce the auxiliary variable

$$\begin{aligned} Z:=x_2X_2+\dots +x_nX_n. \end{aligned}$$

Then, the ODE governing the evolution of the \((X_1,Z)\) couple is

$$\begin{aligned} {\left\{ \begin{array}{ll} X_1'&{}= -ZX_1,\\ Z' &{}= -\frac{x_1}{n-1} X_1Z. \end{array}\right. } \end{aligned}$$

Now, assume \(K>0\). Then,

$$\begin{aligned} \frac{x_1}{n-1}X_1(0)>Z(0), \end{aligned}$$

implies the asymptotic extinction of Z, which bounds each \(X_i\) to the same fate. Since Z(t) is decaying exponentially fast (recall Lemma 2.4), we know that

$$\begin{aligned} \lim _{t\rightarrow +\infty } X_1(t)=X_1(0) \text {exp}\bigg (-\int _0^{+\infty } Z(s)\text {d}s\bigg )>0. \end{aligned}$$

Now, assume \(K<0\). Then,

$$\begin{aligned} \frac{x_1}{n-1}X_1(0)<Z(0), \end{aligned}$$

implies the exponential vanishing of \(X_1\). Lastly, observing that

$$\begin{aligned} \lim _{t\rightarrow +\infty } X_i=X_i(0) \text {exp}\bigg (-\dfrac{x_1}{n-1}\int _0^{+\infty } X_1(s)\text {d}s\bigg )>0, \end{aligned}$$

since \(X_1(t)\) is decaying exponentially fast (recall again Lemma 2.4), we know that all the allied armies will survive the conflict. \(\square\)

Remark 3.2

Modelling combat with a continuous approach can, in some cases, provide counter intuitive results, such as very small armies in (15) not vanishing since they belong to the winning side of the conflict. In a discrete time model, a very small army might vanish entirely even though its allies win the war, as explained, for instance, by Fox (2010).

Remark 3.3

We remark that in the case \(n=3\), the all-fight-all and army-of-one cases completely describe every possible scenario, up to a renumbering of the armies involved.

Remark 3.4

John Lepingwell (1987) proposed a model, in which one of the two armies is split in two different units. The model is similar to the aimed version of (13). However, his model takes into account the proportional distribution of the forces of the single army between the two units of the other two forces (i.e. the single army distributes the attacks in proportion to the number or survivors of the other armies). Even though the idea is interesting, this model is affected by the same problems as the original one, namely solutions eventually assuming negative values.

3.3 Four armies example and general result

Lastly, we provide an interesting example involving four armies, which does not fall in the all-fight-all nor in the army-of-one setting. We remark that the only cases for which, thus far, we have characterized the asymptotic behaviour are the n-dimensional all-fight-all, see Theorem 3.3, and the n-dimensional army-of-one, see Theorem 3.5.

However, for \(n\ge 4\), there can be scenarios which do not fall in either of the two cases described thus far; see Fig. 4. The review of MacKay (2006) contains various examples of multilateral conflicts, in particular the case of mixed forces, in which two or more armies fight from various regions. Moreover, we refer the interested reader to Colegrave and Hyde (1993) for the analysis of the two-against-two based on an approach based on Hamiltonian dynamics, and to Kress et al. (2018) for the study of three-way and multilateral war. In the following, we describe a general scenario of multilateral conflict which, to our knowledge, has not been studied in the literature.

We introduce the following definition, for ease of notation.

Definition 3.6

We define a conflict as a triple \(\kappa =(X(0),{\tilde{C}},D)\), where X(0) is the vector \((X_1(0),\dots ,X_n(0))\) of the initial values of the involved armies, \({\tilde{C}}\) is the adjacency matrix of Definition 3.1, and D is the diagonal matrix such that \(d_{i,i}=x_i\).

A sub-conflict \(\kappa ^j\) is given by a subset \(X^j(0)\) of the armies and the corresponding matrices \({\tilde{C}}^j\) and \(D^j\), stripped of the rows and columns of the armies not considered.

As we already mentioned at the beginning of Sect. 3, to a conflict \(\kappa\) we can associate two systems of ODEs:

$$\begin{aligned} X' = -CDX, \end{aligned}$$

which correspond to the aimed fire system in Lanchester’s approach, and

$$\begin{aligned} X'_i = (-CDX)_i X_i, \end{aligned}$$

which correspond to the unaimed fire one. Here the matrix C is derived from the matrix \({\tilde{C}}\) as in Definition 3.1.

Definition 3.7

A node \(X_i\) is called bridge node if it divides the multilateral conflict \(\kappa\) in \(m\ge 2\) sub-conflicts \(\kappa ^1,\dots ,\kappa ^m\), such that

$$\begin{aligned} \bigcap _{j=1}^m X^j = X_i. \end{aligned}$$
Fig. 4
figure 4

Visual representation of the 4D cases a all-fight-all; b army-of-one; and c a case, neither all-fight-all nor army-of-one, which we discuss in this section

The next example will showcase the strategy for a general n-dimensional proof of Theorem 3.8, which we do not provide, as it would require the use of an unnecessarily cumbersome notation. We hence analyse the scenario depicted in Fig. 4c. In our modelling approach it is described by the following system of ODEs:

$$\begin{aligned} {\left\{ \begin{array}{ll} X_1'&{}= \left( -x_2 X_2 -\frac{x_3}{2}X_3 -\frac{x_4}{2} X_4 \right) X_1,\\ X_2'&{}= -\frac{x_1}{3} X_1 X_2,\\ X_3'&{}= \left( -\frac{x_1}{3} X_1-\frac{x_4}{2} X_4 \right) X_3,\\ X_4'&{}= \left( -\frac{x_1}{3} X_1 -\frac{x_3}{2}X_3 \right) X_4. \end{array}\right. } \end{aligned}$$
(17)

From Definition 3.7 it is easy to see that \(X_1\) is a bridge node for system (17). Intuitively, one would be tempted to “split” the dynamics into two all-fight-all subsystems (namely, \(X_1\) vs. \(X_2\) and \(X_1\) vs. \(X_3\) vs. \(X_4\)), and check the asymptotic behaviour for both of them with Theorem 3.3. The forces of \(X_1\) would then be split in 3, one per each army it is facing.

This would mean finding the maxima in two sets: \(\{x_1X_1(0)/3, x_2X_2(0)\}\) and \(\{x_1X_1(0)/3, x_3X_3(0)/2, x_4X_4(0)/2\}\), as highlighted in the titles of Fig. 5.

There are only two clear-cut cases: the first is when the maximum of both sets is \(x_1X_1(0)/3\), as \(X_1\) will result victorious against all its 3 opponents; see Fig. 5a. The second is when the maximum of each set is \(x_2X_2(0)\) and one among \(x_3X_3(0)/4\) or \(x_4X_4(0)/4\); see Fig. 5b.

The remaining cases are not so simple to interpret at glance, but they can be described analytically. For example if \(X_1\) would win against \(X_3\) and \(X_4\), but \(X_2\) would win against \(X_1\), then one should expect \(X_2\) as the only winner. However, our simulations show that \(X_1\) will eventually perish and the winner among \(X_3\) and \(X_4\) will be determined by their respective strength; see Fig. 5c. Vice versa, if \(X_1\) would win against \(X_2\), but not in the all-fight-all with \(X_3\) and \(X_4\), one would expect to see one between \(X_3\) and \(X_4\) as the only winner, with all the other armies reaching 0. However, in our simulations \(X_1\) perishes (as expected) and \(X_2\) survives; see Fig. 5d.

Assume, for example, that \(X_1\) would lose in the \(X_1\) vs. \(X_2\) sub-conflict, i.e. \(x_1 X_1(0)/3< x_2 X_2(0)\), but would otherwise win in the \(X_1\) vs. \(X_3\) vs. \(X_4\) one, i.e. \(x_1 X_1(0)/3> x_3 X_3(0)/2, x_4 X_4(0)/2\). Then we can focus only on these two variables and bound their respective ODE from above by

$$\begin{aligned} {\left\{ \begin{array}{ll} X_1'&{}= (-x_2 X_2 -\frac{x_3}{2}X_3 -\frac{x_4}{2} X_4)X_1 \le -x_2 X_2 X_1,\\ X_2'&{}= -\frac{x_1}{3} X_1 X_2. \end{array}\right. } \end{aligned}$$

Recalling (6), we know that \(X_1(t)\) is bounded from above by an exponentially vanishing function. Similarly, we can bound the ODEs governing \(X_3\) and \(X_4\) by

$$\begin{aligned} {\left\{ \begin{array}{ll} X_3'&{}= \left( -\frac{x_1}{3} X_1-\frac{x_4}{2} X_4 \right) X_3 \le -\frac{x_4}{2} X_4 X_3,\\ X_4'&{}= \left( -\frac{x_1}{3} X_1 -\frac{x_3}{2}X_3 \right) X_4\le -\frac{x_3}{2}X_3 X_4. \end{array}\right. } \end{aligned}$$

Assume now, without loss of generality, that \(x_3X_3(0)>x_4X_4(0)\). Then, again by (6), we know that \(X_4(t)\) is bounded from above by an exponentially vanishing function. Combining these remarks, we obtain

$$\begin{aligned} \lim _{t\rightarrow +\infty } X_3(t)= X_3(0) \text {exp}\bigg (-\dfrac{x_1}{3}\int _0^{+\infty } X_1(s)\text {d}s -\dfrac{x_4}{2}\int _0^{+\infty } X_4(s)\text {d}s \bigg )>0. \end{aligned}$$

The case in which \(X_1\) would win against \(X_2\), but would otherwise lose in the \(X_1\) vs. \(X_3\) vs. \(X_4\) subconflict, can be analysed similarly.

Fig. 5
figure 5

Evolution in time of the solutions of (17), in four different scenarios. As highlighted by the title, if we looked at the two conflicts separately, we would expect a \(X_1\) to win both; b \(X_2\) and \(X_3\) to win; c \(X_2\) to win against \(X_1\), and \(X_1\) to win against \(X_3\), \(X_4\); and d \(X_1\) to win against \(X_2\), and \(X_4\) to win against \(X_1\), \(X_3\). Notice that only the first two predictions are satisfied

These considerations highlight the contribution of bridge nodes, such as \(X_1\) in our example, which play a crucial role in higher dimensional conflicts. Indeed, if they (asymptotically) vanish, the remaining conflicts become disconnected, and multiple sub-conflict proceed to the respective asymptotic outcome. We summarize the content of this section in the following statement.

Theorem 3.8

Assume \(X_i\) is a bridge node for a conflict \(\kappa\). Assume there exists a sub-conflict in which \(X_i\) (asymptotically) loses. Then, \(X_i\) loses in all other sub-conflicts as well. Moreover, each sub-conflict \(\kappa ^j\) in which \(X_i\) is involved evolves qualitatively as \(\kappa ^j_i\), which is the sub-conflict of \(\kappa ^j\) without \(X_i\). Hence, the winner of each of those sub-conflict is determined comparing the initial force of the remaining armies, i.e. by analyzing the sub-conflict \(\kappa ^j_i\).

We omit the general proof of Theorem 3.8, since it would require a cumbersome notation and would follow the same steps of the discussion which precedes it.

Remark 3.5

The 4D case actually highlights a limitation of our model. Let us consider the example depicted in Fig. 5c: even though \(X_1\) quickly approaches 0, \(X_3\) and \(X_4\) keep attaching each other with only half of their respective forces. This can be interpreted as an initial distributions of forces on multiple fronts, which do not get rearranged even if a front becomes quiet, having defeated the corresponding opponent. However, if we were to replace the homogeneous distribution of forces with a distribution proportional to the remaining enemies in each army (by using the same idea described in Remark 3.4, for example), then we would get the problem we already discussed of some armies becoming negative.

4 Conclusion

Lanchester’s first models have paved the way to more advanced modelling efforts. This paper provides a first step towards a general n-dimensional version of the model, which we believe grants more in-depth explorations. Our formalization may sometimes provide counter intuitive results, which could hopefully be avoided with a discretization of the problem. In this section, we propose various outlooks, which represent the first step of a wider, mathematically consistent body of research.

The research carried out for this paper opened many interesting directions for extending Lanchester’s model of warfare.

First and foremost, the general analysis of system (7), in which we force all the variables to remain non-negative, would be extremely interesting. Qualitatively, one would do so by multiplying the ODE describing the evolution of \(X_i\) by the characteristic function of the set \(\{ X_i \ge 0 \}\). Then, the system would decrease its dimensionality each time a variable reaches zero, meaning the corresponding army loses all its strength. This goes beyond the scope of this paper; however, the authors plan on complementing the results presented here with a similar analysis of (7) in a future work.

From a mathematical point of view, the introduction of the randomness in the model (in terms of stochastic processes) is certainly worthwhile. Some advancements in this direction have already been made, for instance, by Kingman (2002), Osborne (2003). Moreover, the link with game theory, as well as decision theory, is clear and it is well exposed by Chen et al. (2011), Jin-Jiang et al. (2011), and Chen et al. (2012). This perspective connects in particular to the division of forces in a multilateral conflict, which in this first paper we assumed to be homogeneous. In a general scenario, which is the best, possibly heterogeneous and time-varying, division of forces for a given army?

As explained in Sect. 3, a discrete modeling constitutes a stimulating challenge, which may lead to a more realistic idealization of the warfare problem.

Finally, we remark the possibility of developing a new unconventional approach based on graph theory: this view could lay out a more in-depth exploration of higher dimensional cases, depending in particular on the adjacency matrix.

For all this reasons, we believe Lanchester’s type models pose an attractive topic of research in the context of applied mathematics, and the study of the issues arising from their analysis can involve many other areas of mathematics, making this subject still fascinating.