Get error when solving SDEs after last update

Specific Domains Modelling & Simulations
,
1

I tried rerunning my code after the latest 4.2 update to the DifferentialEquations package. I use the reaction_network DSL from DiffEqBiological package, from these I generate and solve SDEs with non-diagonal noise. However I generated an error

DimensionMismatch("array could not be broadcast to match destination")

more detailed (bu not all, quite long):

[1] check_broadcast_shape(::Tuple{Base.OneTo{Int64}}, ::Tuple{Base.OneTo{Int64}}) at ./broadcast.jl:83
 [2] check_broadcast_indices at ./broadcast.jl:89 [inlined]
 [3] broadcast_c! at ./broadcast.jl:210 [inlined]
 [4] broadcast! at ./broadcast.jl:206 [inlined]
 [5] perform_step!(::StochasticDiffEq.SDEIntegrator{StochasticDiffEq.ImplicitEM{0,true

I have been searching some. Trying to figure out what it is, these are my results.
This works fine:

model = @reaction_network rn begin
  (1,1), X ↔ Y
end
prob_sde = SDEProblem(model,[200.,300.],(0.,10.));
sol = solve(prob_sde,ImplicitEM());

while this is not:

model = @reaction_network rn beginbegin
    (1,1), X ↔ Y
    (10,0.1), ∅ ↔ X
end
prob_sde = SDEProblem(model,[200.,300.],(0.,10.));
sol = solve(prob_sde,ImplicitEM());

Solving the corresponding ODEs is no problem.

Next I try implementing it normally:

function f(du,u,p,t) 
    du[1] = u[2] - u[1] + 10 - 0.1*u[1]
    du[2] = u[1] - u[2]
end
function g(du,u,p,t)
  du[1,1] = -sqrt(u[1])
  du[1,2] = sqrt(u[2])
  du[1,3] = sqrt(10)
  du[1,4] = -sqrt(0.1*u[1])
  du[2,1] = sqrt(u[1])
  du[2,2] = -sqrt(u[2])
  du[2,3] = 0
  du[2,4] = 0
end
prob_sde = SDEProblem(f,g,[200.,300.],(0.,10.),noise_rate_prototype=zeros(2,4))
sol = solve(prob_sde,ImplicitEM());
plot(sol)

This does not work. I also get the same error when I copy from the documentation:

f(du,u,p,t) = du .= 1.01u
function g(du,u,p,t)
  du[1,1] = 0.3u[1]
  du[1,2] = 0.6u[1]
  du[1,3] = 0.9u[1]
  du[1,4] = 0.12u[2]
  du[2,1] = 1.2u[1]
  du[2,2] = 0.2u[2]
  du[2,3] = 0.3u[2]
  du[2,4] = 1.8u[2]
end
prob = SDEProblem(f,g,ones(2),(0.0,1.0),noise_rate_prototype=zeros(2,4))
sol = solve(prob_sde,ImplicitEM());
plot(sol)

I am not sure whenever this is an issue with the package and should be raised there, or if it is me doing things in the wrong way, in which it might be good to keep it here for others to see properly.

2

That was a bug in the new adaptivity algorithm with rectangular non-diagonal noise. Pkg.update() to get 4.2.1 and try it now. (BTW, I’d recommend ISSEM() in place of ImplicitEM(), at least in theory though it needs benchmarks).

3

OK, works now.

However I experience a slow down with factor 30 when solving with ISSEM or ImplicitEM(). This increase to a factor 300(!) when including the actually interesting things of the simulations (I have a signal which is turned on and activates the system. If I only simulate at equilibrium until the signal the slowdown is a factor 30, if I include time after the signal is turned on slowdown is increases to a factor 300). This all compares to solving the same system using ImplicitEM() before the update.

I will start looking at this and try to figure out what actually is happening.

4

Nothing about the actual algorithm changed BTW other than adaptivity. So my “guess” of what’s happening (which is a correct guess :wink:) is that it’s taking more steps because it’s adaptive and solving to a much lower tolerance than the dt you chose before, and also doing so more conservatively.

We can fiddle with the conservativeness, and I have an issue to discuss that:

You can set adaptive=false to get the old behavior, or raise the tolerance. But I think that the main issue is that you probably used a dt which gave you a lot higher error than you were expecting (SDEs default to abstol=reltol=1e-2, so “two digits correct” but is conservative so more like 3 or 4). Remember, ImplicitEM converges at a rate of 0.5 and has higher error in places where the derivatives of f and g are high, and the adaptivity is just correcting for these effects.

But try adding beta1=0.4,beta2=0.1 and see if you like that result better. To know how it’s really doing, we’ll need to investigate your SDE with DiffEqDevTools to really estimate the real error.

My guess is that before you were solving it to about 1e-1 error (in terms of strong error), so beta1=0.4,beta2=0.1,abstol=1e-1,reltol=1e-1 would probably be the correct options. But without actually estimating the true solution talking about the global error is really hard, so I invite you to share the problem so we can both learn more about what’s going on (and if you don’t want to share it publicly until after a publication, that’s fine too!)