Числовой интеграл от греха (х) / х

Question

Я хочу вычислить определенный интеграл с квадратурами sin (x) / x в python, используя scipy.С n = 256. Кажется, он не работает хорошо:

 from scipy import integrate

 exact = integrate.quad(lambda x : (np.sin(x))/x, 0, 2*np.pi)[0]
 print("Exact value of integral:", exact)

 # Approx of sin(x)/x by Trapezoidal rule
 x = np.linspace(0, 2*np.pi, 257)
 f = lambda x : (np.sin(x))/x
 approximation = np.trapz(f(x), x)
 print ("Estimated value of trapezoidal O(h^2):", round(approximation, 5), 
   '+', round((2/256)**2, 5))
 print ("real error:", exact - approximation)

 # Approx of sin(x)/x by Simpsons 1/3 rule
 approximation = integrate.simps(f(x), x)
 print("Estimated value of simpsons O(h^4):", round(approximation, 9), 
   '+', round((2/256)**4, 9))
 print ("real error:", exact - approximation)

 plt.figure()
 plt.plot(x, f(x))
 plt.show()

Точный интеграл рассчитан хорошо, но я получаю ошибку с квадратурами ... Что не так?

Exact value of integral: 1.4181515761326284
Estimated value of trapezoidal O(h^2): nan + 6e-05
real error: nan
Estimated value of simpsons O(h^4): nan + 4e-09
real error: nan

Заранее спасибо!

Answer 1

NaN означает «не число».В твоем случае это в основном бесконечность.Когда вы создаете:

x = np.linspace(0, 2*np.pi, 257)

Вы создаете массив со значением 0, затем вы пытаетесь разделить на x, и вы не можете делить на 0 ...

Oneрешение состоит в том, чтобы использовать это:

x = np.linspace(0.1, 2*np.pi, 257)

Что дает вам это:

Exact value of integral: 1.4181515761326284
Estimated value of trapezoidal O(h^2): 1.31822 + 6e-05
real error: 0.099935104987
Estimated value of simpsons O(h^4): 1.318207115 + 4e-09
real error: 0.0999444614012

Чем ближе вы будете к нулю, тем лучше будет приближение!

Answer 2

Вы могли бы заставить функцию возвращать 1 (предел sin (x) / x в 0) вместо NaN для x == 0. Таким образом, вам не придется обманывать и изменять интервал, на которомВы интегрируете, чтобы исключить 0.

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

exact = integrate.quad(lambda x : (np.sin(x))/x, 0, 2*np.pi)[0]
print("Exact value of integral:", exact)


def f(x):
    out = np.sin(x) / x
    # For x == 0, we get nan. We replace it by the 
    # limit of sin(x)/x in 0
    out[np.isnan(out)] = 1
    return out

# Approx of sin(x)/x by Trapezoidal rule

x = np.linspace(0, 2*np.pi, 257)

approximation = np.trapz(f(x), x)
print ("Estimated value of trapezoidal O(h^2):", round(approximation, 5), 
  '+', round((2/256)**2, 5))
print ("real error:", exact - approximation)
 # Approx of sin(x)/x by Simpsons 1/3 rule
approximation = integrate.simps(f(x), x)
print("Estimated value of simpsons O(h^4):", round(approximation, 9), 
  '+', round((2/256)**4, 9))
print ("real error:", exact - approximation)
plt.figure()
plt.plot(x, f(x))
plt.show()

Вывод:

Exact value of integral: 1.4181515761326284
Estimated value of trapezoidal O(h^2): 1.41816 + 6e-05
real error: -7.98955129322e-06
Estimated value of simpsons O(h^4): 1.418151576 + 4e-09
real error: 2.72103006793e-10

Answer 3

В вашем коде есть по крайней мере некоторые проблемы:

1. Ваш пробел начинается с 0, поэтому при оценке функции для интеграции в началетрапецеидальный интеграл, у вас есть: sin(0)/0 = nan.Вы должны использовать числовой ноль вместо точного нуля (в приведенном ниже примере я использовал 1e-12)
2. Когда вы получите первый nan, nan + 1.0 = nan: это означает, что в вашем коде, когда вы суммируете интеграл по интервалам, первый nan запутывает все ваши результаты.
3. только для Python 2 : Деление 2/256 - это деление между двумя целыми числами, и результат равен 0.Вместо этого попробуйте 2.0/256.0 (спасибо @MaxU за указание на это).

Это ваш код исправлен (я запускаю его в python2, это то, что установлено на компьютере, яиспользуя сейчас):

from scipy import integrate
import numpy as np

exact = integrate.quad(lambda x : (np.sin(x))/x, 0, 2*np.pi)[0]
print("Exact value of integral:", exact)

# Approx of sin(x)/x by Trapezoidal rule
x = np.linspace(1e-12, 2*np.pi, 257) # <- 0 has become a numeric 0
f = lambda x : (np.sin(x))/x
approximation = np.trapz(f(x), x)
print ("Estimated value of trapezoidal O(h^2):", round(approximation, 5), 
  '+', round((2.0/256.0)**2, 5))
print ("real error:", exact - approximation)

# Approx of sin(x)/x by Simpsons 1/3 rule
approximation = integrate.simps(f(x), x)
print("Estimated value of simpsons O(h^4):", round(approximation, 9), 
  '+', round((2/256)**4, 9))
print ("real error:", exact - approximation)

с его выводом:

('Exact value of integral:', 1.4181515761326284)
('Estimated value of trapezoidal O(h^2):', 1.41816, '+', 6e-05)
('real error:', -7.9895502944626884e-06)
('Estimated value of simpsons O(h^4):', 1.418151576, '+', 0.0)
('real error:', 2.7310242955991271e-10)

Discalimer Предел sin(x)/x -> 1 for x -> 0, но из-за плавающего округления для sin(1e-12)/1e-13 = 1!