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

• ∂tm − 1

2

• i,j ∂xixj (aij[m](x, t)m) + div(b[m](x, t) m) = 0

in Rd × R+ m(·, 0) = m0(·) in Rd, (1) describes the evolution of the law of the diffusion process X(t)

dX(s′) = b[m](X(s′), s′))ds′ + σ[m](X(s′), s′)dW(s′), X(s) = x, (2) where b[m] : Rd × R+ × R → Rd, given a vector field, depending non locally on m; σ[m] : Rd × R+ × R → Rd×r, (A[m](x, t))i,j = aij[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 m0; m0 ≥ 0 and

• Rd m0(x)dx = 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

• ∂tv − σ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) = −f2(m(x, t))Dv[m](x, t) where v[m] is the solution of a stationary HJB |Dv| = 1 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 (Rd) we have that

d dt

• Rd φ(x)dm(t)(x)

=

• Rd

  1

2

• i,j

aij[m](x, t)∂xixjφ(x) +

• i

bi[m](x, t)φxi(x)   dm(t)(x) =

• Rd
• 1

2tr(A[m](x, t)D2φ(x)) + b[m](x, t)⊤Dφ(x)

• dm(t)(x)

(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

• f 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 tk = 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:

• R

φ(x)dm(tk+1) =

• R

φ(x)dm(tk) (4) + tk+1

tk

• R

[b[m](x, t)Dφ(x) + 1

2σ2[m](x, t)D2φ(x)]dm(t)dt.

8 / 31

Semi-discrete in time

We first approximate (4) as

• R

φ(x)dm(tk+1) =

• R

[φ(x) + ∆tb[mk](x, tk)Dφ(x) + ∆t 2 σ2[mk](x, tk)D2φ(x)]dm(tk+1). Note that the right hand side corresponds to a Taylor expansion. Hence we write

• R

φ(x)dm(tk+1) = 1 2

• R

[φ(x + ∆tb[mk](x, tk) + √ ∆tσ2[mk](x, tk)dm(tk)+ 1 2

• R

[φ(x + ∆tb[mk](x, tk) − √ ∆tσ2[mk](x, tk)]dm(tk).

9 / 31

Fully-discrete in space

Given ∆x > 0 consider a space grid, for i = 1, . . . , M Ei = [xi − 1

2∆x, xi + 1 2∆x],

mi,k :=

• Ei dm(tk).

(5) By the standard rectangular quadrature formula, we get

• j∈Z

φ(xj)mj,k+1 = (6) 1 2

• j∈Z

φ(Φ+

j [mk)])mj,k + 1

2

• j∈Z

φ(Φ−

j [mk)])mj,k,

where, for µ ∈ RM, j ∈ Z, k = 0, . . . , N − 1, we have defined Φ±

j [µ] := xj + ∆t b[µ](xj, tk) ±

√ ∆tσ[µ](xj, tk). (7)

10 / 31

Interpretation of the scheme by means of characteristics

Note that Φ defined in (7) (with µj = m(xj, tk)) can be interpreted as a single Euler step approximation of dX(s) = b[m] (X(s), s)) ds + σ[m] (X(s), s)) dW(s), s ∈ (tk, tk+1), X(tk) = xj, (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 2 (9) the function Φ±

j [mk] corresponds to one realization of

xj + ∆tb[mk](xj, tk) ± √ 2∆tσ[mk](xj, tk)Z.

11 / 31

Fully-discrete in space

• 2
• 1,5
• 1
• 0,5

Figure: P1-basis function, βi

Given i ∈ Z and setting in (6) φ = βi we have for all i ∈ Z      mi,k+1 = G(mk, i, k) =

1 2

• j∈Z
• βi(Φ+

j [mk]) + βi(Φ− j [mk])

• mj,k

mi,0 =

• Ei m0(x)dx

(10) Non-negative : mi,k ≥ 0 for k = 0, . . . , N − 1, i ∈ Z Mass conservative : ∆x

i mi,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 Ω ⊂ R2, we denote by ˆ n(x) the outward normal to Ω at x ∈ ∂Ω. Homogeneous Neumann condition on ΓN ⊆ ∂Ω J[m](x, t) · ˆ n = 0

• n

∂ΓN, (11) where Ji[m](x, t) = bi[m](x, t)m −

k ∂xk(ai,k[m](x, t)m)

Dirichlet boundary condition on ΓD := ∂Ω \ ΓN m = 0

• n

Γ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[mk](xi, tk) ±

• 2dγσℓ[mk](xi, tk) ∈ ΓD} ∧ h,

and redefine Φℓ

i,± as

Φℓ

i,± := xi +

hℓ,±

i,k b[mk](xi, tk) +

• 2d

hℓ

i,±σℓ[mk](xi, tk).

• 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 : Rd → Ω, defined as P(z) :=        z, if z ∈ Ω, 2w∗ − z, if z / ∈ Ω, where w∗ := argmin

w∈Ω

|z − w|. (13)

• 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(xj) = δi,j and

i βi(x) = 1

for each x ∈ Ω. Median dual control volume Ei :=

• T∈T :xi∈∂T

Ei,T , Ei,T := {x ∈ T : βj(x) ≤ βi(x) j = i} . We approximate the solution m by a sequence {mk}k, such that mi,k ≃

• Ei

dm(tk)

16 / 31

Consistency

Let us denote m∆x,∆t a proper extension in Rd × [0, T] of the discrete approximation mi,k and P1 a propability space, metrized by the Kantorowich-Rubinstein distance:

d1(µ1, µ2) := sup

• Rd 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 φ ∈ C∞

0 (R) if (∆xn, ∆tn) is s.t. as n → ∞ (∆xn, ∆tn) → 0 and ∆xn ∆tn → 0, and

m∆xn,∆tn → m in C(0, T; P1(Rd)), then for kn such that tkn → t lim

n→∞

• Rd φ(x) [dm∆xn,∆tn(tkn+1) − dm∆xn,∆tn(0)]

=

• R

t 1

2tr(A[m](x, s)D2φ) + b[m](x, s)⊤Dφ

• dm(s).

17 / 31

Stability

Proposition

Suppose that (∆x)2 = O(∆t), b and σ continuous functions s.t. there exists C > 0 |b(x, t, µ)| ≤ C|x| |σ(x, t, µ)| ≤ C|x|, for any x, t, µ. Then, there exists a constant C > 0 (independent of (∆x, ∆t)) such that:

• Rd |x|2dm∆x,∆t(t) ≤ C

∀ t ∈ [0, T]. (14)

Proposition

Suppose that (∆x)2 = O(∆t). Then, there exists a constant C > 0 (independent of (∆x, ∆t)) such that: d1(m∆x,∆t(t1), m∆x,∆t(t2)) ≤ C|t1 − t2|

1 2

∀ t1, t2 ∈ [0, T]. (15)

18 / 31

Convergence

Theorem

Suppose that (∆x)2 = O(∆t), b and σcontinuous functions Lipschtiz w.r. to x and there exists C > 0 such that for any x, t, µ |b(x, t, µ)| ≤ C|x| |σ(x, t, µ)| ≤ C|x|, then m∆x,∆t → m in C([0, T], P1), where m is solution of (1) and m∆x,∆t is a sequence of solution of (10) (extended to a general dimension d).

19 / 31

Interacting species

Let m1 and m2 be the densities of two interacting species          ∂tm1 = div(m1(∇E1(m1) + ∇U1(m1, m2))) ∂tm2 = div(m2(∇E2(m2) + ∇U2(m1, m2))) b.c. internal energy E1(m) = E2(m) = 3

2m2 ⋆ Kε (repulsive

self-interactions) cross interaction potential U1(m1, m2) = W11 ⋆ m1 + W21 ⋆ m2 U1(m1, m2) = W12 ⋆ m1 + W22 ⋆ m2

• Ref. J.D.Benamou, G. Carlier, M.Laborde

20 / 31

Interacting species

Let us set W11 = 3|x|2

2 ,W21 = 3|x|2 2

W12 = −|x|2

2

,W22 = −x|2

2

we solve the system on Ω = [−1.5, 1.5] × [−1.5, 1.5] and [0, T] = [0, 3], we set h = 0.03 and ∆t = h/2 Homogenous Neumann Condition

21 / 31

Interacting species

−1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5

t = 0 t = 0.6 t = 1.2

−1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5

t = 1.8 t = 2.4 t = 3

22 / 31

Interacting species without self-interactions (E1 = E2 = 0

Figure: t = 3

23 / 31

Interacting species without selfinteractions (E1 = E2 = 0)

−1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5

t = 0 t = 0.6 t = 1.2

−1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5 −1.5 −1 −0.5 0.5 1 1.5

t = 1.8 t = 2.4 t = 3

24 / 31

Interacting species without selfinteractions (E1 = E2 = 0

Figure: t = 3

25 / 31

The Hughes model

In 2002 R. Hughes proposed a macroscopic model for pedestrian dynamics given by    ∂tm(x, t) − div(m(x, t) f2(m(x, t))Du(x, t)) = 0, |Du(x, t)| = 1 f(m(x, t)), (16) where x ∈ Ω , t ∈ (0, T], T ∈ R+, m corresponds to the density of a crowd trying to exit a domain as fast as possible u the weighted shortest distance to a target, for example an exit. Typical choice for f: f(m) = 1 − m where 1 corresponds to the maximum scaled pedestrian density. Few analytic results are available, all of them restricted to spatial dimension one. In this case the dependence of the drift on m is local!, and the convergence analysis does not apply.

26 / 31

A modified version of the Hughes’ model

(Joint work with A.Festa,F.J Silva, M.T. Wolframm) A modified version corresponds to    ∂tm(x, t) − ε∆m(x, t) − div(m(x, t) f2(m(x, t))Du(x, t)) = 0, −ε∆u(x, t) + 1

2|Du(x, t)|2 =

1 2f2(m(x, t)) + δ, (17) in Ω × (0, T). δ > 0 prevents the blow-up of the cost when approaching the maximum density one. ε > 0 represents a small diffusion.

Ref: Di Francesco M, Markowich P A, Pietschmann F P, Wolfram M T (2011)

27 / 31

Exit scenario with two doors

Figure: Density contour lines at time t = 0.3, ε = 0.001. The exits are marked with the letter “E”.

The position of the two exits induces an asymmetric splitting of the crowd

28 / 31

Exit scenario with two doors

Figure: Density contour lines at time t = 0.3 (left) and t = 1.2 (right). M0 = 0.7 and ε = 0.001.

Initially, a large part of the crowd chooses to move to the right exit. After a while this exit gets congested, inducing a part of the population to change objective.

29 / 31

References

Conclusions: we have presented a new SL scheme for a non linear FP main advantage: it is explicit, it allows large time steps, it preserves the positivity and it conserves the mass On going work convergence analysis SL scheme for non linear FP convergence analysis SL scheme for HJB on bounded domain with Neumann end Dirichlet condition Open questions high order SL scheme for FP

30 / 31

References

E.Carlini, F.J. Silva, A Semi-Lagrangian scheme for a degenerate second order Mean Field Game system, DCDS-A, ( 2015) E.Carlini, F.J. Silva, A Semi-Lagrangian scheme for the Fokker-Planck equation, 2nd IFAC Workshop on Control of Systems Governed E.Carlini, A. Festa, F.J. Silva, M.T. Wolfram, A Semi-Lagrangian scheme for a modified version of the Hughes model for pedestrian flow, Dynamic Games and Applications, Springer, (2016)

31 / 31

References

E.Carlini, F.J. Silva, A Semi-Lagrangian scheme for a degenerate second order Mean Field Game system, DCDS-A, ( 2015) E.Carlini, F.J. Silva, A Semi-Lagrangian scheme for the Fokker-Planck equation, 2nd IFAC Workshop on Control of Systems Governed E.Carlini, A. Festa, F.J. Silva, M.T. Wolfram, A Semi-Lagrangian scheme for a modified version of the Hughes model for pedestrian flow, Dynamic Games and Applications, Springer, (2016) THANK YOU

31 / 31