Gauss-Lobatto 积分

                     

贡献者: addis

  • 本文需要更多讲解,便于帮助理解。
预备知识 定积分,换元积分法

  1高斯积分(Gauss quadrature)可以用求和来近似表示积分

\begin{equation} \int_{-1}^1 f(x) \,\mathrm{d}{x} \approx \sum_i w_i f(x_i)~. \end{equation}
对于某区间的一组 $x_i, w_i$($i = 1,\dots,N$),那么当 $f(x)$ 是小于等于某个阶的多项式时(取决于具体的类型),上式取等号。

   我们可以用换元积分法把上述定积分的区间变为任意 $[a,b]$。令

\begin{equation} t = \frac{b-a}{2}x + \frac{b+a}{2}~, \end{equation}
那么换元得(过程略)
\begin{equation} \int_a^b g(t) \,\mathrm{d}{t} \approx \sum_i w'_i g(t_i)~. \end{equation}
其中
\begin{equation} g(t) = \frac{2}{b-a} f \left(\frac{2t - b - a}{b - a} \right) = \frac{2}{b-a} f(x)~, \end{equation}
\begin{equation} t_i = \frac{b-a}{2}x_i + \frac{b+a}{2}~, \end{equation}
\begin{equation} \omega'_i = \frac{b-a}{2}\omega~. \end{equation}

1. Gauss-Lobatto 积分

   Gauss-Lobatto 积分中令采样点(包括两个端点)的个数为 $N$,如果 $f(x)$ 是 $2N-3$(或更低)阶多项式($f(x) = x^{2N-3} + \dots$),则积分没有误差。

   注意 Gauss-Lobatto 积分是对称的

\begin{equation} x_i = -x_{N-i+1}, \qquad w_{i} = w_{N-i+1}~. \end{equation}

   其中 $x_0 = -1, x_N = 1$,剩下的 $x_i$($1 < i < N$)是 $P'_{N-1}(x)$ 的根2,$P_N(x)$ 是勒让德多项式。另见式 10

\begin{equation} w_i = \frac{2}{N(N-1)[P_{N-1}(x_i)]^2}~. \end{equation}

   低阶情况下 $x_i$ 可以表示为带根号的表达式,在 Mathematica 中求解解析式和任意精度数值解如

代码 1:gauss_lobatto.nb
NN = 6;(*Nodes*)
digits = 38;(*digits*)
sol = Solve[D[LegendreP[NN - 1, x], x] == 0, x];
sol = Append[sol, {x -> -1}];
sol = Append[sol, {x -> 1}];
xi = N[x /. sol, digits](*DVR points*)
w = N[2/(NN (NN - 1) (LegendreP[NN - 1, x])^2) /. sol, digits](*DVR weights*)
其中 sol 得到的是 list of rule,/. 用于把这些 rule 作用到前面的表达式上面。

2. 正交归一基底

   每个基底都是 $N-1$ 阶多项式,由于阶数和零点数一样,多项式可以唯一确定,即拉格朗日插值多项式

\begin{equation} \begin{aligned} p_n(x) &= \prod_{i=1}^{n-1} \frac{x-x_i}{x_n-x_i} \prod_{i=n+1}^{N} \frac{x-x_i}{x_n-x_i}\\ &= \frac{x-x_1}{x_n-x_1} \dots \frac{x-x_{n-1}}{x_n-x_{n-1}}\frac{x-x-{n+1}}{x_n-x_{n+1}} \dots \frac{x-x_N}{x_n-x_N}~, \end{aligned} \end{equation}
\begin{equation} p_n(x_{n'}) = \delta_{n, n'}~. \end{equation}

   由于任意两个基底乘积只是 $2N-2$ 阶的多项式,所以用求和代替积分存在误差。可以证明它们之间有近似的正交关系

\begin{equation} \int_{-1}^1 p_i(x) p_j(x) \,\mathrm{d}{x} \approx \sum_k w_k p_i(x_k) p_j(x_k) = w_i \delta_{ij}~. \end{equation}
所以可以定义 $N$ 个近似正交归一的多项式基底
\begin{equation} f_n(x) = \frac{1}{\sqrt{w_n}} p_n(x)~, \end{equation}
满足
\begin{equation} f_i(x_j) = \frac{1}{\sqrt{w_i}} \delta_{ij}~. \end{equation}

   基底的另一种等效的表示方法是利用 $N$ 阶 Gauss-Lobatto 数值积分对应的多项式 $P'_{N-1}(x)$ 来构建。根据定义,其 $N-2$ 个零点分别为 $x_2, x_3\dots x_{N-1}$。给它加入 $x_1=-1$ 与 $x_N=1$ 这两个零点,将其变为 $N$ 阶多项式得

\begin{equation} (1-x^2)P'_{N-1}(x)~. \end{equation}
然而,式 10 要求 $p_n(x_n)=1$,所以我们将式 14 除以它自己在 $x_N$ 处的切线,在 $x=x_N$ 处形成极限类型 $0/0=1$ 即可得到多项式 $p_n(x)$。
\begin{equation} p_n(x) = \frac{(1-x^2)P'_{N-1}(x)}{[(1-x^2)P'_{N-1}(x)]'_{x = x_n}(x-x_n)}~. \end{equation}
该式与式 9 事实上是完全相同的多项式,因为所有具有 $N-1$ 个零点的 $N-1$ 阶多项式都可以因式分解成式 9 的形式乘以一个待定常数。用式 15 便于快速地展开多项式(因为勒让德多项式的系数可以直接通过公式计算)。


1. ^ 参考 Wikipedia 相关页面
2. ^ $P_N(x)$ 在 $[-1,1]$ 有 $l$ 个根。$P'_N(x)$ 有 $N-1$ 个根,所以 $P'_{N-1}(x)$ 有 $N-2$ 个根,加上两个端点就是 $N$ 个。

                     

© 小时科技 保留一切权利