Scipy: Интегриране на функцията на Hermite с квадратурни тегла

Искам да интегрирам продукта от две функции на Hermite, изместени във времето и честотата, като използвам 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)

Как мога да направя това изчисление по-точно? Hermite-функцията от scipy съдържа променлива за тегла, която трябва да се използва за квадратура на Гаус, както е дадено в документацията (http://docs.scipy.org/doc/scipy/reference/special.html#orthogonal-полиноми). Въпреки това не намерих подсказка в документите как да използвам тези тежести.

Надявам се, че можете да помогнете :)

Благодаря, Макс


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