I am trying to solve a 4D chemical Langevin SDE problem with multiplicative noise in the form of independent wiener processes (WPs). Since my system is built upon 4 states and 8 reactions, my drift function F returns a (4x1) Vector and my noise function G a 4x8 matrix. The SDE is dU = f(u,p,t)dt + g(u,p,t)dW
using DifferentialEquations, LinearAlgebra
const stoich_matrix = [1 0 0 0 -1 0 0 0;
0 1 0 0 0 -1 0 0;
0 0 1 0 0 0 -1 0;
0 0 0 1 0 0 0 -1]
function hill_fn(var, theta, eta)
return (1+(var/theta)^eta)^(-1)
end
function propensity_fn(u::Vector, i; atcmax=100, iptgmax=1)
prop_vec = [3.2e-2 + 8.3*hill_fn(u[4]*hill_fn(atcmax*i, 11.65, 2.0), 30.0, 2.00),
1.19e-1 + 2.06*hill_fn(u[3]*hill_fn(iptgmax*(1-i), 9.06e-2, 2.00), 31.94, 2.00),
9.726e-1*u[1], 1.17*u[2], 1.386e-1*u[1], 1.386e-1*u[2], 1.65e-2*u[3], 1.65e-2*u[4]]
return prop_vec
function drift_fn(u, p, t)
controller_input = p(u, t)
propensities = propensity_fn(u[1:4], controller_input)
du = stoich_matrix * propensities
return du
end
function g(u, p, t)
controller_input = p(u, t)
propensities = propensity_fn(u[1:4], controller_input)
diag = Diagonal(sqrt.(propensities))
du = stoich_matrix * diag
return du
end
I’m almost certain this interpretation of the SDE is correct as it has been referenced from a number of papers and researchers. Since by noise term is a 4x8 matrix it must be multiplied by an 8 dimensional noise vector (the 8 WPs) so I created NoiseGrid of 8d vectors to be used as the noise argument in SDEProblem.
u0 = [2.55166047363230, 38.7108543679906, 102.155003051775, 1196.05604522200]
dt = 0.05
tspan = (0.0, 10.0)
ts = Array(tspan[1]:dt:(tspan[2]+dt)) # time array for noise grid
W = sqrt(dt)*randn(length(ts), 8)
Wsum = cumsum(vcat(zeros(1, 8), W[1:end-1, :]), dims = 1)
W_vec = [vec(Wsum[el, :]) for el in 1:1:size(Wsum, 1)]
noisegrid = NoiseGrid(ts, W_vec)
sde = SDEProblem(drift_fn, g, u0, tspan, (x,t)->0, noise = noisegrid)
sol = solve(sde, SOSRA())
Dependent on solver choice I get one of two stacktraces. I believe each is respectively caused by using an Ito or Stratonovich method.
MethodError: no method matching /(::Nothing, ::Float64)
Closest candidates are:
/(::Any, !Matched::ChainRulesCore.AbstractThunk) at C:\Users\ollie\.julia\packages\ChainRulesCore\16PWJ\src\tangent_types\thunks.jl:33
/(!Matched::StridedArray{P}, ::Real) where P<:Dates.Period at C:\Users\ollie\AppData\Local\Programs\Julia-1.7.2\share\julia\stdlib\v1.7\Dates\src\deprecated.jl:44
/(!Matched::Union{SparseArrays.SparseVector{Tv, Ti}, SubArray{Tv, 1, <:SparseArrays.AbstractSparseVector{Tv, Ti}, Tuple{Base.Slice{Base.OneTo{Int64}}}, false}, SubArray{Tv, 1, <:SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, false}} where {Tv, Ti}, ::Number) at C:\Users\ollie\AppData\Local\Programs\Julia-1.7.2\share\julia\stdlib\v1.7\SparseArrays\src\sparsevector.jl:1476
in eval at base\boot.jl:373
in top-level scope at discourse.jl:65
in solve at DiffEqBase\RHAWf\src\solve.jl:747
in #solve#33 at DiffEqBase\RHAWf\src\solve.jl:752
in solve_up at DiffEqBase\RHAWf\src\solve.jl:756
in var"#solve_up#34" at DiffEqBase\RHAWf\src\solve.jl:767
in solve_call at DiffEqBase\RHAWf\src\solve.jl:401
in #solve_call#28 at DiffEqBase\RHAWf\src\solve.jl:429
in __solve at StochasticDiffEq\KQyGf\src\solve.jl:6
in __solve at StochasticDiffEq\KQyGf\src\solve.jl:6
in __solve at StochasticDiffEq\KQyGf\src\solve.jl:6
in __solve at StochasticDiffEq\KQyGf\src\solve.jl:6
in __solve at StochasticDiffEq\KQyGf\src\solve.jl:6
in #__solve#110 at StochasticDiffEq\KQyGf\src\solve.jl:7
in solve! at StochasticDiffEq\KQyGf\src\solve.jl:611
in perform_step! at StochasticDiffEq\KQyGf\src\perform_step\sra.jl:176
DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 4 and 8")
in eval at base\boot.jl:373
in top-level scope at discourse.jl:66
in solve at DiffEqBase\RHAWf\src\solve.jl:747
in #solve#33 at DiffEqBase\RHAWf\src\solve.jl:752
in solve_up at DiffEqBase\RHAWf\src\solve.jl:756
in #solve_up#34 at DiffEqBase\RHAWf\src\solve.jl:771
in solve_call at DiffEqBase\RHAWf\src\solve.jl:401
in #solve_call#28 at DiffEqBase\RHAWf\src\solve.jl:429
in __solve at DiffEqBase\RHAWf\src\solve.jl:997
in #__solve#52 at DiffEqBase\RHAWf\src\solve.jl:1002
in __solve##kw at DifferentialEquations\TZ7Qg\src\default_solve.jl:4
in var"#__solve#1" at DifferentialEquations\TZ7Qg\src\default_solve.jl:9
in at StochasticDiffEq\KQyGf\src\solve.jl:6
in __solve##kw at StochasticDiffEq\KQyGf\src\solve.jl:6
in __solve##kw at StochasticDiffEq\KQyGf\src\solve.jl:6
in __solve##kw at StochasticDiffEq\KQyGf\src\solve.jl:6
in __solve##kw at StochasticDiffEq\KQyGf\src\solve.jl:6
in #__solve#110 at StochasticDiffEq\KQyGf\src\solve.jl:7
in solve! at StochasticDiffEq\KQyGf\src\solve.jl:611
in perform_step! at StochasticDiffEq\KQyGf\src\perform_step\lamba.jl:16
in materialize at base\broadcast.jl:860
in instantiate at base\broadcast.jl:281
in combine_axes at base\broadcast.jl:499
in broadcast_shape at base\broadcast.jl:504
in _bcs at base\broadcast.jl:510
in _bcs1 at base\broadcast.jl:516
I believe the Dimension Mismatch is the more promising of the errors but I’m unsure.
I have also been able to reconstruct my problem as a RODE but it isn’t what I’m after and is also temperamental when it comes to non-negativity constraints. I apologise for the length of the post, couldn’t think of any more concise way to write it. 