Ask Your Question
3

Numerical solution of a system of non linear equations

asked 2011-12-11 22:21:08 -0500

this post is marked as community wiki

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

Thanks in advanced, Francesco

edit retag flag offensive close merge delete

Comments

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

Volker Braun gravatar imageVolker Braun ( 2011-12-12 00:47:35 -0500 )edit

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

Francesco gravatar imageFrancesco ( 2011-12-12 01:11:53 -0500 )edit

4 answers

Sort by ยป oldest newest most voted
2

answered 2011-12-13 05:30:06 -0500

DSM gravatar image

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.

edit flag offensive delete link more
0

answered 2011-12-13 03:50:09 -0500

this post is marked as community wiki

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.

edit flag offensive delete link more
0

answered 2011-12-13 06:10:27 -0500

this post is marked as community wiki

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

Thank you very much! Best, f

edit flag offensive delete link more
0

answered 2017-02-14 04:47:54 -0500

maaaaaaartin gravatar image

updated 2017-02-14 04:48:17 -0500

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

edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2011-12-11 22:21:08 -0500

Seen: 1,693 times

Last updated: Feb 14