crankNicolson2D {sdetorus}

Crank–Nicolson finite difference scheme for the 2D Fokker–Planck equation with periodic boundaries

Description

Implementation of the Crank–Nicolson scheme for solving the Fokker–Planck equation

p(x, y, t)_t = -(p(x, y, t) b_1(x, y))_x -(p(x, y, t) b_2(x, y))_y+

+ \frac{1}{2}(\sigma_1^2(x, y) p(x, y, t))_{xx} + \frac{1}{2}(\sigma_2^2(x, y) p(x, y, t))_{yy} + (\sigma_{12}(x, y) p(x, y, t))_{xy},

where p(x, y, t) is the transition probability density of the toroidal diffusion

dX_t=b_1(X_t,Y_t)dt+\sigma_1(X_t,Y_t)dW^1_t+\sigma_{12}(X_t,Y_t)dW^2_t,

dY_t=b_2(X_t,Y_t)dt+\sigma_{12}(X_t,Y_t)dW^1_t+\sigma_2(X_t,Y_t)dW^2_t.

Usage

crankNicolson2D(u0, bx, by, sigma2x, sigma2y, sigmaxy, N, deltat, Mx, deltax,
  My, deltay, imposePositive = 0L)

Arguments

u0

matrix of size c(Mx * My, 1) giving the initial condition matrix column-wise stored. Typically, the evaluation of a density highly concentrated at a given point. If nt == 1, then u0 can be a matrix c(Mx * My, nu0) containing different starting values in the columns.

bx, by

matrices of size c(Mx, My) containing the evaluation of the drift in the first and second space coordinates, respectively.

sigma2x, sigma2y, sigmaxy

matrices of size c(Mx, My) containing the evaluation of the entries of the diffusion matrix (it has to be positive definite)
rbind(c(sigma2x, sigmaxy), c(sigmaxy, sigma2y)).

N

increasing integer vector of length nt giving the indexes of the times at which the solution is desired. The times of the solution are delta * c(0:max(N))[N + 1].

deltat

time step.

Mx, My

sizes of the equispaced spatial grids in [-\pi,\pi) for each component.

deltax, deltay

space grid discretizations for each component.

imposePositive

flag to indicate whether the solution should be transformed in order to be always larger than a given tolerance. This prevents spurious negative values. The tolerance will be taken as imposePositiveTol if this is different from FALSE or 0.

Details

The function makes use of solvePeriodicTridiag for obtaining implicitly the next step in time of the solution.

If imposePositive = TRUE, the code implicitly assumes that the solution integrates to one at any step. This might b unrealistic if the initial condition is not properly represented in the grid (for example, highly concentrated density in a sparse grid).

Value

• If nt > 1, a matrix of size c(Mx * My, nt) containing the discretized solution at the required times with the c(Mx, My) matrix stored column-wise.

• If nt == 1, a matrix of size c(Mx * My, nu0) containing the discretized solution at a fixed time for different starting values.

References

Thomas, J. W. (1995). Numerical Partial Differential Equations: Finite Difference Methods. Springer, New York. doi:10.1007/978-1-4899-7278-1

Examples

# Parameters
Mx <- 100
My <- 100
N <- 200
x <- seq(-pi, pi, l = Mx + 1)[-c(Mx + 1)]
y <- seq(-pi, pi, l = My + 1)[-c(My + 1)]
m <- c(pi / 2, pi)
p <- c(0, 1)
u0 <- c(outer(dWn1D(x, p[1], 0.5), dWn1D(y, p[2], 0.5)))
bx <- outer(x, y, function(x, y) 5 * sin(m[1] - x))
by <- outer(x, y, function(x, y) 5 * sin(m[2] - y))
sigma2 <- matrix(1, nrow = Mx, ncol = My)
sigmaxy <- matrix(0.5, nrow = Mx, ncol = My)

# Full trajectory of the solution (including initial condition)
u <- crankNicolson2D(u0 = cbind(u0), bx = bx, by = by, sigma2x = sigma2,
                     sigma2y = sigma2, sigmaxy = sigmaxy,
                     N = 0:N, deltat = 1 / N, Mx = Mx, deltax = 2 * pi / Mx,
                     My = My, deltay = 2 * pi / My)

# Mass conservation
colMeans(u) * 4 * pi^2

# Only final time
v <- crankNicolson2D(u0 = cbind(u0), bx = bx, by = by, sigma2x = sigma2,
                     sigma2y = sigma2, sigmaxy = sigmaxy,
                     N = N, deltat = 1 / N, Mx = Mx, deltax = 2 * pi / Mx,
                     My = My, deltay = 2 * pi / My)
sum(abs(u[, N + 1] - v))

## Not run: 
# Visualization of tpd
library(manipulate)
manipulate({
  plotSurface2D(x, y, z = matrix(u[, j + 1], Mx, My),
                main = round(mean(u[, j + 1]) * 4 * pi^2, 4),
                levels = seq(0, 2, l = 21))
  points(p[1], p[2], pch = 16)
  points(m[1], m[2], pch = 16)
}, j = slider(0, N))

## End(Not run)