So, it’s time. I need to put some noise in my nice, clean, simple ODEs. Luckily DifferentialEquations.jl has me covered—the documentation is extensive… But perhaps too extensive: I don’t know where to start.

I *think*, ultimately, I’ll need to construct a Jump Diffusion problem, but need some baby steps and assistance before I get there.

My system is a coupled ODE with three equations. One of those contains the term r \frac{N}{N + H} B_t, where B_t is the dependent variable `u[1]`

, `r`

and `H`

are constants. `N`

is also a constant, but is altered currently via callback at discrete time points.

Ignoring that last point for the moment, I’d like to add some noise to `N`

. I can do this naively, changing `N::Float64`

to `N::Normal(n, sigma)`

, where `n`

is the original value of `N`

, then call `rand(N)`

within my ODE function.

I’ve also tried something like the scalar SDE example, where I set up the problem like:

```
function f(du, u, p, t)
r, N, H = p #N is a Float64 here.
du[1] = r*(N/(N+H))*u[1] + ...
du[2] = ...
du[3] = ...
end
function g(du, u, p, t)
r, N, H = p
du[1] = r*(N/(N+H))*u[1]
du[2] = 0.0
du[3] = 0.0
end
prob = SDEProblem(f,g,u0,tspan,p)
```

Here are the three results for `u[1]`

:

At the moment, I don’t need to be too rigorous about the form of this noise, but I like the look of the SDEProblem’s solution better. However, the interpretation of what I’ve done in `g`

is troubling me. Have I introduced a Weiner process with an amplitude of `r*(N/(N+H))*u[1]`

? If I’m just looking to add some stochasticity to `N`

, is this going about it the wrong way? If I set `du[1]`

in `g`

to just `N`

or `N*u[1]`

the results are non-nonsensical.

So, if I’ve got this part correct, is it a better solution than the `rand(N)`

inside the ODE? As alternatives I tried using Random ODEs and generating my own noise process, but I don’t think either of them make sense for my problem. Basically I’m getting bogged down in the options.

I think I’m on the right path with the SDEProblem though, just need some assistance with a plain language interpretation of what I did there.

Next: assuming I have that right, I’d like to add in discrete jumps via a callback. Say `N=0.7`

at one time point, 10 steps later a jump may change it to `N=1.2`

, then perhaps 20 steps after that, down to `N=0.9`

. I won’t know the size or direction of these jumps before time.

I’d like the Weiner process to take over in the steps between jumps however. This is why I suspect my system will ultimately become a Jump Diffusion problem. Based on what I’ve said here, do you agree that this is the correct path for my process, or am I way off the mark?