DifferentialEquations.jl treat complex different than float

Specific Domains Numerics
, ,
1

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

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

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.

3

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

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

1 Like
4

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.

5

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.

3 Likes
6

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.

1 Like
7

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.

8

Obviously, it was my mistake and the error is not at all that big. I hope that the discussion here will help someone sometime.