Revision history [back]

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.