diff --git a/ext/MTKHomotopyContinuationExt.jl b/ext/MTKHomotopyContinuationExt.jl index e4a9243ff9..ebd8afbfc1 100644 --- a/ext/MTKHomotopyContinuationExt.jl +++ b/ext/MTKHomotopyContinuationExt.jl @@ -196,7 +196,18 @@ function MTK.HomotopyContinuationProblem( error("Equation $eq is not a polynomial in the unknowns. See warnings for further details.") end num, den = handle_rational_polynomials(eq.rhs - eq.lhs, dvs) - push!(denoms, den) + + # make factors different elements, otherwise the nonzero factors artificially + # inflate the error of the zero factor. + if iscall(den) && operation(den) == * + for arg in arguments(den) + # ignore constant factors + symbolic_type(arg) == NotSymbolic() && continue + push!(denoms, abs(arg)) + end + elseif symbolic_type(den) != NotSymbolic() + push!(denoms, abs(den)) + end return 0 ~ num end @@ -232,7 +243,7 @@ Extra keyword arguments: the denominator to be below `denominator_abstol` will be discarded. """ function CommonSolve.solve(prob::MTK.HomotopyContinuationProblem, - alg = nothing; show_progress = false, denominator_abstol = 1e-8, kwargs...) + alg = nothing; show_progress = false, denominator_abstol = 1e-7, kwargs...) sol = HomotopyContinuation.solve( prob.homotopy_continuation_system; show_progress, kwargs...) realsols = HomotopyContinuation.results(sol; only_real = true) diff --git a/test/extensions/homotopy_continuation.jl b/test/extensions/homotopy_continuation.jl index 6d7279899f..9bd64fa7fb 100644 --- a/test/extensions/homotopy_continuation.jl +++ b/test/extensions/homotopy_continuation.jl @@ -98,7 +98,9 @@ end @test sol[x] ≈ 1.0 p = parameter_values(prob) for invalid in [2.0, 3.0] - @test prob.denominator([invalid], p)[1] <= 1e-8 + for err in [-9e-8, 0, 9e-8] + @test any(<=(1e-7), prob.denominator([invalid + err, 2.0], p)) + end end @named sys = NonlinearSystem( @@ -115,14 +117,18 @@ end disallowed_y = [7, 5, 4] @test all(!isapprox(sol[x]; atol = 1e-8), disallowed_x) @test all(!isapprox(sol[y]; atol = 1e-8), disallowed_y) - @test sol[x^2 - 4x + y] >= 1e-8 + @test abs(sol[x^2 - 4x + y]) >= 1e-8 p = parameter_values(prob) for val in disallowed_x - @test any(<=(1e-8), prob.denominator([val, 2.0], p)) + for err in [-9e-8, 0, 9e-8] + @test any(<=(1e-7), prob.denominator([val + err, 2.0], p)) + end end for val in disallowed_y - @test any(<=(1e-8), prob.denominator([2.0, val], p)) + for err in [-9e-8, 0, 9e-8] + @test any(<=(1e-7), prob.denominator([2.0, val + err], p)) + end end @test prob.denominator([2.0, 4.0], p)[1] <= 1e-8 end