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?