julia – MethodError: no method matching setindex! when using DifferentialEquations.jl

I am trying to simulate a multi-dimensional SDE with non-diagonal noise using the DifferentialEquations package in julia. The drift function has the following form:

function mu(dx, x, param, t)
    dx[1] = -param[1]*x[1]*x[2]-param[1]*(x[3]+x[4])*x[1]
    dx[2] = param[1]*x[1]*x[2]-param[3]*x[2]
    dx[3] = param[1]*x[1]*(x[3]+x[4])-param[4]*x[3]
    dx[4] = param[2]*x[5]*(x[3]+x[4])-param[4]*x[4]
    dx[5] = -param[2]*x[5]*(x[3]+x[4])+param[3]*x[2]
    dx[6] = param[4]*(x[3]+x[4])
end

and the diffusion matrix is ​​given by:

function sigma(dx, x, p, t)
    beta = p[1]
    kappa = p[2]
    gamma = p[3]
    zeta = p[4]

    dx[1][1] = -sqrt(beta*x[1]*x[2]/100)
    dx[1][2] = -sqrt(beta*x[1]*(x[3]+x[4])/100)
    dx[2][1] = sqrt(beta*x[1]*x[2]/100)
    dx[2][3] = -sqrt(gamma*x[2]/100)
    dx[3][2] = sqrt(beta*x[1]*(x[3]+x[4])/100)
    dx[3][6] = -sqrt(zeta*x[3]/100)
    dx[4][4] = sqrt(kappa*x[5]*(x[3]+x[4])/100)
    dx[4][5] = -sqrt(zeta*x[4]/100)
    dx[5][4] = sqrt(gamma*x[2]/100)
    dx[5][4] = -sqrt(kappa*x[5]*(x[3]+x[4])/100)
    dx[6][5] = sqrt(zeta*x[4]/100)
    dx[6][6] = sqrt(zeta*x[3]/100)
end

Simulating the ODE specified by the drift function works fine with the following code.

u0 = [0.96; 0.03; 0.01; 0.0; 0.0; 0.0]
tspan = (0.0, 300.0)
p = [0.08; 0.06; 0.04; 0.02]


problem_ODE = ODEProblem(mu, u0, tspan, p)

solution_ODE = solve(problem_ODE, saveat=0.1)

But when I want to solve the SDE:

u0 = [0.96; 0.03; 0.01; 0.0; 0.0; 0.0]
tspan = (0.0, 300.0)
p = [0.08; 0.06; 0.04; 0.02]


problem = SDEProblem(mu, sigma, u0, tspan, p, noise_rate_prototype=zeros(6,6))

solution = solve(problem, EM(), dt=0.01, adaptive=false)

It does not work and gives the following error and Stacktrace:

MethodError: no method matching setindex!(::Float64, ::Float64, ::Int64)

Stacktrace:
  [1] sigma(dx::Matrix{Float64}, x::Vector{Float64}, p::Vector{Float64}, t::Float64)
    @ Main ~/masterthesis/two_variant_model/own_simulation/final_simulation.ipynb:7
  [2] perform_step!(integrator::StochasticDiffEq.SDEIntegrator{EM{true}, true, Vector{Float64}, Float64, Float64, Float64, Vector{Float64}, Float64, Float64, Float64, NoiseProcess{Float64, 2, Float64, Vector{Float64}, Nothing, Nothing, typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE), true, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, true}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, true}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, Nothing, Vector{Float64}, RODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, NoiseProcess{Float64, 2, Float64, Vector{Float64}, Nothing, Nothing, typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE), true, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, true}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, true}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, Nothing, SDEFunction{true, typeof(mu), typeof(sigma), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, typeof(sigma), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Matrix{Float64}}, EM{true}, StochasticDiffEq.LinearInterpolationData{Vector{Vector{Float64}}, Vector{Float64}}, DiffEqBase.DEStats}, StochasticDiffEq.EMCache{Vector{Float64}, Vector{Float64}, Matrix{Float64}}, SDEFunction{true, typeof(mu), typeof(sigma), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, typeof(sigma), Nothing, StochasticDiffEq.SDEOptions{Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Float64, Float64, Float64, Tuple{}, Tuple{}, Tuple{}}, Nothing, Float64, Nothing, Nothing}, cache::StochasticDiffEq.EMCache{Vector{Float64}, Vector{Float64}, Matrix{Float64}})
    @ StochasticDiffEq ~/.julia/packages/StochasticDiffEq/Wl3hr/src/perform_step/low_order.jl:40
  [3] solve!(integrator::StochasticDiffEq.SDEIntegrator{EM{true}, true, Vector{Float64}, Float64, Float64, Float64, Vector{Float64}, Float64, Float64, Float64, NoiseProcess{Float64, 2, Float64, Vector{Float64}, Nothing, Nothing, typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE), true, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, true}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, true}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, Nothing, Vector{Float64}, RODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, NoiseProcess{Float64, 2, Float64, Vector{Float64}, Nothing, Nothing, typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE), true, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, true}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, true}, RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, Nothing, SDEFunction{true, typeof(mu), typeof(sigma), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, typeof(sigma), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Matrix{Float64}}, EM{true}, StochasticDiffEq.LinearInterpolationData{Vector{Vector{Float64}}, Vector{Float64}}, DiffEqBase.DEStats}, StochasticDiffEq.EMCache{Vector{Float64}, Vector{Float64}, Matrix{Float64}}, SDEFunction{true, typeof(mu), typeof(sigma), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, typeof(sigma), Nothing, StochasticDiffEq.SDEOptions{Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Float64, Float64, Float64, Tuple{}, Tuple{}, Tuple{}}, Nothing, Float64, Nothing, Nothing})
    @ StochasticDiffEq ~/.julia/packages/StochasticDiffEq/Wl3hr/src/solve.jl:611
  [4] #__solve#100
    @ ~/.julia/packages/StochasticDiffEq/Wl3hr/src/solve.jl:7 [inlined]
  [5] #solve_call#39
    @ ~/.julia/packages/DiffEqBase/U3LtB/src/solve.jl:155 [inlined]
  [6] #solve_up#41
    @ ~/.julia/packages/DiffEqBase/U3LtB/src/solve.jl:182 [inlined]
  [7] #solve#40
    @ ~/.julia/packages/DiffEqBase/U3LtB/src/solve.jl:168 [inlined]
  [8] top-level scope
    @ ~/masterthesis/two_variant_model/own_simulation/final_simulation.ipynb:1
  [9] eval
    @ ./boot.jl:360 [inlined]
 [10] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1116
 [11] #invokelatest#2
    @ ./essentials.jl:708 [inlined]
 [12] invokelatest
    @ ./essentials.jl:706 [inlined]
 [13] (::VSCodeServer.var"#160#161"{VSCodeServer.NotebookRunCellArguments, String})()
    @ VSCodeServer ~/.vscode-server/extensions/julialang.language-julia-1.6.17/scripts/packages/VSCodeServer/src/serve_notebook.jl:19
 [14] withpath(f::VSCodeServer.var"#160#161"{VSCodeServer.NotebookRunCellArguments, String}, path::String)
    @ VSCodeServer ~/.vscode-server/extensions/julialang.language-julia-1.6.17/scripts/packages/VSCodeServer/src/repl.jl:184
 [15] notebook_runcell_request(conn::VSCodeServer.JSONRPC.JSONRPCEndpoint{Base.PipeEndpoint, Base.PipeEndpoint}, params::VSCodeServer.NotebookRunCellArguments)
    @ VSCodeServer ~/.vscode-server/extensions/julialang.language-julia-1.6.17/scripts/packages/VSCodeServer/src/serve_notebook.jl:13
 [16] dispatch_msg(x::VSCodeServer.JSONRPC.JSONRPCEndpoint{Base.PipeEndpoint, Base.PipeEndpoint}, dispatcher::VSCodeServer.JSONRPC.MsgDispatcher, msg::Dict{String, Any})
    @ VSCodeServer.JSONRPC ~/.vscode-server/extensions/julialang.language-julia-1.6.17/scripts/packages/JSONRPC/src/typed.jl:67
 [17] serve_notebook(pipename::String, outputchannel_logger::Base.CoreLogging.SimpleLogger; crashreporting_pipename::String)
    @ VSCodeServer ~/.vscode-server/extensions/julialang.language-julia-1.6.17/scripts/packages/VSCodeServer/src/serve_notebook.jl:136
 [18] top-level scope
    @ ~/.vscode-server/extensions/julialang.language-julia-1.6.17/scripts/notebook/notebook.jl:32
 [19] include(mod::Module, _path::String)
    @ Base ./Base.jl:384
 [20] exec_options(opts::Base.JLOptions)
    @ Base ./client.jl:285
 [21] _start()
    @ Base ./client.jl:485

Since this is my first time using julia, I would really appreciate any idea how to solve that issue and simulate the SDE.

Another question would be if anyone knows whether there is a way to tell the integration algorithm that the simulated process (the solution of the SDE) shall stay positive. My only idea there was to use the integrator interface to perform integration step by step and check it after every step, but maybe someone knows whether there isd a built-in solution for that.

But the first problem is much more important.

Thanks in advance.

Leave a Comment