9.4: Квадратури Гауса

Раніше ми припускали в нашому аналізі схем числового інтегрування, що точки дискретизації однаково розташовані. Це припущення не є суворо необхідним; наприклад, ми можемо легко змінити формули правила середньої точки та правила трапеції для роботи з неоднаково розташованими точками.

Однак якщо ми вільні вибирати точки дискретизації для виконання числового інтегрування, і ці точки не повинні бути однаково розташовані, то можна використовувати цю свободу для подальшого підвищення точності числового інтегрування. Ідея полягає в оптимізації розміщення точок дискретизації, щоб мінімізувати результуючу числову похибку. Це основна ідея методу інтеграції за гаусовими квадратурами.

Ми не будемо обговорювати деталі цього методу числового інтеграції. Для його використання можна викликати функцію scipy.integrate.quad. Ця функція використовує низькорівневу числову бібліотеку з назвою QUADPACK, яка виконує квадратурну інтеграцію з адаптивними квадратурами, тобто автоматично з'ясовує, скільки точок дискретизації слід використовувати, і де вони повинні бути розташовані, щоб отримати результат з бажана числова точність.

Оскільки QUADPACK з'ясовує точки дискретизації для себе, ви повинні передати quad функцію, що представляє integrand, а не масив integrand значень, як з trapz або simps. Стандартним способом виклику функції є:

t = quad(f, a, b)

який обчислює інтеграл

$$t = \int_a^b f(x) dx.$$

Повертається значення є кортежем виду (t, err), де t - значення інтеграла, а err - оцінка числової похибки. Функція quad також приймає багато інших необов'язкових входів, за допомогою яких можна вказати додаткові входи для передачі до функції integrand, допуску до помилок, кількості використовуваних підінтервалів тощо. Докладніше про це можна дізнатися з документації.

Вибір того, чи виконувати числове інтегрування за допомогою правила Сімпсона (симпсів) або квадратур Гаусса (quad), є ситуативним. Якщо ви вже знаєте значення integrands в заздалегідь визначеному наборі точок дискретизації (наприклад, з результату обчислення скінченної різниці), то використовуйте симпс. Якщо ви можете визначити функцію, яка може швидко обчислити значення integrand в будь-якій точці, використовуйте quad.

Ось приклад використання quad для обчислення інтеграла$\int_{0}^\infty \frac{dx}{x^2 + 1}$:

from scipy import *
from scipy.integrate import quad

def f(x):
    return 1.0 / (x*x+1)

integ, _ = quad(f, 0, 1000)
print(integ)

(Зауважте, що quad повертає два значення; перше - обчислене значення інтеграла, а інше - абсолютна помилка, яка нас не цікавить, тому ми кидаємо його в «викинутої» змінної з ім'ям _. Докладніше дивіться в документації.) Запуск вищевказаної програми виводить результат,$1.569796\dots,$ який відрізняється від точного результату,$\pi/2 = 1.570796\dots,$ відносною помилкою$0.06\%$.

Ось ще один приклад, де integrand приймає додатковий параметр$\int_0^\infty x e^{-\lambda x}\, dx$:

from scipy import *
from scipy.integrate import quad

def f(x, lambd):
    return x * exp(-lambd * x)

integ, _ = quad(f, 0, 100, args=(0.5))
print(integ)

Запущена програма виводить результат$4.0$, який узгоджується з точним результатом$1/\lambda^2$ для обраного значення параметра$\lambda =0.5$.