solve system of non-linear implicit equations numerically
I am attempting to solve for a solution of a system of two non-linear implicit equations using the following code:
x = var('x') y = var('y') P = [(-1,-5), (1,-5), (-5,0), (5,5)] # Defining the function d = sum([sqrt( (x-p)^2 + (y-p)^2 ) for p in P]) show(d) # Differentiate with respect to x and y eqx = d.diff(x) eqy = d.diff(y) # Plot both implicit curves g1 = implicit_plot( eqx==0, (x,-10,10), (y,-10,10), color="blue" ) g2 = implicit_plot( eqy==0, (x,-10,10), (y,-10,10), color="red" ) show(g1 + g2) # note that you can clearly see an intersection of the two curves # Solve for the solution print("Solving...") sol = solve([eqx==0, eqy==0], x, y) # this gets stuck or takes a long time show(sol)
Everything runs, up to the point of the solve function, which continues to run for what appears to be indefinitely. The code show(g1 + g2) shows a graph that clearly shows there exists an intersection for both curves. I tried to use to_poly_solve=True without success. I do not mind an approximate solution, however I was unable to find a numeric solver for a system such as this (find_root afaik only works on one variable) that will work.
Does there exist a numeric solver which is capable of solving a system of this form? What other alternatives are there?