Numerical solution of a system of non linear equations

edit

You could use mpmath's multidimensional solver, but unfortunately it's not well-wrapped either. On the bright side you can do it without any Cython:

import mpmath

var("a b")
eq0 = -2*(a*b/sqrt(-1/2*a^2 + b^2) - 2/(a^2*b^2))*(sqrt(-1/2*a^2 + b^2)*b - 1/(a*b^2)) - 16/(a^5*b^2)==0
eq1 = 4*b^3 + 4*(sqrt(-1/2*a^2 + b^2)*b - 1/(a*b^2))*(b^2/sqrt(-1/2*a^2 + b^2) + sqrt(-1/2*a^2 + b^2) + 2/(a*b^3)) - 8/(a^4*b^3)==0

f = [lambda a,b: eq0.lhs().subs(a=RR(a), b=RR(b)),
     lambda a,b: eq1.lhs().subs(a=RR(a), b=RR(b))]

found_root = mpmath.findroot(f, (2, 2))
found_root = Matrix(RR, found_root.tolist())

print found_root
fa,fb = found_root.list()

print eq0.subs(a=fa,b=fb)
print eq1.subs(a=fa,b=fb)

which gives:

sage: load "rfind.sage"
[1.49174202118713]
[1.10177415033843]
(-2.67785793539588e-13) == 0
(4.93827201353270e-13) == 0

and so we see it's done a decent job. The above was a very manual wrapping, it wouldn't be that much harder to treat the general case.

You can try minimizing sqrt(f1(a,b)**2 + f2(a,b)**2) using sage.minimize.

The included GSL library has a multi-dimensional numerical root finder gsl_multiroot_function() but it seems that Sage does not wrap that functionality. You could write your own wrapper. Should be easy enough using sage/gsl/gsl_roots.pxi as a template.

Alternatively, you can introduce dummy variables for the square roots and clear denominators. This will give you a system of polynomial equations which you can attack using Groebner basis techniques.