I was happily surprised to discovered that homotopy continuation was implemented in Julia. I have been trying the package out, however there are a few non obvious things which I do not understand (this might partially be due to my lack of experience in the area). I am unsure whenever this is the right place to ask my questions, but it was the best I found (if not I apologise).
First, it seems that when there is a singularity present an additional 1.0
is added to the solution array, e.g.
@polyvar x y z w
f = [x*y+z^2, y-3x*z+z^2+2, x^3+y^2+z-1];
sols = solutions(solve(f), only_real=true)
for i = 1:length(sols)
println(real(sols[i].solution))
end
outputs
[3.62126e-46, 1.50418e-25, 1.0, 8.62497e-7]
[-1.00432e-24, -7.17904e-17, 1.0, -7.61496e-7]
[-1.40608, 2.36735, -1.82447]
[-76.8682, 674.107, -227.634]
where both the first and second solution contains a singularity. However, I did not find anything on what the actual meaning of the 1.0
and its position is.
Next I have a few things which all might be regarding my lack of knowledge in the area, I am no expert. However, when I run the part
sols = solutions(solve(f), only_real=true)
for i = 1:length(sols)
println(real(sols[i].solution))
end
several times I do not all time get the previously mentioned solution, sometimes the solution [-76.8682, 674.107, -227.634]
is not found (roughly half of the time). In addition this solution does not actually seem to be a good solution, indeed
sol = [-76.8682, 674.107, -227.634]
for i = 1:length(f)
println(f[i](x=> sol[1] , y=>sol[2], z=>sol[3])," ")
end
returns
-0.15374140000494663
-0.10256040000604116
-1.070417910619625
is the solving method simply inaccurate for certain cases, not running long enough, or am I doing something wrong? The documentations says
The aim of this project is twofold: establishing a fast numerical polynomial solver in Julia and at the same time providing a highly customizable algorithmic environment for researchers for designing and performing individual experiments.
Then the documentations starts with a simple example of solving a polynomial system (what I want to do and also what I am trying here), and then continue with advanced stuff of defining your own homotopies and pathtracking. I am uncertain where the "fast numerical polynomial solver" (what I am interested in) ends and the "highly customizable algorithmic environment for researchers for designing and performing individual experiments" starts.
If the issue lies in my ability to use the methods and that they require some experience, do you maybe have a suggested reading for understanding more?
A further example where I struggle. For 4 variable systems I have been unable to find proper solutions for even simple linear systems like
@polyvar x y z w
f = [w+x-4y, z+3-2x, 3w+x, 10-w+y];
sols = solutions(solve(f), only_real=true)
for i = 1:length(sols)
println(real(sols[i].solution))
end
This gives
[6.66667, -20.0, -3.33333, -43.0]
However,
sol = [6.66667, -20.0, -3.33333, -43.0]
for i = 1:length(f)
println(f[i](x=> sol[1] , y=>sol[2], z=>sol[3], w=>sol[4])," ")
end
outputs
43.666669999999996
-13.66667
-122.33333
33.0
which is distinctively non zero.