SDE simulation error: ERROR: BoundsError: attempt to access 3-element Vector{Float64} at index [1, 2]

Specific Domains Modelling & Simulations
1

I am getting the following error trying to run a simulation
ERROR: BoundsError: attempt to access 3-element Vector{Float64} at index [1, 2]

My code is as follows
using DifferentialEquations
using Plots; plotly()

function β(t)
    0.2*sin(t/200)+0.5
end
function Γ(t)
    0.1*cos(t/200)+0.85
end
function Σ1(t)
    2+0.1*cos(t)
end
function Σ2(t)
    2+0.1*sin(t)
end
function φ12(t)
    0.1+0.001*cos(t)
end
function  φ3(t)
    0.2+0.002*cos(t)
end
function χ(t)
    1.2+0.05*sin(t)
end

function f!(du,u,p,t)
    du[1] = -β(t)*u[1]*u[2]
    du[2] = β(t)*u[1]*u[2]-Γ(t)*u[2]
    du[3] = Γ(t)*u[2]
end

function g!(du,u,p,t)
    du[1,1] = -u[1]*u[2]
    du[1,2] = 0.0
    du[2,1] = u[1]*u[2]
    du[2,2] = u[2]*u[3] 
    du[3,1] = 0.0
    du[3,2] = -u[2]*u[3]
end

u0 = [0.75, 0.15, 0.1]
dt= 1/2^15
W = WienerProcess(0.0,0.0)

prob = SDEProblem(f!, g!, u0,(0.0,1.0), noise = W)
sol = solve(prob, EM(), dt=dt)

plot(sol)

This was previously working fine. If I change g! to

du[1] = ...
du[2] = ..
du[3] = ..

It will run but that isn’t the correct form, I want non-diagonal noise.

EDIT1: If I add noise_rate_prototype with a value it still doesn’t work.
EDIT2: I copied the example here Stochastic Differential Equations · DifferentialEquations.jl And, I then altered it to look identical to what I want and it compiled no issue so I am very lost.

2

Unrelated to the question, but maybe have a look at this article about quoting code.

3

Thanks for pointing that out, I have edited.

1 Like
4

You don’t need to explicitly construct the noise process as the default is to use Wiener processes. You can just do:

prob = SDEProblem(f!, g!, u0, (0.0,1.0); noise_rate_prototype=zeros(3,2))
sol = solve(prob, EM(), dt=dt)
plot(sol)
5

And if you really want to explicitly construct the process even though it isn’t necessary I think this works:

# note the initial value is a vector of length two since you need a
# two-dimensional Brownian motion
W = WienerProcess(0.0, [0.0,0.0])   
prob = SDEProblem(f!, g!, u0, (0.0,1.0); noise = W, noise_rate_prototype=zeros(3,2))
6

Thank you! And I just checked, by placing

W = WienerProcess(0.0, [0.0,0.0])

it does work.

7

Great, glad that works!