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

## Description

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

p(x, t)_t = -(p(x, t) * b(x))_x + 1/2 * (σ^2(x) p(x, t))_{xx},

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

dX_t=b(X_t)dt+σ(X_t)dW_t

.

## Usage

 `1` ```crankNicolson1D(u0, b, sigma2, N, deltat, Mx, deltax, imposePositive = 0L) ```

## Arguments

 `u0` matrix of size `c(Mx, 1)` giving the initial condition. Typically, the evaluation of a density highly concentrated at a given point. If `nt == 1`, then `u0` can be a matrix `c(Mx, nu0)` containing different starting values in the columns. `b` vector of length `Mx` containing the evaluation of the drift. `sigma2` vector of length `Mx` containing the evaluation of the squared diffusion coefficient. `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` size of the equispaced spatial grid in [-π,π). `deltax` space grid discretization. `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, nt)` containing the discretized solution at the required times.

• If `nt == 1`, a matrix of size `c(Mx, 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``` ```# Parameters Mx <- 200 N <- 200 x <- seq(-pi, pi, l = Mx + 1)[-c(Mx + 1)] times <- seq(0, 1, l = N + 1) u0 <- dWn1D(x, pi/2, 0.05) b <- driftWn1D(x, alpha = 1, mu = pi, sigma = 1) sigma2 <- rep(1, Mx) # Full trajectory of the solution (including initial condition) u <- crankNicolson1D(u0 = cbind(u0), b = b, sigma2 = sigma2, N = 0:N, deltat = 1 / N, Mx = Mx, deltax = 2 * pi / Mx) # Mass conservation colMeans(u) * 2 * pi # Visualization of tpd plotSurface2D(times, x, z = t(u), levels = seq(0, 3, l = 50)) # Only final time v <- crankNicolson1D(u0 = cbind(u0), b = b, sigma2 = sigma2, N = N, deltat = 1 / N, Mx = Mx, deltax = 2 * pi / Mx) sum(abs(u[, N + 1] - v)) ```

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