Vectorized operations in NumPy

Question

I have a numpy array of 5 numpy arrays, each with 5 observations. I need to use the trapezoid rule of integration to get 5 results of integrating. I'm trying to use for loop for apply a function to each numpy array and looking for something that has faster implementation.

def trapezoid_rule(a,b,n,f_nparray):
  h = (b-a) / (n-1)
  return (h/2) * (f_nparray[0] + 2 * sum(f_nparray[1:n-1]) + f_nparray[n-1])

simulation_res = np.array([[1,2,3,4,5], [1,3,5,7,9], [1,4,2,5,7], [1,5,2,6,4], [6,2,5,3,4]])

# List of integration results
I = []

for i in simulation_res:
  I.append(trapezoid_rule(0,1,10,i))

Expected output format = [a,b,c,d,e]

try something like (f_nparray[:, [0]] + 2 * sum(f_nparray[:, 1:n-1]) + f_nparray[:, [n-1]]) - applying the trapezoid_rule function to all rows at once. The trickiest part is getting the 3 terms to work together, with the rules of broadcasting. shapes: (n,1), (n,m), (n,1) — hpaulj
– hpaulj, Commented Oct 10, 2022 at 16:12

t.o. · Accepted Answer · 2022-10-10 18:17:13Z

1

You can do this without for loops, but I'll show both. numpy has a built-in trapezoidal rule that works in a vectorized way. If you wanted to use it:

import numpy as np

simulation_res = np.array([[1,2,3,4,5], 
                           [1,3,5,7,9], 
                           [1,4,2,5,7], 
                           [1,5,2,6,4], 
                           [6,2,5,3,4]])
result = np.trapz(simulation_res, axis=1)
print(result) 
# Gives [12.  20.  15.  15.5 15. ]

However if the loop is required for some reason, you could just preallocate some memory in a numpy array and append it

result = np.zeros(5)
for idx, i in enumerate(simulation_res):
    result[idx] = np.trapz(i)
print(result)
# Gives [12.  20.  15.  15.5 15. ]

edited Oct 10, 2022 at 18:17

answered Oct 10, 2022 at 17:30

t.o.

8897 silver badges22 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

ask1504 Over a year ago

How do I provide the limits of integral in the first approach? I'm trying -

t.o. Over a year ago

From what I understand, this numerical quadrature method uses the dx of the array. To specify limits, you would either need a vector function f(x) or know the indices of the start and stop points in the array when you pass it to trapz. If you want to utilize the limits, look at quad_vec. Alternatively @hpaulj gave an answer that works for your question.

hpaulj · Accepted Answer · 2022-10-10 19:45:39Z

With:

def foo(a,b,n,f_nparray):
  h = (b-a) / (n-1)
  return (h/2) * (f_nparray[:,0] + 2 * np.sum(f_nparray[:,1:n-1], axis=1) + f_nparray[:,n-1])

In [20]: foo(0,1,5,simulation_res)
Out[20]: array([3.   , 5.   , 3.75 , 3.875, 3.75 ])

which matches your I.

I had to add axis to np.sum; and in contrast to my comment, just use 0 for the indexing. I just need one value per row of the 2d array.

Compared to np.trapz:

In [29]: np.trapz(simulation_res, axis=1)/4
Out[29]: array([3.   , 5.   , 3.75 , 3.875, 3.75 ])

In [30]: timeit np.trapz(simulation_res, axis=1)/4
34 µs ± 1.66 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [31]: timeit foo(0,1,5,simulation_res)
25.7 µs ± 139 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

Collectives™ on Stack Overflow

Vectorized operations in NumPy

2 Answers 2

2 Comments

Comments

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

2 Comments

Comments

Related