Численное интегрирование с использованием заданных трапециевидных ящиков - PullRequest
Новая
       20

Численное интегрирование с использованием заданных трапециевидных ящиков

0 голосов
/ 11 июля 2019

Я написал функцию, которая выполняет численное интегрирование, используя определенные ячейки для значений x. Другими словами, он использует правило трапеции, используя размер корзины в качестве основы трапеции. Для большинства размеров бина возвращаемое значение является неправильным.

Функция принимает массив x, массив y, минимальное значение x, максимальное значение x и размер корзины. Для очень маленького размера ячейки (<.0001) функция фактически дает правильные интегральные значения, но почти для всего остального значение является неправильным. Например, если вы использовали правило трапеции вручную для y = x ^ 2 от -2 до 1 с размером ячейки 1, вы должны получить 3,5. Мне нужна эта функция для работы со всеми размерами бинов. </p> 

def integrate_bins(f_X, f_Y, bound_min, bound_max, binSize):

    # get bound indices (ASSUMES SORTED)
    bound_lower_idx = f_X.searchsorted(bound_min, 'left')
    bound_upper_idx = f_X.searchsorted(bound_max, 'right')

    # trim arrays
    f_X = f_X [ bound_lower_idx : bound_upper_idx ]
    f_Y = f_Y [ bound_lower_idx : bound_upper_idx ]

    if np.amin(f_X) > bound_min:
        bound_min = np.amin(f_X)
    if np.amax(f_X) < bound_max:
        bound_max = np.amax(f_X)

    # interpolate along binSize
    f_X_new = [i * binSize + bound_min for i in range(0, int((bound_max - bound_min)/binSize))]
    f_Y_new = [np.interp(x, f_X, f_Y) for x in f_X_new]

    # integrate
    return np.trapz(f_Y_new, f_X_new)


#### test case
testx = np.arange(-3,4,.25)
testy = testx**2
print(integrate_bins(testx, testy, -2,1,1))

Только для этого примера результат этого теста равен 3,0, тогда как он должен быть 3,5. Величина ошибки также изменяется в зависимости от размеров ячейки и алгебраической функции, приведенной в контрольном примере.

...