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$阶代数精度

思路： 找一个$2n$次的多项式，让数值积分有误差，即可证明数值公式至多只有$2n-1$阶代数精度。

证明. 对$n+1$个节点$\{x_i\}_{i=0}^n$，记数值积分公式

\[\begin{aligned} 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 \]

则有

\[\int_a^bp(x)dx>0 \]

但

\[I_n(p(x))=\sum_{i=0}^n a_ip(x_i)=\sum_{i=0}^n 0 =0 \]

即$2n+2$次多项式$p(x)$的数值积分有误差。

因此$n+1$个点的数值积分公式$I_n(f)$的代数精度不超过$2n+1$阶。

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

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

进而有

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

即

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

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

思路： 只需证明，格式至少达到$2n-1$阶即可。

证明. 数值积分的误差就是插值多项式误差的积分。 用Newton型的误差表达得到数值积分的误差为

\[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)=c\omega_n(x) \]

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

已经证明，当$f(x)\in P^{2n-1}(x)$为$2n-1$次多项式时， $f[x_1,x_2,\cdots,x_n,x]$为$n-1$次多项式。 由$g_n(x)$的特性，有

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

因而

\[\begin{aligned} E(f)=&\int_a^b f[x_1,x_2,\cdots,x_n,x]\omega_n(x)W(x)dx \\ =&\int_a^b \frac1c g_n(x) f[x_1, x_2, \cdots, x_n, x]W(x)dx=0 \end{aligned} \]

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

定义 1.
称具有最高阶代数精度的数值积分格式为Gauss积分，相应的积分点称为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公式

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

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

其中$x_i$是Legendre多项式

\[P_n(x)=\frac1{2^n n!}\frac{d^n}{dx^n}(x^2-1)^n \]

的根。

\[\omega_i=\frac{2}{(1-x_i^2)P'_n(x_i)} \]

Legendre多项式有递推关系，

\[\begin{aligned} &P_0(x)=1, P_1(x)=x, \\ &P_n(x)=\frac{2n-1}nxP_{n-1}(x)-\frac{n-1}nP_{n-2}(x) \end{aligned} \]

n=2时，$L_2(x) = \frac12(3x^2-1)$

n=3时，$L_3(x) = \frac12(5x^3-3x)$

n=2时，公式

\[\int_{-1}^1 f(x) dx \approx f(-\frac1{\sqrt 3})+f(\frac1{\sqrt 3}) \]

n=3时，有

\[\int_{-1}^1 f(x) dx \approx \frac59f(-\sqrt{\frac35})+\frac89f(0)+\frac59f(\sqrt{\frac35}) \]

例 3. 计算$\displaystyle\int_{-1}^1x^2\cos(x)dx$

解. 利用Gauss公式 

Calculute integration of x^2 cos(x) from -1 to 1
========= 2 points Gauss ==========
0.5586078851299956 Error:  0.08034063127322955
====== 3 points Gauss ========
0.47646879530281677 Error:  -0.0017984585539492781
======== 5 points Gauss ==========
0.47826718331725243 Error:  -7.05395136191278e-08

利用带权$x^2$的Gauss公式 

========= 2 points Gauss with weight x^2==========
0.47646879530281677 Error:  -0.0017984585539492781

对于区间$[a,b]$上的可积函数，做变量代换$x=\frac{a+b}2+\frac{b-a}2t$，得到

\[\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) \]

例 4. 用Gauss-Legendre公式计算$\displaystyle\int_0^1\frac{\sin(x)}xdx$

解. 计算结果 

========= 2 points Gauss ==========
0.9460411368978208, Error: -4.193346936220976e-05
====== 3 points Gauss ========
0.9460831340784723, Error:  6.371128935533932e-08
composite 2 times Gauss:
0.9460831340784723, Error:  9.758444052820892e-10
======== 5 points Gauss ==========
0.9460830703672151, Error:  3.2085445411667024e-14

复化积分计算结果 

======== Composite Trapaezoid ========
   2, Error: -0.0062897855610057896
   4, Error: -0.0015695487017934884
   8, Error: -0.0003922067844818189
  16, Error: -9.804043279693087e-05
  32, Error: -2.4509404415007374e-05
  64, Error: -6.127307119907499e-06
 128, Error: -1.5318240310646658e-06
 256, Error: -3.829558359313978e-07
 512, Error: -9.573894821368611e-08
1024, Error: -2.393473641504329e-08
======= Composite Simpson ========
   2 6.281190640378131e-05
   4 3.863584610797055e-06
   8 2.405212889966535e-07
  16 1.501776447643266e-08
  32 9.383790411376935e-10
  64 5.864508878516972e-11
 128 3.665401315799954e-12
 256 2.291500322826323e-13
 512 1.4210854715202004e-14
1024 9.992007221626409e-16

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) \]

例 5. 计算积分$\displaystyle\int_0^{\infty} e^{-x}\sin(x)dx$

解. 问题有真解$0.5$，数值结果 

calculute \int_{0}^{\infty}e^{-x}\sin(x)dx=0.5
2 points Gauss-Laguerre
0.4324594546798442 "Error: ",-0.0675405453201558
3 points Gauss-Laguerre
0.4960298274805634 "Error: ",-0.003970172519436599
5 points Gauss-Laguerre
0.49890332095606377 "Error: ",-0.0010966790439362328

问题. 对于积分$\displaystyle\int_0^{\infty} f(x)dx$可以写成如下表达式吗？

\[\int_0^{\infty}f(x)dx=\int_0^{\infty} e^{-x} {\color{red}f(x)e^{x}} dx \]

例 6. 求$\displaystyle\int_{0}^{\infty}\frac{\sin(x)}{x}dx=\frac{\pi}{2}$

解. 解得 

2 points Gauss-Laguerre
1.0961083699935568 Error: -0.4746879568013398
3 points Gauss-Laguerre
1.9552766139831168 Error: 0.38448028718822025
5 points Gauss-Laguerre
1.9583762914822787 Error: 0.3875799646873821

例 7. 求$\displaystyle \int_{0}^{\infty}\frac{\cos(x)}{1+x^2}dx=\frac{\pi}{2}e^{-1}$ 

2 points Gauss-Laguerre
0.61258765815713 Error: 0.034723983261669145
3 points Gauss-Laguerre
0.6867235574022316 Error: 0.1088598825067707
5 points Gauss-Laguerre
0.6235180541294388 Error: 0.04565437923397797

例 8. 求$\displaystyle \int_{0}^{\infty}\sin(x^2)dx=\sqrt{\frac{\pi}{8}}$ 

2 points Gauss-Laguerre
-2.996836634033218 Error: -3.623493702690968
3 points Gauss-Laguerre
3.194079602511634 Error: 2.567422533853884
5 points Gauss-Laguerre
5.3082444225535115 Error: 4.681587353895761

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}) \]

其中$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} \]

其中

\[\left\lfloor {\tfrac {n}{2}}\right\rfloor =\begin{cases} \frac{n}2 , & n=2k \\ \frac{n-1}2 , & n=2k+1 \end{cases} \]

有

\[{\displaystyle { \begin{aligned} H_{0}(x)&=1 , H_{1}(x)=2x\\ H_{2}(x)&=4x^{2}-2\\ H_{3}(x)&=8x^{3}-12x\\ \end{aligned} }} \]

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\\ \end{aligned} }} \]
\[{\displaystyle { \begin{aligned} 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) \]

例 9. 求$\displaystyle\int_{-\infty}^{\infty}e^{-x^2}\cos(x)dx$

解. 利用含参变量反常积分的特性可以得到

\[\int_{0}^{\infty}e^{-x^2}\cos(2\beta x)dx=\frac{\sqrt{\pi}}2e^{-\beta^2} \]

因而，问题的真解为$\sqrt{\pi}e^{-\frac14}$。 

calculute \int_{-\infty}^{\infty}e^{-x^2}\cos(x)dx
      True solution: \sqrt{\pi}e^{-(\frac12)^2}
2 points Gauss-Hermite:
1.3474984637168128 Error:  0.032889983326330086
3 points Gauss-Hermite:
1.3820330713880473 Error: -0.0016446243449044218
5 points Gauss-Hermite:
1.3803900759356567 Error: -1.62889251376086e-06

Gauss-Chebyshev公式

区间$[-1,1]$上，权函数$W(x)=\frac1{\sqrt{1-x^2}}$的积分公式，称为Gauss-Chebyshev公式。

\[\int _{-1}^{1}\frac1{\sqrt{1-x^2}}f(x)\,dx\approx \sum _{i=1}^{n}w_{i}f(x_{i}) \]

其中$x_i=\displaystyle \cos\frac{(2k-1)\pi}{2n}$是Chebyshev多项式

\[T_{n}(x)=\cos[n\arccos(x)] \]

的根，积分系数为

\[w_{i}=\frac{\pi}n \]

n=2，有

\[\int_{-1}^1 \frac1{\sqrt{1-x^2}}f(x)dx\approx \frac{\pi}2\left(f(-\frac1{\sqrt 2})+f(\frac1{\sqrt{2}})\right) \]

n=3，有

\[\int_{-1}^1 \frac1{\sqrt{1-x^2}}f(x)dx\approx \frac{\pi}3\left(f(-\frac{\sqrt3}{2})+f(0)+f(\frac{\sqrt3}{2})\right) \]

定理 3.
Gauss积分中的积分系数$A_k>0$。

证明. 设Gauss积分公式

\[G_n(f) = \sum_{k=1}^n A_k f(x_k) \]

其中

\[A_k = \int_a^b W(x) l_k(x)dx \]

$l_k(x)$是n-1次Lagrange插值基函数，$W(x)$为权函数。

  • 注意到Gauss积分具有最高阶代数精度，$2n-1$阶代数精度。
  • 对$i=1,2\cdots,n$，函数$l_i^2(x)$是$2n-2$次多项式，
  • 因而，有
    \[\int_a^b W(x) l_i^2(x)dx = \sum_{k=1}^n A_k l_i^2(x_k)=A_i \]
  • 而积分$\displaystyle \int_a^b W(x) l_i^2(x)dx>0$，从而$A_i>0$

作业： 试分析带权$W(x)$的积分$\displaystyle\int_a^bW(x)f(x)dx$的Gauss积分公式的误差

谢谢

vertical slide 2

ref:

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

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

目录

本节读完

例 10.

[#ex9-1-0].