An operator splitting method for solving a class of Fokker-Planck - - PowerPoint PPT Presentation

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

An operator splitting method for solving a class

• f Fokker-Planck equations

Beatrice Gaviraghi

Institut für Mathematik, Universität Würzburg

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Contents

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Contents

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

How to characterize a stochastic process

We consider a continuous-time stochastic processes X = {Xt}t∈I with I ⊆ R+ and range in Rn.

How to characterize a stochastic process

Stochastic differential equations (SDEs). They describe the evolution of the stochastic processes. Partial(-integro) differential equations (PIDEs). They describe the evolution of the probability density functions (PDFs).

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

How to characterize a stochastic process

We consider a continuous-time stochastic processes X = {Xt}t∈I with I ⊆ R+ and range in Rn.

How to characterize a stochastic process

Stochastic differential equations (SDEs). They describe the evolution of the stochastic processes. Partial(-integro) differential equations (PIDEs). They describe the evolution of the probability density functions (PDFs).

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

How to characterize a stochastic process

We consider a continuous-time stochastic processes X = {Xt}t∈I with I ⊆ R+ and range in Rn.

How to characterize a stochastic process

Stochastic differential equations (SDEs). They describe the evolution of the stochastic processes. Partial(-integro) differential equations (PIDEs). They describe the evolution of the probability density functions (PDFs).

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

How to characterize a stochastic process

We consider a continuous-time stochastic processes X = {Xt}t∈I with I ⊆ R+ and range in Rn.

How to characterize a stochastic process

Stochastic differential equations (SDEs). They describe the evolution of the stochastic processes. Partial(-integro) differential equations (PIDEs). They describe the evolution of the probability density functions (PDFs).

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

An initial value problem

We consider an initial-value problem in the interval [t0, T] with a stochastic differential equation:

• dXt = a(Xt, t)dt + b(Xt, t)dWt

Xt0 = X0 with drift coefficient a, diffusion coefficient b, Wiener process W , and given initial data X0. Under growth and smoothness assumptions on the coefficients, there exists a pathwise-unique solution of this problem.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

An initial value problem

We consider an initial-value problem in the interval [t0, T] with a stochastic differential equation:

• dXt = a(Xt, t)dt + b(Xt, t)dWt

Xt0 = X0 with drift coefficient a, diffusion coefficient b, Wiener process W , and given initial data X0. Under growth and smoothness assumptions on the coefficients, there exists a pathwise-unique solution of this problem.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The discretization of a SDE

A sample path of a stochastic process over the interval [t0, T] can be simulated with The Euler-Maruyama method. Other higher-order methods, which are more computationally expensive due to the evaluation of several stochastic integrals. When a stochastic process is simulated using a numerical method, its values are specified at the points of the discrete time grid t0, t1, ..., tN = T.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Examples of discretized sample paths

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The PDF of the Ornstein-Uhlenbeck process on a periodic domain

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Monte Carlo methods

Disadvantages of Monte Carlo methods:

1 the rate of convergence in distribution is proportional to N− 1 2 ,

where N is the number of the samples;

2 the probability density function of the simulated process is not

available in closed form.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The evolution of the PDF of a stochastic process

Given the model

• dXt = a(Xt, t)dt + b(Xt, t)dWt

Xt0 = X0, the evolution of the PDF f (x, t) is given by the following parabolic PDE ∂tf = −

n

• i=1

∂xi (aif ) + 1 2

n

• i,j=1

∂2

xi,xj

• (bbT)ijf
• which is called Fokker-Planck (FP) equation.

In the FP equation, the space dimension corresponds to number of components of the stochastic process.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The evolution of the PDF of a stochastic process

Given the model

• dXt = a(Xt, t)dt + b(Xt, t)dWt

Xt0 = X0, the evolution of the PDF f (x, t) is given by the following parabolic PDE ∂tf = −

n

• i=1

∂xi (aif ) + 1 2

n

• i,j=1

∂2

xi,xj

• (bbT)ijf
• which is called Fokker-Planck (FP) equation.

In the FP equation, the space dimension corresponds to number of components of the stochastic process.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The evolution of the PDF of a stochastic process

Given the model

• dXt = a(Xt, t)dt + b(Xt, t)dWt

Xt0 = X0, the evolution of the PDF f (x, t) is given by the following parabolic PDE ∂tf = −

n

• i=1

∂xi (aif ) + 1 2

n

• i,j=1

∂2

xi,xj

• (bbT)ijf
• which is called Fokker-Planck (FP) equation.

In the FP equation, the space dimension corresponds to number of components of the stochastic process.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The FP problem

We define the FP problem as follows

• ∂tf = − n

i=1 ∂xi (aif ) + 1 2

n

i,j=1 ∂2 xi,xj (σijf )

f (x, 0) = fX0(x), where σ is defined as σ := bbT and the initial condition is given by the density of the initial random variable. The existence and uniqueness of the solution of this problem follows under growth conditions of the coefficients a and σ and under the definite positivity of the matrix σ.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The FP problem

We define the FP problem as follows

• ∂tf = − n

i=1 ∂xi (aif ) + 1 2

n

i,j=1 ∂2 xi,xj (σijf )

f (x, 0) = fX0(x), where σ is defined as σ := bbT and the initial condition is given by the density of the initial random variable. The existence and uniqueness of the solution of this problem follows under growth conditions of the coefficients a and σ and under the definite positivity of the matrix σ.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Additional conditions on the FP problem

The FP problem differs from a classical parabolic PDE with given initial condition. The solution f must be non negative. The total probability must be preserved at any time t:

• Rd f (x, t) = 1.

The initial condition is given by initial distribution of the underlying process. The boundary conditions (BCs) depend on the underlying process.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Additional conditions on the FP problem

The FP problem differs from a classical parabolic PDE with given initial condition. The solution f must be non negative. The total probability must be preserved at any time t:

• Rd f (x, t) = 1.

The initial condition is given by initial distribution of the underlying process. The boundary conditions (BCs) depend on the underlying process.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Additional conditions on the FP problem

The FP problem differs from a classical parabolic PDE with given initial condition. The solution f must be non negative. The total probability must be preserved at any time t:

• Rd f (x, t) = 1.

The initial condition is given by initial distribution of the underlying process. The boundary conditions (BCs) depend on the underlying process.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Additional conditions on the FP problem

The FP problem differs from a classical parabolic PDE with given initial condition. The solution f must be non negative. The total probability must be preserved at any time t:

• Rd f (x, t) = 1.

The initial condition is given by initial distribution of the underlying process. The boundary conditions (BCs) depend on the underlying process.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Additional conditions on the FP problem

The FP problem differs from a classical parabolic PDE with given initial condition. The solution f must be non negative. The total probability must be preserved at any time t:

• Rd f (x, t) = 1.

The initial condition is given by initial distribution of the underlying process. The boundary conditions (BCs) depend on the underlying process.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Examples of boundary conditions

Examples of boundary conditions: Natural BCs: lim

x→±∞ f (x, t) = 0.

Reflecting BCs (impenetrable wall in x = α). Periodic BCs.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Examples of boundary conditions

Examples of boundary conditions: Natural BCs: lim

x→±∞ f (x, t) = 0.

Reflecting BCs (impenetrable wall in x = α). Periodic BCs.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Examples of boundary conditions

Examples of boundary conditions: Natural BCs: lim

x→±∞ f (x, t) = 0.

Reflecting BCs (impenetrable wall in x = α). Periodic BCs.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Contents

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

A jump-diffusion model

Let us consider the jump-diffusion model

• dXt = a(Xt, t)dt + b(Xt, t)dWt + dPt

Xt0 = X0, with X0, a, b and W as before and P compound Poisson process. A compound Poisson process is a piecewise-continuous process. The discontinuities take place with rate λ > 0 and their amplitudes are independent and identically distributed, with common distribution G.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

A jump-diffusion model

Let us consider the jump-diffusion model

• dXt = a(Xt, t)dt + b(Xt, t)dWt + dPt

Xt0 = X0, with X0, a, b and W as before and P compound Poisson process. A compound Poisson process is a piecewise-continuous process. The discontinuities take place with rate λ > 0 and their amplitudes are independent and identically distributed, with common distribution G.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

A sample path of a jump-diffusion process

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The partial-integro differential equation

The evolution of the probability density function of a jump-diffusion process in R is modelled by the following PIDE ∂tf (x, t) = − ∂x (a(x, t)f (x, t)) + 1 2∂2

x ((σ(x, t)f (x, t))

+ λ

• R

[f (x − y, t) − f (x, t)] g(y)dy where a and σ are defined as before, λ is rate of jumps and g density of the jump distribution. Existence and uniqueness of solutions for PIDEs have been studied in the framework of Hölder spaces and the are related to the definite positivity of the diffusion matrix σ.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The partial-integro differential equation

The evolution of the probability density function of a jump-diffusion process in R is modelled by the following PIDE ∂tf (x, t) = − ∂x (a(x, t)f (x, t)) + 1 2∂2

x ((σ(x, t)f (x, t))

+ λ

• R

[f (x − y, t) − f (x, t)] g(y)dy where a and σ are defined as before, λ is rate of jumps and g density of the jump distribution. Existence and uniqueness of solutions for PIDEs have been studied in the framework of Hölder spaces and the are related to the definite positivity of the diffusion matrix σ.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The initial value problem

Given the PIDE ∂tf (x, t) = − ∂x (a(x, t)f (x, t)) + 1 2∂2

x ((σ(x, t)f (x, t))

+ λ

• R

[f (x − y, t) − f (x, t)] g(y)dy with a, σ, λ and g given, we consider the FP problem with jumps

• ∂tf (x, t) = L(f (x, t)) + I(f (x, t))

f (x, 0) = f0(x) where L is the differential operator, I is the integral operator and the Cauchy data is given by the initial distribution.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The solution of a FP problem with jumps

As in the FP problem without jumps, we have that The solution f must be non negative → maximum principle theorems. The total probability must be preserved at any time t:

• R

f (x, t) = 1. The initial condition is given by initial distribution of the underlying process. The boundary conditions depend on the underlying process. An additional problem consists of the treatment of the integral part.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The solution of a FP problem with jumps

As in the FP problem without jumps, we have that The solution f must be non negative → maximum principle theorems. The total probability must be preserved at any time t:

• R

f (x, t) = 1. The initial condition is given by initial distribution of the underlying process. The boundary conditions depend on the underlying process. An additional problem consists of the treatment of the integral part.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Contents

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Truncation of the spatial domain

An unbounded spatial domain has to be replaced by Ω = [α, β] ⇒ boundary conditions need to be introduced. The best choice turns out to be setting reflecting barriers on ∂Ω, so that we have d dt

• Ω

f (x, t)dx = 0. As a consequence,

• Ω

f (x, t)dx =

• Ω

f0(x)dx = 1 for each t ≥ 0.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Truncation of the spatial domain

An unbounded spatial domain has to be replaced by Ω = [α, β] ⇒ boundary conditions need to be introduced. The best choice turns out to be setting reflecting barriers on ∂Ω, so that we have d dt

• Ω

f (x, t)dx = 0. As a consequence,

• Ω

f (x, t)dx =

• Ω

f0(x)dx = 1 for each t ≥ 0.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Truncation of the spatial domain

An unbounded spatial domain has to be replaced by Ω = [α, β] ⇒ boundary conditions need to be introduced. The best choice turns out to be setting reflecting barriers on ∂Ω, so that we have d dt

• Ω

f (x, t)dx = 0. As a consequence,

• Ω

f (x, t)dx =

• Ω

f0(x)dx = 1 for each t ≥ 0.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The Strang-Marchuk splitting scheme

The splitting time step is defined as δt. Given f n and the time window [tn, tn+1], with width δt, f n+1 is computed as follows. 1. ∂f1(x,t)

∂t

= L(f1(x, t)) f1(x, tn) = f n(x) in the time window [tn, tn+ 1

2 ].

2. ∂f2(x,t)

∂t

= I(f2(x, t)) f2(x, tn) = f1(x, tn+ 1

2 )

in the time window [tn, tn+1]. 3. ∂f3(x,t)

∂t

= L(f3(x, t)) f3(x, tn+ 1

2 ) = f2(x, tn+1)

in the time window [tn+ 1

2 , tn+1].

f n+1(x) = f3(x, tn+1). The Strang-Splitting method is of order O(δt2).

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

We combine the Strang-Marchuk splitting with the Euler method. At each time step n = 0, ..., M − 1, we perform the three required steps as follows:

1 Chang-Cooper method with Euler implicit. 2 trapezoidal rule with Euler explicit. 3 Chang-Cooper method with Euler implicit. Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

We combine the Strang-Marchuk splitting with the Euler method. At each time step n = 0, ..., M − 1, we perform the three required steps as follows:

1 Chang-Cooper method with Euler implicit. 2 trapezoidal rule with Euler explicit. 3 Chang-Cooper method with Euler implicit. Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The flux form of the FP equation

The Fokker-Planck equation ∂tf (x, t) = −∂x (a(x, t)f (x, t)) + 1 2∂2

x ((σ(x, t)f (x, t))

can be rewritten in flux form as follows ∂tf (x, t) = ∂xF(x, t). With zero-flux boundary conditions and initial condition given by a probability density function, the Chang-Cooper scheme is conservative.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The flux form of the FP equation

The Fokker-Planck equation ∂tf (x, t) = −∂x (a(x, t)f (x, t)) + 1 2∂2

x ((σ(x, t)f (x, t))

can be rewritten in flux form as follows ∂tf (x, t) = ∂xF(x, t). With zero-flux boundary conditions and initial condition given by a probability density function, the Chang-Cooper scheme is conservative.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The Chang-Cooper scheme (1)

The right hand term in ∂tf (x, t) = ∂xF(x, t) is approximated at (tn, xj) as follows ∂xF n

j ≈ 1

h(F n

j+ 1

2 − F n

j− 1

2 ),

where for example F n

j+ 1

2 =

• (1 − δn

j )Bn j+ 1

2 + 1

hC n

j+ 1

2

• f n+1

j+1 −

1 hC n

j+ 1

2 − δn

j Bn j+ 1

2

• f n+1

j

.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The Chang-Cooper scheme (2)

The quantity δn

j gives a linear convex combination of f n+1 j

and f n+1

j+1 as follows

f n+1

j+ 1

2 = (1 − δn

j ) f n+1 j+1 + δn j f n+1 j

. When δn

j is chosen as follows

δn

j = 1

wj − 1 exp(wj) − 1, where wj = h Bn

j+ 1

2 /C n

j+ 1

2 , the numerical scheme preserves the

positivity of the solution.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The Chang-Cooper scheme (2)

The quantity δn

j gives a linear convex combination of f n+1 j

and f n+1

j+1 as follows

f n+1

j+ 1

2 = (1 − δn

j ) f n+1 j+1 + δn j f n+1 j

. When δn

j is chosen as follows

δn

j = 1

wj − 1 exp(wj) − 1, where wj = h Bn

j+ 1

2 /C n

j+ 1

2 , the numerical scheme preserves the

positivity of the solution.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Stability and accuracy of the Chang-Cooper scheme

Stability The scheme proposed by Chang and Cooper is stable, in the sense that f n ≤ 2n/2f 0 + δt

n−1

• k=0

2

n−k+1 2

hn+1. Accuracy With first-order backward time differencing, the error has order O(δt + h2). With second-order backward time differencing, the error has

• rder O(δt2 + h2).

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Stability and accuracy of the Chang-Cooper scheme

Stability The scheme proposed by Chang and Cooper is stable, in the sense that f n ≤ 2n/2f 0 + δt

n−1

• k=0

2

n−k+1 2

hn+1. Accuracy With first-order backward time differencing, the error has order O(δt + h2). With second-order backward time differencing, the error has

• rder O(δt2 + h2).

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The trapezoidal rule

The integral part of the PIDE is approximated with the trapezoidal rule as follows I(f

n+ 1

2

j

) = λ

• R
• f (xj − y, tn+ 1

2 ) − f (xj, tn+ 1 2 )

• g(y)dy ≈

λ        h  f

n+ 1

2

1

g(x1 − xj) 2 +

N−1

• i=2

f

n+ 1

2

i

g(xi − xj) + f

n+ 1

2

N

g(xN − xj) 2  

• trapezoidal rule

       − λf

n+ 1

2

j

.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The explicit step

In the time window [tn, tn+1], f n+1 is calculated as follows f n+1 = (1 − λ δt)f n + (λ δt h)G f n, where the matrix G comes from the trapezoidal rule and it is defined as follows Gij =       

g(xi−xj) N

k=1 g(xi−xk)

for i = 1, ..., N and j = 2, ..., N − 1

g(xi−xj) 2N

k=1 g(xi−xk)

for i = 1, ..., N and j = 1, N.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Preservation of the properties of the splitting scheme

1 Positivity of the solution

Chang-Cooper step The positivity is ensured by the appropriate choice of δn at each time step n. Trapezoidal rule A sufficient condition for the positivity of the solution at each time step is the choice of δt as follows δt < 1 λ.

2 Conservativeness of the total probability

Chang-Cooper step Zero flux boundary conditions: F(x, t) = 0 on ∂Ω. Trapezoidal rule Normalization of the matrix G.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Preservation of the properties of the splitting scheme

1 Positivity of the solution

Chang-Cooper step The positivity is ensured by the appropriate choice of δn at each time step n. Trapezoidal rule A sufficient condition for the positivity of the solution at each time step is the choice of δt as follows δt < 1 λ.

2 Conservativeness of the total probability

Chang-Cooper step Zero flux boundary conditions: F(x, t) = 0 on ∂Ω. Trapezoidal rule Normalization of the matrix G.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Preservation of the properties of the splitting scheme

1 Positivity of the solution

Chang-Cooper step The positivity is ensured by the appropriate choice of δn at each time step n. Trapezoidal rule A sufficient condition for the positivity of the solution at each time step is the choice of δt as follows δt < 1 λ.

2 Conservativeness of the total probability

Chang-Cooper step Zero flux boundary conditions: F(x, t) = 0 on ∂Ω. Trapezoidal rule Normalization of the matrix G.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Preservation of the properties of the splitting scheme

1 Positivity of the solution

Chang-Cooper step The positivity is ensured by the appropriate choice of δn at each time step n. Trapezoidal rule A sufficient condition for the positivity of the solution at each time step is the choice of δt as follows δt < 1 λ.

2 Conservativeness of the total probability

Chang-Cooper step Zero flux boundary conditions: F(x, t) = 0 on ∂Ω. Trapezoidal rule Normalization of the matrix G.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Preservation of the properties of the splitting scheme

1 Positivity of the solution

Chang-Cooper step The positivity is ensured by the appropriate choice of δn at each time step n. Trapezoidal rule A sufficient condition for the positivity of the solution at each time step is the choice of δt as follows δt < 1 λ.

2 Conservativeness of the total probability

Chang-Cooper step Zero flux boundary conditions: F(x, t) = 0 on ∂Ω. Trapezoidal rule Normalization of the matrix G.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Preservation of the properties of the splitting scheme

1 Positivity of the solution

Chang-Cooper step The positivity is ensured by the appropriate choice of δn at each time step n. Trapezoidal rule A sufficient condition for the positivity of the solution at each time step is the choice of δt as follows δt < 1 λ.

2 Conservativeness of the total probability

Chang-Cooper step Zero flux boundary conditions: F(x, t) = 0 on ∂Ω. Trapezoidal rule Normalization of the matrix G.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Stability and accuracy of the splitting scheme

Stability The proposed splitting scheme is stable, in the sense that f n ≤2n(G2 + 1)nf 0 + δt

n−1

• k=0

2

n−k+1 2

hk+1+ + δt

n−1

• k=0

2n−k+1(G2 + 1))n−k+1hk+ 1

2 .

Accuracy First order time differencing: f − fh,δtL2

h,δt = O(δt + h2)

Second order time differencing: f − fh,δtL2

h,δt = O(δt2 + h2) Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Stability and accuracy of the splitting scheme

Stability The proposed splitting scheme is stable, in the sense that f n ≤2n(G2 + 1)nf 0 + δt

n−1

• k=0

2

n−k+1 2

hk+1+ + δt

n−1

• k=0

2n−k+1(G2 + 1))n−k+1hk+ 1

2 .

Accuracy First order time differencing: f − fh,δtL2

h,δt = O(δt + h2)

Second order time differencing: f − fh,δtL2

h,δt = O(δt2 + h2) Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The initial guess for the numerical experiments

The initial guess f (x, t) = 1 √ 2πσf exp

• −(x − µf t)2

2σ2

f

• solves the Fokker-Planck problem with a source term
• ∂tf (x, t) = L(f (x, t)) + I(f (x, t)) + h(x, t)

f (0, x) = f0(x). We set the jump amplitude to be Gaussian ⇒ the integral is analytically computable!

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The initial guess for the numerical experiments

The initial guess f (x, t) = 1 √ 2πσf exp

• −(x − µf t)2

2σ2

f

• solves the Fokker-Planck problem with a source term
• ∂tf (x, t) = L(f (x, t)) + I(f (x, t)) + h(x, t)

f (0, x) = f0(x). We set the jump amplitude to be Gaussian ⇒ the integral is analytically computable!

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The Ornstein-Uhlenbeck process with jumps (1)

The process solves the initial-value problem

• dXt = −xdt + CdWt + dPt

X0 = 0, with Gaussian jump amplitude and constant rate of jumps λ. Its PDF is the solution of the initial-value problem ∂tf (x, t) = − ∂x (xf (x, t)) + C∂2

xf (x, t)

+ λ

• R

f (y, t)g(y − x)dy − λf (x, t) f (x, 0) = f0(x), with g Gaussian.

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

The Ornstein-Uhlenbeck process with jumps (2)

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

References

• J. S. Chang and G. Cooper, A Practical Difference Scheme for

Fokker-Planck Equations, Journal of Computational Physics (1970) 6, 1-16.

• M. G. Garroni and J. L. Menaldi, Green Functions for Second-Order

Parabolic Integro-Differential Problems, Longman (1992).

• J. Geiser, Decomposition Methods for Differential Equations: Theory and

Applications, Chapman & Hall, (2009).

• M. Mohammadi and A. Borzì, Analysis of the Chang-Cooper

Discretization Scheme for a Class of Fokker-Planck Equations, J. Numer. Math., to appear (2015).

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck

Stochastic processes and the Fokker-Planck equation The Fokker-Planck partial-integro differential equation The numerical scheme

Thank you for your attention!

Beatrice Gaviraghi An operator splitting method for solving a class of Fokker-Planck