Description Usage Arguments Details Value References Examples

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*

.

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

`u0` |
matrix of size |

`b` |
vector of length |

`sigma2` |
vector of length |

`N` |
increasing integer vector of length |

`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 |

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).

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.

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

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))
``` |

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.