numpy - Sympy 和 Numpy 求解器之间的区别

Question

我在使用 nsolve 时遇到了一些问题很难找到一些函数的解决方案，给出一些初步的猜测。然后我想尝试使用 numpy/scipy 求解器。

这是一个使用 sympy 的程序，可以很好地提供此解决方案:[0.0, -9.05567e-72, 9.42477, 3.14159]

from sympy import *

# Symbols
theta = Symbol('theta')
phi = Symbol('phi')
phi0 = Symbol('phi0')
H0 = Symbol('H0')
# Constants
phi0 = 60*pi.evalf()/180
a = 0.05
t = 100*1e-9
b = 0.05**2/(8*pi.evalf()*1e-7)
c = 0.001/(4*pi.evalf()*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 -(t*a*H0)*k(theta,phi,phi0)+b*t*(cos(theta)**2)+c*t*(sin(2*theta)**2)+t*sin(theta)**4*sin(2*phi)**2
def F_phi(theta,phi,phi0,H0):
    return diff(F(theta,phi,phi0,H0),phi)
def G(phi):
    return F_phi(theta,phi,phi0,H0).subs(theta,pi/2)

H0 = -0.03/(4*pi.evalf()*1e-7)
sol = []
for i in range(5):
    x0=i*pi.evalf()/4
    solution = float(nsolve(G(phi),x0))
    sol.append(solution)
sol = list(set(sol)) # remove duplicate values
print sol

这是同一个程序，但使用 numpy 兼容函数:

from numpy import *
from scipy.optimize import fsolve
# Constants
phi0 = 60*pi/180
a = 0.05
t = 100*1e-9
b = 0.05**2/(8*pi*1e-7)
c = 0.001/(4*pi*1e-7)

def m(theta,phi):
    return array([sin(theta)*cos(phi),sin(theta)*cos(phi),cos(phi)])
def h(phi0):
    return array([cos(phi0),sin(phi0),0])
def k(theta,phi,phi0):
    return dot(m(theta,phi).T,h(phi0))
def F(theta,phi,phi0,H0): 
    return -(t*a*H0)*k(theta,phi,phi0)+b*t*(cos(theta)**2)+c*t*(sin(2*theta)**2)+t*sin(theta)**4*sin(2*phi)**2
def F_phi(theta,phi,phi0,H0):
    return diff(F(theta,phi,phi0,H0),phi)
def G(phi):
    return F_phi(pi/2,phi,phi0,H0)

H0 = -0.03/(4*pi*1e-7)
sol = []
for i in range(5):
    x0=array([i*pi/4]) # x0 as ndarray argument for fsolve
    solution = float(fsolve(G,x0))
    sol.append(solution)
sol = list(set(sol)) # remove duplicate values
print sol

但是当我运行程序时:

Traceback (most recent call last):
File "Test4.py", line 27, in <module>
solution = float(fsolve(G,x0))
File "/usr/lib64/python2.7/site-packages/scipy/optimize/minpack.py", line 127, in fsolve
res = _root_hybr(func, x0, args, jac=fprime, **options)
File "/usr/lib64/python2.7/site-packages/scipy/optimize/minpack.py", line 224, in _root_hybr
raise errors[status][1](errors[status][0])
TypeError: Improper input parameters were entered.

我尝试给 x0 赋值 0，第二个程序(使用 numpy)给出一个接近 0 的数值，但是从 pi/4 开始，它给出了错误消息。我错过了 numpy 中的某些东西吗？

Answer 1

在 numpy 版本函数中 G(array([pi/4]))返回一个空数组:

>> G(array([pi/4]))  
array([], dtype=float64)

问题在于:

return diff(F(theta,phi,phi0,H0),phi)

numpy.diff计算数组的连续元素之间的差异，而 sympy.diff计算导数。您可以自行修改 F_phi函数返回分析计算的导数(如果您知道解决方案)或数值计算。对于数值解，您可以使用:

def F_phi(theta,phi,phi0,H0, eps=1e-12):
    return (F(theta,phi+eps,phi0,H0) - F(theta,phi,phi0,H0))/eps

和解析解(用 sympy 计算):

def F_phi(theta, phi, phi0, H0):
    return -H0*a*t*(-sin(phi)*sin(phi0)*sin(theta) - sin(phi)*sin(theta)*cos(phi0)) + 4*t*sin(2*phi)*sin(theta)**4*cos(2*phi)

请记住，数值解不会像解析那样精确。因此，sympy(分析)和numpy(数值)方法之间可能仍然存在差异。