# Get error when solving SDEs after last update

**URL:** <https://discourse.julialang.org/t/get-error-when-solving-sdes-after-last-update/10137>\
**Category:** Modelling & Simulations\
**Tags:** diffeq, sde\
**Created:** [April 3, 2018, 12:53pm UTC](https://discourse.julialang.org/t/get-error-when-solving-sdes-after-last-update/10137 "2018-04-03T12:53:15Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Torkel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torkel/32/5030_2.png) [@Torkel](https://discourse.julialang.org/u/Torkel)\
**Post date:** [April 3, 2018, 12:53pm UTC](https://discourse.julialang.org/t/get-error-when-solving-sdes-after-last-update/10137/1 "2018-04-03T12:53:15Z")

</div>

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

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

```

more detailed (bu not all, quite long):

```julia
[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:

```julia
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:

```julia
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:

```julia
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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 3, 2018, 7:09pm UTC](https://discourse.julialang.org/t/get-error-when-solving-sdes-after-last-update/10137/2 "2018-04-03T19:09:42Z")

</div>

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).

---

<div class="post-metadata">

**Author:** ![Torkel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torkel/32/5030_2.png) [@Torkel](https://discourse.julialang.org/u/Torkel)\
**Post date:** [April 4, 2018, 10:40am UTC](https://discourse.julialang.org/t/get-error-when-solving-sdes-after-last-update/10137/3 "2018-04-04T10:40:40Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 4, 2018, 1:25pm UTC](https://discourse.julialang.org/t/get-error-when-solving-sdes-after-last-update/10137/4 "2018-04-04T13:25:54Z")

</div>

> [@Torkel](#):
>
> I will start looking at this and try to figure out what actually is happening.

Nothing about the actual algorithm changed BTW other than adaptivity. So my “guess” of what’s happening (which is a correct guess 😉) 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:

> <https://github.com/SciML/StochasticDiffEq.jl/issues/62>
>
> I tend to think that \`beta1\` and \`beta2\` are off. I think that \`beta1=0.4\` and \`…beta2=0.1\` may be better defaults than what we have right now, but this would need extensive testing. That's a very high gain but right now we rarely ever reject, and are usually 10-100x below the error estimate we are aiming for.

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.

> <https://github.com/SciML/SciMLBenchmarks.jl/issues/17>
>
> @onoderat has moved us forward finally in the Non-diagonal SDE department, and s…o we should think about creating some benchmark notebooks in this area.

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!)