Different options for numerical integration
In PyWake we need to integrate over the wind speed and wind directions to compute the AEP and Lifetime DEL. This is currently done by summing over the elements, but there are many other possibilities, as shown in the following snippet
import numpy as np
from scipy.integrate import dblquad
from scipy.interpolate import interp2d
def fun(x, y):
return 5.0 * x ** 2 * np.cos(x * y)
x = np.linspace(-1, 1, 1000)
y = np.linspace(-2, 2, 2500)
dx = np.mean(np.diff(x))
dy = np.mean(np.diff(y))
X, Y = np.meshgrid(x, y)
Z = fun(X, Y)
# Call qagse from QUADPACK. Adaptive method.
int_dblquad, _ = dblquad(fun, -2.0, 2.0, -1.0, 1.0)
# Same as before, but we use linear interpolation. Equivalent to the trapezoidal method.
fun_interp = interp2d(x, y, Z, kind='linear', copy=False)
int_dblquad_interp, _ = dblquad(fun_interp, -2.0, 2.0, -1.0, 1.0)
# 2 ways of applying the trapezoidal method over a non-uniform grid.
int_trapz_1 = np.trapz(np.trapz(Z, x, axis=1), y)
int_trapz_2 = np.trapz(np.trapz(Z, y, axis=0), x)
# 2 ways of applying the trapezoidal method over a uniform grid.
int_trapz_3 = np.trapz(np.trapz(Z, dx=dx, axis=1), dx=dy)
int_trapz_4 = np.trapz(np.trapz(Z, dx=dy, axis=0), dx=dx)
# Full factorial method.
area = 4.0 * 2.0
int_ff = area * np.mean(Z)
I think that it would be nice to provide the user with some more options
cc @mmpe