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$ 个。


致读者: 小时百科一直以来坚持所有内容免费无广告,这导致我们处于严重的亏损状态。 长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。 因此,我们请求广大读者热心打赏 ,使网站得以健康发展。 如果看到这条信息的每位读者能慷慨打赏 20 元,我们一周就能脱离亏损, 并在接下来的一年里向所有读者继续免费提供优质内容。 但遗憾的是只有不到 1% 的读者愿意捐款, 他们的付出帮助了 99% 的读者免费获取知识, 我们在此表示感谢。

                     

友情链接: 超理论坛 | ©小时科技 保留一切权利