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
where is the transition probability density of the circular diffusion
.
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))