Structure-preserving discretization of nonlinear cross-diffusion systems

• Ilaria Perugia (Universität Wien)

Room: Amphitheater, Building B The main challenges in designing numerical methods for approximating nonlinear cross-diffusion systems is that the diffusion matrix may not be symmetric or positive semidefinite, and that a maximum principle may be not available. In this talk, we present a Local Discontinuous Galerkin method for discretizing nonlinear cross-diffusion systems, which is based on the boundedness-by-entropy framework introduced by A. Jüngel in 2015. Motivated by the underlying entropy structure of the PDE system, nonlinear transformations in terms of the entropy variable allow to enforce positivity of approximate solutions. Moreover, by appropriately introducing auxiliary variables, the problem is reformulated so that nonlinearities do not appear within differential operators or interface terms, leading to nonlinear operators that can be naturally evaluated in parallel. The resulting method has the following desirable properties: i) it allows arbitrary degrees of approximation in space; ii) it preserves boundedness of the physical unknowns without requiring postprocessing or slope limiters; iii) nonlinearities do not appear explicitly within differential operators or interface terms, giving the method with a natural parallelizable structure and high efficiency; iv) it respects a discrete version of the entropy stability estimate of the continuous problem. This is a joint work with Sergio Gómez and Ansgar Jüngel.

Computable reliable bounds for Poincaré–Friedrichs constants via Čech–de-Rham complexes

• Martin W. Licht (EPFL)

Room: Amphitheater, Building B We derive computable and reliable upper bounds for Poincaré–Friedrichs constants of classical Sobolev spaces and, more generally, Sobolev de-Rham complexes. The upper bounds are in terms of local Poincaré–Friedrichs constants over subdomains and the smallest singular value of a finite-dimensional operator that is easily assembled from the geometric setting. Thus we reduce the computational effort when computing the Poincaré–Friedrichs constant of finite de-Rham complexes, and we provide computable reliable bounds even for the original Sobolev de-Rham complex. The reduction to a finite-dimensional system uses diagram chasing within a Čech–de-Rham complex. Additionally, we utilize estimates for Poincaré–Friedrichs constants over local finite element patches. Part of this is joint work with Théophile Chaumont-Frelet and Martin Vohralík.

A posteriori goal-oriented error estimators based on equilibrated flux and potential reconstructions

• Emmanuel Creusé (Université Polytechnique Hauts-de-France)

Room: Amphitheater, Building B Nowadays, many engineering problems require computing some quantities of interest, which are usually linear functionals applied to the solution of a partial differential equation. Error estimations of such functionals are called "goal-oriented" error estimations. Such estimations are based on the resolution of an adjoint problem, whose solution is used in the estimator definition, and on the use of some energy-norm error estimators. In this talk, an overview of such techniques in different contexts will be given. We will then provide an upper-bound of the error which can be totally and explicitly computed for various discretization schemes. Finally, the behaviour of such estimators on some numerical benchmarks will be investigated. Two models will be particularly considered : a reaction-diffusion problem, and an eddy-current problem, arising in the context of low-frequency electromagnetism.

Basics for polynomial interpolation on simplices

• Francesca Rapetti (Université Côte d'Azur)

Room: Amphitheater, Building B Polynomial interpolation is a key aspect in numerical analysis, used in very classical settings as for reconstructing a field from measures, computing integrals by quadratures formulas or selecting basis functions in finite element methods. We will review the roles of the three main characters featured in this action, namely, the representation of the domain by a mesh, the polynomial basis, the Vandermonde matrix as the link between the two aspects. We will present a general framework for the interpolation of differential k-forms on simplices that allows to retrieve fundamentals concepts featuring in the classical scalar case. The Lebesgue constant pops up naturally to measure the stability of the interpolation and we will see the Runge phenomenon as a counter-example of this stability. This work is partially done in collaboration with Ana Alonso Rodriguez and Ludovico Bruni Bruno.

Hybrid compatible Finite Element and Finite Volume discretization for viscous and resistive MHD

• Eric Sonnendrücker (Max Planck Institute for Plasma Physics and TU München)

Room: Amphitheater, Building B MHD simulations including small viscous and resistive effects are fundamental for simulations related to magnetic fusion. However, due to the needed long time simulations and the very different wave speeds, implicit or semi-implicit methods are unavoidable. On the other hand, div B = 0 as well as other symmetries and invariants need to be preserved by the numerical algorithm. To this aim, we developed a method based both on Finite Volumes for handling the slow convection and robust for shocks and on the other hand Finite Element Exterior Calculus which enable exact conservation of the main invariant of the system.

The Exterior Calculus Discrete De Rham complex

• Jérôme Droniou (CNRS, Université de Montpellier and Monash University)

Room: Amphitheater, Building B In its standard presentation, the de Rham complex organises the gradient, curl and divergence operator into a sequence that embeds the well-known calculus relations: the image of one operator (e.g. gradient) is included in the kernel of the following one (e.g. curl). The de Rham theorem states that the gaps between these images and kernels, embedded into the cohomology of the complex, is related to the topology of the domain. The importance of this complex and its properties in the stability analysis of models of partial differential equations (such as the Stokes/Navier-Stokes equations, magnetostatic equations, etc.) has been understood for decades. Reproducing the properties of this complex at the discrete level is essential for the design of stable schemes for these models, and is related, e.g., to the design of inf-sup stable methods for saddle point problems. In the last two decades, the Finite Element Exterior Calculus (FEEC) framework has been set up to devise versions of the de Rham complex through the exterior calculus framework, which allows to treat all operators (gradient, curl, divergence) in a unified way as exterior derivatives of differential forms of certain degrees. These discrete complexes are however restricted to particular meshes (mostly made of tetrahedra and hexahedra), which do not easily lends themselves to standard scientific calculus techniques like local mesh refinement or mesh agglomeration (appearing, e.g., in multi-grid methods). In this talk, we will present the Exterior Calculus Discrete De Rham (ECDDR) method. This is a discrete version of the de Rham complex of differential forms, that can be applied on polytopal meshes (made of generic polygons in 2D, generic polyhedra in 3D). As many polytopal methods, its design is based on adopting a higher and systematic view, which not only relaxes the conditions on the meshes, but can also lead to leaner methods than standard Finite Element methods. The design of ECDDR relies on the Stokes formula, which identifies the relevant degrees of freedom, as well as provides expressions for the discrete differential forms and potential reconstructions. We will show that the algebraic properties of the de Rham complex are preserved at the discrete level (including its cohomology), and we will explain how the tools in the ECDDR can be used to design numerical schemes.

Polytopal methods on Riemannian manifolds

• Marien-Lorenzo Hanot (The University of Edinburgh)

Room: Amphitheater, Building B Discretization methods based on differential complexes have many advantageous properties in terms of stability, framework for analysing the discrete formulation, and preservation of important quantities such as the mass, helicity, or the pressure robustness in fluid dynamics. The premise of this kind of approach appeared early on with elements based on the compatibility between the geometry and the differential operator, such as the Nédélec or the Raviart-Thomas elements relating the curl operator to the circulation along edges, and the divergence operator to the flux across surfaces. The connection between the usual differential operators (gradient, curl, and divergence) and the geometry can be better seen through framework of the exterior calculus, where these operators are unified as the exterior derivative applied to differential forms of different degrees. The associated finite element spaces also have a natural description which has been developed into the Finite Element Exterior Calculus (FEEC) framework, leading to an intrinsic definition common to each space and operator in any dimension. Several other methods that replicate the complexes structuring systems of differential equations at the discrete level were then developed to use smoother spaces, different complexes, or more general meshes such as the Discrete De Rham (DDR) method. Although the notions involving exterior calculus are manifestly independent of the underlying metric, most discrete methods must ultimately assume a trivial space for the notion of simplicial/polytopal mesh, and for the notion of polynomial. In this talk, we will present the generalisation of the Exterior Calculus Discrete De Rham (ECDDR) method to general Riamannian manifolds. This construction uses a much more lenient notion of mesh, allowing to consider manifolds described by several charts, and to use potentially any shape for the elements. In particular, it is possible to use curved elements even when working on a flat space. The basis functions are intrinsically defined element-wise. They are based on polynomial spaces of arbitrary order and are adapted to the chart and the metric. We will then present a numerical application of this method to the Maxwell equations on a surface.

Learning based reduction methods in the context of PDE constrained optimization

• Mario Ohlberger (Universität Münster)

Room: Amphitheater, Building B Model order reduction for parameterized partial differential equations is a very active research area that has seen tremendous development in recent years from both theoretical and application perspectives. A particular promising approach is the reduced basis method that relies on the approximation of the solution manifold of a parameterized system by tailored low dimensional approximation spaces that are spanned from suitably selected particular solutions, called snapshots. With speedups that can reach several orders of magnitude, reduced basis methods enable high fidelity real-time simulations for certain problem classes and dramatically reduce the computational costs in many-query applications. While the ”online efficiency” of these model reduction methods is very convincing for problems with a rapid decay of the Kolmogorov n-width, there are still major drawbacks and limitations. Most importantly, the construction of the reduced system in a so called ”offline phase” is extremely CPU-time and memory consuming for large scale systems. For practical applications, it is thus necessary to derive model reduction techniques that do not rely on a classical offline/online splitting but allow for more flexibility in the usage of computational resources. In this talk we focus on learning based reduction methods in the context of PDE constrained optimization and inverse problems and evaluate their overall efficiency. We discuss learning strategies, such as adaptive enrichment as well as a combination of reduced order models with machine learning approaches in the contest of time dependent problems. Concepts of rigorous certification and convergence will be presented, as well as numerical experiments that demonstrate the efficiency of the proposed approaches.

Global space-time low-rank methods for the time-dependent Schrödinger equations

• Virginie Ehrlacher (Ecole des Ponts ParisTech and Inria Paris Centre)

Room: Amphitheater, Building B The aim of this talk is to present novel global space-time methods for the approximation of the time-dependent Schrödinger equation, using Kato theory. The latter can be used in conjunction with low-rank tensor formats (such as Tensor Trains for instance) to derive new variational principles to compute dynamical low-rank approximations of the solution, which are different from the Dirac-Frenkel principle. One significant advantage of this new variational formulation is that the existence of a dynamical low-rank approximation for any finite-time horizon can be proved, whereas dynamical low-rank approximations constructed with the Dirac-Frenkel principle can usually be porved to exist only locally in time. Illustrative numerical results will be presented to highlight the differences between the dynamical low-rank approximations obtained with these different approaches. This is joint work with Clément Guillot and Mi-Song Dupuy.

Runge-Kutta methods are stable

• Eitan Tadmor (University of Maryland)

Room: Amphitheater, Building B The numerical solution of PDEs often ends up with a large system of ODEs, and a canonical choice for the solution of such systems of “method of lines” is the class of Runge-Kutta (RK) methods. Indeed, RK methods are used routinely for integration of large systems of ODEs encountered in various applications. But the standard stability arguments of RK method fail to cover arbitrarily large systems of ODEs. We explain the failure of different approaches, offer a new stability theory and demonstrate a few examples.

Poster session Room: Plenary room, Building A Margherita Castellano (Ecole Polytechnique): A finite volume method for Cahn--Hilliard equations with surfactants Jean Cauvin-Vila (TU Vienna): Structure-preserving finite volume approximation of cross-diffusion systems coupled via a moving interface Farah Chaaban (ENSTA Paris): A volume optimal control-based numerical method for 2D time-harmonic Maxwell's equations with sign-changing coefficients Amélie Dupouy (Inria Lille): Theoretical and numerical analysis of a diffusion problem on a moving domain Youssef Essadaoui (Université Sultan Moulay Slimane): Convergence to equilibrium for a sinc-type model surface growth model Maxime Jonval (Inria Lille, IFPEN): Parametrization and Cartesian representation techniques for robust resolution of chemical equilibria Tino Laidin (Université de Lille): Conservative polynomial approximations and applications to Fokker--Planck equations François Madiot (CEA Saclay): Criticality calculations in neutronics: model order reduction and a posteriori estimators Christina Mahmoud (Université de Montpellier): Uniformly accurate schemes for hyperbolic relaxation systems Ismail Merabet (Kasdi Merbah University Ouargla): Discontinuous finite element method for the contact problem of a linearly elastic shell Julien Moatti (TU Vienna): A finite volume scheme for Maxwell--Stefan systems using Bott--Duffin inverse Jia Jia Qian (Monash University): Discretisations of exterior calculus models from physics Marwa Salah (Université de Montpellier): A serendipity fully discrete div-div complex on polygonal meshes

Modeling Multiphase Multicomponent Porous Flows

• Noel J. Walkington (Carnegie Mellon University)

Room: Amphitheater, Building B This talk will review structural properties of the equations used to model porous flows involving multiple components undergoing phase transitions. These equations only model the gross properties of these problems since a precise description of the physical system is neither available nor computationally tractable. The saddle point structure resulting from the interaction between dissipation and free energy (or entropy) of the fluids will be highlighted. The construction of numerical schemes which are robust in the presence of degeneracy, and solution techniques which exploit the saddle point structure, will be considered.

Asymptotic analysis of a fluid-structure-porous layer coupled model near contact

• Céline Grandmont (Inria Paris Centre and Université Libre de Bruxelles)

Room: Amphitheater, Building B The numerical simulation of systems involving fluid-structure-contact interaction raises many modeling, mathematical and numerical difficulties. It is in particular crucial for numerous biomedical applications such as the simulation of cardiac valve dynamics (native or artificial) for instance. Fluid-structure interaction without contact is already challenging due to the moving geometries and the potential strong coupling between the solid and the fluid subsystems. If contact between solids is to be modeled as well, the complexity increases with additional difficulties among which: - In some configurations and with no-slip boundary conditions, fluid-structure interaction models are unable to predict contact (see, e.g., [1, 4, 5]), this is the so called no collision paradox; - The simple addition of a contact constraint (variational inequality) to a fluid-structure interaction model which allows for contact yields a mechanically inconsistent fluid-structure-contact interaction (see [6, 7]). The first difficulty can be solved by modifying the boundary and interface conditions on the contact walls (see, e.g., [8, 9, 10]). However this is not necessarily enough to obtain a mechanically consistent model. An alternate approach considered in [6, 7] to circumvent these two issues is to consider a poroelastic modeling of the fluid seepage induced by the roughness of the contacting solid. Yet, very little is known on the mathematical foundations of this approache. In this work, we analyze the ability of such a model to encompass contact. We consider a 2D fluid-ball interaction problem and we build on the approach proposed in [4, 11] to evaluate the behaviour of the drag force with respect to the gap between the solid and the wall. The asymptotics of the model with respect to the porous layer parameters are also investigated. Finally, numerical evidence of these theoretical results is provided. REFERENCES [1] Burman, Erik, et al. A mechanically consistent model for fluid-structure interactions with contact including seepage.Computer Methods in Applied Mechanics and Engineering 392 (2022): 114637. [2] Gérard-Varet, D., and Hillairet, M. Computation of the Drag Force on a Sphere Close to a Wall: The Roughness Issue ESAIM: Mathematical Modelling and Numerical Analysis 46.5 (2012): 1201-224. [3] T. Hesla, Collisions of smooth bodies in viscous fluids: A mathematical investigation, PhD thesis: University Minnesota, 2004. [4] M. Hillairet, Lack of collision between solid bodies in a 2D incompressible viscous flow, Commun. Part. Diff. Eq., vol. 32, no. 9, pp. 1345-1371, 2007. [5] M. Hillairet and T. Takahashi, Collisions in three-dimensional fluid structure interaction problems,SIAM J. Math. Anal., vol. 40, pp. 2451-2477, 2009. [6] C. Ager, B. Schott, A.-T. Vuong, A. Popp and W. Wall, A consistent approach for fluid-structure-contact interaction based on a porous flow model for rough surface contact, Internat. J. Numer. Methods Engrg., vol. 119, no. 13, pp. 1345-1378, 2019. [7] S. Frei, F. Gerosa, E. Burman and M. Fernàndez, A mechanically consistent model for fluid-structure interactions with contact including seepage, HAL preprint, 2021. [Online]. Available: https://hal.archives-0uvertes.fr/hal-03174087. [8] L. Balilescu, J. S. Martin and T. Takahashi, Fluid-rigid structure interaction system with Coulomb’s law, SIAM J. Math. Anal, vol. 49, no. 6, pp. 4625-4657, 2017. [9] D. Gerard-Varet, M. Hillairet and C. Wang, The influence of boundary conditions on the contact problem in a 3d Navier-Stokes flow, J. Math. Pures Appl., vol. 103, p. 1- 38, 2015. [10] D. Gérard-Varet and M. Hillairet, Regularity issues in the problem of fluid structure interaction, Arch. Ration. Mech. Anal., vol. 195, p. 375-407, 2010. [11] D. Gérard-Varet, D., M Hillairet, Computation of the Drag Force on a Sphere Close to a Wall: The Roughness Issue. ESAIM: Mathematical Modelling and Numerical Analysis 46.5 (2012): 1201-224

A nonlinear reduced model based on optimal transport for electronic structure calculations

• Geneviève Dusson (CNRS, Université Bourgogne Franche-Comté)

Room: Amphitheater, Building B Electronic structure calculations are widely used to predict the physical properties of molecules and materials. They require to solve nonlinear partial differential and eigenvalue equations. These equations are generally numerically very demanding, especially since they are parameterized by the positions of the nuclei in the molecule and must be solved a large number of times when these positions vary. This is the case for example when simulating the dynamics of a molecule. In this talk, I will present a recent work aimed at efficiently calculating approximate solutions of such parameterized PDEs, with the objective of reducing the overall computational time. For this, I will present a non-linear interpolation method between several solutions, based on optimal transport, and using in particular Wasserstein barycenters. I will illustrate this method with simulations carried out on a 1D toy model.

On the approximation of the von Neumann equation in the semi-classical limit

• Francis Filbet (Université Paul Sabatier)

Room: Amphitheater, Building B We propose a new approach to discretize the von Neumann equation, which is efficient in the semi-classical limit. This method is first based on the so called Weyl’s variables to address the stiffness associated with the equation. Then, by applying a truncated Hermite expansion of the density operator, we successfully handle this stiffness. Additionally, we develop a finite volume approximation for practical implementation and conduct numerical simulations to illustrate the efficiency of our approach. This asymptotic preserving numerical approximation, combined with the use of Hermite polynomials, provides an efficient tool for solving the von Neumann equation in all regimes, near classical or not.

Neural and hybrid methods for elliptic problems

• Emmanuel Franck (Inria Centre at Université de Lorraine)

Room: Amphitheater, Building B First, we'll introduce the neural methods used to solve PDEs, such as PINNs or the Deep Ritz method. It will be shown that these approaches can fit within the framework of classical Galerkin methods, where only the approximation space changes. The advantages and shortcomings of these approaches will be discussed. Next, a "prediction-correction" approach will be proposed, in which neural approaches are used to quickly predict a solution, which is then corrected by a coarse numerical method. This approach will be applied to two problems for which convergence proofs and numerical examples will be given. In the final section we will introduce a recent purely neural approach that partially overcomes the accuracy and convergence shortcomings of the basic approaches.

A multi-dimensional staggered scheme for the diffusive limit in the radiative transfer equation

• Khaled Saleh (Université Claude Bernard, Lyon 1)

Room: Amphitheater, Building B The radiative transfer equation is a kinetic PDE modelling the specific radiation intensity carried by a population of photons described by a statistical description, i.e. a transport equation on the fraction of photons travelling in a given direction. It is well known that as the Knudsen number (which is the ratio of the mean free path length to a representative physical length scale) goes to zero, the radiation intensity tends to a solution of a diffusion problem. In this talk, we present a numerical scheme for the radiative transfer equation that has the asymptotic preserving property: when the Knudsen number is fixed, we prove that the numerical solution converges to a solution of the radiative transport equation for vanishing discretisation parameters. And for a fixed discretisation, the numerical solution converges to the solution of a stable and consistent numerical scheme for the limit diffusion equation. The numerical scheme considered is an extension to the multidimensional setting of the 1D scheme developed by Lemou and Mieussens in their 2008 SIAM paper. It is based on a micro-macro decomposition of the main unknown and on a staggered discretisation: the macroscopic variable is cell-centred while the microscopic variable is face-centred. We show that special consistency problems arise in the multi-dimensional setting due to the fact that composing a consistent discrete divergence with a (weakly) consistent discret gradient does not generally yield a consistent Laplacian operator, even on admissible grids. This is a joint work with Mohamed Ghattassi and Nader Masmoudi.

A structure-preserving semi-implicit IMEX finite volume scheme for ideal magnetohydrodynamics at all Mach and Alfvén numbers

• Andrea Thomann (Inria Centre at Université de Lorraine)

Room: Amphitheater, Building B We present a divergence-free semi-implicit finite volume scheme for the simulation of the ideal magnetohydrodynamics (MHD) equations which is stable for large time steps controlled by the local transport speed at all Mach and Alfvén numbers. An operator splitting technique allows to treat the convective terms explicitly while the hydrodynamic pressure and the magnetic field contributions are integrated implicitly, yielding two decoupled linear implicit systems. The linearity of the implicit part makes the scheme very efficient and is achieved by means of a semi-implicit time linearization. This structure is favorable as second-order accuracy in time can be achieved relying on the class of semi-implicit IMplicit-EXplicit Runge-Kutta (IMEX-RK) methods. In space, implicit cell-centered finite difference operators are designed to discretely preserve the divergence-free property of the magnetic field on three-dimensional Cartesian meshes avoiding a staggering of the mesh. The new scheme is also particularly well suited for low Mach number flows towards the incompressible limit of the MHD equations, since no explicit numerical dissipation is added to the implicitly treated sub-systems and the time step is scale independent. Likewise, highly magnetized flows can benefit from the implicit treatment of the magnetic fluxes, hence improving the computational efficiency of the novel method.

Designing conservative and accurately dissipative numerical integrators in time

• Patrick E. Farrell (University of Oxford)

Room: Amphitheater, Building B Numerical methods for the simulation of transient systems with structure-preserving properties are known to exhibit greater accuracy and physical reliability, in particular over long durations. These schemes are often built on powerful geometric ideas for broad classes of problems, such as Hamiltonian or reversible systems. However, there remain difficulties in devising higher-order- in-time structure-preserving discretizations for nonlinear problems, and in conserving non-polynomial invariants. In this work we propose a new, general framework for the construction of structure-preserving timesteppers via finite elements in time and the systematic introduction of auxiliary variables. The framework reduces to Gauss methods where those are structure-preserving, but extends to generate arbitrary-order structure-preserving schemes for nonlinear problems, and allows for the construction of schemes that conserve multiple higher-order invariants. We demonstrate the ideas by devising novel schemes that exactly conserve all known invariants of the Kepler and Kovalevskaya problems, high-order energy-conserving and entropy-dissipating schemes for the compressible Navier–Stokes equations, and multi-conservative schemes for the Benjamin-Bona-Mahony equation.

A posteriori error control in the max norm for the Monge-Ampère equation

• Dietmar Gallistl (Universität Jena)

Room: Amphitheater, Building B This talk discusses a stability result for the Monge-Ampère operator in a (potentially regularized) Hamilton-Jacobi-Bellman format as a consequence of Alexandrov's classical maximum principle. The main application is guaranteed a posteriori error control in the $L^\infty$ norm for the difference of the Monge-Ampère solution and the convex hull of a fairly arbitrary $C^1$-conforming finite element approximation.

Analysis and numerical approximation of mean field game partial differential inclusions

• Iain Smears (University College London)

Room: Amphitheater, Building B Joint work with Yohance A. P. Osborne Mean field games (MFG) are models for differential games involving large numbers of players, where each player is solving a dynamic optimal control problem that may depend on the overall distribution of players across the state space of the game. In a standard formulation, the Nash equilibria of the game are characterized by the solutions of a coupled system of partial differential equations, involving the Hamilton-Jacobi-Bellman equation for the value function and the Fokker-Planck equation for the density of players over the state space of the game. However, in many realistic applications, the underlying optimal control problems can lead to systems with nondifferentiable Hamiltonians, such as in minimal time problems, problems with bang-bang controls, etc. This leads to the crucial issue that the PDE system is then not well-defined in the usual sense. From a modelling perspective, this corresponds to nonuniquess of optimal controls, and the question of how players choose among the optimal controls and how this determines the resulting density of players. In this talk, we show that a suitable generalization of the problem is provided by relaxing the Fokker-Planck equation to a partial differential inclusion (PDI) involving the subdifferential of the Hamiltonian, which expresses mathematically the idea that, in the nondifferentiable case, the structure of the Nash equilibria can become more complicated since players in the same state may be required to make distinct choices among the various optimal controls. Our analytical contributions include theorems on the existence of solutions of the resulting MFG PDI system under very general conditions on the problem data, allowing for both local/nonlocal and nonsmoothing nonlinear couplings, for both the steady-state and the time-dependent cases in the stochastic setting. We also show that the MFG PDI system conserves uniqueness of the solution for monotone couplings, as a generalization of the result of Lasry and Lions. We also give concrete examples of some nontrivial Nash equilibria that can be modelled by our approach. Regarding the numerical analysis, we also propose and analyse a stabilized finite element method for the PDI system and we present theorems on its well-posedness and its convergence. We also prove that the method is almost quasi-optimal in the sense of near-best approximations for the case of differentiable Hamiltonians, which leads to optimal rates of convergence for solutions with sufficient regularity, and also some theorems illustrating the robustness of the approximation of the value function relative to the density. We present numerical experiments for both steady-state and time-dependent problems.