Расчет площади под математической функцией

Question

У меня есть диапазон данных, которые я аппроксимировал, используя полином степени 2 в Python. Я хочу вычислить площадь под этим полиномом от 0 до 1.

Есть ли исчисление или подобный пакет из numpy, который я могу использовать, или мне просто сделать простую функцию для интеграции этих функций?

Мне немного неясно, каков наилучший подход к определению математических функций.

Спасибо.

Answer 1

Если вы интегрируете только полиномы, вам не нужно представлять общую математическую функцию, используйте numpy.poly1d, который имеет метод integ для интеграции.

>>> import numpy
>>> p = numpy.poly1d([2, 4, 6])
>>> print p
   2
2 x + 4 x + 6
>>> i = p.integ()
>>> i
poly1d([ 0.66666667,  2.        ,  6.        ,  0.        ])
>>> integrand = i(1) - i(0) # Use call notation to evaluate a poly1d
>>> integrand
8.6666666666666661

Для интеграции произвольных числовых функций вы должны использовать scipy.integrate с обычными функциями Python для функций. Для аналитической интеграции функций вы должны использовать sympy. В данном случае это звучит не так, как вы хотите, особенно в последнем.

Answer 2

Смотри, Ма, нет импорта!

>>> coeffs = [2., 4., 6.]
>>> sum(coeff / (i+1) for i, coeff in enumerate(reversed(coeffs)))
8.6666666666666661
>>>

Наша гарантия: работает для полинома любой положительной степени или возврата ваших денег!

Обновление от нашей исследовательской лаборатории: гарантия продлена; с / положительный / неотрицательный /: -)

Обновление Вот версия промышленного уровня, которая является надежной в отношении случайных целых чисел в коэффициентах без вызова функции в цикле и не использует ни enumerate(), ни reversed() в установке :

>>> icoeffs = [2, 4, 6]
>>> tot = 0.0
>>> divisor = float(len(icoeffs))
>>> for coeff in icoeffs:
...     tot += coeff / divisor
...     divisor -= 1.0
...
>>> tot
8.6666666666666661
>>>

Answer 3

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

У вас есть многочлен степени 2: f (x) = ax ² + bx + c

Вы хотите найти область под кривой для x в диапазоне [0,1] .

Антидериватив F (x) = ax ³ / 3 + bx ² / 2 + cx + C

Площадь под кривой от 0 до 1: F (1) - F (0) = a / 3 + b / 2 + c

Так что, если вы рассчитываете площадь только для интервала [0,1] , вы можете рассмотреть используя это простое выражение, а не прибегая к методам общего назначения.

Answer 4

Если кто-то интегрирует квадратичные или кубические полиномы с самого начала, альтернативой получению явных интегральных выражений является использование правила Симпсона; Это глубокий факт, что этот метод точно объединяет полиномы степени 3 и ниже.

Чтобы позаимствовать пример Майка Грэма (я давно не использовал Python; извиняюсь, если код выглядит не очень хорошо):

>>> import numpy
>>> p = numpy.poly1d([2, 4, 6])
>>> print p
   2
2 x + 4 x + 6
>>> integrand = (1 - 0)(p(0) + 4*p((0 + 1)/2) + p(1))/6

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

Конечно, я не упростил выражение для integrand, чтобы указать, что 0 и 1 можно заменить произвольными значениями u и v, и код все равно будет работать для поиска интеграл функции от u до v.

Answer 5

' quad ' в scipy.integrate - это метод общего назначения для интегрирования функций одной переменной за определенный интервал. В простом случае (например, описанном в вашем вопросе) вы передаете свою функцию и нижний и верхний пределы соответственно. 'quad' возвращает кортеж, состоящий из интегрального результата и верхней границы члена ошибки.

from scipy import integrate as TG

fnx = lambda x: 3*x**2 + 9*x    # some polynomial of degree two
aoc, err = TG.quad(fnx, 0, 1)

[Примечание: после того, как я опубликовал это, я опубликовал ответ перед моим, который представляет полиномы с использованием 'poly1d' в Numpy. Мой скриптлет чуть выше также может принимать полином в такой форме:

import numpy as NP

px = NP.poly1d([2,4,6])
aoc, err = TG.quad(px, 0, 1)
# returns (8.6666666666666661, 9.6219328800846896e-14)