2012-11-14 3 views
2

nsolve에 약간의 문제가있어서 초기 추측을 제공하는 일부 기능에 대한 해결책을 찾기가 어려웠습니다. 나는 그때 수줍음/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에서 시작하여 오류 메시지를 표시합니다. 나는 수치스럽고 무언가를 놓쳤는가?

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

문제는 일치한다 : sympy.diff이 유도체를 산출하는 반면

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

numpy.diff는, 어레이의 연속 된 요소 사이의 차이를 계산 빈 배열을 반환 G(array([pi/4])) NumPy와 버전 함수에서

+0

귀하의 문제에 대한 해결책을 제공하는 답변을 반드시 수락하십시오 (답변 옆의 체크 표시 사용). – btel

답변

3

. 자신의 F_phi 함수를 수정하여 분석적으로 계산 된 미분을 (솔루션을 알고있는 경우) 또는 수치 적으로 반환 할 수 있습니다. 수치 솔루션을 위해 당신은 사용할 수 있습니다 : (sympy 계산)

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

및 분석 솔루션 :

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 (수치 적) 접근법 간에는 여전히 차이가있을 수 있습니다.

+0

예, 두 라이브러리에서 모두 동일하다고 생각되는 diff 함수에서 문제가 발생했습니다. 실제로, numpy의 해는 sympy의 해와는 다른 것으로 보인다. (이 마지막 경우에서, 해는 nsolve를 사용하여 수치 적으로 계산되었거나,이 함수의 범위를 오해 했는가?). 게다가, eps의 가치에 따라, 나는 다른 해결책을 가지고있다 : 당신이 내게 준 값으로 RuntimeWarning과 관련된 경고가 있었지만, 해결책의 세트는'eps <1e-18'에 대해 안정화되는 것처럼 보인다. sympy와 함께 위의 솔루션. – aymenbh

+0

더 정확한 분석 솔루션을 추가했습니다. – btel

+0

'phi'가 주기적이기 때문에'solution % (2 * pi)'를 계산할 수도 있습니다. 이 수정으로'sympy'와'numpy' 솔루션이 매우 근접합니다! – btel