Блог

У меня возникли проблемы с пониманием np.trapz. Я должен сам написать трапециевидное правило, а затем сравнить его с np.trapz Однако есть одна загвоздка. Скажем, интеграл равен от a до b. Я должен найти интеграл для a = 1 b = 1, a = 0 b = 2, a = 0 b = 3… a = 0 b = 10 и постройте график этих значений. Вот мой код для трапециевидной функции, которую я создал:

 ## trapezoidal ##
import numpy as np
import matplotlib.pyplot as plt

def f(x):
    o1 = 0.3 ## matter density of universe
    o2 = 0.7 ## density of dark enerrgy
    o3 = 0 ## density of radiation
    c = 3*10**3 ## constant
    return c/((o1*(1 x)**3   o3*(1 x)**2   o2)**(1/2))

data = [] ## data array

for b in range (0, 10):## for loop going through z2 as b

    ## definitions 
    a = 0
    N = 100
    h = (b-a)/N
    integral = 0

    integral = 0.5*(f(a)   f(b)) ## endpoints

    for i in range(1,N):
        integral  = f(a   (i*h))
    integral = integral*h
    integral = integral/(1 b)
    data.append(integral) ## appending each itteration into data array

 print(data)

plt.plot([1,2,3,4,5,6,7,8,9,10],data)
plt.xlabel('z')
plt.ylabel('DA')
  

и вот что я попытался для np.trapz, хотя я думаю, что я безнадежно ошибаюсь.
## np.trapz ##
импортировать numpy как np
импортировать matplotlib.pyplot как plt

 def f(x):
    o1 = 0.3 ## matter density of universe
    o2 = 0.7 ## density of dark enerrgy
    o3 = 0 ## density of radiation
    c = 3*10**3 ## constant
    return c/((o1*(1 x)**3   o3*(1 x)**2   o2)**(1/2))

x = np.arange(1,10,1)
##area = []

for i in range (0,10):
    x = [0,i]
    y = f/(1 i)
    area.append(np.trapz(y,x))

print(area)

plt.plot([1,2,3,4,5,6,7,8,9,10],area)
plt.xlabel('z')
plt.ylabel('DA')
  

Я попытаюсь показать, как реализовать правило трапеции, и сравнить его с его реализацией numpy . Итак, у вас есть функция def f(x) , которую вы хотите интегрировать из a в b .

Вы ничего не указали о принятии ошибок / допусков, но при сравнении различных методов численной оценки интегралов (.., уравнений или аналогичных) часто часто сравнивается масштаб их ошибок при увеличении размера системы, количества шагов итерации и так далее.

Правило трапеции является распространенным методом учебника, и существует множество бесплатной литературы (например, в Википедии). Следовательно, вы легко проверяете правильность своей реализации, потому что ее поведение уже хорошо известно. Кроме того, всегда полезно протестировать свой метод на задаче, которую можно решить аналитически (с помощью ручки и бумаги). Это быстро выявит ошибки в вашем коде / реализации.

Итак, давайте сделаем это! Рассмотрим тривиальную функцию g(x) = x . Интеграция g(x) из 0 в 10 легко выполняется и приводит к (1/2) * (10^2 - 0^2) = 50 . Вот g(x) в python,

 def g(x):
  return x
  

Теперь мы реализуем функцию для вычисления интеграла от любой функции func , используя правило трапеции, как определено в Википедии (ссылка выше),

 def my_trapz(func, a, b, n_steps):
  X = np.linspace(a,b,n_steps)
  integral = (func(a) func(b))/2.0   sum([func(x) for x in X])
  return integral * (b-a)/n_steps
  

Здесь func находится функция python ( def ), которая принимает ровно один аргумент. a является нижней и b верхней границей интеграла, в то время n_steps как количество значений x для оценки в пределах диапазона [a,b] . Чем больше n_steps , тем точнее будет интеграл. numpy.linspace(a,b,n) Функция создает массив n линейно расположенных чисел между a и b . теперь мы можем вызвать,

 a = 0.0
b = 10.0
n_steps = 100

integral = my_trapz(g, a, b, n_steps)
print(integral) # outputs 50.4999999..
  

Вы заметите, что чем выше значение, которое вы вводите n_steps , тем ближе результат 50 , как и ожидалось. Теперь мы хотим сравнить нашу реализацию с той numpy , что есть. Прочитав его документацию, мы видим, что интеграл вычисляется по,

 X = np.linspace(a,b,n_steps)

integral = np.trapz([g(x) for x in X], dx=(b-a)/n_steps)
print(integral) # outputs 49.9023437.., slightly better than out implementation
  

Здесь мы использовали понимание списка непосредственно в вызове функции. В качестве альтернативы, мы могли бы вместо этого сделать

 g_values = []
for x in :
  g_values.append(g(x))
integral = np.trapz(g_values, dx=(b-a)/n_steps)
  

Мы можем сравнить их, согласовав набор различных n_steps и исследуя, как работают два метода: для каждого значения n_steps в списке n_steps_list вычислите соответствующий интеграл, используя оба numpy и наш собственный метод. Затем постройте интегральный результат как функцию n_steps . Это,

 a = 0.0
b = 10.0
n_steps_list = [2**n for n in range(3,10)] # =[8, 16, 32, ..., 1024]

integral_np = []
integral_my = []

for n_steps in n_steps_list:
  X = np.linspace(a,b,n_steps)
  dx = (b-a)/n_steps
  integral_np.append(np.trapz([g(x) for x in X], dx=dx))
  integral_my.append(my_trapz(g, a, b, n_steps))

plt.scatter(n_steps_list, integral_np, color='g', label='np')
plt.scatter(n_steps_list, integral_my, color='r', label='my')
plt.xlabel('number of steps (resolution)')
plt.ylabel('Integral result')
plt.legend()
plt.show()
  

Результирующий график показан ниже,

Мы видим, что оба метода сходятся экспоненциально к ожидаемому значению, 50 . Это согласуется с анализом ошибок, представленным в статье Википедии.

Теперь я предлагаю вам попробовать повторить этот процесс, но заменив мою тривиальную g(x) на вашу исходную функцию f(x) и выяснить, что происходит.

 def f(x):
    o1 = 0.3 ## matter density of universe
    o2 = 0.7 ## density of dark enerrgy
    o3 = 0 ## density of radiation
    c = 3*10**3 ## constant
    return c/((o1*(1 x)**3   o3*(1 x)**2   o2)**(1/2))