计算定积分的插值方法

hbghlyj

矩形公式与梯形公式
把定积分考虑作面积,解释作和的极限$$\int_a^bf(x)~\mathrm dx=\lim\sum_{i=1}^nf(\xi_i)(x_i-x_{i-1})$$以下我们总是把区间$(a,b)$分为$n$等份,把每一部分的长度记作$h$,于是$$h=\frac{b-a}n,x_i=a+ih\quad(x_0=a,x_n=a+nh=b)$$当$x=x_i$时,被积函数$y=f(x)$的值记作$y_i(i=0,1,2,\cdots,n)$,有$y_i=f(x_i)=f(a+ih)$这些量我们算作已知的,若函数$f(x)$是由解析式给定的,它们可以直接由计算得到;若是由图示法给定的,可以由图上量出来.令$\xi_i=x_{i-1}$或$x_i$,我们得到两个求近似值的矩形公式$$\int_a^bf(x)~\mathrm dx\sim\frac{b-a}n(y_0+y_1+\cdots+y_{n-1})$$$$\int_a^bf(x)~\mathrm dx\sim\frac{b-a}n(y_1+y_2+\cdots+y_n)$$其中记号"~"表示近似相等.
数目$n$越大时,就是$h$越小时,这两个公式越准确,当$n\to\infty$,而$h\to0$时,取极限,就达到这定积分的正确的量.如此,当纵坐标的数目渐增时,这两个矩形公式的误差趋向零.
当给定的函数$f(x)$在区间$(a,b)$上是单调的时,由纵坐标的数目可以很简单的确定出误差的一个上限.在这情形下,由图显然可以直接看出,两个矩形公式中任何一个的误差不超过有斜线的矩形的面积和,就是不超过以$\frac{b-a}{n}=h$为底,以$y_{n}-y_{0}$为高的有斜线的矩形的面积,也就是不超过量$\frac{b-a}{n}(y_n-y_0)$.矩形公式用矩形周界中一部分的横竖线段组成的阶形折线来代替曲线$y=f(x)$,就得到所要的面积的近似表达式.

若用与给定的曲线所差很少的其他的线来代替这样的阶形折线,就得到另外的近似表达式.这辅助线与曲线$y=f(x)$越近,用这辅助线下的面积来代替所要的面积时,误差越小.例如,若用给定的曲线的内接折线,当$x=x_i$时,它的纵坐标与曲线的纵坐标全同,换句话说,就是用这些内接梯形的面积和来代替所要的面积,就得到近似梯形公式$$\int_{a}^{b} f(x)~\mathrm dx\sim h\left[\frac{y_{0}+y_{1}}{2}+\frac{y_{1}+y_{2}}{2}+\cdots+\frac{y_{n-1}+y_{n}}{2}\right]=\frac{b-a}{2 n}\left[y_{0}+2 y_{1}+2 y_{2}+\cdots+2 y_{n-1}+y_{n}\right]$$
切线公式与波塞恩公式
现在我们把分份的数目增多一倍,就是把每一份再分成两半.这样作就得到$2n$份$x_{0}, x_{\frac{1}{2}}=a+\frac{h}{2}, x_{1}=a+h, \cdots, x_{i}=a+i h,x_{i+\frac{1}{2}}=a+\left(i+\frac{1}{2}\right) h, \cdots, x_{n}=b$各对应于纵坐标$y_{0}, y_{\frac{1}{2}}, y_{1}, \cdots, y_{i}, y_{i+\frac{1}{2}}, \cdots, y_{n}$.我们把$y_{0}, y_{1}, \cdots, y_{n}$叫做偶纵坐标,把$y_{\frac{1}{2}}, y_{\frac{3}{2}}, \cdots, y_{n-\frac{1}{2}}$叫做奇纵坐标.

过每一个奇纵坐标的端点作切线,与相邻的两个偶纵坐标相交,用这样作成的诸梯形的面积和来代替所要的面积.如此得到的近似公式叫作切线公式.$$\int_{a}^{b} f(x)~\mathrm{d} x \sim \frac{b-a}{n}\left[y_{\frac{1}{2}}+y_{\frac{3}{2}}+\cdots+y_{n-\frac{1}{2}}\right]=:\sigma_{1}$$再考虑联结诸相邻奇纵坐标端点的弦得到的内接梯形,加上联结纵坐标$y_0$与$y_{\frac12},y_{n-\frac12}$与$y_n$端点的两个弦在两端作成的梯形.这样得到的梯形的面积和记作$$\sigma_2:=\frac{b-a}{2 n}\left[\frac{y_{0}+y_{n}}{2}-\frac{y_{\frac{1}{2}}+y_{n-\frac{1}{2}}}{2}+2 y \frac{1}{2}+2 y \frac{3}{2}+\cdots+2 y_{n-\frac{1}{2}}\right]$$若曲线$y=f(x)$在区间$(a,b)$上没有扭转点,就是单纯的向上凹或单纯的向下凹,则曲线下面积$S$在$\sigma_1$与$\sigma_2$之间,所以自然可以取等差中项$\frac{\sigma_1+\sigma_2}2$作$S$的近似表达式,这是波恩塞公式$$\int_a^bf(x)~\mathrm dx\sim\frac{b-a}{2 n}\left[\frac{y_{0}+y_{n}}{4}-\frac{y_{\frac{1}{2}}+y_{n-\frac{1}{2}}}{4}+2 y_{\frac{1}{2}}+2 y_{\frac{3}{2}}+\cdots+2 y_{n-\frac{1}{2}}\right]$$不难看出,在上述关于曲线的假定下,这个公式的误差不超过$\frac{\sigma_{1}-\sigma_{2}}{2}=\left(\frac{y_{\frac{1}{2}}+y_{n-\frac{1}{2}}}{2}-\frac{y_{0}+y_{n}}{2}\right) \frac{b-a}{4 n}$的绝对值,由梯形中线的性质,不难证明,上式中括号内的表达式等于联结两个极端纵坐标与联结两个极端奇纵坐标的弦在正中纵坐标上截下的一段之长.
辛普森公式
如上,分作双数部分,过每三个纵坐标$y_{0}, y_{\frac{1}{2}}, y_{1} ; y_{1}, y_{\frac{3}{2}}, y_{2} ; \cdots ; y_{n-1}, y_{n-\frac{1}{2}}, y_{n}$的端点各作二次抛物线,用这一组抛物线来代替给定的曲线.我们先推导一个公式:
点$A_i(x_i,y_i),i=1,2,3$在抛物线$y=ax^2+bx+c$上,且$x_1,x_2,x_3$成公差为$d$的等差数列,则抛物线与$x$轴与$x=x_1$与$x=x_3$围成的曲边梯形的面积为$\int_{x_1}^{x_3}(ax^2+bx+c)~\mathrm dx=\frac13a(x_3^3-x_1^3)+\frac12b(x_3^2-x_1^2)+c(x_3-x_1)=\frac13a(x_3-x_1)(x_3^2+x_3x_1+x_1^2)+\frac12b(x_3-x_1)(x_3+x_1)+c(x_3-x_1)$
$=\frac d3(2a(x_3^2+x_3x_1+x_1^2)+3b(x_3+x_1)+6c)=\frac d3(a(x_3^2+(x_3+x_1)^2+x_1^2)+b(x_3+x_1)+4b\frac{x_3+x_1}2+6c)=\frac d3(a(x_3^2+4x_2^2+x_1^2)+b(x_3+4x_2+x_1)+6c)=\frac d3(y_1+4y_2+y_3)$.
由这个公式算出这样作成的每一条抛物线下的面积,我们就得到辛普森公式$$\int_a^bf(x)~\mathrm dx\sim\frac{b-a}{6 n}\left[y_{0}+4 y \frac{1}{2}+2 y_{1}+4 y_{\frac{3}{2}}+2 y_{2}+\cdots+2 y_{n-1}+4 y_{n-\frac{1}{2}}+y_{n}\right]$$关于上面的作法,我们用到一个事实:总可以选择适当的$a,b$与$c$,使抛物线$ax^2+bx+c$通过给定的平面上横坐标不同的三个点.这个由线性方程组的理论很容易证明.
实用上,曲线很平时,会得到很准确的结果,在曲线骤然改变形式的点的附近,应当计算得比较精细,那就需要把每份分得更小.
当作近似计算时,安排分划的计划是最重要的.为此,并且为比较上述几个不同的近似公式的准确性,我们介绍下面几个例子.$$S=\int_{0}^{\frac{\pi}{2}} \sin x \mathrm{~d} x=1$$$n=10, \frac{b-a}{n}=0.15707963, \frac{b-a}{2 n}=0.07853981, \frac{b-a}{6 n}=0.02617994$
矩形公式(短) S≈0.9194080
矩形公式(长) S≈1.0765828
切线公式       S≈1.0010290
梯形公式       S≈0.9979430
波恩塞公式   S≈0.9995487
辛普森公式   S≈1.0000000$$S=\int_{0}^{1} \frac{\ln (1+x)}{1+x^{2}} \mathrm{~d} x=\frac{\pi}{8} \ln 2=0.2721982613 \cdots$$$n=10, \frac{b-a}{2 n}=\frac{1}{20}, \frac{b-a}{6 n}=\frac{1}{60}$
波恩塞公式   S≈0.2719918
辛普森公式   S≈0.27219846$$S=\int_{0}^{1} \frac{\mathrm{d} x}{1+x}=\ln 2=0.69314718 \cdots$$$n=20, \frac{b-a}{2 n}=\frac{1}{40}, \frac{b-a}{6 n}=\frac{1}{20}$
梯形公式      S≈0.69330333
波恩塞公式  S≈0.69318196
辛普森公式  S≈0.69314716

---

拟柱体的体积公式
拟柱体是指所有的顶点都在两个平行平面中的多面体.拟柱体的侧面可能是三角形,梯形或平行四边形.一般的柱体,棱台,台塔,球台等都属于拟柱体.
下面给出一般拟柱体(不包括台塔)体积的计算公式:$$V=\frac h6(S_{h}+S_{0}+4S_{\frac {h}{2}})$$其中,$h$为高,${\displaystyle S_{\frac {h}{2}}}$在高度${\displaystyle {\frac {h}{2}}}$平行于底面的截面积;${\displaystyle S_{h}}$,高度$h$,就是顶面;$S_{0}$,高度0,就是底面.其来源为对不超过三次的多项式,以辛普森积分法求定积分之结果,‎对于拟柱体而言,截面积是距离底面的高度‎‎的不超过二次的多项式‎‎‎.

---

Newton–Cotes quadrature
It is assumed that the value of a function $f$ defined on $[a,b]$ is known at $n+1$ equally spaced points $a \leq x_0 < x_1 < \dots < x_n \leq b$. There are two classes of Newton–Cotes quadrature: they are called "closed" when $x_0 = a$ and $x_n = b$, i.e., they use the function values at the interval endpoints, and "open" when $x_0 > a$ and $x_n < b$, i.e., they do not use the function values at the endpoints. Newton–Cotes formulas using $n+1$ points can be defined (for both classes) as$$\int_a^b f(x) \,dx \approx \sum_{i=0}^n w_i\, f(x_i)$$where

• for a closed formula, $x_i = a + i h$, with $h = \frac{b-a}{n}$,
• for an open formula, $x_i = a + (i+1) h$, with $h = \frac{b-a}{n+2}$.

The number $h$ is called ''step size'', $w_i$ are called ''weights''.
The weights can be computed as the integral of Lagrange basis polynomials. They depend only on $x_i$ and not on the function $f$.
Let $L(x)$ be the interpolation polynomial in the Lagrange form for the given data points $(x_0, f(x_0)), \dots ,(x_n, f(x_n))$, then$$\int_a^b f(x) \,dx \approx \int_a^b L(x)\,dx = \int_a^b \left( \sum_{i=0}^n f(x_i)\, l_i(x) \right) \, dx
= \sum_{i=0}^n f(x_i) \underbrace{\int_a^b l_i(x)\, dx}_{w_i}.$$Closed Newton–Cotes formulas
This table lists some of the Newton–Cotes formulas of the closed type. For $0\leq i\leq n$, let $x_{i}=a+i{\tfrac {b-a}{n}}=a+ih$, and the notation $f_{i}$ be a shorthand for $f(x_{i})$.

$n$	Step size $h$	Common name	Formula	Error term
1	$b-a$	Trapezoidal rule	$\displaystyle {\frac {h}{2}}(f_{0}+f_{1})$	${\displaystyle -{\frac {1}{12}}h^{3}f^{(2)}(\xi )}$
2	${\frac  {b-a}{2}}$	Simpson's rule	${\frac {h}{3}}(f_{0}+4f_{1}+f_{2})$	${\displaystyle -{\frac {1}{90}}h^{5}f^{(4)}(\xi )}$
3	${\displaystyle {\frac {b-a}{3}}}{\frac  {b-a}{3}}$	Simpson's 3/8 rule	${\displaystyle {\frac {3h}{8}}(f_{0}+3f_{1}+3f_{2}+f_{3})}$	${\displaystyle -{\frac {3}{80}}h^{5}f^{(4)}(\xi )}$
4	${\displaystyle {\frac {b-a}{4}}}$	Boole's rule	${\displaystyle {\frac {2h}{45}}(7f_{0}+32f_{1}+12f_{2}+32f_{3}+7f_{4})}$	${\displaystyle -{\frac {8}{945}}h^{7}f^{(6)}(\xi )}$

Boole's rule is sometimes mistakenly called Bode's rule, as a result of the propagation of a typographical error in Abramowitz and Stegun, an early reference book.
The exponent of the segment size $h$ in the error term shows the rate at which the approximation error decreases. The order of the derivative of $f$ in the error term gives the lowest degree of a polynomial which can no longer be integrated exactly (i.e., with error equal to zero) with this rule. The number $\xi$  must be taken from the interval $(a,b)$.
Open Newton–Cotes formulas
This table lists some of the Newton–Cotes formulas of the open type. Again, $f_{i}$ is a shorthand for $f(x_{i})$, with ${\displaystyle x_{i}=a+i{\frac {b-a}{n+2}}}$.

n	Step size h	Common name	Formula	Error term
0	${\displaystyle {\frac {b-a}{2}}}$	Rectangle rule, or midpoint rule	${\displaystyle 2hf_{1}\,}$	${\displaystyle {\frac {1}{3}}h^{3}f^{(2)}(\xi )}$
1	${\displaystyle {\frac {b-a}{3}}}$	Trapezoid method	${\displaystyle {\frac {3}{2}}h(f_{1}+f_{2})}$	${\displaystyle {\frac {1}{4}}h^{3}f^{(2)}(\xi )}$
2	${\displaystyle {\frac {b-a}{4}}}$	Milne's rule	${\displaystyle {\frac {4}{3}}h(2f_{1}-f_{2}+2f_{3})}$	${\displaystyle {\frac {14}{45}}h^{5}f^{(4)}(\xi )}$
3	${\displaystyle {\frac {b-a}{5}}}$		${\displaystyle {\frac {5}{24}}h(11f_{1}+f_{2}+f_{3}+11f_{4})}$	${\displaystyle {\frac {95}{144}}h^{5}f^{(4)}(\xi )}$

战巡

........
这年头没什么人用这些方法求定积分的，要估算定积分用的也是用蒙特卡洛法

TripleIntegrals

战巡 发表于 2022-12-16 10:48
........
这年头没什么人用这些方法求定积分的，要估算定积分用的也是用蒙特卡洛法 ...

其实是有的。
计算机领域估算定积分一般会用到“自适应辛普森法”，就是对此的一个扩展。尤其是信息竞赛会用到。