Getting started

General workflow

The general workflow for using HighDimPDE.jl is as follows:

  • Define a Partial Integro Differential Equation problem
HighDimPDE.HighDimPDEModule

Build Status

HighDimPDE.jl

HighDimPDE.jl is a Julia package to solve Highly Dimensional non-local, non-linear PDEs of the form

$

\begin{aligned} (\partialt u)(t,x) &= \int{\Omega} f\big(t,x,{\bf x}, u(t,x),u(t,{\bf x}), ( \nablax u )(t,x ),( \nablax u )(t,{\bf x} ) \big) d{\bf x} \ & \quad + \big\langle \mu(t,x), ( \nablax u )( t,x ) \big\rangle + \tfrac{1}{2} \text{Trace} \big(\sigma(t,x) [ \sigma(t,x) ]^* ( \text{Hess}x u)(t, x ) \big). \end{aligned} $

where $u \colon [0,T] \times \Omega \to \mathbb{R}, \Omega \subseteq \mathbb{R}^{d}$ is subject to initial and boundary conditions, and where $d$ is large.

Documentation

Installation

Open Julia and type the following

using Pkg;
Pkg.add("HighDimPDE.jl")

This will download the latest version from the git repo and download all dependencies.

Getting started

See documentation and test folders.

Reference

  • Boussange, V., Becker, S., Jentzen, A., Kuckuck, B., Pellissier, L., Deep learning approximations for non-local nonlinear PDEs with Neumann boundary conditions. arXiv (2022)

<!– - Becker, S., Braunwarth, R., Hutzenthaler, M., Jentzen, A., von Wurstemberger, P., Numerical simulations for full history recursive multilevel Picard approximations for systems of high-dimensional partial differential equations. arXiv (2020)

  • Beck, C., Becker, S., Cheridito, P., Jentzen, A., Neufeld, A., Deep splitting method for parabolic PDEs. arXiv (2019)
  • Han, J., Jentzen, A., E, W., Solving high-dimensional partial differential equations using deep learning. arXiv (2018) –>

<!– ## Acknowledgements HighDimPDE.jl is inspired from Sebastian Becker's scripts in Python, TensorFlow and C++. Pr. Arnulf Jentzen largely contributed to the theoretical developments of the solver algorithms implemented. –>

HighDimPDE.PIDEProblemType
PIDEProblem(g, f, μ, σ, x, tspan, p = nothing, x0_sample=nothing, neumann_bc=nothing)

Defines a Partial Integro Differential Problem, of the form

\[\begin{aligned} \frac{du}{dt} &= \tfrac{1}{2} \text{Tr}(\sigma \sigma^T) \Delta u(x, t) + \mu \nabla u(x, t) \\ &\quad + \int f(x, y, u(x, t), u(y, t), ( \nabla_x u )(x, t), ( \nabla_x u )(y, t), p, t) dy, \end{aligned}\]

with $u(x,0) = g(x)$.

Arguments

  • g : initial condition, of the form g(x, p, t).
  • f : nonlinear function, of the form f(x, y, u(x, t), u(y, t), ∇u(x, t), ∇u(x, t), p, t).
  • μ : drift function, of the form μ(x, p, t).
  • σ : diffusion function σ(x, p, t).
  • x: point where u(x,t) is approximated. Is required even in the case where x0_sample is provided.
  • tspan: timespan of the problem.
  • p: the parameter vector.
  • x0_sample : sampling method for x0. Can be UniformSampling(a,b), NormalSampling(σ_sampling, shifted), or NoSampling (by default). If NoSampling, only solution at the single point x is evaluated.
  • neumann_bc: if provided, neumann boundary conditions on the hypercube neumann_bc[1] × neumann_bc[2].
  • Select a solver algorithm
  • Solve the problem

Examples

Let's illustrate that with some examples.

MLP

Local PDE

Let's solve the Fisher KPP PDE in dimension 10 with MLP.

\[\partial_t u = u (1 - u) + \frac{1}{2}\sigma^2\Delta_xu \tag{1}\]

using HighDimPDE

## Definition of the problem
d = 10 # dimension of the problem
tspan = (0.0,0.5) # time horizon
x0 = fill(0.,d)  # initial point
g(x) = exp(- sum(x.^2) ) # initial condition
μ(x, p, t) = 0.0 # advection coefficients
σ(x, p, t) = 0.1 # diffusion coefficients
f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = max(0.0, v_x) * (1 -  max(0.0, v_x)) # nonlocal nonlinear part of the
prob = PIDEProblem(g, f, μ, σ, x0, tspan) # defining the problem

## Definition of the algorithm
alg = MLP() # defining the algorithm. We use the Multi Level Picard algorithm

## Solving with multiple threads 
sol = solve(prob, alg, multithreading=true)

Non local PDE with Neumann boundary conditions

Let's include in the previous equation non local competition, i.e.

\[\partial_t u = u (1 - \int_\Omega u(t,y)dy) + \frac{1}{2}\sigma^2\Delta_xu \tag{2}\]

where $\Omega = [-1/2, 1/2]^d$, and let's assume Neumann Boundary condition on $\Omega$.

using HighDimPDE

## Definition of the problem
d = 10 # dimension of the problem
tspan = (0.0,0.5) # time horizon
x0 = fill(0.,d)  # initial point
g(x) = exp( -sum(x.^2) ) # initial condition
μ(x, p, t) = 0.0 # advection coefficients
σ(x, p, t) = 0.1 # diffusion coefficients
mc_sample = UniformSampling(fill(-5f-1, d), fill(5f-1, d))
f(x, y, v_x, v_y, ∇v_x, ∇v_y, t) = max(0.0, v_x) * (1 -  max(0.0, v_y)) 
prob = PIDEProblem(g, f, μ, 
                    σ, x0, tspan) # defining x0_sample is sufficient to implement Neumann boundary conditions

## Definition of the algorithm
alg = MLP(mc_sample = mc_sample ) 

sol = solve(prob, alg, multithreading=true)

DeepSplitting

Let's solve the previous equation with DeepSplitting.

using HighDimPDE

## Definition of the problem
d = 10 # dimension of the problem
tspan = (0.0, 0.5) # time horizon
x0 = fill(0f0,d)  # initial point
g(x) = exp.(- sum(x.^2, dims=1) ) # initial condition
μ(x, p, t) = 0.0f0 # advection coefficients
σ(x, p, t) = 0.1f0 # diffusion coefficients
x0_sample = UniformSampling(fill(-5f-1, d), fill(5f-1, d))
f(x, y, v_x, v_y, ∇v_x, ∇v_y, p, t) = v_x .* (1f0 .- v_y)
prob = PIDEProblem(g, f, μ, 
                    σ, x0, tspan, 
                    x0_sample = x0_sample)

## Definition of the neural network to use
using Flux # needed to define the neural network

hls = d + 50 #hidden layer size

nn = Flux.Chain(Dense(d, hls, tanh),
        Dense(hls, hls, tanh),
        Dense(hls, 1)) # neural network used by the scheme

opt = ADAM(1e-2)

## Definition of the algorithm
alg = DeepSplitting(nn,
                    opt = opt,
                    mc_sample = x0_sample)

sol = solve(prob, 
            alg, 
            0.1, 
            verbose = true, 
            abstol = 2e-3,
            maxiters = 1000,
            batch_size = 1000)

Solving on the GPU

DeepSplitting can run on the GPU for (much) improved performance. To do so, just set use_cuda = true.

sol = solve(prob, 
            alg, 
            0.1, 
            verbose = true, 
            abstol = 2e-3,
            maxiters = 1000,
            batch_size = 1000,
            use_cuda=true)