SciPy - 积分

SciPy - 数值积分 (`scipy.integrate`)

当函数无法进行解析积分（analytical integration），或者解析解过于复杂时，数值积分（Numerical integration，也称为求积，quadrature）变得至关重要。SciPy 的 scipy.integrate 模块提供了一套强大且通用的数值积分函数。

scipy.integrate 中的关键函数包括：

函数	描述
`quad`	通用单变量积分（自适应求积）。
`dblquad`	通用二重积分。
`tplquad`	通用三重积分。
`nquad`	通用 N 维多重积分。
`fixed_quad`	使用 n 阶高斯求积法积分 func(x)。
`quadrature`	使用高斯求积法积分到指定容差。
`romberg`	龙贝格积分。
`trapz` or `trapezoid`	对离散数据点使用梯形法则。
`simps` or `simpson`	对离散数据点使用辛普森法则。
`solve_ivp`	求解常微分方程组 (ODEs) 的初值问题 (Initial Value Problems)。
`solve_bvp`	求解常微分方程组 (ODEs) 的边值问题 (Boundary Value Problems)。

（注意：trapz 和 simps 现在也分别可用作 trapezoid 和 simpson，以提高清晰度。polyint 和 poly1d 是 NumPy 函数，常与积分结合使用。）

使用 `quad` 进行单变量积分

quad 函数是计算单变量函数定积分的主要工具： $\int_{a}^{b} f(x)dx$。

基本语法是 scipy.integrate.quad(func, a, b, args=(), ...)，其中：

func 是要积分的 Python 函数或可调用对象（callable）。
a 是积分下限。
b 是积分上限。
args 是一个可选元组，用于传递给 func 的额外参数。

quad 返回一个元组 (integral_value, absolute_error_estimate)，分别表示积分值和绝对误差估计。

让我们对高斯函数 $f(x) = e^{-x^2}$ 从 0 积分到 1：

import numpy as np
import scipy.integrate

# Define the function to integrate
def gaussian_function(x):
    return np.exp(-x**2)

# Alternatively, using a lambda function:
# gaussian_function_lambda = lambda x: np.exp(-x**2)

# Set integration limits
lower_limit = 0
upper_limit = 1

# Perform the integration
integral_value, error_estimate = scipy.integrate.quad(gaussian_function, lower_limit, upper_limit)

print(f"Integral of e^(-x^2) from {lower_limit} to {upper_limit}:")
print(f"Value: {integral_value}")
print(f"Estimated Absolute Error: {error_estimate}")

输出：

Integral of e^(-x^2) from 0 to 1:
Value: 0.7468241328124271
Estimated Absolute Error: 8.291413475940725e-15

quad 函数可以使用 np.inf 和 -np.inf 处理无穷积分限。它还可以直接积分许多标准的 NumPy 函数，只要这些函数接受单个数组参数并返回一个数组。

使用 args 参数传递参数的示例：

import numpy as np
import scipy.integrate

# Define a function f(x, a, b) = a*x + b
def linear_func(x, a, b):
    return a * x + b

# Integrate f(x, 2, 1) from 0 to 1
integral_value, _ = scipy.integrate.quad(linear_func, 0, 1, args=(2, 1)) # a=2, b=1
print(f"Integral of 2x + 1 from 0 to 1: {integral_value}") # Expected: [x^2 + x]_0^1 = 2

输出：

Integral of 2x + 1 from 0 to 1: 2.0

多重积分

对于多重积分，SciPy 提供了 dblquad、tplquad 和 nquad。

使用 `dblquad` 进行二重积分

dblquad 函数计算二重积分： $\int_{a}^{b} \int_{g(x)}^{h(x)} f(x,y) ,dy ,dx$。

语法： scipy.integrate.dblquad(func, a, b, gfun, hfun, args=())

func(y, x, ...): 要积分的函数（注意参数顺序：y 在前，然后是 x）。
a, b: x 的积分下限和上限（外层积分）。
gfun(x), hfun(x): 定义 y 积分下限和上限的函数（内层积分），它们可以依赖于 x。

让我们计算 $\int_{0}^{1/2} \int_{0}^{\sqrt{1-4y^2}} 16xy ,dx ,dy$。等等，原始示例是 $\int_{0}^{1/2} dy \int_{0}^{\sqrt{1-4y^2}} 16xy :dx$。这意味着 $x$ 是内层积分变量，$y$ 是外层积分变量。SciPy 的 dblquad 在 func(y,x) 时，积分 $dy dx$，其中 $y$ 是内层变量，$x$ 是外层变量。如果问题是 $\int_{y_low}^{y_high} dy \int_{x_low(y)}^{x_high(y)} f(x,y) ,dx$，我们需要匹配 dblquad 的预期或交换变量。原始示例 f = lambda x, y: 16 * x * y 使用 dblquad(f, 0, 0.5, g, h)，这意味着 f 是 f(内层变量, 外层变量)。考虑到限制 g 和 h 是针对内层变量的，而 0, 0.5 是针对外层变量的。如果 g 和 h 依赖于 y（如 h = lambda y: sqrt(1-4*y**2)），那么 y 必须是 dblquad 的外层变量。让我们重新验证 dblquad 的 func 和积分限参数顺序。

dblquad 的文档说明是：func(y, x) 的积分范围是 $y$ 从 gfun(x) 到 hfun(x)，而 $x$ 从 a 到 b。所以 $x$ 是外层变量，$y$ 是内层变量。我们的积分 $\int_{0}^{1/2} dy \int_{0}^{\sqrt{1-4y^2}} 16xy :dx$ 需要重新表述，如果 $y$ 是外层变量的话。让我们坚持原始的积分形式作为示例，假设用户期望 $y$ 是外层变量，$x$ 是内层变量。这意味着 $x$ 的积分限可以依赖于 $y$。dblquad 要求 func(内层变量, 外层变量)。因此，如果 $x$ 是内层，$y$ 是外层，则函数是 func(x,y)。$x$ 的积分限是 $g(y)$ 和 $h(y)$，$y$ 的积分限是 $a$ 和 $b$。这意味着：dblquad(lambda x, y: 16*x*y, y_lower, y_upper, lambda y_val: x_lower_limit(y_val), lambda y_val: x_upper_limit(y_val))。

让我们计算 $\int_{y=0}^{y=0.5} \int_{x=0}^{x=\sqrt{1-4y^2}} 16xy ,dx ,dy$。这里，$y$ 是外层变量，$x$ 是内层变量。

import numpy as np
import scipy.integrate

# Function to integrate: f(x, y) = 16xy
# dblquad expects func(inner_variable, outer_variable)
# Here, x is inner, y is outer.
integrand = lambda x, y: 16 * x * y

# Limits for the outer variable y
y_lower = 0
y_upper = 0.5

# Limits for the inner variable x (can be functions of y)
x_lower_limit_func = lambda y_val: 0
x_upper_limit_func = lambda y_val: np.sqrt(1 - 4 * y_val**2)

# Perform double integration
integral_value, error_estimate = scipy.integrate.dblquad(
    integrand,
    y_lower,          # Outer variable (y) lower limit
    y_upper,          # Outer variable (y) upper limit
    x_lower_limit_func, # Inner variable (x) lower limit function (of y)
    x_upper_limit_func  # Inner variable (x) upper limit function (of y)
)

print(f"Double integral value: {integral_value}")
print(f"Estimated error: {error_estimate}")

输出：

Double integral value: 0.5000000000000001
Estimated error: 1.7092350012594845e-14

`tplquad` 和 `nquad`

tplquad 以类似方式用于三重积分，而 nquad 则将其推广到 N 维积分。对于 nquad，您需要为每个维度提供一个积分限函数列表或范围列表。这些函数功能强大，但在高维情况下计算量可能很大。

积分采样数据

如果您拥有数据点 $(x_i, y_i)$ 而不是函数，则可以使用梯形法则（trapz 或 trapezoid）或辛普森法则（simps 或 simpson）等方法。

import numpy as np
from scipy.integrate import trapezoid, simpson

# Sample data points
x_samples = np.linspace(0, np.pi, 10)
y_samples = np.sin(x_samples)

# Integrate using trapezoidal rule
integral_trapz = trapezoid(y_samples, x_samples)
print(f"Integral of sin(x) from 0 to pi (Trapezoidal with 10 points): {integral_trapz}")

# Integrate using Simpson's rule (requires an odd number of samples, or an even number of intervals)
x_samples_simps = np.linspace(0, np.pi, 11) # 11 points = 10 intervals
y_samples_simps = np.sin(x_samples_simps)
integral_simps = simpson(y_samples_simps, x_samples_simps)
print(f"Integral of sin(x) from 0 to pi (Simpson with 11 points): {integral_simps}")

# Analytical result: integral of sin(x) from 0 to pi is [-cos(x)]_0^pi = -(-1) - (-1) = 2
print(f"Analytical result: {2.0}")

输出：

Integral of sin(x) from 0 to pi (Trapezoidal with 10 points): 1.9796508989590044
Integral of sin(x) from 0 to pi (Simpson with 11 points): 2.0001095173150043
Analytical result: 2.0

scipy.integrate 模块是 Python 中科学计算的基石。对于更复杂的问题，例如求解微分方程，请参考 solve_ivp（初值问题）和 solve_bvp（边值问题）等函数。有关详细用法和高级选项，请务必查阅 SciPy 官方文档。

SciPy - 积分

SciPy - 数值积分 (scipy.integrate)

使用 quad 进行单变量积分