Python Riemann Sum не дает нуля для равных положительных и отрицательных площадей

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

Следующий код отображает нечетную функцию f(x) = x^3 от -1 до 1, поэтому площадь должна быть равна нулю. Вместо этого мой код приближает его к 1,68065561477562 e^-15.

Чем это вызвано? Это ошибка округления в delta_x, x или y? Я знаю, что могу просто округлить значение до нуля, но мне интересно, есть ли другая проблема или способ решить эту проблему.

Я попытался использовать класс Decimal.decimal для delta_x, но получил еще меньшее число.

Код Python:

import matplotlib.pyplot as plt
import numpy as np

# Approximates and graphs integral using Riemann Sum


# example function: f(x) = x^3
def f_x(x):
    return x**3

# integration range from a to b with n rectangles
a, b, n = -1, 1, 1000

# calculate delta x, list of x-values, list of y-values, and approximate area under curve
delta_x = (b - a) / n

x = np.arange(a, b+delta_x, delta_x)

y = [f_x(i) for i in x]

area = sum(y) * delta_x

# graph using matplotlib
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x, y)
ax.bar(x, y, delta_x, alpha=.5)
plt.title('a={}, b={}, n={}'.format(a, b, n))
plt.xlabel('A = {}'.format(area))
plt.show()

python matplotlib integral calculus
person halfmpty    schedule 17.11.2017    source источник

Ответы (3)


arrow_upward
1
arrow_downward

Вы должны знать, что то, что вы вычисляете, не является интегралом Римана в первоначальном смысле. Вы делите интервал на n интервалов, а затем суммируете по n+1 интервалам (здесь n = 1000, но len(x) == 1001). Таким образом, результат может быть близок к тому, что вы ожидаете, но это, безусловно, не лучший способ его получить.

Используя сумму Римана, вы разделите свой интервал на n интервалов, а затем просуммируете значения эти n контейнеров. У вас есть выбор: вычислять левую сумму Римана, правую сумму Римана или, возможно, брать средние точки.

import numpy as np

def f_x(x):
    return x**3

# integration range from a to b with n rectangles
a, b, n = -1, 1, 1000

delta_x = (b - a) / float(n)

x_left = np.arange(a, b, delta_x)
x_right = np.arange(a+delta_x, b+delta_x, delta_x)
x_mid = np.arange(a+delta_x/2., b+delta_x/2., delta_x)

print len(x_left),  len(x_right), len(x_mid)          ### 1000 1000 1000


area_left = f_x(x_left).sum() * delta_x
area_right = f_x(x_right).sum() * delta_x
area_mid = f_x(x_mid).sum() * delta_x

print area_left     # -0.002
print area_right    #  0.002
print area_mid      # 1.81898940355e-15

В то время как средняя сумма уже дает хороший результат, для симметричных функций было бы неплохо выбрать n четное и взять среднее значение левой и правой суммы,

print 0.5*(area_right+area_left)   # 1.76204537072e-15

Это одинаково близко к 0.

Теперь стоит отметить, что numpy.arange сам по себе выдает некоторые ошибки. Лучшим выбором было бы использование numpy.linspace

x_left = np.linspace(a, b-delta_x, n)
x_right = np.linspace(a+delta_x, b, n)
x_mid = np.linspace(a+delta_x/2., b-delta_x/2., n)

уступающий

print area_left     # -0.002
print area_right    #  0.002
print area_mid      # 8.52651282912e-17
print 0.5*(area_right+area_left)   # 5.68121938382e-17

5.68121938382e-17 довольно близко к 0. Причина, по которой оно не совсем равно 0, действительно заключается в неточностях с плавающей запятой. .

Известным примером этого может быть 0.1 + 0.2 - 0.3, который дает 5.5e-17 вместо 0. Это должно показать, что эта простая операция вносит ту же ошибку порядка 1e-17, что и интегрирование по Риману.

person ImportanceOfBeingErnest    schedule 17.11.2017
comment
Да, как ни странно, если вы когда-либо используете только нечетное количество подинтервалов, то взятие левой суммы Римана от -1 до 1+delta_x эквивалентно использованию аппроксимации средней точки для интервала от -1-delta_x/2 до 1+delta_x/2 и, таким образом, остается симметричным. - person eugenhu; 17.11.2017

arrow_upward
1
arrow_downward

Да, это просто из-за неточности с плавающей запятой.

Для построения раздела np.linspace() может быть более подходящим, поскольку np.arange() может включать или не включать конечную точку в зависимости от того, как она округляется при использовании нецелочисленных размеров шага.

От numpy.arange() документация:

При использовании нецелочисленного шага, например 0,1, результаты часто не будут согласовываться. Для этих случаев лучше использовать linspace.

person eugenhu    schedule 17.11.2017
comment
В принципе это правильно, хотя код из вопроса имеет еще один недостаток, более серьезный, чем обсуждаемая здесь точность. Смотрите мой ответ. - person ImportanceOfBeingErnest; 17.11.2017

arrow_upward
0
arrow_downward

Ваш код кажется мне действительным. Он суммирует y, возвращенный f_x, и учитывает ошибку аппроксимации наличия только 1000 подинтервалов, добавляя delta_x ко второму аргументу arange. К сожалению, я думаю, что здесь играет роль ошибка округления.

person MattS    schedule 17.11.2017