How to Simulation Continuous Time White Noise Autocorrelation Dirac Kronecker

White noise can be implemented as a "forcing function". A random component within the model function does not work, because it confuses the solver.

In the following examples, white noise is added at a fixed interval of 1. It can technically be implemented as interpolation with approxfun or just as a lookup in a vector. The lookup is faster, but I use approxfun as the more general approach. The function W can then be passed to the model function as additional named argument.

The non-continuous noise signal is a certain challenge for the solver, as it violates the continuity assumption of an ODE, but the lsoda solver is quite robust and can handle this well at the cost of small internal time steps.

For immediate impulses (that may be described as Dirac's delta), different approaches are possible.

Approach 1: Approximation by a smooth function (integral = 1)

Here an example, where we approximate delta by a Gaussian function (dnorm = normal density function). The bandwidth of the signal can be changed with parameter sigma (cf. Wikipedia). A rectangular signal would also be possible but a smooth differentiable signal is preferred.

          library("deSolve")  model <- function(t, state, parameters, W) {   with(as.list(c(state, parameters)), {          # approximate delta with a Gaussian that has integral 1     delta <- dnorm(t, sd = sigma)     dv <- -2 * alpha * v - beta ^ 2 * p - (kappa2 - kappa3) * delta +           (beta - 2 * alpha) * W(t)     dp <- v - kappa1 * delta + W(t)     # return the rate of change     list(c(dv, dp), delta = delta, W = W(t))   }) }  state <- c(v = 0.8, p = 0.5)  parameters <- c(alpha = 1,                 beta = 2,                 kappa1 = 1,                 kappa2 = 1,                 kappa3 =-1,                 sigma = 0.1, # width of the delta approximation (> 0)                 mean = 0     # location of delta, shift parameter  )   # start with - 100 to show the influence of delta at time t=0 times <- seq(-100, 100, by = 1)  ## predefined white noise in **discrete** time ## sd specifies amplitude of noise as standard deviation W <- approxfun(times, rnorm(length(times), sd = 0.1), method="constant")  out <- ode(y = state, times = times, func = model, parms = parameters, W = W) head(out)  par(oma = c(0, 0, 3, 0)) plot(out, xlab = "time", ylab = "-") mtext(outer = TRUE, side = 3, "ODE", cex = 1.5)                  

Approach 2: Use of an event function

The idea is to stop the solver at the time point(s) of delta, then to modify the states in an "event function", and then continue numerical integration. This can either be done by chaining calls of ode or with the event mechanism implemented in ode.

The following code shows an example with several events (impulses) occurring at time steps timedelta. A parameter delta is used for better visibility in the code, but as it is 1.0 by definition, it can also be omitted.

          library("deSolve")  model <- function(t, state, parameters, W) {   with(as.list(c(state, parameters)), {     dv <- -2 * alpha * v - beta ^ 2 * p  + (beta - 2 * alpha) * W(t)     dp <- v  + W(t)     list(c(dv, dp), W = W(t))   }) }  ## the event function returns instantly modified states eventfunc <- function(t, state, parameters, ...) {   with(as.list(c(state, parameters)), {     c(v - (kappa2 - kappa3) * delta,       p - kappa1 * delta     )   }) }  parameters <- c(alpha = 1, beta = 2, kappa1 = 1, kappa2 = 1, kappa3 =-1, delta = 1) state      <- c(v = 0.8, p = 0.5) times      <- seq(0, 100, by = 0.1) timedelta  <- seq(0, 50, by = 10) W <- approxfun(times, rnorm(length(times), sd = 0.1), method="constant")  out <- ode(y = state, times = times, func = model,    parms = parameters,    events = list(func = eventfunc, time = timedelta), W = W)  plot(out, xlab = "time", ylab = "-", mfrow=c(1, 3), which = c("W", "v", "p"))                  

fletcheraffix1945.blogspot.com

Source: https://stackoverflow.com/questions/73138215/solving-ode-with-white-noise-and-dirac-delta-function