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.