Numerical approximation of Nash equilibria for a class of non-cooperative differential games

In this paper we propose a numerical method to obtain an approximation of Nash equilibria for multi-player non-cooperative games with a special structure. We consider the infinite horizon problem in a case which leads to a system of Hamilton-Jacobi equations. The numerical method is based on the Dynamic Programming Principle for every equation and on a global fixed point iteration. We present the numerical solutions of some two-player games in one and two dimensions. The paper has an experimental nature, but some features and properties of the approximation scheme are discussed.

💡 Research Summary

The paper addresses the problem of computing Nash equilibria in non‑cooperative differential games with multiple players, focusing on the infinite‑horizon, exponentially discounted setting. Under standard regularity assumptions, the value functions of such games satisfy a coupled system of first‑order Hamilton–Jacobi (HJ) equations: for each player i, λ_i u_i(x) = H_i(x,∇u_1,…,∇u_m), where the Hamiltonian H_i depends on the optimal feedback controls of all players, which in turn are functions of the gradients of the value functions. This system is generally ill‑posed, especially in dimensions larger than one, and analytical results are scarce.

To obtain numerical approximations, the authors extend semi‑Lagrangian dynamic‑programming schemes—originally developed for optimal control and zero‑sum differential games—to the coupled HJ system arising from Nash equilibria. They restrict the exposition to two‑player games with scalar controls and consider dynamics in one or two spatial dimensions. The discount rates are taken equal to one for simplicity.

The discretisation proceeds as follows. A fictitious time step h>0 is introduced, leading to a semi‑discrete Bellman‑type recursion:
u_i(x) = min_{a_i∈A_i} { (1/(1+h)) u_i(x + h f(x,a_i,a_{−i}^)) + (h/(1+h)) ψ_i(x,a_i,a_{−i}^) }.
Because the forward point x + h f(·) does not generally lie on the spatial grid, a linear interpolation operator Λ(a) is constructed; it maps the vector of nodal values to the interpolated values at the off‑grid points. The spatial domain Ω is covered by a uniform grid G of N nodes, and each admissible control set A_i is replaced by a finite set A_i^# (e.g., a uniform sampling).

The fully discrete system can be written compactly as
U_i = min_{a_i∈A_i^#} { (1/(1+h)) Λ(a_i,a_{−i}^) U_i + (h/(1+h)) Ψ_i(a_i,a_{−i}^) },
where U_i∈ℝ^N collects the nodal values of u_i and Ψ_i contains the discretised running costs. The two equations are coupled through the optimal control pair a^* = (a_1^,a_2^), which itself depends on the current approximation of the value functions.

The authors propose a global fixed‑point iteration: starting from an initial guess U^{(0)}, for each grid node they exhaustively search the finite control product A_1^# × A_2^# to find the pair a^* that minimises the right‑hand side of the Bellman recursion for both players simultaneously. Once a^* is identified, the update (13) is applied to obtain U^{(k+1)}. The iteration stops when the sup‑norm change falls below prescribed tolerances ε_1, ε_2.

A key theoretical observation is that, unlike the single‑player or zero‑sum min‑max case where the fixed‑point operator is a strict contraction with factor 1/(1+h), the Nash case lacks a global contraction proof because the optimal control mapping a^(U) is piecewise constant and may jump when U varies. Nevertheless, on regions where a^ remains unchanged, the Jacobian of the operator satisfies ‖J_F‖_∞ = 1/(1+h) < 1, so the map is “piecewise contractive”. The paper defines this notion formally and suggests checking the Jacobian numerically to assess local convergence.

Numerical experiments illustrate the method. In a one‑dimensional test, the dynamics are f(x,a_1,a_2)=a_1+a_2 and the running costs are ψ_i(x,a_1,a_2)=h_i(x)+a_i^2 with h_i≡0. The associated HJ system admits three formal solutions; only one satisfies the admissibility conditions (absolute continuity, sub‑linear growth, and a one‑sided derivative condition) and is recovered by the scheme on a grid of 51 points over