Как выполнить интеграцию (трапецию) в Python с x ^ 2?

Question

Моя задача - сделать сначала интеграцию, а затем интеграцию трапеции с Python из f (x) = x ^ 2

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(-10,10)   
y = x**2

l=plt.plot(x,y)
plt.show(l)

Now I want to integrate this function to get this: F(x)=(1/3)x^3 with the picture:

This should be the output in the end:

введите описание изображения здесь

Может ли кто-нибудь объяснить мне, как получить первообразную F (x) от f (x) = x ^ 2 с python? Я хочу сделать это с помощью нормальной интеграции и интеграции трапеции. Для трапецеидального интегрирования от (-10 до 10) и размера шага 0,01 (ширина трапеций). В итоге я хочу получить функцию F (x) = (1/3) x ^ 3 в обоих случаях. Как я могу этого достичь?

Спасибо за помощь

Answer 1

Есть два ключевых наблюдения:

• правило трапеции относится к числовому c интегрированию, выход которого не является интегральной функцией, а числовой.
• интегрирование зависит от произвольная константа, которая не включена в ваше определение F(x)

Имея это в виду, вы можете использовать scipy.integrate.trapz() для определения интегральной функции:

import numpy as np
from scipy.integrate import trapz


def numeric_integral(x, f, c=0):
    return np.array([sp.integrate.trapz(f(x[:i]), x[:i]) for i in range(len(x))]) + c

или, что более эффективно, используя scipy.integrate.cumtrapz() (который выполняет вычисления сверху):

import numpy as np
from scipy.integrate import cumtrapz


def numeric_integral(x, f, c=0):
    return cumtrapz(f(x), x, initial=c) 

Это графики, как показано ниже:

import matplotlib.pyplot as plt


def func(x):
    return x ** 2


x = np.arange(-10, 10, 0.01)
y = func(x)
Y = numeric_integral(x, func)

plt.plot(x, y, label='f(x) = x²')
plt.plot(x, Y, label='F(x) = x³/3 + c')
plt.plot(x, x ** 3 / 3, label='F(x) = x³/3')
plt.legend()

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

участок

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

Answer 2

Функция cumtrapz от scipy предоставит первообразную с использованием трапециевидной интеграции:

from scipy.integrate import cumtrapz
yy = cumtrapz(y, x, initial=0)

# make yy==0 around x==0 (optional)
i_x0 = np.where(x >= 0)[0][0]
yy -= yy[i_x0]

Answer 3

Трапецеидальное интегрирование

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(-10, 10, 0.1)
f = x**2

F = [-333.35]
for i in range(1, len(x) - 1):
    F.append((f[i] + f[i - 1])*(x[i] - x[i - 1])/2 + F[i - 1])
F = np.array(F)

fig, ax = plt.subplots()

ax.plot(x, f)
ax.plot(x[1:], F)

plt.show()

Здесь я применил теоретическую формулу (f[i] + f[i - 1])*(x[i] - x[i - 1])/2 + F[i - 1], а интегрирование выполняется в блоке:

F = [-333.35]
for i in range(1, len(x) - 1):
    F.append((f[i] + f[i - 1])*(x[i] - x[i - 1])/2 + F[i - 1])
F = np.array(F)

Обратите внимание на , что для построения x и F они должны иметь одинаковое количество элементов; поэтому я игнорирую первый элемент x, поэтому они оба имеют элемент 199. Это результат метода трапеции d: если вы интегрируете массив f из n элементов, вы получите массив F из n-1 элементов. Более того, я установил начальное значение F на -333.35 на x = -10, это произвольная константа из процесса интегрирования, я решил это значение, чтобы передать функцию около начала координат.

Аналитическая интеграция

import sympy as sy
import numpy as np
import matplotlib.pyplot as plt

x = sy.symbols('x')
f = x**2
F = sy.integrate(f, x)

xv = np.arange(-10, 10, 0.1)
fv = sy.lambdify(x, f)(xv)
Fv = sy.lambdify(x, F)(xv)

fig, ax = plt.subplots()

ax.plot(xv, fv)
ax.plot(xv, Fv)

plt.show()

Здесь я использую модуль symboli c math через sympy. Интеграция выполняется в блоке:

F = sy.integrate(f, x)

Обратите внимание , что в этом случае F и x уже имеют одинаковое количество элементов. Тем более что код проще.