2.1 `scipy.integrate`：数值积分与微分方程求解

2.1 scipy.integrate：数值积分与微分方程求解

SciPy 核心模块详解：2.1 scipy.integrate - 数值积分与微分方程求解

2.1.1 数值积分

数值积分，也称为求积，是计算定积分近似值的过程。当无法找到积分的解析解，或者积分函数过于复杂时，数值积分就变得至关重要。scipy.integrate 提供了多种数值积分方法，主要分为以下几类：

• quad: 通用积分器，适用于单变量函数的定积分。

• dblquad: 用于二重积分。

• tplquad: 用于三重积分。

• nquad: 用于多重积分 (n 维)。

• fixed_quad: 使用固定阶数的高斯求积。

• quadrature: 使用自适应高斯求积。

• romberg: 使用 Romberg 积分法。

• trapz: 使用梯形法则进行积分。

• simps: 使用辛普森法则进行积分。

2.1.1.1 quad 函数

quad 函数是 scipy.integrate 中最常用的积分函数。它的基本语法如下：

result, error = quad(func, a, b, args=(), ...)
​

• func: 要积分的函数，必须是单变量函数。

• a: 积分下限。

• b: 积分上限。

• args: 传递给 func 的额外参数，以元组形式传递。

• result: 积分的近似值。

• error: 积分的绝对误差估计。

示例：计算 \int_{0}^{1} x^2 dx

import numpy as np
from scipy import integrate
def f(x):
  return x**2
result, error = integrate.quad(f, 0, 1)
print("积分结果:", result)
print("误差估计:", error)
​

示例：传递额外参数

def f(x, a):
  return a * x**2
a = 2
result, error = integrate.quad(f, 0, 1, args=(a,))
print("积分结果:", result)
print("误差估计:", error)
​

示例：处理无穷积分

def f(x):
  return np.exp(-x)
result, error = integrate.quad(f, 0, np.inf)
print("积分结果:", result)
print("误差估计:", error)
​

2.1.1.2 多重积分 (dblquad, tplquad, nquad)

对于多重积分，scipy.integrate 提供了 dblquad (二重积分), tplquad (三重积分) 和 nquad (n 重积分) 函数。

示例：计算 \int_{0}^{1} \int_{0}^{2} xy \, dx \, dy

from scipy import integrate
def f(x, y):
  return x * y
result, error = integrate.dblquad(f, 0, 1, lambda x: 0, lambda x: 2)
print("积分结果:", result)
print("误差估计:", error)
​

注意：dblquad 和 tplquad 的积分限可以是函数，这使得可以处理更复杂的积分区域。nquad 提供了更通用的多重积分方法，可以处理任意维度的积分。

2.1.1.3 基于样本的积分 (trapz, simps)

当函数的值在离散点上给出时，可以使用 trapz (梯形法则) 和 simps (辛普森法则) 进行数值积分。

示例：使用梯形法则积分

import numpy as np
from scipy import integrate
x = np.linspace(0, np.pi, 100)
y = np.sin(x)
result = integrate.trapz(y, x)
print("积分结果:", result)
​

示例：使用辛普森法则积分

import numpy as np
from scipy import integrate
x = np.linspace(0, np.pi, 100)
y = np.sin(x)
result = integrate.simps(y, x)
print("积分结果:", result)
​

2.1.2 常微分方程 (ODE) 求解

scipy.integrate 提供了多种求解常微分方程 (ODE) 的方法。主要函数是 solve_ivp。

2.1.2.1 solve_ivp 函数

solve_ivp 函数是一个通用的 ODE 求解器，可以处理各种类型的 ODE 系统。它的基本语法如下：

result = solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, ...)
​

• fun: 定义 ODE 系统的函数。它必须接受时间 t 和状态向量 y 作为输入，并返回 y 的导数。

• t_span: 积分的时间区间，例如 (t0, tf)。

• y0: 初始状态向量。

• method: 使用的求解器方法。常见的选项包括 'RK45' (默认), 'RK23', 'DOP853', 'BDF', 'LSODA' 等。

• t_eval: 可选参数，指定要计算解的时间点。

• result: 一个包含求解结果的对象。

示例：求解一阶 ODE: \frac{dy}{dt} = -y, y(0) = 1

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
def dydt(t, y):
  return -y
t_span = (0, 5)
y0 = [1]  # 初始条件
result = solve_ivp(dydt, t_span, y0, dense_output=True)
print(result)
t = np.linspace(0, 5, 100)
y = result.sol(t)[0]
plt.plot(t, y)
plt.xlabel('t')
plt.ylabel('y(t)')
plt.title('Solution of dy/dt = -y')
plt.grid(True)
plt.show()
​

示例：求解二阶 ODE: \frac{d^2y}{dt^2} + y = 0, y(0) = 1, \frac{dy}{dt}(0) = 0 (简谐振动)

为了使用 solve_ivp 求解二阶 ODE，需要将其转换为一阶 ODE 系统。令 y_1 = y 和 y_2 = \frac{dy}{dt}，则有：

\frac{dy_1}{dt} = y_2

\frac{dy_2}{dt} = -y_1

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
def d2ydt2(t, y):
  y1, y2 = y
  return [y2, -y1]
t_span = (0, 10)
y0 = [1, 0]  # 初始条件 [y(0), dy/dt(0)]
result = solve_ivp(d2ydt2, t_span, y0, dense_output=True)
t = np.linspace(0, 10, 100)
y1 = result.sol(t)[0]
y2 = result.sol(t)[1]
plt.plot(t, y1, label='y(t)')
plt.plot(t, y2, label='dy/dt(t)')
plt.xlabel('t')
plt.ylabel('Amplitude')
plt.title('Solution of d^2y/dt^2 + y = 0')
plt.grid(True)
plt.legend()
plt.show()
​

示例：使用 t_eval 指定输出时间点

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
def dydt(t, y):
  return -y
t_span = (0, 5)
y0 = [1]
t_eval = np.linspace(0, 5, 10) # 指定输出时间点
result = solve_ivp(dydt, t_span, y0, t_eval=t_eval)
print("时间点:", result.t)
print("解:", result.y)
plt.plot(result.t, result.y[0])
plt.xlabel('t')
plt.ylabel('y(t)')
plt.title('Solution of dy/dt = -y')
plt.grid(True)
plt.show()
​

2.1.2.2 选择合适的求解器方法

solve_ivp 提供了多种求解器方法，选择合适的求解器取决于 ODE 系统的性质。

• 'RK45': 默认方法，适用于大多数非刚性 ODE 系统。

• 'RK23': 低阶 Runge-Kutta 方法，适用于精度要求不高的场合。

• 'DOP853': 高阶 Runge-Kutta 方法，适用于高精度计算。

• 'BDF': 隐式方法，适用于刚性 ODE 系统。

• 'LSODA': 自动切换方法，可以自动选择刚性或非刚性求解器。

刚性 ODE 系统是指那些解的变化速度非常不同的系统。刚性系统通常需要隐式方法才能有效地求解。

2.1.3 scipy.integrate 的高级应用

scipy.integrate 还可以用于更高级的应用，例如：

• 参数估计: 通过求解 ODE 并将结果与实验数据进行比较，可以估计 ODE 模型中的参数。

• 灵敏度分析: 可以分析 ODE 模型的解对参数变化的敏感程度。

• 最优控制: 可以使用 scipy.integrate 求解最优控制问题，找到使系统性能最优的控制策略。

2.1.4 总结

scipy.integrate 模块是 SciPy 库中用于数值积分和求解常微分方程的强大工具。它提供了多种积分器和 ODE 求解器，可以处理各种类型的函数和 ODE 系统。通过学习和掌握 scipy.integrate 的使用方法，可以解决科学计算和工程应用中的许多实际问题。