Development of Lanchester-Type Spatial Models with Obtaining Localized Solutions for the Interaction of Two Groups

Nikita D. Borisov

Phystech School of Applied Mathematics and Informatics, Moscow Institute of Physics and Technology (National Research University), Dolgoprudny, Russia.

DOI: 10.4236/jamp.2025.137132                

Abstract

This paper presents an advanced mathematical framework for modeling combat dynamics between two opposing forces using nonlinear reaction-diffusion equations. Extending classical Lanchester models, we incorporate spatially dependent diffusion coefficients to capture modern battlefield complexities. A robust numerical scheme based on the Alternating Direction Implicit (ADI) method is developed, ensuring stability and second-order accuracy in spatiotemporal discretization. The model integrates logistic growth, combat attrition, and tactical diffusion processes, validated through analytical benchmarks. Simulations reveal intricate pattern formation, transient dynamics, and boundary adherence, demonstrating applicability to military strategy optimization. The complete formulation and numerical implementation are thoroughly discussed, providing insights into nonlinear system behavior under varying tactical conditions.

Keywords

Combat Dynamics, Lanchester Models, Nonlinear Systems, Numerical Stability, Alternating Direction Implicit Method

1. Introduction

This paper develops a mathematical model for combat operations between two military groups using nonlinear diffusion-reaction equations. The system extends classical Lanchester-type models [1]-[3] by incorporating spatial dynamics through nonlinear diffusion terms dependent on troop concentrations:

$\frac{\partial U}{\partial t}=\text{Combat terms}+\nabla \cdot \left(k\left(U,V\right)\nabla U\right)$ (1)

where diffusion coefficients $k\left(U,V\right)$ depend nonlinearly on troop density $U\left(x,y,t\right)$ and $V\left(x,y,t\right)$ . Key challenges include handling nonlinear diffusion and ensuring numerical stability—issues resolved through our ADI-based approach [4] [5].

The model captures essential battlefield phenomena: troop movements, engagement dynamics, and force depletion [6] [7]. Key challenges include handling the nonlinear diffusion coefficients and developing stable numerical solutions.

We present an efficient numerical solver for reaction-diffusion systems using the Alternating Direction Implicit (ADI) method, achieving second-order accuracy in both space and time [8]. This approach maintains computational tractability while capturing complex nonlinear interactions.

2. Mathematical Model

The combat dynamics between two military groups is described by a system of nonlinear diffusion-type partial differential equations [9] [10]:

$\{\begin{array}{l}\frac{\partial U}{\partial t}=\alpha U\left(a_1-b_{1}U-c_{1}V\right)-d_{1}UV+{\mu }_1\nabla \cdot \left(k_1\left(U,V\right)\nabla U\right),\hfill \\ \frac{\partial V}{\partial t}=\beta V\left(a_2-b_{2}V-c_{2}U\right)-d_{2}VU+{\mu }_2\nabla \cdot \left(k_2\left(U,V\right)\nabla V\right),\hfill \end{array}$ (2)

where $U=U\left(x,y,t\right)$ and $V=V\left(x,y,t\right)$ represent the troop density of the first and second groups at point $\left(x,y\right)$ and time $t$ . The system is solved in the spatial-temporal domain:

$D=D_{xy}\times \left(0<t<T_{\max }\right)=\left\{\left(X_0<x<X_1\right)\times \left(Y_0<y<Y_1\right)\right\}\times \left(0<t<T_{\max }\right),$ (3)

which can be simplified to a unit square without loss of generality.

The diffusion operator has the standard divergence form:

$\nabla \cdot \left(k_s\nabla F\right)=\frac{\partial }{\partial x}\left(k_s\frac{\partial F}{\partial x}\right)+\frac{\partial }{\partial y}\left(k_s\frac{\partial F}{\partial y}\right),s=1,2.$ (4)

The nonlinearity of the system arises from the dependence of diffusion coefficients $k_1\left(U,V\right)$ and $k_2\left(U,V\right)$ on the unknown functions, which may be specified by power-law relations [11]:

$k_1\left(U,V\right)=k_1^{\left(0\right)}{U}^{{\sigma }_1}{V}^{{\delta }_1},k_2\left(U,V\right)=k_2^{\left(0\right)}{U}^{{\sigma }_2}{V}^{{\delta }_2},$ (5)

where $k_1^{\left(0\right)},k_2^{\left(0\right)},{\sigma }_{1,2},{\delta }_{1,2}$ are appropriately chosen constants.

The system is complemented by initial conditions:

$\{\begin{array}{l}U\left(x,y,0\right)=U_0\left(x,y\right),\hfill \\ V\left(x,y,0\right)=V_0\left(x,y\right),\hfill \end{array}\left(x,y\right)\in \left(0,1\right)\times \left(0,1\right),$ (6)

We consider a reaction-diffusion system of equations in a square domain $\text{Ω}=\left[0,L\right]\times \left[0,L\right]$ with boundary $\partial \Omega$ . For the components $U\left(x,y,t\right)$ and $V\left(x,y,t\right)$ , homogeneous Dirichlet boundary conditions are imposed, where zero values are maintained on all domain boundaries [12]:

$\{\begin{array}{ll}U\left(0,y,t\right)=U\left(L,y,t\right)=0,\hfill & y\in \left[0,L\right],t\ge 0\hfill \\ U\left(x,0,t\right)=U\left(x,L,t\right)=0,\hfill & x\in \left[0,L\right],t\ge 0\hfill \\ V\left(0,y,t\right)=V\left(L,y,t\right)=0,\hfill & y\in \left[0,L\right],t\ge 0\hfill \\ V\left(x,0,t\right)=V\left(x,L,t\right)=0,\hfill & x\in \left[0,L\right],t\ge 0\hfill \end{array}$ (7)

3. Numerical Methodology

The ADI method decomposes the spatial operator into directional components [5]:

$ℒ={ℒ}_x+{ℒ}_y$ (8)

where the directional operators for group U are defined as:

${ℒ}_{x}U=\frac{\partial }{\partial x}\left(k_1\left(U,V\right)\frac{\partial U}{\partial x}\right)\approx \frac{k_{i+1/2,j}\left(U_{i+1,j}^n-U_{i,j}^n\right)-k_{i-1/2,j}\left(U_{i,j}^n-U_{i-1,j}^n\right)}{h_x^2}$ (9)

${ℒ}_{y}U=\frac{\partial }{\partial y}\left(k_1\left(U,V\right)\frac{\partial U}{\partial y}\right)\approx \frac{k_{i,j+1/2}\left(U_{i,j+1}^n-U_{i,j}^n\right)-k_{i,j-1/2}\left(U_{i,j}^n-U_{i,j-1}^n\right)}{h_y^2}$ (10)

The diffusion coefficients at half-points are averaged using:

$k_{i+1/2,j}=\frac{1}{2}\left(k_1\left(U_{i+1,j}^n,V_{i+1,j}^n\right)+k_1\left(U_{i,j}^n,V_{i,j}^n\right)\right)$ (11)

3.1. Temporal Discretization

The ADI scheme employs the following fractional steps [8]:

$\frac{U^{n+1/2}-U^n}{\tau /2}={ℒ}_x{U}^{n+1/2}+{ℒ}_y{U}^n+ℛ\left(U^n,V^n\right)$ (12)

$\frac{U^{n+1}-U^{n+1/2}}{\tau /2}={ℒ}_x{U}^{n+1/2}+{ℒ}_y{U}^{n+1}+ℛ\left(U^{n+1/2},V^n\right)$ (13)

where $ℛ$ contains all reaction terms. The matrix form for Step 12 becomes:

$\left(I-\frac{\tau }{2}{A}_x\right)U^{n+1/2}=\left(I+\frac{\tau }{2}{A}_y\right)U^n+\frac{\tau }{2}{ℛ}^n$ (14)

with $A_x$ , $A_y$ being tridiagonal matrices representing ${ℒ}_x$ , ${ℒ}_y$ respectively.

3.2. Nonlinear Treatment

For nonlinear terms, we employ Picard iteration [4]:

1: Initialize $U^{\left(0\right)}=U^n$ ;

2: for $k=1,2,\cdots$ until $\Vert U^{\left(k\right)}-U^{\left(k-1\right)}\Vert <{10}^{-6}$ do;

3: Update diffusion coefficients $k_1^{\left(k\right)}=k_1\left(U^{\left(k-1\right)},V^n\right)$ ;

4: Solve tridiagonal systems for $U^{\left(k\right)}$ ;

5: end for.

3.3. Stability Analysis

The scheme is unconditionally stable when linearized. For constant $k_1$ , the amplification factor G satisfies:

$\left|G\right|={\left|1-\frac{\tau }{2}\left({\lambda }_x+{\lambda }_y\right)+\frac{{\tau }^2}{4}{\lambda }_x{\lambda }_y\right|}^{-1}\le 1$ (15)

where ${\lambda }_x$ , ${\lambda }_y$ are eigenvalues of $A_x$ , $A_y$ respectively. The stability region is shown in Figure 1. The local truncation error is [8]:

${ℰ}_{local}=\frac{{\tau }^2}{4}\left({ℒ}_x{ℒ}_y-{ℒ}_y{ℒ}_x\right)U+\underset{\_}{O}\left({\tau }^3+h^2\right)$ (16)

Figure 1. Stability region of the ADI scheme (shaded) compared to explicit Euler (dashed line).

4. Model Verification

4.1. Analytical Benchmarks

To validate our numerical implementation, we conduct verification against two classes of exact solutions:

4.1.1. Traveling Wave Solution (Fisher-KPP Type)

For constant diffusion $k_1={\mu }_1$ and $V\equiv 0$ , the system reduces to [11]:

$\frac{\partial U}{\partial t}=\underset{\text{Logistic growth}}{\underbrace{\alpha U\left(1-U\right)}}+\underset{\text{Diffusion}}{\underbrace{{\mu }_1{\nabla }^{2}U}}$ (17)

The exact traveling wave solution $U\left(x,t\right)=\varphi \left(x-ct\right)$ satisfies:

$-{\mu }_1{\varphi }^{″}+c{\varphi }^{\prime }+\alpha \varphi \left(1-\varphi \right)=0$ (18)

For wave speed $c=\frac{5}{\sqrt{6}}$ , the solution has implicit form:

${\displaystyle {\int }_{{\varphi }_0}^{\varphi }\frac{\text{d}\psi }{\psi \sqrt{1-\psi +\frac{6{\mu }_1}{\alpha c^2}\ln \psi }}}=\xi -{\xi }_0$ (19)

Numerical implementation uses the asymptotic approximation:

$\varphi \left(\xi \right)\approx {\left(1+{\text{e}}^{\xi /\sqrt{6{\mu }_1/\alpha }}\right)}^{-1},\xi =x-ct$ (20)

4.1.2. Self-Similar Blow-Up Solution

For power-law diffusion $k_1=U^{\sigma }$ with $\sigma =2$ and cubic reaction $-U^3$ , the solution exhibits finite-time blow-up [9] [13]:

$U\left(r,t\right)={\left(T-t\right)}^{-1/2}\text{exp}\left[-\frac{r^2}{8\left(T-t\right)\text{ln}\left(1/\left(T-t\right)\right)}\left(1+\underset{\_}{O}\left(\frac{1}{\text{ln}\left(1/\left(T-t\right)\right)}\right)\right)\right]$ (21)

The scaling law for $L^{\infty }$ norm is:

${\Vert U\left(\cdot ,t\right)\Vert }_{L^{\infty }}~{\left(T-t\right)}^{-1/2},t\to T^{-}$ (22)

Results are shown in Figure 2.

Figure 2. Analytical benchmark verification.

4.2. Grid Convergence Study

Temporal ($h=0.01,0.005,0.0025$ ) and spatial ($\tau =0.1,0.05,0.025$ ) refinements confirm second-order convergence [8]. The Richardson extrapolation error estimate satisfies (Figure 3):

${\Vert U_{h,\tau }-U_{\text{exact}}\Vert }_2\propto h^2+{\tau }^2.$ (23)

Figure 3. Verification results: (a) Analytical vs numerical wavefronts; (b) L2-error vs grid resolution.

5. Results of Numerical Modeling

The numerical simulations of the reaction-diffusion system reveal complex spatiotemporal dynamics between the two interacting components $U\left(x,y,t\right)$ and $V\left(x,y,t\right)$ . The implemented ADI (Alternating Direction Implicit) method with $\underset{\_}{O}\left({\tau }^2+h^2\right)$ accuracy successfully captures the evolution of the system, demonstrating stable behavior even for the strongly nonlinear case [5] [8].

Key observations from the simulations include:

• Formation of distinct spatial patterns in U with amplitude modulation reaching $\left|U_{\max }\right|\approx 2.0$ , while V shows more subdued dynamics with $\left|V_{\max }\right|\approx 0.5$ [12];

• Emergence of transient structures during the initial phase ($t<0.2T_{\max }$ ) followed by stabilization [10];

• Effective temporal modulation of reaction coefficients through the $\alpha \left(1-\sin \left(\pi t\right)\right)\left(t/T_{\max }+\beta \right)$ term [14];

• Preservation of boundary conditions (${U|}_{\text{Γ}}={V|}_{\text{Γ}}=0$ ) throughout the simulation [12].

Figure 4. Evolution of $U\left(x,y,t\right)$ showing pattern formation stages.

The visualizations in Figure 4 clearly show the competition between diffusion processes (controlled by ${\mu }_1=0.1$ , ${\mu }_2=0.2$ ) and nonlinear reaction terms, resulting in complex but stable pattern formation [9]. The numerical scheme proves robust for this class of problems, handling both the stiff diffusion terms and nonlinear couplings effectively [8].

6. Conclusions and Suggestions

The numerical implementation of the coupled reaction-diffusion system using the ADI method with temporal parameter modulation has yielded several significant computational insights. The solver effectively captures the evolving spatiotemporal patterns of the interacting components $U\left(x,y,t\right)$ and $V\left(x,y,t\right)$ , with $U$ showing amplified patterns (peak magnitude ≈2.0) compared to $V$ (peak magnitude ≈1.0), consistent with the designed initial conditions and parameter choices (${\alpha }_1=1.0$ , ${\beta }_1=-0.3$ ) [12].

Figure 5 demonstrates the localized spatial patterns of force concentrations, where subfigure (a) shows the primary group’s density distribution with distinct peak formations, while subfigure (b) reveals the opposing group’s more dispersed configuration. These visualizations confirm the model’s capability to capture both concentrated and diffuse combat scenarios.

The power-law form of diffusion coefficients in Equation (5) is justified by the following military considerations:

Nonlinear Mobility Effects

• ${\sigma }_s>0$ : Enhanced troop mobility at high friendly density (e.g., concentrated logistics support).

• ${\delta }_s<0$ : Restricted movement due to enemy presence (area denial effects).

• ${\sigma }_s={\delta }_s=0$ : Reduces to classical Lanchester models.

Figure 5. Spatial force distribution showing localization effects.

Tactical Scenarios Modeled

Blitzkrieg:${\sigma }_1\gg 1,{\delta }_1\approx 0$ (Fast advance at high concentration);

Defensive:${\sigma }_1\approx 0,{\delta }_1<0$ (Enemy slows movement);

Guerilla:${\sigma }_1<0$ (Dispersion at high density).

Mathematical Advantages

• Preserves scaling invariance for analytical solutions;

• Accommodates bifurcation analysis ($\sigma ,\delta$ as critical parameters);

• Captures threshold effects (e.g., minimal $U$ for offensive).

This formulation provides a parsimonious yet flexible representation of modern combat dynamics where troop mobility depends nonlinearly on both friendly and enemy force distributions.

The temporal modulation scheme $\alpha \left(1-\sin \left(\pi t\right)\right)\left(t/T_{\max }+\beta \right)$ successfully introduces controlled non-stationary behavior while preserving numerical stability, with default modulation parameters ${\alpha }_{\text{mod}}=0.5$ and ${\beta }_{\text{mod}}=0.2$ [14]. The method maintains strict Dirichlet boundary conditions ($U=V=0$ at boundaries) throughout simulations and handles the nonlinear coupling terms ($a_{12}UV$ and $a_{21}UV$ ) robustly, with diffusion coefficients ${\mu }_1=0.05$ and ${\mu }_2=0.1$ producing physically meaningful gradient evolution [9].

The implementation shows particular promise for modeling systems with:

• Effective handling of asymmetric initial conditions (2:1 amplitude ratio between U and V) [12].

• Robust temporal parameter modulation without stability constraints [14].

• Accurate treatment of nonlinear reaction terms (quadratic self-interaction and cross-component coupling) [10].

• Efficient solution of implicit systems through sparse matrix algebra [8].

The simulated spatiotemporal dynamics of troop concentrations $U\left(x,y,t\right)$ and $V\left(x,y,t\right)$ exhibit mathematically predictable pattern formation with direct military applications. The characteristic length scale of emerging combat patterns is determined by the critical engagement distance:

${\lambda }_c\approx \sqrt{\frac{{\mu }_1}{\alpha }}$

where ${\mu }_1$ represents force mobility and $\alpha$ the engagement intensity. This fundamental scale suggests that in typical combat scenarios with ${\mu }_1\approx 0.1$ km²/hr and $\alpha \approx 0.05$ hr⁻¹, the natural operational separation between units should not exceed ${\lambda }_c\approx 1.4$ km to maintain effective mutual support. Regions exhibiting high Laplacian values ${\nabla }^{2}U$ correspond to optimal kill zones for indirect fire systems, while gradient fields $\nabla V$ reveal probable enemy lines of advance.

Future extensions could explore adaptive time-stepping and three-dimensional implementations while maintaining the current scheme’s $\underset{\_}{O}\left({\tau }^2+h^2\right)$ accuracy and stability properties [15].

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Atkinson, M.P., Kress, M. and MacKay, N.J. (2020) Targeting, Deployment, and Loss-Tolerance in Lanchester Engagements. Operations Research, 69, 71-81. [Google Scholar] [CrossRef]
[2]	Kostić, M. and Jovanović, A. (2020) Lanchester’s Differential Equations as Operational Command Decision Making Tools. Military Technical Courier, 68, 523-540.
[3]	McCartney, M. (2023) Battling with Lanchester’s Equations in the Classroom. International Journal of Mathematical Education in Science and Technology, 54, 451-461. [Google Scholar] [CrossRef]
[4]	Fryazinov, I.V. and Bakirova, M.I. (1972) On Economical Difference Schemes for Heat Equation in Curvilinear Coordinates. USSR Computational Mathematics and Mathematical Physics, 12, 87-100.
[5]	Marchuk, G.I. (1988) Splitting Methods, Nauka. (In Russian)
[6]	Korepanov, V.O., Chkhartishvili, A.G. and Shumov, V.V. (2023) Basic Models of Combat Operations. Urban and Building Science, 103, 40-77. (In Russian)
[7]	Liu, Y., Zhang, X., Du, H., Wang, G. and Zeng, D. (2021) Construction and Simulation of Lanchester Battle Equations Based on Space-Based Information Support. Journal of Physics: Conference Series, 2384, Article 012011. [Google Scholar] [CrossRef]
[8]	Samarskii, A.A. (2001) The Theory of Difference Schemes. Marcel Dekker.
[9]	Kurdyumov, S.P., Kurkina, E.S. and Tel’kovskaya, O.V. (1989) Blow-Up Regimes in Two-Component Media. Applied Mathematical Modelling, 35, 34-50 (In Russian).
[10]	Zang, W.B. (1999) Synergetic Economics: Time and Change in Nonlinear Economic Theory. Springer.
[11]	Galaktionov, V.A., Dorodnitsyn, V.A., Elenin, G.G., Kurdyumov, S.P., and Samarskii, A.A. (1986) Quasilinear Heat Equation with Source: Blow-Up, Localization and Exact Solutions. Mathematical Modeling.
[12]	Inovenkov, I.N., Nefedov, V.V. and Tikhomirov, V.V. (2022) Computer Modeling of Urban Population Dynamics. Modern Information Technologies and IT-Education, 18, 300-309.
[13]	Belavin, V.A. and Kurdyumov, S.P. (2000) Blow-Up Regimes in Demographic Systems. Computational Mathematics, 40, 227-239.
[14]	Fradkov, A. (2021) Dynamics and Stability under Iterated Sanctions and Counter-sanctions. Automation and Remote Control, 82, 401-415.
[15]	Yoo, B.J. (2020) Statistical Review and Explanation for Lanchester Model. International Journal of Applied Mathematics and Statistics, 59, 123-135.