fvade2d: Compute generator for a 2D advection-diffusion equation

In Uffe-H-Thygesen/SDEtools: Utilities for Stochastic Differential Equations

Compute generator for a 2D advection-diffusion equation

Description

fvade2d discretizes the advection-diffusion equation

dC/dt = -div ( u C - D grad C)

on a rectangular domain using the finite-volume method. Here, u=(ux,uy) and D=diag(Dx,Dy)

Usage

fvade2d(
  ux,
  uy,
  Dx,
  Dy,
  xgrid,
  ygrid,
  bc = list(N = "r", E = "r", S = "r", W = "r"),
  Dxy = function(x, y) 0
)

Arguments

ux	function mapping state (x,y) to advective term (numeric scalar)
uy	function mapping state (x,y) to advective term (numeric scalar)
Dx	function mapping state (x,y) to diffusivity (numeric scalar)
Dy	function mapping state (x,y) to diffusivity (numeric scalar)
xgrid	The numerical grid. Numeric vector of increasing values, giving cell boundaries
ygrid	The numerical grid. Numeric vector of increasing values, giving cell boundaries
bc	Specification of boundary conditions. See details.

Details

Boundary conditions: bc is a list with elements N,E,S,W. Each element is either "r": Reflective boundary "a": Absorbing boundary: Assume an absorbing boundary cell, which is not included "e": Extend to include an absorbing boundary cell "p": Periodic. When hitting this boundary, the state is immediately transferred to the opposite boundary, e.g. N->S.

Return value: The function fvade returns a generator (or sub-generator) G of a continuous-time Markov Chain. This chain jumps between cells defined by xgrid and ygrid. When using the generator to solve the Kolmogorov equations, note that G operates on probabilities of each cell, not on the probability density in each cell. The distinction is particularly important when the grid is non-uniform.

Value

a quadratic matrix, the generator of the approximating continuous-time Markov chain, with (length(xgrid)-1)*(length(ygrid)-1) columns

Examples

# Generator of a predator-prey model
xi <- seq(0,1.5,0.02)
yi <- seq(0,1.6,0.02)
xc <- 0.5*(utils::head(xi,-1)+utils::tail(xi,-1))
yc <- 0.5*(utils::head(yi,-1)+utils::tail(yi,-1))

ux <- function(x,y) x*(1-x)-y*x/(1+x)
uy <- function(x,y) y*x/(1+x)-y/3
D <- function(x,y) 0.01

G <- fvade2d(ux,uy,Dx=D,Dy=D,xi,yi)

phiv <- StationaryDistribution(G)
phim <- unpack.field(phiv,length(xc),length(yc))
image(xi,yi,t(phim))

Uffe-H-Thygesen/SDEtools documentation built on Jan. 15, 2025, 6:41 p.m.