Scipy.integrate.quad: примеры (PYTHON)

Простое использование функции quad

Пример вычисления интеграла от полинома.

import scipy.integrate as spi
import numpy as np

def f(x):
    return x**2 + 2*x + 1

result, error = spi.quad(f, 0, 2)
print(f"Интеграл: {result}")
print(f"Ошибка: {error}")

Интеграл: 8.666666666666666
Ошибка: 9.621932880084691e-14

Использование аргумента args для передачи параметров.

def g(x, a, b):
    return a * np.exp(-b * x)

result, error = spi.quad(g, 0, np.inf, args=(2.0, 1.0))
print(f"Интеграл: {result}")

Интеграл: 2.0

Работа с бесконечными пределами интегрирования.

result, error = spi.quad(lambda x: np.exp(-x**2), -np.inf, np.inf)
print(f"Интеграл Гаусса: {result:.10f}")

Интеграл Гаусса: 1.7724538509

Использование параметра full_output.

result, error, info = spi.quad(np.sin, 0, np.pi, full_output=1)
print(f"Результат: {result}")
print(f"Количество вычислений функции: {info['neval']}")

Результат: 2.0
Количество вычислений функции: 21

Другие функции интегрирования в SciPy

В модуле scipy.integrate присутствуют и другие функции для численного интегрирования, каждая со своей областью применения.

• quad_vec: Векторизованная версия quad, предназначенная для интегрирования функций, возвращающих векторы (массивы). Полезна для одновременного вычисления нескольких интегралов или интегрирования вектор-функций.
• fixed_quad: Вычисляет интеграл, используя фиксированный порядок квадратуры Гаусса. Быстрее, чем quad, для гладких функций, когда известна необходимая точность, но не адаптируется автоматически.
• quadrature и romberg: Альтернативные адаптивные методы интегрирования. romberg использует экстраполяцию Ричардсона и может быть эффективнее для очень гладких функций без особенностей.
• trapz, simps, romb: Правила интегрирования для уже заданных отсчетов функции (данных), а не callable-объекта. Используются, когда функция известна только в дискретных точках.
• nquad: Функция для вычисления кратных интегралов (двойных, тройных и т.д.). Внутри себя часто использует quad рекурсивно.
• odeint: Решает системы обыкновенных дифференциальных уравнений (ОДУ), что является родственной, но отличной от интегрирования функциональности.

Функция quad является основным универсальным выбором для интегрирования скалярной функции одной переменной, особенно когда функция может быть вычислена в любой точке интервала. Для векторных функций предпочтительнее quad_vec. Если данные представлены в виде массива, используют trapz или simps.

Аналоги функции quad в других языках

Численное интегрирование поддерживается во многих языках программирования, часто через специализированные библиотеки.

JavaScript (библиотека math.js)

const math = require('mathjs');
function f(x) { return Math.exp(-x * x); }
// В math.js нет прямой аналогии quad, часто используют простые методы.
// Пример реализации правила Симпсона:
function integrate(a, b, n) {
    let h = (b - a) / n;
    let sum = f(a) + f(b);
    for (let i = 1; i < n; i++) {
        sum += (i % 2 === 0 ? 2 : 4) * f(a + i * h);
    }
    return sum * h / 3;
}
console.log(integrate(0, 1, 100)); // Приближенное значение интеграла exp(-x^2)

0.7468241328124273

Java (Apache Commons Math)

import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.integration.*;

UnivariateFunction func = new UnivariateFunction() {
    public double value(double x) { return Math.sin(x); }
};

Integrator integrator = new SimpsonIntegrator(); // Есть и адаптивный интегратор
double result = integrator.integrate(1000, func, 0, Math.PI); // 1000 - максимальное количество итераций
System.out.println(result); // Должно быть близко к 2.0

C# (MathNet.Numerics)

using MathNet.Numerics.Integration;

Func f = x => Math.Sin(x);
double result = NewtonCotesTrapeziumRule.IntegrateAdaptive(f, 0, Math.PI, 1e-10);
Console.WriteLine(result); // ~2.0

Golang (без стандартной библиотеки, пример пакета gonum)

// Пример с использованием gonum.org/v1/gonum/integrate/quad
// Код упрощен
// import "gonum.org/v1/gonum/integrate/quad"
// f := func(x float64) float64 { return math.Sin(x) }
// result, err := quad.Fixed(f, 0, math.Pi, 100, nil, 0) // Неадаптивный метод
// Адаптивный метод в gonum: quad.Fixed -> менее автоматизирован, чем scipy.quad

PHP (не имеет широко распространенной стандартной библиотеки для научных вычислений)

Обычно используют пользовательские реализации или расширения вроде PECL's stats.

SQL (не предназначен для таких вычислений)

Численное интегрирование в чистом SQL непрактично и не является стандартной функцией.

Lua

Нет встроенной поддержки, требуются внешние библиотеки, например, через SciLua или GSL-Lua.

Kotlin (может использовать Java-библиотеки, такие как Apache Commons Math)

Синтаксис и подход аналогичны Java.

Основное отличие scipy.integrate.quad от многих аналогов - это высокая степень автоматизации, надежность адаптивного алгоритма и удобная интеграция с экосистемой SciPy для научных расчетов. Многие реализации на других языках требуют явного указания количества разбиений или менее удобны в использовании.

Распространенные ошибки при работе с quad

Ошибки часто возникают из-за непонимания поведения функции или особенностей численных методов.

1. Интегрирование функции, возвращающей массив вместо скаляра.

import scipy.integrate as spi
import numpy as np

def wrong_func(x):
    return np.array([x, x**2]) # Возвращает массив

try:
    result, error = spi.quad(wrong_func, 0, 1)
except Exception as e:
    print(f"Ошибка: {type(e).__name__}: {e}")

Ошибка: TypeError: объект 'numpy.ndarray' не может быть интерпретирован как целое число

2. Интегрирование функции с особенностью (сингулярностью) внутри интервала без указания точки разрыва в параметре points.

def f_with_singularity(x):
    return 1.0 / (np.abs(x - 0.5) ** 0.5) # Особенность в x=0.5

# Без указания points результат может быть неверным или будет ошибка.
result, error = spi.quad(f_with_singularity, 0, 1, points=[0.5], limit=1000)
print(f"Результат с учетом особенности: {result:.6f}")

Результат с учетом особенности: 3.999999

3. Слишком высокая требуемая точность (epsabs, epsrel) для сложной функции, приводящая к большому количеству вычислений и возможному превышению limit.

def oscillatory(x):
    return np.sin(1/(x+0.01))

try:
    result, error = spi.quad(oscillatory, 0, 1, epsabs=1e-15, epsrel=1e-15, limit=50)
except Exception as e:
    print(f"Ошибка: {e}")

Ошибка: Достигнут предел количества подынтервалов (50).

4. Использование целочисленных аргументов, которые могут привести к целочисленному делению внутри функции.

def func_with_division(x, a, b):
    # Если a и b целые, деление в Python 2.x давало целочисленный результат.
    # В Python 3.x это не проблема.
    return a / b * x

# Правильно передавать числа с плавающей точкой.
result, error = spi.quad(func_with_division, 0, 1, args=(2.0, 3.0))
print(result)

0.33333333333333337

5. Неправильная работа с бесконечными пределами, если функция убывает недостаточно быстро.

def slow_decay(x):
    return 1 / (x + 1) # Убывает как 1/x

try:
    result, error = spi.quad(slow_decay, 0, np.inf)
except Exception as e:
    print(f"Интеграл расходится: {e}")

Интеграл расходится: Результат, вероятно, расходится или медленно сходится.

Изменения в функции quad

Функция scipy.integrate.quad является зрелой и стабильной, поэтому значительные изменения в ее интерфейсе происходят редко. Основные изменения связаны с улучшением алгоритма и исправлением ошибок.

• SciPy 1.8.0 (2022): Улучшена обработка некоторых крайних случаев для весовых функций и интегралов с параметрами.
• SciPy 1.6.0 (2020): Улучшена производительность и точность для некоторых типов интегралов, особенно с бесконечными пределами. В документации были уточнены описания некоторых параметров.
• Более ранние версии: Основные изменения касались внутренней реализации алгоритма QUADPACK (перевод с Fortran на C), что повысило надежность и переносимость. Добавлена поддержка интегралов с весовыми функциями типа Коши ('cauchy').

Пользователям рекомендуется обращаться к официальным release notes SciPy для получения точной информации об изменениях в конкретных версиях. Для большинства пользователей, использующих функцию в стандартных сценариях, изменения остаются прозрачными.

Расширенные примеры использования quad

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

import scipy.integrate as spi
import numpy as np

# Функция с логарифмической особенностью в 0
result, error = spi.quad(lambda x: np.log(x), 0, 1, points=[0], weight='alg', wvar=(0,0))
print(f"Интеграл от ln(x) от 0 до 1: {result:.6f}")
print(f"Точное значение: -1")

Интеграл от ln(x) от 0 до 1: -1.000000
Точное значение: -1

Интегрирование комплекснозначной функции (возвращающей действительную и мнимую части по отдельности).

def complex_func(x):
    return np.exp(1j * x) # e^(i*x)

# Интегрируем действительную и мнимую части по отдельности
real_part, real_err = spi.quad(lambda x: np.real(complex_func(x)), 0, np.pi)
imag_part, imag_err = spi.quad(lambda x: np.imag(complex_func(x)), 0, np.pi)
print(f"Результат: {real_part:.6f} + {imag_part:.6f}j")
print(f"Точное значение: 2.0 + 0.0j")

Результат: 2.000000 + 0.000000j
Точное значение: 2.0 + 0.0j

Использование для вычисления функции ошибок (erf) путем интегрирования.

def erf_integral(z):
    """Вычисление erf(z) через интеграл от exp(-t^2)."""
    prefactor = 2 / np.sqrt(np.pi)
    result, error = spi.quad(lambda t: np.exp(-t**2), 0, z)
    return prefactor * result

print(f"erf(1) приближенно: {erf_integral(1):.10f}")
print(f"erf(1) из scipy.special: {spi.special.erf(1):.10f}")

erf(1) приближенно: 0.8427007929
erf(1) из scipy.special: 0.8427007929

Интеграл с параметром, зависящим от переменной интегрирования в верхнем пределе (как пример составной задачи).

def integral_with_variable_limit(a, b):
    """Вычисляет ∫_a^b ( ∫_0^x cos(t) dt ) dx."""
    # Внутренний интеграл как функция от x
    def inner_integral(x):
        res, _ = spi.quad(np.cos, 0, x)
        return res
    # Внешний интеграл от результата внутреннего
    result, error = spi.quad(inner_integral, a, b)
    return result

print(f"Результат для a=0, b=pi: {integral_with_variable_limit(0, np.pi):.6f}")
# Аналитически: ∫ sin(x) dx от 0 до pi = 2.0
print(f"Аналитический результат: 2.0")

Результат для a=0, b=pi: 2.000000
Аналитический результат: 2.0

Обработка интеграла с разрывом первого рода (скачком) с помощью параметра points.

def piecewise_func(x):
    if x < 0.3:
        return x
    elif x < 0.7:
        return x**2
    else:
        return np.sqrt(x)

# Указываем точки разрыва производной (скачки функции нет, но есть изломы)
result, error = spi.quad(piecewise_func, 0, 1, points=[0.3, 0.7])
print(f"Интеграл кусочно-заданной функции: {result:.8f}")

Интеграл кусочно-заданной функции: 0.45817761