|
1 | 1 | module MTKHomotopyContinuationExt |
2 | 2 |
|
3 | 3 | using ModelingToolkit |
| 4 | +using ModelingToolkit.SciMLBase |
4 | 5 | using ModelingToolkit.Symbolics: unwrap |
5 | 6 | using ModelingToolkit.SymbolicIndexingInterface |
6 | 7 | using HomotopyContinuation |
@@ -98,17 +99,26 @@ function ModelingToolkit.HomotopyContinuationProblem( |
98 | 99 |
|
99 | 100 | mtkhsys = MTKHomotopySystem(nlfn.f, p, nlfn.jac, hvars, length(eqs)) |
100 | 101 |
|
101 | | - prob = ModelingToolkit.HomotopyContinuationProblem{typeof(mtkhsys), typeof(u0)}( |
102 | | - sys, mtkhsys, u0) |
| 102 | + return ModelingToolkit.HomotopyContinuationProblem(u0, mtkhsys, sys) |
103 | 103 | end |
104 | 104 |
|
105 | 105 | function CommonSolve.solve(prob::ModelingToolkit.HomotopyContinuationProblem; kwargs...) |
106 | | - sol = HomotopyContinuation.solve(prob.hcsys; kwargs...) |
107 | | - rsols = HomotopyContinuation.real_solutions(sol) |
108 | | - rsol = findmin(rsols) do val |
109 | | - norm(prob.u0 - val) |
| 106 | + sol = HomotopyContinuation.solve(prob.homotopy_continuation_system; kwargs...) |
| 107 | + realsols = HomotopyContinuation.results(sol; only_real = true) |
| 108 | + if isempty(realsols) |
| 109 | + u = state_values(prob) |
| 110 | + resid = prob.homotopy_continuation_system(u) |
| 111 | + retcode = SciMLBase.ReturnCode.ConvergenceFailure |
| 112 | + else |
| 113 | + distance, idx = findmin(realsols) do result |
| 114 | + norm(result.solution - state_values(prob)) |
| 115 | + end |
| 116 | + u = real.(realsols[idx].solution) |
| 117 | + resid = prob.homotopy_continuation_system(u) |
| 118 | + retcode = SciMLBase.ReturnCode.Success |
110 | 119 | end |
111 | | - return rsol |
| 120 | + |
| 121 | + return SciMLBase.build_solution(prob, :HomotopyContinuation, u, resid; retcode) |
112 | 122 | end |
113 | 123 |
|
114 | 124 | end |
0 commit comments