crankNicolson2D {sdetorus} | R Documentation |
Crank–Nicolson finite difference scheme for the 2D 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 toroidal diffusion
Usage
crankNicolson2D(u0, bx, by, sigma2x, sigma2y, sigmaxy, N, deltat, Mx, deltax,
My, deltay, imposePositive = 0L)
Arguments
u0 |
matrix of size |
bx , by |
matrices of size |
sigma2x , sigma2y , sigmaxy |
matrices of size |
N |
increasing integer vector of length |
deltat |
time step. |
Mx , My |
sizes of the equispaced spatial grids in |
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 |
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 * My, nt)
containing the discretized solution at the required times with thec(Mx, My)
matrix stored column-wise.If
nt == 1
, a matrix of sizec(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
# 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[1], 0.5), dWn1D(y, p[2], 0.5)))
bx <- outer(x, y, function(x, y) 5 * sin(m[1] - x))
by <- outer(x, y, function(x, y) 5 * sin(m[2] - 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[1], p[2], pch = 16)
points(m[1], m[2], pch = 16)
}, j = slider(0, N))
## End(Not run)