2

我遇到了一些问题，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 中错过了什么吗？

4

1 回答 1

3

在 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（数值）方法之间可能仍然存在差异。

于 2012-11-14T15:32:41.700 回答