数值积分

Gauss积分

张瑞
中国科学技术大学数学科学学院

Gauss积分

在Newton-Cotes积分中，可以看到$n$为偶数时，代数精度为$n+1$。


问题. $n$个点的数值积分公式，最多可以具有几阶代数精度？

例 1. 求$a_0$, $a_1$, $x_0$, $x_1$，使得数值积分公式

\[\int_{-1}^1f(x)dx\approx a_0f(x_0)+a_1f(x_1) \]

具有尽可能高的代数精度

解. 依定义，有4个未知量，列出4个方程

\[\begin{cases} a_0\cdot 1+a_1\cdot 1=\int_{-1}^1 1dx=2 \\ a_0\cdot x_0+a_1\cdot x_1=\int_{-1}^1 xdx=0 \\ a_0\cdot x_0^2+a_1\cdot x_1^2=\int_{-1}^1 x^2dx=\frac23 \\ a_0\cdot x_0^3+a_1\cdot x_1^3=\int_{-1}^1 x^3dx=0 \\ \end{cases} \]

由(2), $a_0x_0=-a_1x_1$，代入(4)

\[a_0x_0^3=-a_1x_1^3=a_0x_0x_1^2 \]

即有$x_0^2=x_1^2$，则有$x_0=-x_1$，代入(2)后，有

\[a_0x_0+a_1(-x_0)=0 \]

即有$a_0=a_1$，代入（1）后，可得

\[a_0=a_1=1 \]

联合$x_0=-x_1$，代入(3)，可得

\[x_0=-x_1=\frac{1}{\sqrt3} \]

可得到

\[\begin{cases} a_0=1, a_1=1 \\ x_0=\frac1{\sqrt3}, x_1=-\frac1{\sqrt3} \end{cases} \]

可以验证

\[a_0x_0^4+a_1x_1^4=\frac29\neq\int_{-1}^1x^4dx=\frac25 \]

即，格式的代数精度为$3$

问题. 有可能更高吗？

定理 1.
$n$个积分点的数值积分公式，至多具有$2n-1$阶代数精度

证明. 记

\[\begin{aligned} I(f)=&\int_a^bf(x)dx \\ I_n(f)=&\sum_{i=0}^na_if(x_i) , \quad a_i=\int_a^bl_i(x)dx \end{aligned} \]

取

\[p(x)=(x-x_0)^2\cdots(x-x_n)^2 \]

则有

\[I(p(x))>0, \quad I_n(p(x))=0 \]

即$2n+2$次多项式$p(x)$的数值积分有误差，因此数值积分公式$I_n(f)$的代数精度不超过$2n+1$阶

问题. 如何构造最高阶精度的公式？

更一般地，考虑如下的带权积分

\[I(f)=\int_a^bW(x)f(x)dx, \quad W(x)\geq 0 \]

$W(x)$称为权函数。定义两个可积函数的内积为

\[(f,g)=\int_a^bW(x)f(x)g(x)dx \]

在线性代数中学过，可以用Schmidt正交化过程，将一组基函数变为正交基函数。

利用Schmidt正交化过程

\[\begin{cases} g_0(x)=f_0(x) \\ \displaystyle g_i(x)=f_i(x)-\sum_{j=0}^{i-1}\frac{(f_i(x),g_j(x))}{(g_j(x),g_j(x))}g_j(x), j=1,2,\cdots,n \end{cases} \]

可以将$n$次多项式函数空间的一组基函数$\{1,x,x^2,\cdots,x^n\}$，变为正交基函数$\{g_0(x),\cdots,g_n(x)\}$，且有

\[\begin{aligned} &g_n(x)\in P^n(x) \\ &(g_n(x),g_j(x))=0 , j=0,1,\cdots,n-1 \end{aligned} \]

可以得到

\[\int_a^b g_n(x)p(x)dx=0, \forall p(x)\in P^{n-1}(x) \]

定理 2.
以$n$阶正交多项式$g_n(x)$的$n$个零点为积分点的数值积分公式具有$2n-1$阶的代数精度

证明. 只需证明，格式至少达到$2n-1$阶即可。由数值积分的误差公式，

\[E(f)=\int_a^b f[x_1,x_2,\cdots,x_n,x]\omega_n(x)W(x)dx \]

其中

\[\omega_n(x)=(x-x_1)\cdots(x-x_n) \]

注意到$x_i$是$g_n(x)$的零点，则有

\[g_n(x)=a(x-x_1)\cdots(x-x_n)=a\omega_n(x) \]

其中$a$为某非0实数。

另外，当$f(x)\in P^{2n-1}(x)$为$2n-1$次多项式时，$f[x_1,x_2,\cdots,x_n,x]$为$n-1$次多项式。

利用$g_n(x)$的正交性，可得$E(f)=0$。

即，格式具有$2n-1$阶代数精度。

定义 1.
称具有最高阶代数精度的数值积分格式为Gauss积分，相应的积分点称为Gauss点


Gauss积分的构造方法

  1. 求出区间$[a,b]$上权函数为$W(x)$的正交多项式$g_n(x)$
  2. 求出$g_n(x)$的所有零点
  3. 以这些零点作为积分点，构造的数值积分公式即为Gauss公式

例 2. 求积分$\displaystyle \int_{-1}^1x^2f(x)dx$的2点Gauss公式

解. 按Schmidt正交化过程，有

\[\begin{aligned} g_0(x)&=1 \\ g_1(x)&=x-\frac{(x,g_0(x))}{(g_0(x),g_0(x))}g_0(x) \\ g_2(x)&=x^2--\frac{(x^2,g_0(x))}{(g_0(x),g_0(x))}g_0(x)-\frac{(x^2,g_1(x))}{(g_1(x),g_1(x))}g_1(x) \\ &=x^2-\frac35 \end{aligned} \]

可以得到积分点为$\pm\sqrt{\frac35}$。相应的积分系数为

\[\begin{aligned} a_1=\int_{-1}^{1} x^2l_1(x)dx=\int_{-1}^1 x^2\frac{x-x_2}{x_1-x_2}dx=\frac13 \\ a_2=\int_{-1}^{1} x^2l_1(x)dx=\int_{-1}^1 x^2\frac{x-x_1}{x_2-x_1}dx=\frac13 \\ \end{aligned} \]

有2点Gauss公式

\[\int_{-1}^1 x^2f(x)dx\approx \frac13(f(-\sqrt{\frac35})+f(\sqrt{\frac35})) \]

Gauss-Legendre公式

区间$[-1,1]$上，权函数为$W(x)=1$的Gauss型公式称为Gauss-Legendre公式。

\[\int_{-1}^1f(x)dx\approx\sum_{i=1}^n\omega_if(x_i) \]

由

\[\int_a^bf(x)dx=\frac{b-a}2\int_{-1}^1f(\frac{a+b}2+\frac{b-a}2t)dt \]

可以得到区间$[a,b]$上的Gauss求积公式

\[\int_a^bf(x)dx\approx\frac{b-a}2\sum_{i=1}^n\omega_if(\frac{a+b}2+\frac{b-a}2x_i) \]

Gauss-Laguerre公式

区间$[0,\infty)$上，权函数$W(x)=e^{-x}$的积分公式，称为Gauss-Laguerre公式。

\[\int _{0}^{+\infty }e^{-x}f(x)\,dx\approx \sum _{i=1}^{n}w_{i}f(x_{i}) \]

其中$x_i$是Laguerre多项式

\[L_{n}(x)={\frac {e^{x}}{n!}}{\frac {d^{n}}{dx^{n}}}\left(e^{-x}x^{n}\right) \]

的根，积分系数为

\[w_{i}={\frac {x_{i}}{(n+1)^{2}[L_{n+1}(x_{i})]^{2}}} \]

Laguerre多项式满足公式

\[\begin{aligned} L_0(x)&=1 \\ L_1(x)&=1-x \\ L_{k+1}(x)&={\frac {(2k+1-x)L_{k}(x)-kL_{k-1}(x)}{k+1}} \end{aligned} \]

及边界值

\[L_k(0)=1, L'_k(0)=-k \]

可以表达为

\[L_n(x)=\sum _{k=0}^{n}{\binom {n}{k}}{\frac {(-1)^{k}}{k!}}x^{k} \]
\[\begin{aligned} L_0(x)=&1 \\ L_1(x)=&-x+1 \\ L_2(x)=&\frac {1}{2}(x^{2}-4x+2) \\ L_3(x)=&\frac {1}{6}(-x^{3}+9x^{2}-18x+6) \\ L_4(x)=&\frac {1}{24}(x^{4}-16x^{3}+72x^{2}-96x+24) \\ L_5(x)=&\frac {1}{120}(-x^{5}+25x^{4}-200x^{3}+600x^{2}-600x+120) \\ L_6(x)=&\frac {1}{720}(x^{6}-36x^{5}+450x^{4}-2400x^{3}+5400x^{2}-4320x+720) \end{aligned} \]

n=2时2点公式

\[\int_0^{+\infty}e^{-x}f(x)dx\approx \frac{2+\sqrt2}4f(2-\sqrt2)+\frac{2-\sqrt2}4f(2+\sqrt2) \]

Gauss-Hermite公式

区间$(-\infty,\infty)$上，权函数$W(x)=e^{-x^2}$的积分公式，称为Gauss-Hermite公式。

\[ \int_{-\infty }^{+\infty }e^{-x^{2}}f(x)\,dx\approx \sum_{i=1}^{n}w_{i}f(x_{i}) \]

其中$n$是积分点的个数。$x_i$是Hermite多项式

\[H_n(x)=(-1)^{n}e^{x^{2}}{\frac {\mathrm {d} ^{n}}{\mathrm {d} x^{n}}}e^{-x^{2}} \]

的根，积分系数为

\[ w_{i}={\frac {2^{n-1}n!{\sqrt {\pi }}}{n^{2}[H_{n-1}(x_{i})]^{2}}}. \]

Hermite多项式为

\[H_{n}(x)=\sum_{k=0}^{n}a_{n,k}x^{k} \]

它满足

\[H_{n+1}(x)=2xH_{n}(x)-H_{n}'(x)=2xH_{n}(x)-2nH_{n-1}(x) \]

则有系数关系式为

\[a_{n+1,k}={\begin{cases} -a_{n,k+1}&k=0\\ 2a_{n,k-1}-(k+1)a_{n,k+1}&k>0 \end{cases}} \]

和$a_{0,0}=1$, $a_{1,0}=0$, $a_{1,1}=2$。

Hermite多项式有如下的显式表达

\[H_{n}(x)=n!\sum _{m=0}^{\left\lfloor {\tfrac {n}{2}}\right\rfloor }{\frac {(-1)^{m}}{m!(n-2m)!}}(2x)^{n-2m} \]

Hermite多项式的前11个的值为

\[{\displaystyle { \begin{aligned} H_{0}(x)&=1\\ H_{1}(x)&=2x\\ H_{2}(x)&=4x^{2}-2\\ H_{3}(x)&=8x^{3}-12x\\ H_{4}(x)&=16x^{4}-48x^{2}+12\\ H_{5}(x)&=32x^{5}-160x^{3}+120x\\ H_{6}(x)&=64x^{6}-480x^{4}+720x^{2}-120\\ H_{7}(x)&=128x^{7}-1344x^{5}+3360x^{3}-1680x\\ H_{8}(x)&=256x^{8}-3584x^{6}+13440x^{4}-13440x^{2}+1680\\ H_{9}(x)&=512x^{9}-9216x^{7}+48384x^{5}-80640x^{3}+30240x\\ H_{10}(x)&=1024x^{10}-23040x^{8}+161280x^{6}-403200x^{4}+302400x^{2}-30240 \end{aligned} }} \]

n=2时，2点公式

\[\int_{-\infty}^{+\infty}e^{-x^2}f(x)dx\approx \frac{\sqrt{\pi}}2f(-\frac{\sqrt2}2) +\frac{\sqrt{\pi}}2f(\frac{\sqrt2}2) \]

n=3时，3点公式

\[\int_{-\infty}^{+\infty}e^{-x^2}f(x)dx\approx \frac{\sqrt{\pi}}6f(-\frac{\sqrt6}2) +\frac{2\sqrt\pi}3f(0) +\frac{\sqrt{\pi}}6f(\frac{\sqrt6}2) \]

谢谢

ref:

https://en.wikipedia.org/wiki/Laguerre_polynomials

https://en.wikipedia.org/wiki/Laguerre_polynomials

vertical slide 2

目录

本节读完

例 3.

[#ex9-1-0].