Двойной интеграл по массивам с scipy.trapz (без функции analyti c) - PullRequest
Новая
       97

Двойной интеграл по массивам с scipy.trapz (без функции analyti c)

0 голосов
/ 15 января 2020

Я должен объединить два массива (f (x) и g (x + r)) и функцию wfg(r) = tri angular function со следующим интегралом:

enter image description here

с аналитическими c формами f (x) и g (x) неизвестны.

Моя первоначальная попытка состояла в том, чтобы интегрировать f (x) g (x + r) по dx, используя scipy.trapz, а затем умножить этот результат на интеграл wfg(r) относительно dr: 

for ri in r:
    if np.abs(ri) < l:
        term1 = lambda ri, l : (1 - np.abs(ri)/l) 
        tmp1 += integrate.quad(term1, r[-1],r[0] ,args = (l))[0]
    else:
        tmp1 += 0

return np.abs(tmp1 * integrate.trapz(gpeaks1*gpeaks2_r,x_fg))

где интеграл от wfg(r) = tmp1, f(x) = gpeaks1 и g(x+r) = gpeaks2.

Но математически это неверно. Как бы я go об этом?

1 Ответ

1 голос
/ 15 января 2020

Хорошо, я понял следующее (поправьте меня, если я ошибаюсь):

Что у вас есть:

  • Массив x_fg с формой (x_size, 1)
  • Массив r_fg с формой (1, r_size)
  • Массив gpeaks1[i] = f(x_fg[i]) с формой (x_size, 1)
  • Массив gpeaks2[i, j] = g(x_fg[i] + r_fg[j]) с формой (x_size, r_size)

В качестве промежуточного результата у вас также есть:

  • Массив h[j] = ∫_i gpeaks1[i] * gpeaks2[i, j] с формой (1, r_size)
    • Вычислено с h = integrate.trapz(gpeaks1 * gpeaks2, x_fg, axis=-2)

О w_fg: после вашего кода к букве w_fg - это число , а не массив . Но ваш вопрос, по-видимому, указывает на то, что w_fg - это функция веса, и поэтому она должна быть представлена ​​в виде массива.

Итак, мой первый вопрос: "что такое w_fg?"

Предполагая:

  • Массив w_fg[j] = w(r[j]) с формой (1, r_size)

Мой второй вопрос: «Что вы хотите рассчитать?»

Если это ∫_j w_fg[j] h[j], тогда просто сделайте integrate.trapz(w_fg * h).

Предыдущий ответ, который фактически не отвечает на вопрос

Итак, учитывая два одномерных массива f и g представляющих функции, вам нужен массив h, представляющий h(r)=∫f(x)g(x+r)dx.

. Вы можете полностью сделать это с помощью scipy.trapz, но интегралы вида h(r)=∫f(x)g(x+r)dx могут быть намного легче вычислены с помощью np.correlate или эквивалент в scipy.

Наиболее эффективным способом является использование БПФ. Сейчас у меня нет времени, чтобы написать объяснение, но я добавлю его позже.

Предлагаю вам прочитать о свертках и взаимных корреляциях .

...