3

After fixing the imports in the previous question (Python AttributeError:cos) and changing a little bit the functions using the sympy ones:

from sympy import *
from sympy import Symbol
from sympy.solvers import nsolve

# Symbols
theta = Symbol('theta')
phi = Symbol('phi')
phi0 = Symbol('phi0')
H0 = Symbol('H0')
# Constants
a = 0.05 
b = 0.05**2/(8*pi*1e-7)
c = 0.001/(4*pi*1e-7)
phi0 = 60*pi/180 
H0 = -0.03/(4*pi*1e-7)
def m(theta,phi):
    return Matrix([[sin(theta)*cos(phi), sin(theta)*cos(phi), cos(phi)]])
def h(phi0):
    return Matrix([[cos(phi0), sin(phi0), 0]])
def k(theta,phi,phi0):
    return m(theta,phi).dot(h(phi0))
def F(theta,phi,phi0,H0): 
    return -(a*H0)*k(theta,phi,phi0)+b*(cos(theta)**2)+c*(sin(2*theta)**2)+sin(theta)**4*sin(2*phi)**2
def F_phi(theta,phi,phi0,H0):
    return simplify(diff(F(theta,phi,phi0,H0),phi))
def G(phi):
    return F_phi(pi/2,phi,phi0,H0)
solution = nsolve(G(phi), phi)
print solution

I get this Traceback:

Traceback (most recent call last):
File "Test.py", line 31, in <module>
solution = nsolve(G(phi), phi)
File "/usr/lib64/python2.7/site-packages/sympy/solvers/solvers.py", line 2050, in nsolve
return findroot(f, x0, **kwargs)
File "/usr/lib64/python2.7/site-packages/mpmath/calculus/optimization.py", line 908, in findroot
x0 = [ctx.convert(x0)]
File "/usr/lib64/python2.7/site-packages/mpmath/ctx_mp_python.py", line 662, in convert
return ctx._convert_fallback(x, strings)
File "/usr/lib64/python2.7/site-packages/mpmath/ctx_mp.py", line 561, in _convert_fallback
raise TypeError("cannot create mpf from " + repr(x))
TypeError: cannot create mpf from phi
Community
  • 1
  • 1
aymenbh
  • 151
  • 4
  • 10

1 Answers1

3

The second argument to nsolve should be an initial guess, not a Symbol. From the docstring

Solve a nonlinear equation system numerically::

    nsolve(f, [args,] x0, modules=['mpmath'], **kwargs)

f is a vector function of symbolic expressions representing the system.
args are the variables. If there is only one variable, this argument can
be omitted.
x0 is a starting vector close to a solution.

So you want something like nsolve(G(phi), 0) or nsolve(G(phi), 3), etc. (depending on what solution you want).

asmeurer
  • 86,894
  • 26
  • 169
  • 240
  • Thanks! I didn't check the nsolve function properly. Actually I was transposing a Mathematica program in Python. But it seems that this function gives you one solution depending on your initial guess. Is it possible to get the list of all possible solutions? (in this case, it's 0 and pi). – aymenbh Nov 12 '12 at 09:55
  • I just used a for loop against a small range. Unfortunately, it raises an exception for values that are too far from the solution, so you'll have to wrap it in a try block. In theory, SymPy should be able to solve this sort of thing symbolically, but the algorithms aren't implemented yet. – asmeurer Nov 12 '12 at 17:33
  • I actually tried it with a for loop and it works fine. Hope to see the algorithms to solve symbolically this kind of equations implemented soon! – aymenbh Nov 12 '12 at 18:30
  • 1
    Actually, it can find one of the solutions (0) if you first apply `expand(trig=True)`. SymPy is still a long way off from solving trig expressions (for one thing, it needs a good way to express and manipulate infinitely many solutions). – asmeurer Nov 13 '12 at 01:24
  • I used a for loop on uniformly distributed values (as initial guesses) over [0, pi] (i*pi/4, i in range(5)) and I kept track of the solutions in a list. If you try it you can see that you may have the same solution for different i. What I find strange though is that I get 9.42 as a solution for i=1 but if you keep only values in [0,pi] you end up by 2 solutions. – aymenbh Nov 13 '12 at 10:26
  • Probably `n*pi` is a solution for any `n`. nsolve is not guaranteed to find any solution in particular, even the one closest to the initial guess. – asmeurer Nov 13 '12 at 19:08
  • Yes, I'm experiencing problems with nsolve in some situations right now. And the program seems to take much more time to be executed than the equivalent in Mathematica. Is it normal? Do you know if Numpy/Scipy are better than Sympy? – aymenbh Nov 14 '12 at 00:24
  • SymPy's focus is on symbolics. If you are going to use numerics, then something like numpy/scipy will be better, yes. – asmeurer Nov 14 '12 at 01:45