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

Question

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.

Answer 1

I got an answer to my problem, which was just me being annoyingly unaware of the different syntaxes to access elements in matrices in julia and python. Changing the diffusion function to:

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

solved the problem.

But then the second question got more important since I now get domain errors, when there appear some negative values. Solved this as well by using integrator intferace to check after each step. There is a built-in isoutofdomain option or one could use callback-function PositiveDomain, but this just ensures non-negativity up to an set tolerance, which still gives domain errors.