crankNicolson1D: Crank-Nicolson finite difference scheme for the 1D...

In egarpor/sdetorus: Statistical Tools for Toroidal Diffusions

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 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 [-\pi,\pi).
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. \Sexpr[results=rd]{tools:::Rd_expr_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))

egarpor/sdetorus documentation built on March 4, 2024, 1:23 a.m.