DiffEqDevTools.jl
DiffEqDevTools.jl copied to clipboard
Fixed L2 error for complex valued problems
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:
For
vec2 = rand(Float64, 1000)
, I found the following using benchmark tools:
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.