# DifferentialEquations.jl treat complex different than float

**URL:** <https://discourse.julialang.org/t/differentialequations-jl-treat-complex-different-than-float/78412>\
**Category:** Numerics\
**Tags:** float, differentialequation, complex-numbers\
**Created:** [March 24, 2022, 4:07pm UTC](https://discourse.julialang.org/t/differentialequations-jl-treat-complex-different-than-float/78412 "2022-03-24T16:07:36Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Elyco](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elyco/32/18418_2.png) [@Elyco](https://discourse.julialang.org/u/Elyco)\
**Post date:** [March 24, 2022, 4:07pm UTC](https://discourse.julialang.org/t/differentialequations-jl-treat-complex-different-than-float/78412/1 "2022-03-24T16:07:37Z")

</div>

Hi,

I encountered something weird when I was trying to solve a system of ODEs in two equivalent (if I don’t miss something very basic) ways.

Here is two ODE, that should be mathematically the same

```julia
C = -0.01+2im
function complex_func(du::Vector{ComplexF64}, u::Vector{ComplexF64}, p, t::Float64)
    du[1] = C*u[1]
end
function real_func(du::Vector{Float64},u::Vector{Float64},p,t::Float64)
    du[1] = real(C)*u[1]-imag(C)*u[2]
    du[2] = imag(C)*u[1]+real(C)*u[2]
end

```

Solving them should be identical. But to my surprise there is a difference between the solution that the solver (`Tsit5`) gives me. The difference is indeed small `~1e-7` for `tspan = (0.0, 1.0)` but I expected the model to be solved internally the same.  
I want to emphasize that I don’t care that they differ from the exact solution (which is easy to calculate), I care about the difference between each other.

```julia
julia> exact = u0c[1]exp(C*tf); abs(exact) # this just for understanding the scale
0.8935496450117586

julia> abs(sol_r[end][1] + 1im*sol_r[end][2] - exact[1]) / abs(exact)
2.00716203862083e-6 # great precision for me, I don't about this.

julia> abs((sol_r[end][1] + 1im*sol_r[end][2] - sol_c[end][1])/(sol_r[end][1] + 1im*sol_r[end][2] + sol_c[end][1]))
9.5661426125702e-8 # but why is that? shouldn't they just be the same?

```

Now, `1e-7` is not a big deal for me, but I have a different ODEs to solve that may be affected by that.  
In my work I have two equivalent ways and the difference of the solution is about 60%  
When I decrease `dtmax` the two ways converge to each other but it’s still weird for me.

Thanks!

> **More code for details**
>
> ```julia
> using DifferentialEquations
> C = -0.01+2im
> function complex_func(du::Vector{ComplexF64}, u::Vector{ComplexF64}, p, t::Float64)
> du[1] = C*u[1]
> end
> function real_func(du::Vector{Float64},u::Vector{Float64},p,t::Float64)
> du[1] = real(C)*u[1]-imag(C)*u[2]
> du[2] = imag(C)*u[1]+real(C)*u[2]
> end
> u_r = rand(2)
> u_c = u_r[1]+1im*u_r[2]
> du_c = [0.0im]
> complex_func(du_c, [u_c], 0.0, 0.0)
> du_r = [0.0, 0.0]
> real_func(du_r, u_r, 0.0, 0.0)
> du_r[1] + 1im*du_r[2] == du_c[1]
> 
> tf = 1.0
> tspan = (0.0, tf)
> u0r = rand(2)
> u0c = [u0r[1]+1im*u0r[2]]
> prob_c = ODEProblem(complex_func, u0c, tspan)
> prob_r = ODEProblem(real_func, u0r, tspan)
> 
> integrator_c = init(prob_c, Tsit5(), abstol = 1e-10)
> integrator_r = init(prob_r, Tsit5(), abstol = 1e-10)
> 
> sol_c = solve!(integrator_c)
> sol_r = solve!(integrator_r)
> 
> sol_r[end][1] + 1im*sol_r[end][2] - sol_c[end][1]
> (sol_r[end][1] + 1im*sol_r[end][2] - sol_c[end][1])/(sol_r[end][1] + 1im*sol_r[end][2] + sol_c[end][1])
> exact = u0c[1]exp(C*tf)
> # =#
> 
> exact = u0c[1]exp(C*tf); abs(exact)
> abs(sol_r[end][1] + 1im*sol_r[end][2] - exact[1]) / abs(exact)
> abs((sol_r[end][1] + 1im*sol_r[end][2] - sol_c[end][1])/(sol_r[end][1] + 1im*sol_r[end][2] + sol_c[end][1]))
> 
> ```

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [March 24, 2022, 4:16pm UTC](https://discourse.julialang.org/t/differentialequations-jl-treat-complex-different-than-float/78412/2 "2022-03-24T16:16:31Z")

</div>

I suspect the difference might result from the difference in real vs complex norms used calculate tolerances for adaptive integrations.

Or even `C*u[1]` using a different order of floating point operations than the unrolled version in the `real_func` case.

---

<div class="post-metadata">

**Author:** ![Elyco](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elyco/32/18418_2.png) [@Elyco](https://discourse.julialang.org/u/Elyco)\
**Post date:** [March 24, 2022, 7:55pm UTC](https://discourse.julialang.org/t/differentialequations-jl-treat-complex-different-than-float/78412/3 "2022-03-24T19:55:22Z")

</div>

I can understand what you saying about real vs complex norm, but do you thing that can be result in difference of ~50%?? 😯

About the order of floating point operations, I think that the error resulting from that should be much less than `1e-8`.

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [March 24, 2022, 8:52pm UTC](https://discourse.julialang.org/t/differentialequations-jl-treat-complex-different-than-float/78412/4 "2022-03-24T20:52:46Z")

</div>

Yeah, informed people (which I’m not but @oscar_smith is) would be probably counting ulp here. So if I were interested (not sure if I am), I would try this with the simplest solver available first.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [March 24, 2022, 8:56pm UTC](https://discourse.julialang.org/t/differentialequations-jl-treat-complex-different-than-float/78412/5 "2022-03-24T20:56:16Z")

</div>

This isn’t an ulp issue. The issue is that you are solving to a different tolerance between the real and complex version since the norm is different. This means that pretty much everything ends up happening differently.

---

<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:** [March 24, 2022, 10:23pm UTC](https://discourse.julialang.org/t/differentialequations-jl-treat-complex-different-than-float/78412/6 "2022-03-24T22:23:42Z")

</div>

Remember that solver tolerances are local, not global. So the true error is usually an order of magnitude or two more than the chosen tolerance. The relative tolerance is not changed here at all, so I’m not too surprised that the global difference is 1e-7ish.

---

<div class="post-metadata">

**Author:** ![Elyco](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elyco/32/18418_2.png) [@Elyco](https://discourse.julialang.org/u/Elyco)\
**Post date:** [March 25, 2022, 12:57pm UTC](https://discourse.julialang.org/t/differentialequations-jl-treat-complex-different-than-float/78412/7 "2022-03-25T12:57:15Z")

</div>

Thank you all for your explanations!

So I guess the question is now why in my case the error is ~50%, I hoped it’s okay to bring my code as is. The structure is pretty much the same as the original post.

…

I was trying to make a MWE of my original code but the difference between the two methods became significantly smaller. So I need to dig into it before continuing here.

---

<div class="post-metadata">

**Author:** ![Elyco](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elyco/32/18418_2.png) [@Elyco](https://discourse.julialang.org/u/Elyco)\
**Post date:** [March 28, 2022, 7:06am UTC](https://discourse.julialang.org/t/differentialequations-jl-treat-complex-different-than-float/78412/8 "2022-03-28T07:06:17Z")

</div>

Obviously, it was my [mistake](https://diffeq.sciml.ai/stable/basics/faq/#My-model-is-reporting-unstable-results.-What-can-I-do?) and the error is not at all that big. I hope that the discussion here will help someone sometime.