# Numerical solution of a system of non linear equations

This post is a wiki. Anyone with karma >750 is welcome to improve it.

Hi everybody! I would like to know whether there exists a simple way to solve in a numerical way a system of non linear equations (the system cannot be solved analytically). In particular, the system (in the variables a,b) I'd like to solve is the following:

-2(ab/sqrt(-1/2a^2 + b^2) - 2/(a^2b^2))(sqrt(-1/2a^2 + b^2)b - 1/(ab^2)) - 16/(a^5*b^2)=0

4b^3 + 4(sqrt(-1/2a^2 + b^2)b - 1/(ab^2))(b^2/sqrt(-1/2a^2 + b^2) + sqrt(-1/2a^2 + b^2) + 2/(ab^3)) - 8/(a^4b^3)=0

edit retag close merge delete

Can you be a bit more vague, your description contains too many details. Seriously, what kind of equations? algebraic, differential, polynomial, transcendental, ...

( 2011-12-12 07:47:35 +0200 )edit

Yes, I corrected my original post with the equations I need to solve numerically

( 2011-12-12 08:11:53 +0200 )edit

Sort by » oldest newest most voted

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.

more

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

more

This post is a wiki. Anyone with karma >750 is welcome to improve it.

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.

more

This post is a wiki. Anyone with karma >750 is welcome to improve it.

Thank you very much! Best, f

more