Scipy: интегрирование функции Эрмита с квадратурными весами

Я хочу интегрировать продукт двух функций Эрмита со сдвигом по времени и частоте, используя scipy.integrate.quad.

Однако, поскольку включены полиномы большого порядка, возникают числовые ошибки. Вот мой код:

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 содержит переменную веса, которую следует использовать для гауссовой квадратуры, как указано в документации (http://docs.scipy.org/doc/scipy/reference/special.html#orthogonal-polynomials). Однако я не нашел в документах подсказки, как использовать эти веса.

Я надеюсь, что вы можете помочь :)

Спасибо, Макс


python numpy scipy signal-processing numerical-integration
person Maximilian Matthé    schedule 20.12.2012    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Ответ заключается в том, что полученный вами результат численно максимально близок к нулю. Я не думаю, что действительно возможно получить гораздо лучшие результаты, если вы работаете с числами с плавающей запятой - вы сталкиваетесь с общей проблемой численного интегрирования.

Учти это:

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

Таким образом, значение подынтегральной функции было вычислено с точностью, сравнимой с эпсилон с плавающей запятой. Из-за ошибки округления при вычитании значений колеблющегося подынтегрального выражения большой величины на самом деле невозможно получить лучшие результаты.

Итак, как действовать? По моему опыту, сейчас вам нужно подойти к проблеме не численно, а аналитически. Важно отметить, что преобразование Фурье полиномов Эрмита, умноженное на весовую функцию, известно, так что здесь вы можете постоянно работать в пространстве Фурье.

person pv.    schedule 20.12.2012
comment
Я думаю, что это также вопрос относительной ошибки. Я работал с нормализованными ортонормированными полиномами Эрмита, и, насколько я помню, мои ошибки были намного меньше в абсолютном выражении, но, думаю, не в относительном. - person Josef; 21.12.2012