# crankNicolson2D: Crank-Nicolson finite difference scheme for the 2D... In sdetorus: Statistical Tools for Toroidal Diffusions

## 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+

p(x, y, t)_t = -(p(x, y, t) * b_1(x, y))_x -(p(x, y, t) * b_2(x, y))_y + 1/2 * (σ_1^2(x, y) *p(x, y, t))_{xx} + 1/2 * (σ_2^2(x, y) p(x, y, t))_{yy} + (σ_{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+σ_1(X_t,Y_t)dW^1_t+σ_{12}(X_t,Y_t)dW^2_t,

dY_t=b_2(X_t,Y_t)dt+σ_{12}(X_t,Y_t)dW^1_t+σ_2(X_t,Y_t)dW^2_t.

## Usage

 ```1 2``` ```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 [-π,π) 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

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42``` ```# 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, 0.5), dWn1D(y, p, 0.5))) bx <- outer(x, y, function(x, y) 5 * sin(m - x)) by <- outer(x, y, function(x, y) 5 * sin(m - 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, p, pch = 16) points(m, m, pch = 16) }, j = slider(0, N)) ## End(Not run) ```

sdetorus documentation built on Aug. 19, 2021, 9:06 a.m.