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 を指定しようとしましたが、2 番目のプログラム (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 に答える