SDE `EnsembleProblem` unhappy with `trajectories>=2`?

Specific Domains Modelling & Simulations
,
1

The title is a bit generic, but I have trouble seeing a more specific version of my SDE problem. Consider the following problem:

using DifferentialEquations
using LinearAlgebra

function make_problem_2()
    m1 = rand(Complex{Float64}, 2, 2)
    m1 ./= norm(m1)
    m2 = rand(Complex{Float64}, 2, 2)
    m2 ./= norm(m2)

    u0 = rand(Complex{Float64}, 2, 2)

    W = WienerProcess(0.0,0.0)
    tspan = (0,10.)

    f = let
        m1 = m1
        f(du, u, p, t) = du .= m1 * u
    end
    g = let
        m2 = m2
        g(du, u, p, t) = du .= m2 * u
    end
    SDEProblem(f,g,u0,tspan,noise=W)
end
##
for i in 1:100
    prob = make_problem_2()
    eprob = EnsembleProblem(prob)
    @time esol = solve(eprob,
        #ImplicitRKMil(autodiff=false),
        #saveat=0.01,
        #dtmax=0.01,
        trajectories=2)
    @assert esol[1].retcode == :Success
end

If I set trajectories=1 it works fine. If I set it to >=2 then it breaks with

┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase ~/.julia/packages/SciMLBase/YVSzo/src/integrator_interface.jl:345
  0.026924 seconds (67.83 k allocations: 4.377 MiB)
┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/YVSzo/src/integrator_interface.jl:351
┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/YVSzo/src/integrator_interface.jl:351
  0.007925 seconds (26.74 k allocations: 1.684 MiB)
ERROR: AssertionError: (esol[1]).retcode == :Success

Why is the instability showing up only with more than one trajectories? Seen both in 1.7.2 and in 1.9.0-DEV.700 Commit 53338ca3424 (2022-06-01 14:43 UTC).

] st gives

  [0c46a032] DifferentialEquations v7.1.0
  [f6369f11] ForwardDiff v0.10.30
  [91a5bcdd] Plots v1.29.0
  [6e0679c1] QuantumOptics v1.0.4
2

You’re reusing the same noise problem with both trajectories, so it’s not thread-safe. You want a noise per problem.

3

Thanks! Could you provide the idiomatic way to create a new noise per problem? I have a few silly ways to do it, but I would like to make a PR to add the correct one to the documentation.

Currently the documentation makes you think this is automatically done because deepcopy is supposed to be used by default due to safetycopy=true: Parallel Ensemble Simulations · DifferentialEquations.jl

Similarly, the tutorial makes it look as if nothing special is needed to copy a new noise instance: Stochastic Differential Equations · DifferentialEquations.jl

4

Make it as part of the prob func and add noise=.to the created prob there

5

Minor documentation update posted as a PR here EnsembleProblem safetycopy clarification by Krastanov · Pull Request #561 · SciML/DiffEqDocs.jl · GitHub

1 Like