DiffEqDevTools.jl icon indicating copy to clipboard operation
DiffEqDevTools.jl copied to clipboard

Fixed L2 error for complex valued problems

Open henhen724 opened this issue 5 months ago • 7 comments

Currently, L2 error calculations in DiffEqDevTools.jl use sqrt(recursive_mean(vecvecapply((x) -> float(x) .^ 2, sol - timeseries_analytic))) the issue is if the solution is complex this returns the sum of the square of the complex differences instead of the sum of the absolute value squared (i.e. the L2 error should never be complex). For example,

wp = WorkPrecisionSet(prob,abstols,reltols,setups,1//2^(10);numruns=5,names=names,maxiters=1e7,error_estimate=:l2,appxsol_setup=Dict(:alg=>RKMilGeneral(;ii_approx=IICommutative())))
plot(wp)

leads to the error: isless(::ComplexF64, ::ComplexF64)

I corrected this by replacing x .^2 with conj.(x).*x: errors[:l2] = sqrt(recursive_mean(vecvecapply((x) -> conj.(float(x)) .* float(x), sol - timeseries_analytic)))

I checked using benchmark tools that this has approximately the same runtime as the line replaced.

julia> function f1!(x)
           return float(x) .^ 2
       end
f1! (generic function with 1 method)

julia> function f2!(x)
           return float(abs.(x)) .^ 2
       end
f2! (generic function with 1 method)

julia> function f3!(x)
           return conj.(float(x)).*float(x)
       end
f3! (generic function with 1 method)

For vec = rand(ComplexF64, 1000), I found the following using benchmark tools: complexvecbenchmark For vec2 = rand(Float64, 1000), I found the following using benchmark tools: float64vecbenchmark

Checklist

  • [ ] Appropriate tests were added (I think this change is too minor to necessitate new unit test.)
  • [x] Any code changes were done in a way that does not break public API
  • [x] All documentation related to code changes were updated
  • [x] The new code follows the contributor guidelines, in particular the SciML Style Guide and COLPRAC.
  • [x] Any new documentation only uses public API

Additional context

Add any other context about the problem here.

henhen724 avatar Sep 10 '24 07:09 henhen724