I’m a bit confused about what your code is doing here. What is the purpose of
It’s not used. And why is this here
Why do you write two different functions for ordinary floats and intervals instead of a single one that handles both?
I would have expected that
all(fun1(x, n, l) .∈ fun2(ix, n, l))
to be true, but it’s not.
If I write a single function, like this, I get something more understandable:
mypow(x, n::Integer) = x^n
mypow(x::Interval, n::Integer) = pow(x, n)
function fun0((x::AbstractVector, l::Integer)
F = diagm(0 => x)
@inbounds for j in eachindex(x) # dangerous to use 1:n where n is an input, since you have @inbounds
xj2 = x[j] * x[j]
Threads.@threads for i in eachindex(x)
xi2 = x[i] * x[i]
xixj = x[i] * x[j]
for k in 0:l-1
θ = 2π * k / l # I did not make this an interval. Why should it be an interval?
cosθ = cos(θ)
cos2θ = cos(2θ)
if i != j
d = 1 / mypow(sqrt(xi2 + xj2 - 2xixj * cosθ), 5)
F[i,i] += -(4xi2 + xj2 + xj2 * cos2θ) * d
F[i,j] += 4 * (xi2 + xj2) * cosθ * d
elseif k > 0
F[i, j] -= 1 / (sqrt(1 - cosθ) * mypow(x[i], 3))
Then call it with:
julia> x = float.(1:n); # make vector not matrix!
julia> xi = interval.(x); # much simpler notation (also, use vector)
julia> all(fun0(x, 20) .∈ fun0(xi, 20)) # this is what I would have expected
julia> @btime fun0(x, 20);
4.941 ms (68758 allocations: 1.03 MiB)
julia> @btime fun0(xi, 20); # pretty fast
6.383 ms (80461 allocations: 2.33 MiB)
Edit: Forgot to interpolate
@btime, but it makes no difference.
Edit2: For comparison with the old implementation for non-intervals:
julia> @btime fun1($x, 20, 20);
6.335 ms (80782 allocations: 1.22 MiB)