I have a function here:
\frac{dx}{dt} = kx(N-x) - lx(1-\alpha\frac{x}{N}) + \sigma k x \xi(t)
\xi(t) is Random perturbation term, \sigma is the variance of it,
k,l,N, \alpha is a const.
using Plots
using DifferentialEquations
k, l, N, α = 0.02, 0.01, 1.0, 0.2
f(u,p,t)=k*u*(N-u) - l*u * (1 - α * u/N)
u₀ = 0.001
tspan = (0.0,100.0)
prob = ODEProblem(f,u₀,tspan)
res = solve(prob,Tsit5(), reltol=1e-8, abstol=1e-8)
plot(res)
how I add the \xi(t) in this function and solve?
That’s actually a stochastic differential equation and not an ordinary differential equation. Take a look at Stochastic Differential Equations · DifferentialEquations.jl for an example of how to solve them with DifferentialEquations.jl.
1 Like
As @dawbarton said, its an SDE–basically, you define a second function for the stochastic term and construct an SDEProblem
instead of an ODEProblem
. If I’m understanding your equation correctly it should look like this:
f(u, p, t) = k*u*(N-u) - l*u * (1 - α * u/N)
g(u, p, t) = k * σ * u
prob = SDEProblem(f, g, u₀, tspan)
res = solve(prob)
As a side note, your equations will solve faster if you avoid global variables by passing in your parameters as a named tuple:
params = (k=0.02, l=0.01, N=1.0, α=0.2, σ=1.0)
f(u, p, t) = p.k*u*(p.N-u) - p.l*u * (1 - p.α * u/p.N)
g(u, p, t) = p.k * p.σ * u
prob = SDEProblem(f, g, u₀, tspan, params)
@time res = solve(prob)
3 Likes
em, so g(u, p, t) = k*σ*u
is σkx
, dW = \frac{d(\xi(t))}{dt}, is that right?
if I want change the variance of (\xi(t)), how to do that?
In DifferentialEquations.jl, the full equation is du = f(u,p,t)dt + g(u,p,t)dW
, such that g(u,p,t)
is the prefactor of the Wiener process dW
as you correctly say. This prefactor directly sets the variance of the Wiener process, in your case this is \sigma k x according to your first equation. So if you say that for your problem \sigma is the variance of \xi than you should simply change \sigma to change the variance.
However, formally the Wiener process is not differentiable such that dW=\frac{d\xi}{dt} is not 100% correct but heuristically okay, see also the wikipedia site for a nice explanation.