inexacterror for SDEs with complex numbers

[WouterJRV]

When I implement example 4 from the SDE tutorial, i.e.

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

and solve, everything works fine. However, by adding a complex contribution, i.e. changing the first line to f(du,u,p,t) = du .= 1.01u+0.001im*u i get an error. One may object that this is because u is initialized as real, but even when replacing the initial ones(2) by [complex(1.0,0.0);complex(1.0,0.0)], even though the SDEproblem is defined as Juno.LazyTree(SDEProblem with uType Array{Complex{Float64},1} and tType Float64. In-place: true, getfield(Atom, Symbol("##35#36")){SDEProblem{Array{Complex{Float64},1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,typeof(f3),typeof(g3),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(g3),Nothing,Array{Float64,2}}}(SDEProblem with uType Array{Complex{Float64},1} and tType Float64. In-place: true timespan: (0.0, 1.0) u0: Complex{Float64}[1.0+0.0im, 1.0+0.0im])) ,

when calling solve I get

`InexactError: Float64(Float64, 0.30024598961268734 + 3.314181761768757e-11im) Type at complex.jl:37 [inlined] convert at number.jl:7 [inlined] setindex! at array.jl:771 [inlined] g3(::Array{Float64,2}, ::Array{Complex{Float64},1}, ::Nothing, ::Float64) at untitled-2b361671bdb2e6b89235821cf0b22d2b:40 perform_step!(::StochasticDiffEq.SDEIntegrator{LambaEM{true},true,Array{Complex{Float64},1},Complex{Float64},Float64,Nothing,Float64,Float64,Complex{Float64},NoiseProcess{Complex{Float64},2,Float64,Array{Complex{Float64},1},Nothing,Nothing,typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST),typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE),true,DataStructures.Stack{Tuple{Float64,Array{Complex{Float64},1},Nothing}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Complex{Float64},1},Nothing},true},RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},Array{Complex{Float64},1},RODESolution{Complex{Float64},2,Array{Array{Complex{Float64},1},1},Nothing,Nothing,Array{Float64,1},NoiseProcess{Complex{Float64},2,Float64,Array{Complex{Float64},1},Nothing,Nothing,typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST),typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE),true,DataStructures.Stack{Tuple{Float64,Array{Complex{Float64},1},Nothing}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Complex{Float64},1},Nothing},true},RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},SDEProblem{Array{Complex{Float64},1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,typeof(f3),typeof(g3),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(g3),Nothing,Array{Float64,2}},LambaEM{true},StochasticDiffEq.LinearInterpolationData{Array{Array{Complex{Float64},1},1},Array{Float64,1}}},StochasticDiffEq.LambaEMCache{Array{Complex{Float64},1},Array{Complex{Float64},1},Array{Float64,2},Array{Complex{Float64},1}},SDEFunction{true,typeof(f3),typeof(g3),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(g3),StochasticDiffEq.SDEOptions{Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Float64,Float64,Complex{Float64},Array{Float64,1},Array{Float64,1},Array{Float64,1}},Nothing,Complex{Float64}}, ::StochasticDiffEq.LambaEMCache{Array{Complex{Float64},1},Array{Complex{Float64},1},Array{Float64,2},Array{Complex{Float64},1}}, ::Function) at lamba.jl:64 perform_step! at lamba.jl:34 [inlined] solve!(::StochasticDiffEq.SDEIntegrator{LambaEM{true},true,Array{Complex{Float64},1},Complex{Float64},Float64,Nothing,Float64,Float64,Complex{Float64},NoiseProcess{Complex{Float64},2,Float64,Array{Complex{Float64},1},Nothing,Nothing,typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST),typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE),true,DataStructures.Stack{Tuple{Float64,Array{Complex{Float64},1},Nothing}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Complex{Float64},1},Nothing},true},RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},Array{Complex{Float64},1},RODESolution{Complex{Float64},2,Array{Array{Complex{Float64},1},1},Nothing,Nothing,Array{Float64,1},NoiseProcess{Complex{Float64},2,Float64,Array{Complex{Float64},1},Nothing,Nothing,typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST),typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE),true,DataStructures.Stack{Tuple{Float64,Array{Complex{Float64},1},Nothing}},ResettableStacks.ResettableStack{Tuple{Float64,Array{Complex{Float64},1},Nothing},true},RSWM{:RSwM3,Float64},RandomNumbers.Xorshifts.Xoroshiro128Plus},SDEProblem{Array{Complex{Float64},1},Tuple{Float64,Float64},true,Nothing,Nothing,SDEFunction{true,typeof(f3),typeof(g3),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(g3),Nothing,Array{Float64,2}},LambaEM{true},StochasticDiffEq.LinearInterpolationData{Array{Array{Complex{Float64},1},1},Array{Float64,1}}},StochasticDiffEq.LambaEMCache{Array{Complex{Float64},1},Array{Complex{Float64},1},Array{Float64,2},Array{Complex{Float64},1}},SDEFunction{true,typeof(f3),typeof(g3),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(g3),StochasticDiffEq.SDEOptions{Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Float64,Float64,Complex{Float64},Array{Float64,1},Array{Float64,1},Array{Float64,1}},Nothing,Complex{Float64}}) at solve.jl:392

__solve#41(::Base.Iterators.Pairs{Symbol,Bool,Tuple{Symbol},NamedTuple{(:de...`

Note that this is not a problem for all SDE problems, eg. adding an imaginary contribution to the Lorentz attractor example seemed to work without error.

[ChrisRackauckas]

ones(2). It looks like you used real numbers instead of complex numbers for your initial condition? You need to make your initial condition match the state type that you want.

[WouterJRV]

Hi, thanks for your answer. As I mentioned, I changed that initial condition to complex before without solving the issue. Now I just saw that I didn't define the noise_rate_prototype as complex, after doing this the solver runs without error. I'll still do some sanity checks on the solutions, but probably now it will work fine.

[ChrisRackauckas]

Okay cool. One thing to note is that multiplicative noise when complex is actually naturally commutative noise, so be careful but the higher order methods still provably converge just not at 1.5 (but sometimes at strong order 1.0).

[WouterJRV]

Okay cool. One thing to note is that multiplicative noise when complex is actually naturally commutative noise, so be careful but the higher order methods still provably converge just not at 1.5 (but sometimes at strong order 1.0).

The complex results seem to work fine to me :) .

Are you saying that using a higher order method intended for noise conditions that are not satisfied, causes losing some order of convergence; but are in general not "wrong" in the sense that they would give a bias to the solutions?

[ChrisRackauckas]

If you multiply two complex numbers in the noise term for multiplicative noise and then look at the real interpretation, you will notice it's not actually multiplicative noise anymore. If the g function is analytic and diagonal, it will still actually satisfy commutative noise properties though. This doesn't bias the solution because it all still converges, it just doesn't converge at a higher rate because the higher order terms of the stochastic Taylor series no longer cancel.

[WouterJRV]

This doesn't bias the solution because it all still converges, it just doesn't converge at a higher rate because the higher order terms of the stochastic Taylor series no longer cancel.

Is this generically true for methods for which the noise requirements are not satisfied? That would have some interesting implications, if e.g. noise is 'almost diagonal', it may be often good to try a diagonal method anyway without too much risk?

[ChrisRackauckas]

The diagonal methods won't allow it in other cases because they specialize on the size of the noise being the same for broadcast, but mathematically yes. Just like how with ODEs, a "bad method" accidentally converges with order 1.0 subject to same very basic constraints (sum(a_ij) = c_i and sum(b_i) = 1 IIRC), those same basic conditions work for order 0.5 in a stochastic RK method. However, you won't get cancellation of higher order terms so you do a bunch of extra work while still only being order 0.5, but maybe it still has less error by reducing the truncation error coefficients. In fact, the semi-implicit Trapezoid method is a common example of this. As an SRK method, it's only strong order 0.5, but it usually gives quite a significant constant decrease in the error over semi-implicit Euler-Maruyama, so that's why it's still generally used even though order-wise they are the same. In many cases, the size of the truncation error coefficients matters more than order, since dt is not infinitely small.

[WouterJRV]

Now I just saw that I didn't define the noise_rate_prototype as complex, after doing this the solver runs without error.

When I think about it, this was perhaps rather because using complex variables makes dW want to be complex. I guess the way to keep using Re[dW] and Im[dW] independently is still by manually redefining the noise process? Alternatively these can also be obtained by having access to both dW and dW* but probably this is harder? I think it would be handy if there were a standard option to choose whether noise is complex :)

[ChrisRackauckas]

http://docs.juliadiffeq.org/latest/features/noise_process.html#Wiener-Process-(White-Noise)-1

You can use this to make a real noise process for a complex-valued SDE.

[WouterJRV]

http://docs.juliadiffeq.org/latest/features/noise_process.html#Wiener-Process-(White-Noise)-1

You can use this to make a real noise process for a complex-valued SDE.

Changing wiener_randn by randn in WHITE_NOISE_DIST ?

[WouterJRV]

Naïvely, I would think that in a situation with multiple noise processes, where both the noise_rate_prototype and noise arguments are present, noise should be either a vector of noise-processes W=[WienerProcess(0.0,0.0);WienerProcess(0.0,0.0)] or the noise process with a vector input W=[WienerProcess(0.0,[0.0;0.0]),

but none of it works. If it would, real noise is simply obtained (if it not already becomes real from the real W0) by taking Abs(WienerProcess)