CODEWITHSUNDEEP
pythonmathscipynumerical-integration
Slayre
Slayre

Reputation: 25

Double integral over arrays w/ scipy.trapz (no analytic function)

I must integrate two arrays (f(x) and g(x+r)) and a function wfg(r) = triangular function with the following integral:

enter image description here

with the analytic forms of f(x) and g(x) unknown.

My initial attempt was to integrate f(x)g(x+r) w.r.t dx using scipy.trapz and then multiply this result by the integral of wfg(r) w.r.t 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))

where the integral of wfg(r) = tmp1, f(x) = gpeaks1, and g(x+r) = gpeaks2.

But this doesn't seem correct mathematically. How would I go about this?

Upvotes: 0

Views: 424

Answers (1)

Maxpxt
Maxpxt

Reputation: 189

Ok, so what I understood is this (correct me if I'm wrong):

What you have:

  • Array x_fg with shape (x_size, 1)
  • Array r_fg with shape (1, r_size)
  • Array gpeaks1[i] = f(x_fg[i]) with shape (x_size, 1)
  • Array gpeaks2[i, j] = g(x_fg[i] + r_fg[j]) with shape (x_size, r_size)

As an intermediate result, you also have:

  • Array h[j] = ∫_i gpeaks1[i] * gpeaks2[i, j] with shape (1, r_size)
    • Computed with h = integrate.trapz(gpeaks1 * gpeaks2, x_fg, axis=-2)

About w_fg: following your code to the letter, w_fg is a number, not an array. But your question seems to indicate w_fg is a weight function and, as such, should be represented by an array.

So, my fist question is "what is w_fg?"

Assuming:

  • Array w_fg[j] = w(r[j]) with shape (1, r_size)

My second question is "what do you want to calculate?"

If it's ∫_j w_fg[j] h[j] then just do integrate.trapz(w_fg * h).

Previous answer, which doesn't actually answer the question

So, given two 1D arrays f and g representing functions, you want the array h which represents h(r)=∫f(x)g(x+r)dx.

You can totally do this with scipy.trapz, but integrals of the form h(r)=∫f(x)g(x+r)dx can be much more easily computed with np.correlate or an equivalent in scipy.

The most efficient way is by using FFTs. I don't have time right now to write an explanation, but I'll add one later.

I suggest you read about convolutions and cross-correlations.

Upvotes: 1

Related Questions