2

scipy.integrate.quad를 사용하여 두 개의 시간 및 주파수 편이 된 Hermite 함수의 곱을 통합하려고합니다.Scipy : 직각 가중치를 이용한 Hermite 함수의 통합

그러나 큰 순서 다항식이 포함되므로 숫자 오류가 발생합니다. 여기 내 코드입니다 :

import numpy as np 
import scipy.integrate 
import scipy.special as sp 
from math import pi 


def makeFuncs(): 
    # Create the 0th, 4th, 8th, 12th and 16th order hermite function 
    return [lambda t, n=n: np.exp(-0.5*t**2)*sp.hermite(n)(t) for n in np.arange(5)*4] 

def ambgfun(funcs, i, k, tau, f): 
    # Integrate f1(t)*f2(t+tau)*exp(-j2pift) over t from -inf to inf 
    f1 = funcs[i] 
    f2 = funcs[k] 
    func = lambda t: np.real(f1(t) * f2(t+tau) * np.exp(-1j*(2*pi)*f*t)) 
    return scipy.integrate.quad(func, -np.inf, np.inf) 

def main(): 
    f = makeFuncs() 

    print "A00(0,0):", ambgfun(f, 0, 0, 0, 0) 
    print "A01(0,0):", ambgfun(f, 0, 1, 0, 0) 
    print "A34(0,0):", ambgfun(f, 3, 4, 0, 0) 

if __name__ == '__main__': 
    main() 

미트 기능은 직교, 따라서 모든 적분은 제로 같아야한다. 그러나 출력은 다음과 같이 표시되지 않습니다.

A00(0,0): (1.7724538509055159, 1.4202636805184462e-08) 
A01(0,0): (8.465450562766819e-16, 8.862237123626351e-09) 
A34(0,0): (-10.1875, 26.317246925873935) 

이 계산을보다 정확하게 만들려면 어떻게해야합니까? scipy로부터의 hermite-function은 문서 (http://docs.scipy.org/doc/scipy/reference/special.html#orthogonal-polynomials)에 나와있는 것처럼 Gaussian Quadrature에 사용해야하는 가중치 변수를 포함합니다. 그러나, 나는 문서에서 이러한 가중치를 사용하는 방법에 대한 힌트를 찾지 못했습니다.

나는 대답은 당신이 얻을 결과가 수치가수록 제로에 가까운 점이다 최대

답변

1

, 당신이 :)

감사합니다 도움이 될 수 있기를 바랍니다. 부동 소수점 숫자로 작업하면 더 좋은 결과를 얻는 것이 실제로 가능하지 않다고 생각합니다. 수치 적 통합에 일반적인 문제가 있습니다.

import numpy as np 
from scipy import integrate, special 
f = lambda t: np.exp(-t**2) * special.eval_hermite(12, t) * special.eval_hermite(16, t) 

abs_ig, abs_err = integrate.quad(lambda t: abs(f(t)), -np.inf, np.inf) 
ig, err = integrate.quad(f, -np.inf, np.inf) 

print ig 
# -10.203125 
print abs_ig 
# 2.22488114805e+15 
print ig/abs_ig, err/abs_ig 
# -4.58591912155e-15 1.18053770382e-14 

적분의 값은 따라서, 부동 소수점 비교 엡실론 정확도로 계산되었다 :

이 고려 해보자. 큰 크기로 진동하는 피고 산의 값을 뺄 때의 반올림 오류로 인해 더 나은 결과를 얻는 것은 실제로 가능하지 않습니다.

어떻게 진행하나요? 내 경험상, 지금해야 할 일은 수치 적이 아니라 분석적으로 문제에 접근하는 것입니다. 중요한 것은 Hermite 다항식의 푸리에 변환 곱하기 함수가 알려져 있으므로 여기에서 항상 푸리에 공간에서 작업 할 수 있습니다.

+0

나는 또한 상대적인 오류라고 생각합니다. 나는 정규화되고 직교 정규 형 다항식으로 작업하고 있었고, 내가 기억하는 한 내 오류는 절대적인 용어로는 훨씬 작았지만 상대적인 의미는 아닌 것 같다. – user333700

관련 문제