A Semi-Lagrangian discretization of non linear fokker Planck - PowerPoint PPT Presentation

A Semi-Lagrangian discretization of non linear fokker Planck equations E. Carlini Universit` a Sapienza di Roma joint works with F.J. Silva RICAM, Linz November 21-25, 2016 1 / 31

Outline 1 A Semi-Lagrangian scheme for a nonlinear Fokker-Planck equation 2 Convergence Analysis 3 Numerical simulation 2 / 31

A nonlinear Fokker-Planck equation The nonlinear FP equation � � in R d × R + ∂ t m − 1 i,j ∂ x i x j ( a ij [ m ]( x, t ) m ) + div ( b [ m ]( x, t ) m ) = 0 2 in R d , m ( · , 0) = m 0 ( · ) (1) describes the evolution of the law of the diffusion process X ( t ) b [ m ]( X ( s ′ ) , s ′ ))d s ′ + σ [ m ]( X ( s ′ ) , s ′ )d W ( s ′ ) , d X ( s ′ ) = (2) X ( s ) = x, where b [ m ] : R d × R + × R → R d , given a vector field, depending non locally on m ; σ [ m ] : R d × R + × R → R d × r , ( A [ m ]( x, t )) i,j = a ij [ m ]( x, t ) := ( σ [ m ]( x, t ) σ [ m ]( x, t ) ⊤ ) ij , is the diffusion matrix (possible degenerate),depending non locally on m ; the density of the initial law is given by m 0 ; � m 0 ≥ 0 and R d m 0 ( x )d x = 1 . 3 / 31

Some applications Non local interactions due to collective phenomena (biophysics, social behavior) Mean Filed Games: b [ m ]( x, t ) = − DH ( Dv [ m ]( x, t )) where v [ m ] is the solution of a backward HJB � ∂ t v − σ 2 2 ∆ v + H ( Dv ) = f ( x, m ( t )) v ( x, T ) = g ( x, m ( T ) Since b [ m ]( x, t ) depends on m ( s ) in all time 0 ≤ s ≤ T we gets an Implicit scheme: Fixed Point Iterations needed to solve it Hughes model b [ m ]( x, t ) = − f 2 ( m ( x, t )) Dv [ m ]( x, t ) where v [ m ] is the solution of a stationary HJB 1 | Dv | = f ( m ( x, t )) Since b [ m ]( x, t ) depends on m ( s ) only in time past 0 ≤ s ≤ t we gets an Explicit scheme 4 / 31

Weak solution For many problems of interest, (1) has only a formal meaning, henceforth we mean m is a weak solution of (1) if for all φ ∈ C ∞ c ( R d ) we have that � d R d φ ( x )d m ( t )( x ) d t   � � �  1  d m ( t )( x ) = a ij [ m ]( x, t ) ∂ x i x j φ ( x ) + b i [ m ]( x, t ) φ x i ( x ) 2 R d i,j i � � � 1 2 tr( A [ m ]( x, t ) D 2 φ ( x )) + b [ m ]( x, t ) ⊤ Dφ ( x ) = d m ( t )( x ) R d (3) 5 / 31

Numerical approximation of FP Some references Linear Chang and Cooper (1970): finite difference scheme s.t. preserves positivity, equilibrium states and mass of the distribution (explicit version is stable under a parabolic type CFL condition) Kushner (1976): finite difference via probabilistic method Naess and Johnsen (1993) : Path Integration Method (time integration of the evolution probability density function, where the transition probability density is approximated by a Gaussian) Jordan, Kinderlehrer,Otto (1998): variational scheme for FokkerPlanck equations for which the drift term is given by the gradient of a potential (JKO) (It preserves positivity, and mass of the distribution, costly in high dimension). Achdou,Camilli,Capuzzo Dolcetta (2012): implicit finite difference scheme s.t. preserves positivity, and mass of the distribution. Non Linear Drozdov, Morillo (1995): finite difference scheme s.t. preserves equilibrium states and mass of the distribution, (high order). Benamou, Carlier, Laborde (2015): semi implicit variant of JKO 6 / 31

Numerical approximation based on SL for FP We propose a Semi-Lagrangian (SL) scheme for non linear FP equation s.t. it is first order accurate (numerically) it is explicit and it allows for large time steps it preserves the positivity of the density and conserves its integral equals to 1 Ref. Silva, Camilli (2013), Carlini and Silva (2014,2015) The scheme has been applied to numerically compute the solution of Mean Field Games model (with F. Silva) a regularized Hughes model for pedestrian flow (with A.Festa, F.Silva, M.T. Wolfram) Hughes model for pedestrian flow with different congestion penalty function (with A.Festa, F.Silva) 7 / 31

A Semi-discrete in time SL for a nonlinear 1d Fokker-Planck equation Given ∆ t , we define t k = k ∆ t , k = 0 , . . . , N . We multiply the first equation in (1) by a smooth test function φ with compact support and integrate by parts to get: � � φ ( x )d m ( t k +1 ) = φ ( x )d m ( t k ) (4) R R � t k +1 � [ b [ m ]( x, t ) Dφ ( x ) + 1 2 σ 2 [ m ]( x, t ) D 2 φ ( x )]d m ( t )d t. + t k R 8 / 31

Semi-discrete in time We first approximate (4) as � φ ( x )d m ( t k +1 ) = R � [ φ ( x ) + ∆ tb [ m k ]( x, t k ) Dφ ( x ) + ∆ t 2 σ 2 [ m k ]( x, t k ) D 2 φ ( x )]d m ( t k +1 ) . R Note that the right hand side corresponds to a Taylor expansion. Hence we write � φ ( x )d m ( t k +1 ) = R � √ 1 ∆ tσ 2 [ m k ]( x, t k )d m ( t k )+ [ φ ( x + ∆ tb [ m k ]( x, t k ) + 2 R � √ 1 ∆ tσ 2 [ m k ]( x, t k )]d m ( t k ) . [ φ ( x + ∆ tb [ m k ]( x, t k ) − 2 R 9 / 31

Fully-discrete in space Given ∆ x > 0 consider a space grid, for i = 1 , . . . , M E i = [ x i − 1 2 ∆ x, x i + 1 2 ∆ x ] , � (5) m i,k := E i d m ( t k ) . By the standard rectangular quadrature formula, we get � φ ( x j ) m j,k +1 = (6) j ∈ Z � � 1 j [ m k )]) m j,k + 1 φ (Φ + φ (Φ − j [ m k )]) m j,k , 2 2 j ∈ Z j ∈ Z where, for µ ∈ R M , j ∈ Z , k = 0 , . . . , N − 1 , we have defined √ Φ ± j [ µ ] := x j + ∆ t b [ µ ]( x j , t k ) ± ∆ tσ [ µ ]( x j , t k ) . (7) 10 / 31

Interpretation of the scheme by means of characteristics Note that Φ defined in (7) (with µ j = m ( x j , t k ) ) can be interpreted as a single Euler step approximation of d X ( s ) = b [ m ] ( X ( s ) , s )) d s + σ [ m ] ( X ( s ) , s )) d W ( s ) , s ∈ ( t k , t k +1 ) , X ( t k ) = x j , (8) with a random walk discretization of the Brownian motion W ( · ) . Indeed, considering a random value Z in R such that f P ( Z = 1) = P ( Z = − 1) = 1 (9) 2 the function Φ ± j [ m k ] corresponds to one realization of √ x j + ∆ tb [ m k ]( x j , t k ) ± 2∆ tσ [ m k ]( x j , t k ) Z. 11 / 31

Fully-discrete in space -2 -1,5 -1 -0,5 0 0,5 1 1,5 2 Figure: P1-basis function, β i Given i ∈ Z and setting in (6) φ = β i we have for all i ∈ Z  m i,k +1 = G ( m k , i, k )   � � � β i (Φ + j [ m k ]) + β i (Φ − 1 = j [ m k ]) m j,k (10) j ∈ Z 2  �  m i, 0 = E i m 0 ( x )d x Non-negative : m i,k ≥ 0 for k = 0 , . . . , N − 1 , i ∈ Z Mass conservative : ∆ x � i m i,k = 1 for k = 0 , . . . , N − 1 Generalizable to any dimension Generalizable to handle Dirichlet and Neumann Boundary conditions Generalizable to handle degeneracy of the diffusion matrix 12 / 31

Boundary Condition Let us consider (1) in Ω ⊂ R 2 , we denote by ˆ n ( x ) the outward normal to Ω at x ∈ ∂ Ω . Homogeneous Neumann condition on Γ N ⊆ ∂ Ω J [ m ]( x, t ) · ˆ n = 0 on ∂ Γ N , (11) where J i [ m ]( x, t ) = b i [ m ]( x, t ) m − � k ∂ x k ( a i,k [ m ]( x, t ) m ) Dirichlet boundary condition on Γ D := ∂ Ω \ Γ N m = 0 on Γ D . (12) The idea is to reflect the discrete characteristics when they intersect Γ N and to cut them at Γ D . 13 / 31

Boundary Condition: Dirichlet, d = 2 Let us define, for l = 1 , 2 � � h ℓ i, ± = inf { γ > 0 ; x + γb [ m k ]( x i , t k ) ± 2 dγσ ℓ [ m k ]( x i , t k ) ∈ Γ D } ∧ h, and redefine Φ ℓ i, ± as � i, ± := x i + � h ℓ, ± 2 d � Φ ℓ h ℓ i,k b [ m k ]( x i , t k ) + i, ± σ ℓ [ m k ]( x i , t k ) . Ref. Falcone, Ferretti(2013) 14 / 31

Boundary Condition: Neumann, d = 2 In order to consider the reflection when exiting at Γ N , we use a symmetrized Euler scheme: we project the discrete characteristics on Ω using the operator P : R d → Ω , defined as   z, if z ∈ Ω ,   2 w ∗ − z, if z / ∈ Ω , P ( z ) := (13)  where w ∗ := argmin   | z − w | . w ∈ Ω Ref. M. Bossy and E. Gobet and D.Talay, (2004) 15 / 31

Unstructured grids Let T be a triangulation of Ω , with ∆ x the maximum diameter of the triangles. We call { β i ; i = 1 , ..., M } the set of affine shape functions related to the triangular mesh, such that β i ( x j ) = δ i,j and � i β i ( x ) = 1 for each x ∈ Ω . Median dual control volume � E i := E i,T , T ∈T : x i ∈ ∂T E i,T := { x ∈ T : β j ( x ) ≤ β i ( x ) j � = i } . We approximate the solution m by a sequence { m k } k , such that � m i,k ≃ d m ( t k ) E i 16 / 31

Consistency Let us denote m ∆ x, ∆ t a proper extension in R d × [0 , T ] of the discrete approximation m i,k and P 1 a propability space, metrized by the Kantorowich-Rubinstein distance: �� � d 1 ( µ 1 , µ 2 ) := sup R d f ( x )d( µ 1 ( x ) − µ 2 ( x )) ; � Df � ∞ ≤ 1 . Proposition (Weak consistency) Assuming that b and σ are continuous and Lipschitz w.r. to x , for every 0 ( R ) if (∆ x n , ∆ t n ) is s.t. as n → ∞ (∆ x n , ∆ t n ) → 0 and ∆ x n φ ∈ C ∞ ∆ t n → 0 , and C (0 , T ; P 1 ( R d )) , m ∆ x n , ∆ t n → m in then for k n such that t k n → t � lim R d φ ( x ) [ dm ∆ x n , ∆ t n ( t k n +1 ) − dm ∆ x n , ∆ t n (0)] n →∞ � � t � 1 � 2 tr( A [ m ]( x, s ) D 2 φ ) + b [ m ]( x, s ) ⊤ Dφ = d m ( s ) . 0 R 17 / 31