几种数值积分方法

发表于 | 分类于 | 评论数：

Gaussian Quadrature

适用于积分区间内平滑的积分.

Computes the sample points and weights for Gauss-Legendre quadrature. The sample points are the roots of the n-th degree Legendre polynomial . These sample points and weights correctly integrate polynomials of degree or less over the interval with weight function .

\begin{align} \int_{-1}^1 f(x)\mathrm{d}x = \sum_{i=1}^n w_i f(x_i) \end{align}

其中 \(x_i\) 是 \(n\) 阶 Legendre 多项式的根. \(x_i\) 和权重 \(w_i\) 可以由 scipy.special.roots_ legendre 生成:

Computes the sample points and weights for Gauss-Legendre quadrature. The sample points are the roots of the n-th degree Legendre polynomial . These sample points and weights correctly integrate polynomials of degree or less over the interval with weight function .

积分区间不在 \([-1, 1]\) , 可以通过换换元来转换

\begin{align} \int_a^b f(x) \mathrm{d}x = \frac{b-a}{2}\int_{-1}^1 f(\frac{b-a}{2}t + \frac{b+a}{2})\mathrm{d}t \end{align} 
from scipy.special import roots_legendre as leg
def gauquad(f, a, b, n=50):
    '''
    定义 Gaussian quadrature 积分
    函数 f 的积分区间为 [a,b]
    取 n 个 Legendre 的根
    def Gaussian quadrature integration
    integrate function f from a to b
    take n Legendre roots
    '''
    ft = lambda t: f( (b-a)*t/2 +(a+b)/2 ) * (b-a)/2
    x, w = leg(n)
    I = 0
    for i in range(n):
	I = I + w[i]*ft(x[i])
    err = 0                     # 为了与 scipy.integrate 中积分函数的输出一致, 多一个 err 参数.
    return I,err

Cauchy Principal Value (CPV) Gaussian Quadrature

CPV 用 Gaussian quadrature 时不能给出正确的结果. 这时可以将函数的积部 分和偶部分分开算.

计算积分

\begin{align} \mathcal{P}\int_{-a}^af(x) \mathrm{d}x = \lim_{r\to 0^+}\left[ \int_{-a}^{-r} f(x)\mathrm{d}x +\int_{r}^af(x)\mathrm{d}x \right] \end{align}

其中 \(0\) 为奇点.

关于奇点的 Odd 部分

\begin{align} g(x) = \frac{1}{2}\left[ f(x)-f(-x) \right] \end{align}

Even 部分

\begin{align} h(x) = \frac{1}{2}\left[ f(x) + f(-x)\right] \end{align}

总的积分

\begin{align} \mathcal{P}\int_{-a}^a f(x) \mathrm{d}x = \lim_{r\to 0^+} \int_r^a \left[ f(x) + f(-x)\right] \mathrm{d}x \end{align}

如果奇点不是 \(0\) , 也可以通过换元转化成上式的形式.

code

from scipy import integrate
from scipy.integrate import fixed_quad
def sgq(f, a, b, sp, n=10):
    """带有主值积分的积分.
    Singular Gaussian quadrature."""
    diffa = sp - a
    diffb = b - sp
    def ff(t):
	ff = f(t+sp) + f(-t+sp)
	return ff
    if diffa<diffb:
	sgq1, err = gauquad(ff, 0, sp-a, n=rootNum)
	sgq2, err = gauquad(f, 2*sp-a, b, n=rootNum)
    else:
	sgq1, err = gauquad(ff, 0, b-sp, n=rootNum)
	sgq2, err = gauquad(f, a, 2*sp-b, n=rootNum)
    sgq = sgq1 + sgq2
    return sgq, err

Tanh-Sinh Quadrature

对于端点发散的情况, 可以用 Tanh-Sinh quadrature (Double Exponential) 方法. 对于积分

\begin{align} \int_{-1}^1 f(x)\mathrm{d}x \end{align}

通过换元

\begin{align} x = \tanh \left( \frac{1}{2}\pi \sinh t \right) \end{align}

将积分区间 \(x\in [-1, 1]\) 换为 \(t\in (-\infty, +\infty)\) .

\begin{align} \int_{-1}^1 f(x)\mathrm{d}x = \int_{-\infty}^{+\infty}f[x(t)] w(t)\mathrm{d}t \end{align}

其中

\begin{align} w(t) = \frac{\mathrm{d}}{\mathrm{d}t}x(t) =\frac{ \frac{1}{2}\pi \cosh t}{\cosh^2\left( \frac{1}{2}\pi \sinh t \right)} \end{align}

权重 \(w(t)\) 是一个以原点为中心的钟形函数, 因此可以抵消原来在端点的发散.

然后将积分离散化, 将积分区间作适当截断, 做数值计算. 

def ts(f, a, b, n=51):
    """Tanh-sinh quadrature 方法.
    适用于端点发散的情况."""
    up = 4
    h = 2*up / (n-1)
    t = np.linspace(-up, up, n, endpoint=True)
    x = np.tanh(1/2*np.pi*np.sinh(t))
    w = 1/2*h*np.pi*np.cosh(t)
    w = w/(np.cosh(1/2*np.pi*np.sinh(t))**2)
    gc = 0
    for i in range(n):
	p = (x[i]*(b-a) + a + b)/2
	gc = gc + f(p)*w[i]
    err = 0
    gc = gc * (b-a)/2
    return gc, err