Skip to content

Commit

Permalink
feat: use HomotopyContinuationProblem in NonlinearProblem if poss…
Browse files Browse the repository at this point in the history
…ible
  • Loading branch information
AayushSabharwal committed Nov 30, 2024
1 parent fecca74 commit 55fd1cc
Show file tree
Hide file tree
Showing 2 changed files with 21 additions and 2 deletions.
6 changes: 5 additions & 1 deletion src/systems/nonlinear/nonlinearsystem.jl
Original file line number Diff line number Diff line change
Expand Up @@ -492,10 +492,14 @@ end

function DiffEqBase.NonlinearProblem{iip}(sys::NonlinearSystem, u0map,
parammap = DiffEqBase.NullParameters();
check_length = true, kwargs...) where {iip}
check_length = true, use_homotopy_continuation = true, kwargs...) where {iip}
if !iscomplete(sys)
error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `NonlinearProblem`")
end
prob = safe_HomotopyContinuationProblem(sys, u0map, parammap; check_length, kwargs...)
if prob isa HomotopyContinuationProblem
return prob
end
f, u0, p = process_SciMLProblem(NonlinearFunction{iip}, sys, u0map, parammap;
check_length, kwargs...)
pt = something(get_metadata(sys), StandardNonlinearProblem())
Expand Down
17 changes: 16 additions & 1 deletion test/extensions/homotopy_continuation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,12 +22,19 @@ import HomotopyContinuation
0 ~ x^2 + 4x + 4
0 ~ y * z + 4x^2]
@mtkbuild sys = NonlinearSystem(eqs)
prob = HomotopyContinuationProblem(sys, [x => 1.0, y => 1.0, z => 1.0], [])
u0 = [x => 1.0, y => 1.0, z => 1.0]
prob = HomotopyContinuationProblem(sys, u0)
@test prob[x] == prob[y] == prob[z] == 1.0
@test prob[x + y] == 2.0
sol = solve(prob; threading = false)
@test SciMLBase.successful_retcode(sol)
@test norm(sol.resid)0.0 atol=1e-10

prob2 = NonlinearProblem(sys, u0)
@test prob2 isa HomotopyContinuationProblem
sol = solve(prob2; threading = false)
@test SciMLBase.successful_retcode(sol)
@test norm(sol.resid)0.0 atol=1e-10
end

struct Wrapper
Expand Down Expand Up @@ -92,36 +99,44 @@ end
"Exponent", "not an integer", "not a polynomial"] HomotopyContinuationProblem(
sys, [])
@test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError
@test NonlinearProblem(sys, []) isa NonlinearProblem

@mtkbuild sys = NonlinearSystem([x^x - x ~ 0])
@test_throws ["Cannot convert", "Unable", "symbolically solve",
"Exponent", "unknowns", "not a polynomial"] HomotopyContinuationProblem(
sys, [])
@test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError
@test NonlinearProblem(sys, []) isa NonlinearProblem
@mtkbuild sys = NonlinearSystem([((x^2) / sin(x))^2 + x ~ 0])
@test_throws ["Cannot convert", "both polynomial", "non-polynomial",
"recognized", "sin", "not a polynomial"] HomotopyContinuationProblem(
sys, [])
@test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError
@test NonlinearProblem(sys, []) isa NonlinearProblem

@variables y = 2.0
@mtkbuild sys = NonlinearSystem([x^2 + y^2 + 2 ~ 0, y ~ sin(x)])
@test_throws ["Cannot convert", "recognized", "sin", "not a polynomial"] HomotopyContinuationProblem(
sys, [])
@test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError
@test NonlinearProblem(sys, []) isa NonlinearProblem

@mtkbuild sys = NonlinearSystem([x^2 + y^2 - 2 ~ 0, sin(x + y) ~ 0])
@test_throws ["Cannot convert", "function of multiple unknowns"] HomotopyContinuationProblem(
sys, [])
@test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError
@test NonlinearProblem(sys, []) isa NonlinearProblem

@mtkbuild sys = NonlinearSystem([sin(x)^2 + 1 ~ 0, cos(y) - cos(x) - 1 ~ 0])
@test_throws ["Cannot convert", "multiple non-polynomial terms", "same unknown"] HomotopyContinuationProblem(
sys, [])
@test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError
@test NonlinearProblem(sys, []) isa NonlinearProblem

@mtkbuild sys = NonlinearSystem([sin(x^2)^2 + sin(x^2) - 1 ~ 0])
@test_throws ["import Nemo"] HomotopyContinuationProblem(sys, [])
@test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError
@test NonlinearProblem(sys, []) isa NonlinearProblem
end

import Nemo
Expand Down

0 comments on commit 55fd1cc

Please sign in to comment.