+++
date = ‘2024-12-21T13:49:00+08:00’
draft = true
title = ‘数值计算方法(2) 数值积分方法’
+++
初次发布于我的个人文档
上一期讲了插值方法，这一次自然是要运用一下插值方法了。所以这一期的主题是用插值方法计算定积分。
机械求积方法
下面我们来介绍一下怎么用插值法来得到数值方法计算函数积分。
对于大部分函数，我们其实都是很难求其积分的，甚至很多函数例如 e x 2 e^{x^2} ex2这样的函数压根就没有初等原函数。所以我们需要寻找数值方法来计算他们的积分。
对于大多数函数，我们不仅很难求积分，也许微分也很难求，所以我们希望的数值方法一定是不能包含这俩的。
为了简化问题，我们最终的求积公式就用函数值的线性组合试试看了。
这种最简单的求积公式就是机械求积。
也就是找到某些点的函数值，用这些点的函数值的某个线性组合来近似积分值。
即机械求积公式认为 ∫ a b f ( x ) d x ≈ ∑ k = 0 n A k y k \int_a^b f(x)dx \approx \sum_{k=0}^n A_k y_k ∫abf(x)dx≈∑k=0nAkyk
所以关键点就是点的选取和组合系数 A k A_k Ak的选择了。
一种可能的思路就是用上一期提的插值多项式来近似替代这个函数。而且我们是已知了差值余项的，余项的积分就是我们数值积分方法的误差。
而这些都是多项式啊，都很好计算的。
这种求积分的数值方法得到的积分公式我们都叫做差值型的求积公式。
而在插值方法里，点的选择是任意的，那我们就暂且也让大家任意选择。
这样只需要知道组合系数 A k A_k Ak就知道怎么求积了。
而对于插值型求积公式，我们很容易就能得到 A k A_k Ak。
插值型求积公式就是用f(x)的插值多项式近似f(x)，从而有
∫ a b f ( x ) d x ≈ ∫ a b p n ( x ) d x \int_a^bf(x)dx \approx \int_a^bp_n(x)dx ∫abf(x)dx≈∫abpn(x)dx
再把拉格朗日插值公式代入得
∫ a b f ( x ) d x ≈ ∫ a b ∑ k = 1 n y k l k ( x ) d x \int_a^bf(x)dx \approx \int_a^b \sum_{k=1}^n y_k l_k(x)dx ∫abf(x)dx≈∫ab∑k=1nyklk(x)dx
= ∑ k = 1 n y k ∫ a b l k ( x ) d x =\sum_{k=1}^n y_k\int_a^bl_k(x)dx =∑k=1nyk∫ablk(x)dx
而我们的机械求积公式是 ∫ a b f ( x ) d x = ∑ k = 0 n A k y k \int_a^b f(x)dx=\sum_{k=0}^n A_k y_k ∫abf(x)dx=∑k=0nAkyk，从而
A k = ∫ a b l k ( x ) d x A_k=\int_a^b l_k(x)dx Ak=∫ablk(x)dx
代数精度
你也许会觉得这个办法非常粗糙，确实！我们可以来分析一下他的精度。余项么，自然就是拉格朗日插值余项的积分了，这个没啥好说的。
我们在机械求积里其实最关注代数精度了。
我们拿多项式函数来表彰一个求积公式的精度。
取多项式 p ( x ) = a n x n + a n − 1 x n − 1 + . . . + a 0 p(x)=a_nx^n+a_{n-1}x^{n-1}+...+a_0 p(x)=anxn+an−1xn−1+...+a0
由于多项式函数总是很容易求出积分的准确值的。
我们就来看看你的求积公式是从几次多项式开始才变得不准确的。这个次数就是求积公式的代数精度。
事实上，我们前面那种求积方式的代数精度是很容易看出来的。
你用n次多项式插值f(x)而我们代入验证的f(x)是m次多项式。
如果n=m，那么你的插值多项式就是我们代入验证的m次多项式。
如果n比m还大呢？求积公式的插值多项式会退化为m次多项式，所以积分自然也还是准确的。
因此，n次插值多项式得到的求积公式至少具有n次代数精度。
我这里是用的任意n次多项式来进行验证，不过你在外面看到的定义可能是说用 x n x^n xn代入验证，事实上我们俩说的是等价的。
其他人的说法是你代入 1 , x , x 2 , . . . , x n 1,x,x^2,...,x^n 1,x,x2,...,xn，一直到 x n − 1 x^{n-1} xn−1次方，你都发现求积公式是准确的，但n次方开始不准确了，这时说求积公式有n-1次代数精度。
显然，他们的定义是我说的定义的特殊情况，所以满足我的定义的，自然满足他们的定义。
而满足他们定义的求积公式，对 x n − 1 x^{n-1} xn−1都准确。
则我取任意的n-1次多项式， p ( x ) = a n − 1 x n − 1 + . . . + a 0 p(x)=a_{n-1}x^{n-1}+...+a_0 p(x)=an−1xn−1+...+a0
求积分 ∫ a b p ( x ) d x = a n − 1 ∫ a b x n − 1 d x + a n − 2 ∫ a b x n − 2 d x + . . . + a 0 ∫ a b 1 d x \int_a^bp(x)dx=a_{n-1}\int_a^bx^{n-1}dx+a_{n-2}\int_a^bx^{n-2}dx+...+a_0\int_a^b1dx ∫abp(x)dx=an−1∫abxn−1dx+an−2∫abxn−2dx+...+a0∫ab1dx
1 , x , x 2 , . . . , x n − 1 1,x,x^2,...,x^{n-1} 1,x,x2,...,xn−1的积分准确，那么p(x)的积分也就是准确的了。
所以我们俩的代数精度的定义是等价的，只不过我这个看起来要强一些而已。
牛顿-科特斯公式
前面我们给的机械求积公式略显粗糙，但也能取得至少n次代数精度，接下来我们会对求积公式进行加强，对它有更高的要求。
组合系数 A k A_k Ak是死的了，所以下面的加强都是从插值点入手。
例如牛顿-科特斯公式就是进一步要求，插值点的间距相等（步长为h），然后得到一个特殊的结论。
我们竟然将区间长度提出来了，
∫ a b f ( x ) d x = ( b − a ) ∑ k = 0 n C k f ( x k ) \int_a^b f(x)dx=(b-a)\sum_{k=0}^n C_k f(x_k) ∫abf(x)dx=(b−a)∑k=0nCkf(xk)
并且， C k C_k Ck竟然还是常数。
我们称 C k C_k Ck为科特斯系数，可以直接查科特斯系数表得到，与具体的函数f(x)无关，只和你插值多项式的次数n有关!
这就是牛顿-科特斯公式。
其中n=1得到的求积公式被称为梯形公式，n=2的是辛普森公式。
一般我们也就只用到辛普森公式为止，这是因为前面我们说过插值多项式有龙格现象，也就是会过拟合。我们这里也是一样，你用的求积公式次数太高，插值多项式在区间端点附近疯狂震荡，而且这时候你稍微变一下计算用的数值最后的计算结果差距会非常大（这种情况我们称为数值稳定性差）。
所以记住常用的求积公式，三个矩形公式，也就是
左矩形公式：把f(x)直接当做常函数f(a)求矩形的面积 ∫ a b f ( x ) d x = ( b − a ) f ( a ) \int_a^bf(x)dx=(b-a)f(a) ∫abf(x)dx=(b−a)f(a)
还有右矩形公式： ∫ a b f ( x ) d x = ( b − a ) f ( b ) \int_a^bf(x)dx=(b-a)f(b) ∫abf(x)dx=(b−a)f(b)
以及中矩形公式： ∫ a b f ( x ) d x = ( b − a ) f ( a + b 2 ) \int_a^bf(x)dx=(b-a)f(\frac{a+b}{2}) ∫abf(x)dx=(b−a)f(2a+b)
梯形公式，把f(x)看作(a,f(a)),(b,f(b))的直线求梯形面积,当然也是什么的n=1时的牛顿科特斯公式 ∫ a b f ( x ) d x = ( b − a ) 1 ∗ f ( b ) + 1 ∗ f ( a ) 1 + 1 \int_a^bf(x)dx=(b-a)\frac{1*f(b)+1*f(a)}{1+1} ∫abf(x)dx=(b−a)1+11∗f(b)+1∗f(a)
既然梯形公式是把f(x)看作(a,f(a)),(b,f(b))的直线求梯形面积你其实大可不必用上面我写的复杂的式子，直接上底+下底再乘高除2就完事了。
我也这样是为了让你记住辛普森公式。
辛普森公式涉及三个点，然后是等距的，这里假定求a,c的积分，b就是a和c的中点，从而
∫ a c f ( x ) d x = ( c − a ) 1 ∗ f ( a ) + 4 ∗ f ( b ) + 1 ∗ f ( c ) 1 + 4 + 1 \int_a^cf(x)dx=(c-a)\frac{1*f(a)+4*f(b)+1*f(c)}{1+4+1} ∫acf(x)dx=(c−a)1+4+11∗f(a)+4∗f(b)+1∗f(c)
重点在系数1,4,1，记住这个就可以了。
这也是上面的牛顿-科特斯系数表不把 4 6 \frac{4}{6} 64化简的原因。
这里也点出了牛顿-科特斯系数表的性质即每一行的和为1。
好了，总之这样你应该就记住了这几个公式了。
但是问题来了，这玩意儿啊都是用的低阶插值，所以代数精度不也就低了。
左右矩形公式只用了1个点，所以是0阶插值的结果，有0阶代数精度。
梯形公式两个点，是线性插值的结果，只有1阶代数精度。
辛普森公式3个点，是抛物插值的结果，2阶代数精度，这不就完全垮了啊！
你应该能想到应该怎么做，分段呗，继续分段。
复化方法
这种分段的方法在前面分段低次插值和样条插值里你已经见过了，只不过在这里叫复化方法。
对应的求积公式叫复化求积公式。
我们将区间[a,b]n等分，然后用低阶求积公式求每一段的积分值，然后累加就是整个区间的积分值了，复化方法就是这么简单。
但是我一写公式你又头晕了。
记住上面那一句话的本质即可，公式你对照着本质看就能看懂了。
我们将[a,b]按步长h等分成n份，分点就是 x k = a + k h , k = 0 , 1 , 2 , . . . , n x_k=a+kh,k=0,1,2,...,n xk=a+kh,k=0,1,2,...,n，没毛病吧。区间有n个，所以要n+1个分点，这里下标从0开始，所以k最多取n。
然后求出这n段的积分 I k I_k Ik，注意，分点是0到n，有n+1个分点，区间有n个，如果下标还是从0开始那么k最多到n-1。
然后将这n个积分求和就是 ∫ a b f ( x ) d x = ∑ k = 0 n − 1 I k \int_a^b f(x)dx=\sum_{k=0}^{n-1}I_k ∫abf(x)dx=∑k=0n−1Ik。
到这里是最一般的复化求积公式，没错哈。
然后我们再用低阶求积公式求具体的 I k I_k Ik。
例如，用梯形公式求区间 [ x k , x k + 1 ] [x_k,x_{k+1}] [xk,xk+1]的积分。注意哈k从0开始，而 x 0 = a x_0=a x0=a，而这边是 [ a , a + h ] , [ a + h , a + 2 h ] [a,a+h],[a+h,a+2h] [a,a+h],[a+h,a+2h]这样分区间的，所以区间是k到k+1。
好，来用梯形公式得到 I k = h ∗ f ( x k ) + f ( x k + 1 ) 2 I_k=h*\frac{f(x_k)+f(x_{k+1})}{2} Ik=h∗2f(xk)+f(xk+1)
再试试用辛普森公式求解， I k = h ∗ f ( x k ) + 4 ∗ f ( x k + 1 2 ) + f ( x k + 1 ) 6 I_k=h*\frac{f(x_k)+4*f(x_{k+\frac{1}{2}})+f(x_{k+1})}{6} Ik=h∗6f(xk)+4∗f(xk+21)+f(xk+1)
这里 x k + 1 2 x_{k+\frac{1}{2}} xk+21指的是 x k x_k xk和 x k + 1 x_{k+1} xk+1的中点。
接下来要把 I k I_k Ik代入复化求积公式了哈。
这样就得到复化梯形公式 ∫ a b f ( x ) d x = ∑ k = 0 n − 1 h ∗ f ( x k ) + f ( x k + 1 ) 2 \int_a^b f(x)dx=\sum_{k=0}^{n-1}h*\frac{f(x_k)+f(x_{k+1})}{2} ∫abf(x)dx=∑k=0n−1h∗2f(xk)+f(xk+1)
复化辛普森公式 ∫ a b f ( x ) d x = ∑ k = 0 n − 1 h ∗ f ( x k ) + 4 ∗ f ( x k + 1 2 ) + f ( x k + 1 ) 6 \int_a^b f(x)dx=\sum_{k=0}^{n-1}h*\frac{f(x_k)+4*f(x_{k+\frac{1}{2}})+f(x_{k+1})}{6} ∫abf(x)dx=∑k=0n−1h∗6f(xk)+4∗f(xk+21)+f(xk+1)
不过你在外面看到的式子可能不是这样的，我们来化简一下。
首先把常数提出去。
复化梯形公式 ∫ a b f ( x ) d x = h 2 ∑ k = 0 n − 1 f ( x k ) + f ( x k + 1 ) \int_a^b f(x)dx=\frac{h}{2}\sum_{k=0}^{n-1}f(x_k)+f(x_{k+1}) ∫abf(x)dx=2h∑k=0n−1f(xk)+f(xk+1)
复化辛普森公式 ∫ a b f ( x ) d x = h 6 ∑ k = 0 n − 1 f ( x k ) + 4 ∗ f ( x k + 1 2 ) + f ( x k + 1 ) \int_a^b f(x)dx=\frac{h}{6}\sum_{k=0}^{n-1}f(x_k)+4*f(x_{k+\frac{1}{2}})+f(x_{k+1}) ∫abf(x)dx=6h∑k=0n−1f(xk)+4∗f(xk+21)+f(xk+1)
接下来，注意一下哈， x k , x k + 1 x_k,x_{k+1} xk,xk+1都是区间的分点，所以我在累加的时候除了两端的分点只加了一次，其他分点都加了两次，而 f ( x k + 1 2 ) f(x_{k+\frac{1}{2}}) f(xk+21)是区间分点的中点也只加了一次。
所以式子可以再化简一下这样就是常见的两个复化求积公式了。
复化梯形公式 ∫ a b f ( x ) d x = h 2 [ f ( a ) + 2 ∑ k = 1 n − 1 f ( x k ) + f ( b ) ] \int_a^b f(x)dx=\frac{h}{2}[f(a)+2\sum_{k=1}^{n-1}f(x_k)+f(b)] ∫abf(x)dx=2h[f(a)+2∑k=1n−1f(xk)+f(b)]
复化辛普森公式 ∫ a b f ( x ) d x = h 6 [ f ( a ) + 2 ∑ k = 1 n − 1 f ( x k ) + ∑ k = 0 n − 1 4 ∗ f ( x k + 1 2 ) + f ( b ) ] \int_a^b f(x)dx=\frac{h}{6}[f(a)+2\sum_{k=1}^{n-1}f(x_k)+\sum_{k=0}^{n-1}4*f(x_{k+\frac{1}{2}})+f(b)] ∫abf(x)dx=6h[f(a)+2∑k=1n−1f(xk)+∑k=0n−14∗f(xk+21)+f(b)]
再来解析一下，f(a)和f(b)是区间的端点，只加了一次，所以系数是1，而 f ( x k + 1 2 ) f(x_{k+\frac{1}{2}}) f(xk+21)是区间分点的中点也只加了一次，所以系数也不变。
式子里的 2 ∑ k = 1 n − 1 f ( x k ) 2\sum_{k=1}^{n-1}f(x_k) 2∑k=1n−1f(xk)是区间的内分点，都加了两次所以前面系数带了2。
总之，复化方法就是这样。但是复化方法有一个不好的地方，就是他的步长h依然是我们随便选取的，我们能不能让他自己找到一个合适的步长h？
梯形公式递推化
基于这样的想法就有了递推化梯形公式。
原本我们将区间[a,b]进行n等分，有n+1个分点，可以得到近似公式 T n T_n Tn，如果 T n T_n Tn精度不够，我们就把步长折半，也就是把区间2n等分，得到 T 2 n T_{2n} T2n。
梯形公式递推化就是简化了 T 2 n T_{2n} T2n的计算。
我们是可以按复化梯形公式重新计算 T 2 n T_{2n} T2n的，但是这样比较麻烦。
梯形公式递推化注意到了我们把区间2n等分时其实新增的分点是原来分点的中点，而原来的n+1个分点可以继续利用。由此得到了一个 T n T_n Tn到 T 2 n T_{2n} T2n的递推公式。
T 2 n = T n 2 + b − a 2 n ∑ k = 1 n f [ a + ( 2 k − 1 ) b − a 2 n ] T_{2n}=\frac{T_n}{2}+\frac{b-a}{2n}\sum_{k=1}^nf[a+(2k-1)\frac{b-a}{2n}] T2n=2Tn+2nb−a∑k=1nf[a+(2k−1)2nb−a]
额呵呵，这公式一看就不是什么人话。
我们说人话，简单点。
来个例子就秒懂了。
计算 ∫ 0 1 sin x x d x \int_0^1\frac{\sin x}{x}dx ∫01xsinxdx
首先用梯形公式不等分计算 T 1 = f ( 0 ) + f ( 1 ) 2 = 0.9207355 T_1=\frac{f(0)+f(1)}{2}=0.9207355 T1=2f(0)+f(1)=0.9207355
然后用递推公式计算 T 2 = T 1 2 + 1 2 f ( 1 2 ) T_2=\frac{T_1}{2}+\frac{1}{2}f(\frac{1}{2}) T2=2T1+21f(21)
好，观察一下前面的式子，前面的 T 1 2 \frac{T_1}{2} 2T1在任何时候都是一样的，它永远是上一个值的一半。
然后后面的项 1 2 f ( 1 2 ) \frac{1}{2}f(\frac{1}{2}) 21f(21)你先别管 1 2 \frac{1}{2} 21。[0,1]二等分增加了1个分点 1 2 \frac{1}{2} 21，所以后面是新增加分点 f ( 1 2 ) f(\frac{1}{2}) f(21)的和。
那么下一次呢，把 [ 0 , 1 2 ] , [ 1 2 , 1 ] [0,\frac{1}{2}],[\frac{1}{2},1] [0,21],[21,1]再二等分增加分点 1 4 , 3 4 \frac{1}{4},\frac{3}{4} 41,43，所以之后就是 f ( 1 4 ) + f ( 3 4 ) f(\frac{1}{4})+f(\frac{3}{4}) f(41)+f(43)
那这个系数 1 2 \frac{1}{2} 21是怎么递推的？
这个系数是 b − a 2 n \frac{b-a}{2n} 2nb−a，分子是积分区间长度b-a，分母是2,4,6,8,10以此类推。
这下秒懂了吧。
但是递推化梯形公式还是有问题的，就是梯形公式收敛速度太慢了。
龙贝格算法
这就要龙贝格算法解决了。
龙贝格算法其实也是一种事后估计方法。但我想说的是，龙贝格算法更是一种理查德森外推。
这个理查德森外推是一种通用的数值方法，龙贝格算法是他的特例，之后你还会遇见的。
他是这样的，如果你有一个数值公式 T ( h ) T(h) T(h)，h是某个参数，不一定是步长，任何参数都可以。我们这个方法是一个非常通用的方法。
你对 T ( h ) T(h) T(h)进行泰勒展开得到误差估计
T ( h ) − I = α 1 h + α 2 h 2 + . . . + α n h n + . . . T(h)-I=\alpha_1 h+\alpha_2h^2+...+\alpha_nh^n+... T(h)−I=α1h+α2h2+...+αnhn+...
接下来我想提高这个数值公式的精度为 h 2 h^2 h2应该怎么做？
答案是把h折半，也就是把 h 2 \frac{h}{2} 2h代入公式得到
T ( h 2 ) − I = α 1 h 2 + α 2 h 2 4 + . . . + α n h n 2 n + . . . T(\frac{h}{2})-I=\alpha_1 \frac{h}{2}+\alpha_2\frac{h^2}{4}+...+\alpha_n\frac{h^n}{2^n}+... T(2h)−I=α12h+α24h2+...+αn2nhn+...
然后我们把h的一次方消掉，也就二式乘以2减去1式得到
2 T ( h 2 ) − T ( h ) 2 − 1 − I = − 1 2 α 2 h 3 − . . . \frac{2T(\frac{h}{2})-T(h)}{2-1}-I=-\frac{1}{2}\alpha_2 h^3-... 2−12T(2h)−T(h)−I=−21α2h3−...
也就是说，左边那个 2 T ( h 2 ) − T ( h ) 2 − 1 \frac{2T(\frac{h}{2})-T(h)}{2-1} 2−12T(2h)−T(h)与精确值I的误差是 h 2 h^2 h2阶的，所以 2 T ( h 2 ) − T ( h ) 2 − 1 \frac{2T(\frac{h}{2})-T(h)}{2-1} 2−12T(2h)−T(h)是精度更高的近似值，我们记作 T 1 ( h ) T_1(h) T1(h)
然后如法炮制可以把精度依次提高到h的三次方、四次方、五次方…
而且不难推算出通用的递推公式
T m ( h ) = 2 m T m − 1 ( h 2 ) − T m − 1 ( h ) ) 2 m − 1 T_m(h)=\frac{2^mT_{m-1}(\frac{h}{2})-T_{m-1}(h))}{2^m-1} Tm(h)=2m−12mTm−1(2h)−Tm−1(h))
它的误差是h的m+1次方阶的。
龙贝格算法就是理查德森外推法在求积公式里的特例或者说运用。
只不过在求积公式里 T ( h ) T(h) T(h)的泰勒展开是
T ( h ) − I = α 1 h 2 + α 2 h 4 + . . . + α n h 2 n + . . . T(h)-I=\alpha_1 h^2+\alpha_2h^4+...+\alpha_nh^{2n}+... T(h)−I=α1h2+α2h4+...+αnh2n+...
也就是没有奇数阶，所以我要升阶的话折半一次是
T ( h 2 ) − I = α 1 h 2 4 + α 2 h 4 16 + . . . + α n h 2 n 4 n + . . . T(\frac{h}{2})-I=\alpha_1 \frac{h^2}{4}+\alpha_2\frac{h^4}{16}+...+\alpha_n\frac{h^{2n}}{4^n}+... T(2h)−I=α14h2+α216h4+...+αn4nh2n+...
那么每次都是下式乘以4减去上式，
导致递推公式是
T m ( h ) = 4 m T m − 1 ( h 2 ) − T m − 1 ( h ) ) 4 m − 1 T_m(h)=\frac{4^mT_{m-1}(\frac{h}{2})-T_{m-1}(h))}{4^m-1} Tm(h)=4m−14mTm−1(2h)−Tm−1(h))
并且误差阶每次加2而不是加1。
当然，这是我的讲法，更常见的讲法是按误差的事后估计法来推导的，也就是下面这样：
我们先来计算一下误差，有误差才能事后估计嘛。
从简单的梯形公式开始。
按照拉格朗日插值余项，梯形公式的误差 I − T n = ∑ k = 0 n − 1 − h 3 f ′ ′ ( α k ) 12 I-T_n=\sum_{k=0}^{n-1}-\frac{h^3f''(\alpha_k)}{12} I−Tn=∑k=0n−1−12h3f′′(αk)
而这里为了消掉中值，我们认为h很小，从而 h f ′ ′ ( α k ) ≈ f ′ ( x k + 1 ) − f ′ ( x k ) hf''(\alpha_k) \approx f'(x_{k+1})-f'(x_k) hf′′(αk)≈f′(xk+1)−f′(xk)
也因此梯形公式的误差 I − T n = ∑ k = 0 n − 1 − h 3 f ′ ′ ( α k ) 12 ≈ − h 2 12 [ f ′ ( b ) − f ′ ( a ) ] I-T_n=\sum_{k=0}^{n-1}-\frac{h^3f''(\alpha_k)}{12}\approx \frac{-h^2}{12} [f'(b)-f'(a)] I−Tn=∑k=0n−1−12h3f′′(αk)≈12−h2[f′(b)−f′(a)]
这个误差和 h 2 h^2 h2是一个数量级的。
类似地，可以推出辛普森公式 S n S_n Sn的误差和 h 4 h^4 h4一个数量级，而在牛顿-科特斯公式里n=3的时候得到的三次插值公式对应的求积公式我们称为科特斯公式是 h 6 h^6 h6数量级的。（科特斯公式前面没让你用，但是这里确实需要了，你想记的话也可以回去记啦）
基于此我们进行误差事后估计。
对了，在这里 T n , S n , C n T_n,S_n,C_n Tn,Sn,Cn的角标n指的是复化公式里区间的段的数量。
我们这里说要变步长其实就是变区间的段数嘛。
复化梯形公式误差和 h 2 h^2 h2是一个数量级的，所以步长变为原来的一半也就是区间分段分为原来的2倍，则误差大概变成了原来的四分之一。
所以 I − T 2 n I − T n ≈ 1 4 \frac{I-T_{2n}}{I-T_n}\approx\frac{1}{4} I−TnI−T2n≈41
熟悉的式子。
我们又可以得到更精确的 T 2 n T_{2n} T2n的估计值 I = 4 3 T 2 n − 1 3 T n I=\frac{4}{3}T_{2n}-\frac{1}{3}T_n I=34T2n−31Tn了。
不过这里我要把式子写成
S n = 4 1 4 1 − 1 T 2 n − 1 4 1 − 1 T n S_n=\frac{4^1}{4^1-1}T_{2n}-\frac{1}{4^1-1}T_n Sn=41−141T2n−41−11Tn
没错，熟悉的结果，对复化梯形公式进行事后估计得到了复化辛普森公式。
接下来你应该都猜得到结果了。
C n = 4 2 4 2 − 1 S 2 n − 1 4 2 − 1 S n C_n=\frac{4^2}{4^2-1}S_{2n}-\frac{1}{4^2-1}S_n Cn=42−142S2n−42−11Sn
R n = 4 3 4 3 − 1 C 2 n − 1 4 3 − 1 C n R_n=\frac{4^3}{4^3-1}C_{2n}-\frac{1}{4^3-1}C_n Rn=43−143C2n−43−11Cn
额，这里 R n R_n Rn是复化龙贝格公式啦。
总之，现在按这样的方法进行加工，就可以逐步提高求积精度，这就是龙贝格算法。
你应该也能想到我要掏出这样的图了。
现在这样步长h是可以随着计算折半的，其实我也觉得到这里就很完美了。
高斯求积公式
谁知道高斯求积公式出现了。
高斯求积公式彻底解决了插值点的选取问题，用这个方法选取的一定是最优插值点，而且代数精度是最高的。
选取适当的 x 0 , x 1 , . . . , x n x_0,x_1,...,x_n x0,x1,...,xn这n+1个点，高斯求积公式可以达到至少2n+1的代数精度。
这里选择的这些插值点被称为高斯点。
如果选取了n个高斯点得到的高斯公式称为n点高斯公式。（注意我前面说的是n+1个点有2n+1阶代数精度，那n个点就有2n-1阶代数精度）
1点高斯公式就是中矩形公式。
我们来试着求求看两点公式高斯。
注意，在高斯公式里，我们只能求[-1,1]的积分。
如果你想求解[a,b]的积分，可以先换元。
做一个简单的变换，令 x = b − a 2 t + a + b 2 x=\frac{b-a}{2}t+\frac{a+b}{2} x=2b−at+2a+b，就可以把积分区间变成[-1,1]了。
现在想求两点高斯公式就是要找组合系数 A 1 , A 2 A_1,A_2 A1,A2和合适的分点 x 1 , x 2 x_1,x_2 x1,x2使得
∫ − 1 1 f ( x ) d x = A 1 f ( x 1 ) + A 2 f ( x 2 ) \int_{-1}^1f(x)dx=A_1f(x_1)+A_2f(x_2) ∫−11f(x)dx=A1f(x1)+A2f(x2)
并且两点高斯公式具有2*2-1阶代数精度，所以他对 1 , x , x 2 , x 3 1,x,x^2,x^3 1,x,x2,x3都精确成立。
这样你会得到一个四元三次方程组，根本求解不来啊亲。
还好，有数学家已经“注意到”高斯点的求解了。
n点高斯公式的n个高斯点，就是n次勒让德多项式的零点。
所以真正实用的求积方法是，想用n点高斯公式先写出n次勒让德多项式，再求出其零点，然后零点就是求积点。现在已知了求积节点，不就回到了最开始的机械求积方法了？用那里的方法求出 A k A_k Ak就可以了。
那么这些结论都是怎么来的？
我们等下再说，因为下面我们要直接来一个加强版的。
反正到高斯求积公式这里我们的结论就是高斯点是n次勒让德多项式的零点，有了高斯点就可以用机械求积的方法得出组合系数 A k A_k Ak了。
带权高斯求积公式
所谓高斯求积公式的加强版就是增加了权函数的机制。
这是因为现实里有些函数啊他可能长这样：
∫ a b x e x 2 d x \int_{a}^bxe^{x^2}dx ∫abxex2dx
如果直接用高斯求积公式是可以求解的，但是效率不高。
于是又提出了带权高斯求积公式。
这里我们把x看做权函数 ρ ( x ) \rho(x) ρ(x)， e x 2 e^{x^2} ex2看做要求的积分 f ( x ) f(x) f(x)，于是我们要求的就相当于是带权函数的积分 ∫ a b ρ ( x ) f ( x ) d x \int_{a}^b\rho (x)f(x)dx ∫abρ(x)f(x)dx了。
我们前面说的高斯求积公式其实就是权函数为1，区间为[-1,1]的特殊情况，所以这确实是一个加强版的公式。
那么这个加权的想法是怎么来的呢？
其实还是从线性代数的角度想出来的。
在[快乐数学]傅里叶变换和拉普拉斯变换，链接https://www.bilibili.com/opus/953797131335893047?spm_id_from=333.1387.0.0里我们提到，区间上的连续函数的全体构成线性空间，并且我们可以定义内积结构使之变为内积空间，在那里我们说：
可以定义内积 ∫ a b f ( x ) g ( x ) d x \int_a^bf(x)g(x)dx ∫abf(x)g(x)dx为f(x)和g(X)的内积。
但其实，这个内积条件是可以加强的，我们可以定义带权内积 ∫ a b ρ ( x ) f ( x ) g ( x ) d x \int_a^b\rho (x)f(x)g(x)dx ∫abρ(x)f(x)g(x)dx，这样也能得到一个内积空间。
基于此就可以想到给高斯求积公式也整个带权版。
为啥呢？这是由高斯求积公式的性质决定的。事实上高斯求积公式和其对应的内积空间密切相关。
我们设 x i x_i xi是n点带权高斯公式 ∫ a b ρ ( x ) f ( x ) d x \int_{a}^b\rho (x)f(x)dx ∫abρ(x)f(x)dx的高斯点，那么可以定义一个n次多项式 ω n = ( x − x 1 ) ( x − x 2 ) . . . ( x − x n ) \omega_n=(x-x_1)(x-x_2)...(x-x_n) ωn=(x−x1)(x−x2)...(x−xn)，这个多项式就是由带权内积 ∫ a b ρ ( x ) f ( x ) g ( x ) d x \int_a^b\rho (x)f(x)g(x)dx ∫abρ(x)f(x)g(x)dx诱导的内积空间的一个正交多项式。
也就是说:
x 0 , x 1 , x 2 , . . . x n x_0,x_1,x_2,...x_n x0,x1,x2,...xn这n+1个点为高斯点当且仅当 ω n \omega_n ωn与任意次数小于等于n的多项式 P ( x ) P(x) P(x)带权正交。
我们来证明一下。
先假设 x 0 , x 1 , x 2 , . . , x n x_0,x_1,x_2,..,x_n x0,x1,x2,..,xn这n+1个点是带权高斯公式 ∫ a b ρ ( x ) f ( x ) d x \int_{a}^b\rho (x)f(x)dx ∫abρ(x)f(x)dx的高斯点，那么带权高斯公式 ∫ − 1 1 ρ ( x ) f ( x ) d x = ∑ k = 0 n A k f ( x k ) \int_{-1}^1\rho (x)f(x)dx=\sum_{k=0}^nA_kf(x_k) ∫−11ρ(x)f(x)dx=∑k=0nAkf(xk)就应该有2n+1次代数精度。
这是带权高斯公式能达到的最低代数精度嘛，没问题的。
这也就意味着，这个带权高斯公式对任意次数不小于2n+1的多项式都严格成立，没有任何误差。
那我就取 f ( x ) = P m ( x ) f(x)=P_m(x) f(x)=Pm(x)了，这里的 P m ( x ) P_m(x) Pm(x)就是定理里说的任意次数小于等于n的多项式。
既然 P m ( x ) P_m(x) Pm(x)次数小于等于n，而 ω ( x ) \omega(x) ω(x)次数是n+1(这里有n+1个高斯点，所以次数是n+1而不是前面说的n)，那么 ω ( x ) P m ( x ) \omega(x)P_m(x) ω(x)Pm(x)次数就小于等于2n+1，刚好在带权高斯公式严格成立的范围内。
现在要证明的是 ω ( x ) \omega(x) ω(x)与之带权正交，所以就是要算一下内积呗。
由于带权高斯公式严格成立没有误差所以 ∫ a b ρ ( x ) ω ( x ) P m ( x ) d x = ∑ k = 0 n A k ω ( x k ) P m ( x k ) \int_{a}^b\rho (x)\omega(x)P_m(x)dx=\sum_{k=0}^nA_k\omega(x_k)P_m(x_k) ∫abρ(x)ω(x)Pm(x)dx=∑k=0nAkω(xk)Pm(xk)
由于 ω ( x ) \omega(x) ω(x)的零点刚好就是 x k x_k xk，所以等式右边是0，这也就意味着 ω n \omega_n ωn与任意次数小于等于n的多项式 P ( x ) P(x) P(x)带权正交。
那反过来，如果我找到了若干点 x k x_k xk把他拼成 ω ( x ) \omega(x) ω(x)，并且满足 ω n \omega_n ωn与任意次数小于等于n的多项式 P ( x ) P(x) P(x)带权正交，只要能证明求积公式 ∫ a b ρ ( x ) P m ( x ) d x = ∑ k = 0 n A k P m ( x k ) \int_{a}^b\rho (x)P_m(x)dx=\sum_{k=0}^nA_kP_m(x_k) ∫abρ(x)Pm(x)dx=∑k=0nAkPm(xk)有至少2n+1次代数精度就说明 x k x_k xk是高斯点了。
这是可以证明的。
我们用 ω ( x ) \omega(x) ω(x)除以 P m ( x ) P_m(x) Pm(x)，得到 P m ( x ) = ω ( x ) q ( x ) + r ( x ) P_m(x)=\omega(x)q(x)+r(x) Pm(x)=ω(x)q(x)+r(x)，由于r(x)是余项所以他的次数小于除式 ω ( x ) \omega(x) ω(x)的次数n+1也就是小于等于n，而被除式 P m P_m Pm次数小于等于2n+1，所以商式的次数小于等于2n+1-除式次数(n+1)=n。
那么我们把这个带余除法的式子代入就有
∫ a b ρ ( x ) P m ( x ) d x = ∫ a b ρ ( x ) ω ( x ) q ( x ) d x + ∫ a b ρ ( x ) r ( x ) d x \int_{a}^b\rho (x)P_m(x)dx=\int_{a}^b\rho(x)\omega(x)q(x)dx+\int_{a}^b\rho(x)r(x)dx ∫abρ(x)Pm(x)dx=∫abρ(x)ω(x)q(x)dx+∫abρ(x)r(x)dx
我们条件里说 ω n \omega_n ωn与任意次数小于等于n的多项式 P ( x ) P(x) P(x)带权正交，而q(x)次数小于等于n并且前一个积分恰好就是 ω ( x ) \omega_(x) ω(x)与q(x)的带权内积，所以前一个积分是0。
接着看后一个积分，由于高斯求积公式还是插值型的求积公式，所以后一个积分仍然至少有n次代数精度，而后一个积分的被积函数次数小于等于n所以他的求积公式仍然是准确的，因此有 ∫ a b ρ ( x ) r ( x ) d x = ∑ k = 0 n A k r ( x k ) \int_{a}^b\rho(x)r(x)dx=\sum_{k=0}^nA_kr(x_k) ∫abρ(x)r(x)dx=∑k=0nAkr(xk)
注意哦，你看我们的带余除法的式子：
P m ( x ) = ω ( x ) q ( x ) + r ( x ) P_m(x)=\omega(x)q(x)+r(x) Pm(x)=ω(x)q(x)+r(x)
你把 x k x_k xk代入 ω ( x k ) \omega(x_k) ω(xk)是0，所以 P m ( x k ) = r ( x k ) P_m(x_k)=r(x_k) Pm(xk)=r(xk)因此
∫ a b ρ ( x ) P m ( x ) d x = ∫ a b ρ ( x ) ω ( x ) q ( x ) d x + ∫ a b ρ ( x ) r ( x ) d x = ∑ k = 0 n A k r ( x k ) = ∑ k = 0 n A k P m ( x k ) \int_{a}^b\rho (x)P_m(x)dx=\int_{a}^b\rho(x)\omega(x)q(x)dx+\int_{a}^b\rho(x)r(x)dx=\sum_{k=0}^nA_kr(x_k)=\sum_{k=0}^nA_kP_m(x_k) ∫abρ(x)Pm(x)dx=∫abρ(x)ω(x)q(x)dx+∫abρ(x)r(x)dx=∑k=0nAkr(xk)=∑k=0nAkPm(xk)
也就是说，现在我们的求积公式有至少2n+1阶代数精度，所以 x k x_k xk就是高斯点。
有这层关系在我们才有了求高斯点的方法。
在带权内积 ∫ a b ρ ( x ) f ( x ) g ( x ) d x \int_a^b\rho (x)f(x)g(x)dx ∫abρ(x)f(x)g(x)dx诱导的内积空间中总是有基函数 1 , x , x 2 , . . . , x n , . . . 1,x,x^2,...,x^n,... 1,x,x2,...,xn,...。
虽然这组基不是正交基但是我们在高等代数里学过如何用一组基得出正交基，甚至得出标准正交基。
这个方法就是施密特正交化。
如果我们取权函数 ρ ( x ) = 1 \rho(x)=1 ρ(x)=1，区间[a,b]为[-1,1]，这就是最开始说的高斯求积公式嘛，对基函数 1 , x , x 2 , . . . , x n , . . . 1,x,x^2,...,x^n,... 1,x,x2,...,xn,...进行施密特正交化就能得到新的一组正交基，然后再进行标准化就能得到标准正交基，而这组标准正交基就是所谓的勒让德多项式了。
对于线性空间的一组正交基来说，它与任意其他基之间都是正交的。
而你从中选一个n次多项式 P n ( x ) P_n(x) Pn(x)出来，再从内积空间中随便选一个次数小于n的多项式g(x)，它肯定不是由这个 P n ( x ) P_n(x) Pn(x)线性表出的嘛，所以他一定与 P n ( x ) P_n(x) Pn(x)线性无关，从而 P n ( x ) P_n(x) Pn(x)就与任意次数小于n的多项式正交了。
OK，你现在再回去看定理。
x 0 , x 1 , x 2 , . . . x n x_0,x_1,x_2,...x_n x0,x1,x2,...xn这n+1个点为高斯点当且仅当 ω n \omega_n ωn与任意次数小于等于n的多项式 P ( x ) P(x) P(x)带权正交。
而内积空间中的正交基中的任意一个n+1次多项式都和次数小于n+1也就是次数小于等于n的多项式带权正交。
由于 ω ( x ) = ( x − x 1 ) ( x − x 2 ) . . . ( x − x n ) \omega(x)=(x-x_1)(x-x_2)...(x-x_n) ω(x)=(x−x1)(x−x2)...(x−xn)的最高次系数是1，所以这个 ω ( x ) \omega(x) ω(x)是啥啊？
不就是你从内积空间中的标准正交基里次数为n+1次的多项式吗？
终于结束了。
也就是说，你想求n+1点带权高斯公式 ∫ a b ρ ( x ) f ( x ) d x \int_{a}^b\rho (x)f(x)dx ∫abρ(x)f(x)dx的高斯点，就是要找 ω ( x ) \omega(x) ω(x)的零点，而 ω ( x ) \omega(x) ω(x)就是由带权内积 ∫ a b ρ ( x ) f ( x ) g ( x ) d x \int_a^b\rho (x)f(x)g(x)dx ∫abρ(x)f(x)g(x)dx诱导的内积空间的标准正交基里的n+1次多项式 P ( x ) P(x) P(x)，所以要找高斯点就是找 P ( x ) P_(x) P(x) 的零点。
对于前面说的 ρ ( x ) = 1 \rho(x)=1 ρ(x)=1这样的特例带权内积 ∫ a b ρ ( x ) f ( x ) g ( x ) d x \int_a^b\rho (x)f(x)g(x)dx ∫abρ(x)f(x)g(x)dx诱导的内积空间的标准正交基就是勒让德多项式，我从勒让德多项式里取次数为n+1 的多项式，求他的零点，就得到高斯点了。
那如果是最一般的带权高斯公式呢？
方法一样的，我还是先找到带权内积 ∫ a b ρ ( x ) f ( x ) g ( x ) d x \int_a^b\rho (x)f(x)g(x)dx ∫abρ(x)f(x)g(x)dx诱导的内积空间的标准正交基。那怎么找呢？
高等代数里已经学习过了，将基 1 , x , x 2 , . . . , x n , . . . 1,x,x^2,...,x^n,... 1,x,x2,...,xn,...进行施密特正交化得到一组正交基。
然后呢，我再进行标准化，就得到标准正交基了。
然后再取标准正交基里次数为n+1的多项式，再求他的零点即可。
但是你仔细想想，真的需要标准化吗？
其实不用，所谓的标准化就是在原来正交基的基础上给基除以其长度嘛。
但是我只要求零点就可以了，你基缩放个常数倍影响零点吗？不影响的。
所以其实根本就不用什么标准化，只需要正交化然后算零点就可以了。
那反过来 ρ ( x ) = 1 \rho(x)=1 ρ(x)=1的情况我一定要去记整个勒让德多项式吗？
其实也不用的。
勒让德多项式是什么 n ! ( 2 n ) ! d n d x n [ ( x 2 − 1 ) n ] \frac{n!}{(2n)!}\frac{d^n}{dx^n}[(x^2-1)^n] (2n)!n!dxndn[(x2−1)n]
前面有一个系数 n ! ( 2 n ) ! \frac{n!}{(2n)!} (2n)!n!，这玩意存不存在影响零点吗？不影响的，这是标准化添加的额外的系数，我只求零点的话不要也罢。
所以如果我想求n点高斯公式的高斯点，只要计算多项式 d n d x n [ ( x 2 − 1 ) n ] = 0 \frac{d^n}{dx^n}[(x^2-1)^n]=0 dxndn[(x2−1)n]=0的零点就可以了。
这样记忆量会下降一点。
但是话又说回来了，真的每次都要这么麻烦吗？
其实也不用。
其实啥也不用记。
我们可以直接利用刚刚的性质求解高斯点。
例如我想找 ∫ 0 1 x f ( x ) d x \int_0^1\sqrt xf(x)dx ∫01xf(x)dx的两点高斯求积公式。
我们是要构造 ω ( x ) = ( x − x 1 ) ( x − x 2 ) \omega(x)=(x-x_1)(x-x_2) ω(x)=(x−x1)(x−x2)求零点的嘛，显然 ω ( x ) \omega(x) ω(x)是一个首项系数为1的二次多项式，我们直接设为 x 2 + b x + c x^2+bx+c x2+bx+c。
然后回去看看性质。
x 0 , x 1 , x 2 , . . . x n x_0,x_1,x_2,...x_n x0,x1,x2,...xn这n+1个点为高斯点当且仅当 ω n \omega_n ωn与任意次数小于等于n的多项式 P ( x ) P(x) P(x)带权正交。
所以 ω ( x ) \omega(x) ω(x)与多项式1，x带权正交。
而这里权函数是 x \sqrt x x，既然 ω ( x ) \omega(x) ω(x)与1，x带权正交的话就有：
∫ 0 1 x ω ( x ) 1 d x = 0 \int_0^1\sqrt x \omega(x)1dx=0 ∫01xω(x)1dx=0
∫ 0 1 x ω ( x ) x d x = 0 \int_0^1\sqrt x \omega(x)xdx=0 ∫01xω(x)xdx=0
这不是个含有两个未知量的两个方程吗？我直接解方程不就完事了。
解得 b = − 10 9 , c = 5 21 b=-\frac{10}{9},c=\frac{5}{21} b=−910,c=215
然后求 x 2 − 10 9 x + 5 21 = 0 x^2-\frac{10}{9}x+\frac{5}{21}=0 x2−910x+215=0的根就得到两个高斯点为 0.821162 , 0.289949 0.821162,0.289949 0.821162,0.289949了。
接下来你可以用前面的方法求解系数 A 1 , A 2 A_1,A_2 A1,A2了。
扫一扫
4322
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
