crankNicolson1D {sdetorus} | R Documentation |
Crank–Nicolson finite difference scheme for the 1D Fokker–Planck equation with periodic boundaries
Description
Implementation of the Crank–Nicolson scheme for solving the Fokker–Planck equation
p(x, t)_t = -(p(x, t) b(x))_x + \frac{1}{2}(\sigma^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+\sigma(X_t)dW_t
.
Usage
crankNicolson1D(u0, b, sigma2, N, deltat, Mx, deltax, imposePositive = 0L)
Arguments
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 |
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 sizec(Mx, nt)
containing the discretized solution at the required times.If
nt == 1
, a matrix of sizec(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
# 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))