2014-11-09 3 views
-1

나는 함수를 통합하려고한다. 그래서 이런 짓을 확인하기 위해, 내가 무엇을 기대 내가지고있어 결과가 없습니다Scipy 적분이 잘못 나온다

def integrand(a, b): 
    return scipy.integrate.quad(function, a, b) 

: 나는 그것을 통합하고있어

이제
def function(x): 
    something = ... 
    something_else = ... 
    return exp(something)/sqrt(something_else) 

:이 기능은 비 부정적 보장

for x in range(0, 10000): 
    if integrand(0,x+1) < integrand(0,x): 
     raise ValueError("Weird!") 

물론, '이상한'예외가 있습니다. 어떻게 그렇게 될수 있니?

+1

문제를 나타내는 실제 예제를 표시 할 수 있습니까? – BrenBarn

+0

실제 예제는 10,000 포인트 C14 데이트 보정 곡선이 포함되어 있기 때문에 어렵습니다. – zmbq

+1

작은 테스트 케이스 없이는 문제를 디버그하기가 어렵습니다. 문제를 작은 테스트 케이스로 줄이기 위해 노력하십시오. – BrenBarn

답변

1

예제가 없으면 코드를 디버그하는 것이 어렵지만 여기에 제시된 내용을 기반으로하는 몇 가지 제안 사항이 있습니다.

quad은 두 값 (적분의 추정값 및 결과의 절대 오차 추정값)이 포함 된 튜플을 반환합니다. 코드에서 첫 번째 값만 사용하는 것처럼 보이지는 않습니다. quad에 대한 호출 후

def integrand(a, b): 
    return scipy.integrate.quad(function, a, b)[0] 

참고 [0]integrand을 변경해보십시오.

다음으로 "이상한"사례를 발견하면 integrand(0, x)integrand(0, x+1)의 값을 인쇄해야합니다. 그것들이 대략 동일하다면, 수치 적 부정확성은 당신이 이론적으로 기대하는 순서를 파괴 할 수 있습니다. quad에 의해 반환 된 두 번째 값을 살펴보십시오. 이 값이 integrand(0, x)integrand(0, x+1) 사이의 차이보다 크면 예상 증가량이 계산 오류보다 작기 때문에 실제적으로 추정치가 단조 증가한다고 기대할 수 없습니다. 이것이 문제인 경우 epsabs 및/또는 epsrel을 기본값 인 1.49e-8보다 작은 값으로 설정해 볼 수도 있지만 지금까지만 가능합니다. 결국에는 부동 소수점 계산의 한계에 부딪치게됩니다.

관련 문제