数值积分是工程师和科学家经常使用的基本工具，用来计算无法解析求解的定积分的近似解。 
 如：$\Phi(x)=\int_0^x\frac{t^3}{e^t-1}dt$不存在$\Phi(x)$的解析解，要求$\Phi(5)$。 
 那么我们就要通过数值积分的方法来计算，数值积分的目的是，通过在有限个采样点上计算$f (x)$的值来逼近 $f (x)$在区间$[a,b]$上的定积分. 
 设$a=x_0<x_1<…<x_M=b$. 称形如 
 

 
 

 且具有性质

$\int_a^bf(x)dx=Q[f]+E[f]$的公式为数值积分或
面积公式。项

$E[ f ]$称为积分的截断误差，值

$\{x_k\}^M_{k=0}$称为
面积节点， 

$\{w_k\}^M_{k=0}$称为
权。
 
下面介绍几种常用的数值积分方法。
 
基于多项式差值的面积公式
 
通过M+1个等距点$\{(x_k,f(x_k))\}^M_{k=0}$存在唯一的次数小于等于M的多项式$P_M(x)$。当用该多项式来近似$[a,b]$上的$f(x)$时，$P_M(x)$的积分就近似等于$f(x)$的积分，这类公式称为牛顿－科特斯公式。当使用采样点$x_0=a和x_M=b$时，称为闭型牛顿－科特斯公式。
 
 
左/中/右矩形公式、梯形公式
 
左/中/右矩形公式 
 

 
 

 梯形公式 

 

 
 

 图形如下 

 

 

 
 
辛普森公式
 
推导过程
 
若$f(x)$在$[a,b]$上有定义，将区间$[a,b]$分割成$n$等分（取$n$为偶数），则有$a=x_0<x_1<…<x_M=b$，其中 
 
 $$x_i=a+i\Delta x（\forall i=0,1,\dots,n,\Delta x=\frac{b-a}{n}为步长）$$  

 这里我们想通过

$(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2))$三点抛物线

$g(x)=\alpha x^2+\beta x + \gamma$来取代

$f(x)$在

$[x_0,x_2]$的定义，今儿求出它的近似积分值(如图)，最后用累加的方式求得

$f(x)$在

$[a,b]$上的近似积分。 

 


 
 

 由假设我们有 

 

  
 
 
 
 

 所以可得

 $$\int_a^bf(x)dx=\int_{x_0}^{x_2}+\int_{x_2}^{x_4}+\dots+\int_{x_{n-2}}^{x_n}f(x)dx$$  

 

 

 
误差估计
 
若令$S_n=\frac{\Delta x}{3}[f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)\dots+2f(x_{n-2})+4f(x_{n-1})+f(x_n)]$且$f^{(4)}(x)在[a,b]$上连续，则我们可以估计出辛普森公式的误差值为 
 

 

 
例题1
 
试用辛普森公式估计$\int^1_0e^{-x^2}dx$，取$n=6$ 
 解： 
 令$f(x)=e^{-x^2},\Delta x={1\over 6}$则 
 

 

 
 
拉格朗日插值
 
分段线性插值
 
例如：函数$f(x)={1\over 1+x^2}$如果在区间$[-5，5]$上取$11$个等距节点:$x_k=-5+k (k=0,1,2,...,10)$,由Lagrange插值公式可得到$f(x)$的$10$次$L_{10}(x)$。如图所示： 
 
 
 
 

 

$L_{10}(x)$仅在区间的中部能较好的逼近函数

$f(x)$,在其它部位差异较大，而且越接近端点，逼近效果越差。 

 可以证明，当节点无限加密时，

$L_n(x)$也只能在很小的范围内收敛，这一现象称为
Runge现象。它表明通过增加节点来提高逼近程度是不适宜的，因而不采用高次多项式插值。
 
推导过程
 
已知函数$f(x)$在区间$[x_k ,x_{k+1} ]$的端点上的函数值$y_k =f(x_k ), y_{k+1} = f(x_{k+1} )$,求一个一次函数$y=P_1 (x)$使得$y_k =f(x_k ),y_{k+1} =f(x_{k+1} )$, 其几何意义是已知平面上两点$(x_k ,y_k ),(x_{k+1} ,y_{k+1} )$,求一条直线过该已知两点。 
 由直线的点斜式公式可知： 
 
 $$P_1(x)=yk+\frac{y_{k+1}-y_k}{x_{k+1}-x_k}(x-x_k)$$  

 把此式按照

$y_k 和y_{k+1}$写成两项: 

 

 $$P_1(x)=\frac{x-x_{k+1}}{x_k-x_{k+1}}y_k+\frac{x-x_{k}}{x_{k+1}-x_{k}}y_{k+1}$$  

 记 

 

 $$l_k=\frac{x-x_{k+1}}{x_k-x_{k+1}}\ \ \ \ \ l_{k+1}=\frac{x-x_{k}}{x_{k+1}-x_{k}}$$  

 并称它们为一次插值基函数。该基函数的特点如下： 

 

 $$l_i(x_k)=
\begin{cases}
1& \text{k=i}\\
0& \text{k!=i}
\end{cases}$$  

 对于

$i=0$， 

 

 $$l_0(x)=\frac{x-x_{1}}{x_0-x_{1}},x\in[x_0,x_1]$$  

 其它点上，

$l_0(x)=0$; 

 对于

$i=1,2,…,n-1$， 

 

 $$l_i(x)=
\begin{cases}
\frac{x-x_{i-1}}{x_i-x_{i-1}}& {x\in[x_{i-1},x_i]}\\
\frac{x-x_{i+1}}{x_i-x_{i+1}}& {x\in[x_{i},x_{i+1}]}
\end{cases}$$  

 其它点上，

$l_i(x)=0$; 对于

$i=n$， 

 

 $$l_n(x)=\frac{x-x_{n-1}}{x_n-x_{n-1}},x\in[x_{n-1},x_n]$$  

 其它点上，

$l_n(x)=0$. 于是， 

 

 $$P_1(x)=\sum_{k=0}^ny_kl_k(x)$$  

 此表达式与前面的表达式是相同的,这是因为在区间

$[x_k ,x_{k+1}]$上, 只有

$l_k(x),l_{k+1}(x)$是非零的，其它基函数均为零。 

 从而 

 

 $$P_1 (x) = y_k l_k (x) + y_{k+1} l_{k+1} (x)$$ 
 
此形式称之为拉格朗日型插值多项式。其中, 插值基函数与$y_k 、y_{k+1}$ 无关，而由插值结点$x_k 、x_{k+1}$ 所决定。一次插值多项式是插值基函数的线性组合, 相应的组合系数是该点的函数值$y_k 、y_{k+1}$ . 
 

 

 
误差估计
 
根据拉格朗日一次插值函数的余项,可以得到分段线性插值函数的插值误差估计: 对$x∈[a,b]$,当$x∈[x_k ,x_{k+1} ]$时, 
 
 $$R(x)=\frac{1}{2}f^{\prime\prime}(\xi)(x-x_k)(x-x_{k+1})$$  

 则

$R(x)\leq {h^2\over 8}M$，其中 

 

 $$h=\max_{0\leq k\leq n-1}|x_{k+1}-x_k|,M=\max_{x\in (a,b)}f^{\prime\prime}(x)$$  

 于是有下面的定理： 

 如果

$f(x)$在

$[a,b]$上二阶连续可微，则分段连续函数

$φ(x)$的余项有以下误差估计： 

 

 $$|R(x)|=|f(x)-\phi(x)|\leq {h^2\over 8}M$$  

 其中 

 

 $$h=\max_{0\leq k\leq n-1}|x_{k+1}-x_k|,M=\max_{x\in (a,b)}f^{\prime\prime}(x)$$  

 于是可以加密插值结点, 缩小插值区间, 使

$h$减小, 从而减小插值误差。




例题1

已知函数$y=f(x)={1\over 1+x^2}$在区间

$[0,5]$上取等距插值节点(如下表)， 

 求区间上分段线性插值函数,并利用它求出

$f(4.5)$近似值。 

 


 
 

 解: 

 在每个分段区间

$[k,k+1]$上， 

 

 $$P_1(x)=\frac{x-(k+1)}{k-(k+1)}y_k+\frac{x-(k)}{k+1-(k)}y_{k+1}=-y_k(x-k-1)+y_{k+1}(x-k)$$  

 于是 

 

 $$P_1(4.5)=-0.05882×(4.5-5)+0.03846×(4.5-4)=0.04864$$ 
 
拉格朗日型二次插值多项式
 
推导过程
 
已知函数$y=f(x)$在点$x_{k-1} ,x_k ,x_{k+1}$上的函数值$y_{k-1} =f(x_{k-1} ),y_k =f(x_k ), y_{k+1} =f(x_{k+1} )$, 求一个次数不超过二次的多项式$P_2 (x)$, 使其满足, 
 
 $$P_2 (x_{k-1} )=y_{k-1} , P_2 (x_k )=y_k , P_2 (x_{k+1} )=y_{k+1}$$  

 其几何意义为: 

 已知平面上的三个点 

 

$(x_{k-1} ,y_{k-1} ),(x_k ,y_k ),(x_{k+1} ,y_{k+1} )$, 

 求一个二次抛物线, 使得该抛物线经过这三点。
 
插值基本多项式
 
有三个插值结点$x_{k-1},x_k ,x_{k+1}$构造三个插值基本多项式，要求满足： 
 (1) 基本多项式为二次多项式； 
 (2) 它们的函数值满足下表： 
 

 
 

 因为

$l_{k-1} (x_k )= 0,l_{k-1} (x_{k+1} )=0$, 故有因子

$(x-x_k )(x-x_{k+1} )$, 而其已经是一个二次多项式, 仅相差一个常数倍, 可设 

 

 $$l_{k-1} (x)=a(x-x_k )(x-x_{k+1} )$$  

 又因为 

 

 $$l_{k-1} (x_{k-1})=1 \Rightarrow a(x_{k-1} -x_k )(x_{k-1} -x_{k+1})=1$$  

 得　 

 

 $$a={1\over (x_{k-1}-x_k)(x_{k-1}-x_{k+1})}$$  

 从而 

 

 $$l_{k-1}(x)={(x-x_k)(x-x_{k+1})\over (x_{k-1}-x_k)(x_{k-1}-x_{k+1})}$$  

 同理得 

 

 $$l_{k}(x)={(x-x_{k-1})(x-x_{k+1})\over (x_{k}-x_{k-1})(x_{k}-x_{k+1})}$$  

 

 $$l_{k+1}(x)={(x-x_{k-1})(x-x_{k})\over (x_{k+1}-x_{k-1})(x_{k+1}-x_{k})}$$ 
 
拉格朗日型二次插值多项式
 
由前述, 拉格朗日型二次插值多项式: 
 
 $$P_2 (x)=y_{k-1} l_{k-1} (x)+y_k l_k (x)+y_{k+1} l_{k+1} (x)$$  

 

$P_2 (x)$是三个二次插值多项式的线性组合,因而其是次数不超过二次的多项式,且满足: 

 

 $$P_2 (x_i )=y_i , (i=k-1,k,k+1)$$ 
 
例题2
 
已知：
 
$x_i$
10
15
20
$y_i=lgx_i$
1
1.1761
1.3010
 
利用此三值的二次插值多项式求$lg(12)$的近似值。 
 解：设$x_0 =10,x_1 =15,x_2 =20$,则: 
 
 
 
 
 
 

 故: 

 
 
 
 
  
 
 
  
 
 
 

 所以 

 
 
 
 
 
 
 
误差估计
 
我们在$[a,b]$上用多项式$P_n (x)$ 来近似代替函数$f(x)$, 其截断误差记作$R_n (x)=f(x)-P_n (x)$ 
 当$x$在插值结点$x_i$上时$R_n (x_i )=f(x_i )-P_n(x_i )=0$,下面来估计截断误差：
 
定理1： 
 设函数$y=f(x)$的$n$阶导数$y(n) =f^{(n)} (x)$在$[a,b]$上连续， 
 $y^{(n+1)} = f^{(n+1)} (x)$在$(a,b)$上存在； 
 插值结点为:$a≤x_0 <x_1 <…<x_n ≤b$, 
 $P_n (x)$是$n$次拉格朗日插值多项式； 

则对任意$x∈[a,b]$有： 
 
 $$R_n(x)={1\over (n+1)!}f^{(n+1)}(\xi)\omega_{n+1}(x)$$  

 其中

$ξ∈(a,b)$, 

$ξ$依赖于x：$ω_{n+1} (x)=(x-x_0 )(x-x_1 )…(x-x_n )$ 

 设

$\max_{a\leq x\leq b}|f^{n+1}(x)|\leq M_{n+1}$，则 

 

 $$|R_n(x)|\leq {1\over (n+1)!}M_{n+1}\omega_{n+1}(x)$$  

 易知，线性插值的截断误差为: 

 

 $$R_1(x)={1\over 2}f^{\prime \prime}(\xi)(x-x_0)(x-x_1)$$  

 二次插值的截断误差为: 

 

 $$R_2(x)={1\over 6}f^{(3)}(\xi)(x-x_0)(x-x_1)(x-x_2)$$ 
 
例题3
 
在例2中，用$lg10,lg15和lg20$计算$lg12$. 
 $P_2(12)=1.0766$, 
 $e = |1.0792-1.0766|=0.0026$ 
 估计误差： 
 
 
 
 
 

 
 
三次样条插值法
 
三次样条曲线原理
 

 $$x:a=x_0<x_1<…<x_M=b\\ y:y_0\ \ \ y_1\ \ \ \dots\ \ \ y_n$$ 
 
定义
 
样条曲线$S(x)$是一个分段定义的公式。给定$n+1$个数据点，共有$n$个区间，三次样条方程满足以下条件：


  
a. 在每个分段区间$[x_i,x_{i+1}]（i = 0, 1, …, n-1，x$递增$），S(x)=S_i(x)$ 都是一个三次多项式。
 
b. 满足$S(x_i)=y_i （i = 0, 1, …, n ）$
 
c.$S(x)$，导数$S^{\prime}(x)$，二阶导数$S^{\prime\prime}(x)$在$[a, b]$区间都是连续的，即$S(x)$曲线是光滑的
  
所以n个三次多项式分段可以写作： 
 
 $$S_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^3,i = 0, 1, …, n-1$$  

 其中

$a_i, b_i, c_i, d_i$代表

$4n$个未知系数。
 
求解
 
已知: 
 a. $n+1$个数据点$[x_i, y_i], i = 0, 1, …, n$ 
 b. 每一分段都是三次多项式函数曲线进行逼近 
 c. 节点达到二阶连续 
 d. 左右两端点处特性（自然边界，固定边界，非节点边界）
 
根据定点，求出每段样条曲线方程中的系数，即可得到每段曲线的具体表达式。
 
插值和连续性： 
 , 其中 $i = 0, 1, …, n-1$ 
 微分连续性： 
 , 其中$i = 0, 1, …, n-2$ 
 样条曲线的微分式： 
  
 将步长$h_i=x_{i+1}-x_i$带入样条曲线的条件： 
 a. 由$S_i(x_i)=y_i (i = 0, 1, …, n-1)$推出 
 
 $$a_i=y_i$$  

 b. 由

$S_i(x_{i+1})=y_{i+1} (i = 0, 1, …, n-1)$推出 

 

 $$a_i+h_ib_i+h_i^2c_i+h_i^3d_i=y_{i+1}$$  

 c. 由

$S^\prime_i(x_{i+1})=S^\prime_{i+1}(x_{i+1})(i = 0, 1, …, n-2)$推出 

 

 $$S_i^\prime(x_{i+1})=b_i+2c_i(x_{i+1}-x_i)+3d_i(x_{i+1}-x_i)^2=b_i+2c_ih+3d_ih^2\\ S_{i+1}^\prime(x_{i+1})=b_{i+1}+2c_i(x_{i+1}-x_{i+1})+3d_i(x_{i+1}-x_{i+1})^2=b_{i+1}$$  

 由此可得： 

 

 $$b_i+2h_ic_i+3h_i^2d_i-b_{i+1}=0$$  

 d. 由

$S_i^{\prime\prime}(x_{i+1})=S_{i+1}^{\prime\prime}(x_{i+1}) (i = 0, 1, …, n-2)$推出 

 

 $$2c_i+6h_id_i-2c_{i+1}=0$$ 
 
设$m_i=S_i^{\prime\prime}(x_i)=2c_i$，则 
 a. $2c_i+6h_id_i-c_{i+1}=0$可写为： 
 $m_i+6h_id_i-m_{i+1}=0$，推出 
 
 $$d_i={m_{i+1}-m_i\over6h_i}$$  

 b. 将

$c_i, d_i$带入 

$y_i+h_ib_i+h_i^2c_i+h_i^3d_i=y_{i+1}$可得： 

 

 $$b_i={y_{i+1}-y_i\over h_i}-{h_i\over 2}m_i-{h_i\over 6}(m_{i+1}-m_i)$$  

 c. 将

$b_i, c_i, d_i$带入

$b_i+2h_ic_i+3h_id_i=b_{i+1} (i = 0, 1, …, n-2)$可得： 

 

 $$h_im_i+2(h_i+h_{i+1})m_{i+1}+h_{i+1}m_{i+2}=6[{y_{i+2}-y_{i+1}\over h_{i+1}}-{y_{i+1}-y_i\over h_i}]$$ 
 
端点条件
 
由i的取值范围可知，共有n-1个公式， 但却有n+1个未知量m 。要想求解该方程组，还需另外两个式子。所以需要对两端点$x_0$和$x_n$的微分加些限制。 选择不是唯一的，3种比较常用的限制如下。
 
a. 自由边界(Natural)
 
首尾两端没有受到任何让它们弯曲的力，即$S^{\prime\prime}=0$。具体表示为$m_0=0$和 $m_n=0$ 
 则要求解的方程组可写为： 
 

  
 
 
  
 
 

 
b. 固定边界(Clamped)
 
首尾两端点的微分值是被指定的，这里分别定为A和B。则可以推出 
 
 
 
 
  
 
 
  
 
 
 

 将上述两个公式带入方程组，新的方程组左侧为 

 

 

 
c. 非节点边界(Not-A-Knot)
 
指定样条曲线的三次微分匹配，即 
 
 $$S_o^{\prime\prime\prime}(x_1)=S_1^{\prime\prime\prime}(x_1)\\ S	_{n-2}^{\prime\prime\prime}(x_{n-1})=S_{n-1}^{\prime\prime\prime}(x_{n-1})$$  

 根据

$S_i^{\prime\prime\prime}(x)=6d_i$和

$d_i={m_{i+1}-m_i\over 6h_i}$，则上述条件变为 

 

 $$h_1(m_1-m_0)=h_0(m_2-m_1)\\ h_{n-1}(m_{n-1}-m_{n-2})=h_{n-2}(m_n-m_{n-1})$$  

 新的方程组系数矩阵可写为： 

 


 

 
右下图可以看出不同的端点边界对样条曲线的影响： 
 
 
 
 
 

 
算法总结
 
假定有n+1个数据节点 
 
 $$(x_0,y_0),(x_1,y_1),(x_2,y_2),\dots,(x_n,y_n),$$  

 a. 计算步长

$h_i=x_{i+1}-x_i(i = 0, 1, …, n-1)$ 

 b. 将数据节点和指定的首位端点条件带入矩阵方程 

 c. 解矩阵方程，求得二次微分值

$m_i$。 

 d. 计算样条曲线的系数： 

 


 
 

 其中

$i = 0, 1, …, n-1$ 

 e. 在每个子区间

$x_i\leq x\leq x_{i+1}$中，创建方程 

 

 $$g_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^3$$ 

 

 

  
 
 
数值积分的基本思想
 
由积分中值定理可知，在积分区间 
 

   内存在一点 
 
 

   ，成立
 
 

 

  
 
 

  
 
 

 

   式的几何意义即为：底为 
 
 

   而高为 
 
 

   的矩形的面积恰等于所求曲边梯形的面积 
 
 

   。因此，要想求出 
 
 

   式左端积分，我们只需要知道三个值： 
 
 

   即可，这里 
 
 

   是显然的，问题在于 
 
 

   的具体位置一般是不清楚的，从而 
 
 

   未知。我们暂且将 
 
 

   称为区间 
 
 

   上的平均高度，我们的目标就是寻求一种求出平均高度 
 
 

   的算法，这样我们的问题也就解决了。
 
 

 
首先我们想到的是用区间两端点的“高度” 
 

   与 
 
 

   的算术平均值作为平均高度 
 
 

   的近似值，从而导出下述求积公式
 
 

 

  
 
 

  
 
 

 

   式即为我们我们熟知的
 
 梯形公式。写出MATLAB函数为
 

  
function
 
 
而若用区间中点 
 

   的“高度” 
 
 

   近似取代平均高度 
 
 

   ，则可导出
 
 中矩形公式（简称
 矩形公式）
 

 

  
 
 

  
 
 
其MATLAB函数为
 

  
function
 
 
总结：更一般地，我们可以在区间 
 

   上适当选取某些节点 
 
 

   ，然后用 
 
 

   的加权平均得到平均高度 
 
 

   的近似值，这样构造出来的求积公式具有下列通式
 
 

 

  
 
 

  
 
 
其中 
 

   称为
 
 求积节点， 
 

   称为
 
 求积系数，也称作伴随节点 
 

   的
 
 权。
 
这种数值积分方法通常称为机械求积，其优势在于将积分求值问题归结为被积函数值的计算，很适合在计算机上使用。
 
代数精度
 
由于数值积分方法是近似方法，为保证精度，自然希望我们的求积公式能够对“尽可能多”的函数准确成立，这便提出了代数精度的概念。
 

  如果某个求积公式对于次数不超过 
  

    的多项式均能够准确成立，但对于 
  
  

    次多项式就不准确成立，则称该求积公式具有
  
  
  

    次代数精度
  
  （或
  代数精确度）。
 
 
前面提到的梯形公式 
 

   和矩形公式 
 
 

   均具有一次代数精度。
 
 
求积公式的收敛性与稳定性
 

  在求积公式 
  

    中，若 
  
  

   
  
  

    其中 
  
  

    ，则称求积公式 
  
  

    是
  
  收敛的。
  
对任给 
  

    ，若 
  
  

    ，只要 
  
  

    就有 
  
  

   
  
  

    成立，其中 
  
  

    为计算 
  
  

    时产生误差 
  
  

    后实际得到的值，即 
  
  

    ，则称求积公式 
  
  

    是
  
  稳定的。
 
 
定理 
 

  
 
 若求积公式 
 

  
 
 中系数 
 

  
 
 ，则此求积公式是稳定的。
 
牛顿-柯特斯(Newton-Cotes)公式
 
将积分区间 
 

   划分为 
 
 

   等份，此时步长为 
 
 

   ，选取等距节点 
 
 

   构造出的插值型求积公式
 
 

 

  
 
 

  
 
 
称为牛顿-柯特斯(Newton-Cotes)公式，其中 
 

   称为
 
 柯特斯系数。可以通过下式确定
 

 

  
 
 

  
 
 
当 
   

     时， 
   
   

     这时的求积公式便是之前的梯形公式 
   
   

     。
   
当 
   

     时， 
   
   

    
   
   

     此时对应的求积公式便是
   
辛普森(Simpson)公式： 
   

    
   
   

     写成MATLAB函数为
   
 

  
function
 
 
当 
   

     时的牛顿-柯特斯公式则特别称为
   
柯特斯公式，其形式为 
   

    
   
   

     这里 
   
   

     。
   
 
其MATLAB代码为
 

  
function
 
 
复合求积公式
 
由于牛顿-柯特斯公式在 
 

   时不具有稳定性，故不可能再通过提高阶的方法来提高求积精度。
 
 复合求积法便是通过把积分区间分成若干个子区间（通常是等分），再在每个子区间上使用低阶求积公式的方法来提高精度的。
 
其实细心的同学可以发现，前面我所展示的MATLAB代码使用的便是对每个求积公式使用复合法完成的。常用的有复合梯形公式和复合辛普森公式。
 
龙贝格(Romberg)求积公式
 
前面介绍的复合求积方法可提高求积精度，如若精度仍不够，则可通过将步长逐次减半的方式来提高精度。如对复合梯形公式可导出其递推公式
 

 

  
 
 

  
 
 
其中 
 

   表示在 
 
 

   基础上步长 
 
 

   减半后的积分值。
 
 
定理 
 

  
 
  设 
 

   ，则有
 
 
 

 

  
 
 
 

   其中
 
  
 

   ，
 
 系数 
 

   与 
 
 
 

   无关。
 
 
 
从 
 

   式可以看出， 
 
 

   是 
 
 

   阶，若用 
 
 

   代替 
 
 

   ，有
 
 

  
 
 

  
 
 
再用 
 

   乘 
 
 

   式再减去 
 
 

   式后再除以 
 
 

   得
 
 

 

  
 
 

  
 
 
这里 
 

   是与 
 
 

   无关的系数。
 
 
从 
 

   式可看出，用近似积分值 
 
 

   ，其误差阶为 
 
 

   ，显然误差阶是提高了。类似这种将计算 
 
 

   的近似值的误差阶由 
 
 

   提高到 
 
 

   的方法称为
 
 外推算法，也称为
 理查德森(Richardson)外推算法。
 

  这是“数值分析”中一个重要的技巧，只要真值与近似值的误差能够表示成 
  

    的幂级数，如 
  
  

    式所示，都可以使用外推算法，提高精度。
  
 
 
龙贝格(Romberg)算法
 

 

   次外推加速为
 
 

 

  
 
 

  
 
 
余项为
 

 

  
 
 

  
 
 
此方法常称为理查德森外推加速方法。
 
设用 
 

   表示二分 
 
 

   次后求得的梯形值，且以 
 
 

   表示序列 
 
 

   的 
 
 

   次加速值，则由递推公式 
 
 

   可得
 
 

 

  
 
 

  
 
 
上次则称为龙贝格求积算法，计算过程如下：
 
取 
   

     ，求 
   
   

     。令 
   
   

     ( 
   
   

     记区间 
   
   

     的二分次数)。
   
求梯形值 
   

     ，即按递推公式 
   
   

     计算 
   
   

     。
   
求加速值。
若 
   

     (预先给定的精度)，则终止计算，并且 
   
   

     ；否则令 
   
   

     转 
   
   

     继续计算。
   
 
MATLAB代码为
 

  
function
 
 
自适应积分方法
 
复合求积方法通常适用于被积函数变化不太大的积分，如果在求积区间中被积函数变化很大，有的部分函数值变化剧烈，另一部分却变化平缓。这时如果统一将区间等分再用复合求积公式计算积分将会导致计算量很大，我们想实现的是在满足误差要求的前提下，对变化剧烈部分将区间细分，而平缓部分则可使用大步长，也即针对被积函数在区间上不同情形采用不同的步长，使得再满足精度前提下积分计算工作量尽可能小。其算法技巧是在不同区间上预测被积函数变化的剧烈程度确定相应步长。这就是自适应积分方法的基本思想。
 
下面是一个基于复合辛普森法的自适应积分算法的MATLAB代码：
 

  
function
 
 
下面我们计算积分 
 

   ，可输入下列语句调用上述MATLAB函数进行计算
 
 

  
clear

 
 
结果为
 

  
I 
 
 
高斯求积公式
 
下面研究带权积分 
 

   （ 
 
 

   为权函数）的求积公式
 
 

 

  
 
 

  
 
 
下面看定义
 

  如果求积公式 
  

    具有 
  
  

    次代数精度，则称其节点 
  
  

    为
  
  高斯点，相应公式 
  

    称为
  
  高斯型求积公式。
 
 
这里主要包括四种高斯型求积公式，即
 
高斯-勒让德求积公式
高斯-切比雪夫求积公式
高斯-拉盖尔求积公式
高斯-埃尔米特求积公式
 
多重积分
 
这里主要说一下多重积分计算的思路，归结为一句话，就是：多重积分化累次积分，再由里到外使用数值积分公式进行求积计算。

 
各种数值积分方法总结（龙贝格积分、自适应积分、高斯积分等）
 
1 (一重积分)常用的数值积分方法
1.1 牛顿-科茨（Newton-Cotes）积分公式
1.1.1 梯形法则（2点积分）
1.1.2-3 辛普森法则（3点积分和4点积分）
1.1.2 辛普森1/3法则（3点积分）
1.1.3 辛普森3/8法则（4点积分）
    
1.1.4 布尔法则（5点积分）
1.1.5牛顿-科茨积分公式总结
   
（穿插）积分改进思路1
1.2 复合积分（复合梯形积分、复合辛普森积分等）
1.2.1 复合梯形积分
1.2.2 复合辛普森积分
   
（穿插）积分改进思路2
1.3常用的高精度积分
1.3.1 龙贝格积分
（穿插）积分改进思路3
1.3.2 自适应积分
1.3.3 高斯积分
高斯-勒让德公式
高斯积分误差：
高斯-切比雪夫公式
广义积分：高斯-拉盖尔公式
广义积分：高斯-埃尔米特公式
   
  
  
2 多重积分方法
二重积分的辛普森公式
二重积分的蒙特卡洛方法
多重（高维）积分的蒙特卡洛方法
 

 
本文整理了各种类型的数值积分，从简单的梯形积分、辛普森积分到高精度的自适应积分、龙贝格积分和高斯积分；对于多重积分（高维积分），本文先简单介绍一下蒙特卡洛方法（适用于积分维数大于4的积分），后续会发帖专门针对多重积分方法做介绍。
 （如有疏漏，欢迎指正，谢谢~）
 
1 (一重积分)常用的数值积分方法
 
常用积分方法包括梯形积分、辛普森积分、自适应积分、龙贝格积分和高斯积分等；其中自适应积分、龙贝格积分和高斯积分属于高精度积分。
 
1.1 牛顿-科茨（Newton-Cotes）积分公式
 
牛顿-科茨（Newton-Cotes）积分公式：
 
 包含了梯形法则、辛普森法则、布尔法则等。
 
1.1.1 梯形法则（2点积分）
 
梯形法则积分计算公式：
 
其中，
 
 于是得：
 
 梯形公式要计算两个点的函数值f(a)和f(b)。所以说是两点积分
 
 梯形法则的误差为：
 
 
1.1.2-3 辛普森法则（3点积分和4点积分）
 
辛普森法则分为辛普森1/3法则和辛普森3/8法则，就是用二次/三次插值曲线来代替直线进行积分。
 
 
1.1.2 辛普森1/3法则（3点积分）
 
辛普森1/3法则法则积分计算公式：
 
 
1.1.3 辛普森3/8法则（4点积分）
 
辛普森3/8法则法则积分计算公式：
 
 
1.1.4 布尔法则（5点积分）
 
布尔法则与辛普森法则类似，只是更高阶的多项式积分，其公式为：
 
 
1.1.5牛顿-科茨积分公式总结
 

 注意：1.对于3点积分与4点积分具有相同的误差阶；5点积分与6点积分具有相同的误差阶。对于点数更多的积分也一样。因此优先使用奇数个点的积分。
 2.对于点数大于等于9的牛顿-科茨积分，求积公式的稳定性得不到保证（具体原因可以找数值积分的数看看，本文源自《现代数值积分》第5章 数值积分与数值微分），积分不能保证收敛，因此实际计算中一般不采用高阶牛顿-科茨积分公式。（处于效率考虑，一般也很少用超过5个点的牛顿-科茨积分公式）。
 
（穿插）积分改进思路1
 
牛顿-科茨积分公式思想：根据积分点构造近似的插值多项式，进而利用多项式积分表示原积分。
 对于高次的多项式插值，会出现龙格现象（对于高次的多项式插值，插值多项式会出现不收敛的现象，成为龙格现象，想深入了解可以去查询相关资料。并且次数越高，越容易出现不收敛现象）。
 
针对龙格现象，更好的选择切比雪夫节点来进行插值（想深入了解可以去查询相关资料，本文部分源自《现代数值积分》第3章 多项式插值与样条插值）。
 

 由于切比雪夫节点的特殊性，对于高阶插值的可以收到较好的效果，因此，对于“牛顿-科茨积分公式总结”中提到的点数大于等于9的牛顿-科茨积分，可以采用切比雪夫计算节点，然后构造插值多项式，再计算积分（此方法实现和推广比较麻烦，此处不做扩展，感兴趣的可以自己下去推导实现）。
 
1.2 复合积分（复合梯形积分、复合辛普森积分等）
 
上述主要介绍的是单个区间的积分。为了计算更准确，我们可以将区间划分为若干个小的区间，然后对每个小区间进行积分，然后对各个小区间的积分求和。由于对于单个小区间的积分可分为梯形法则，辛普森法则，和布尔法则，因此复合积分的每个小区间的积分可以是这些方法中的任意一种。
 下面主要讲一下 复合梯形积分和复合辛普森积分。
 
1.2.1 复合梯形积分
 
小区间的划分又分为等长区间和不等长区间两种划分方法。
 不等长区间，复合梯形积分计算公式：
 
 区间不等长情况应用场景很有限。
 
等长区间，复合梯形积分计算公式：
 
 等长区间，复合梯形积分计算误差：
 
 
1.2.2 复合辛普森积分
 
同样的，复合辛普森积分区间的划分也可以分为等长区间和不等长区间两种划分方法。
 由于不等长区间应用很少，因此直接介绍等长区间，复合辛普森1/3积分计算公式：
 

 等长区间，复合辛普森1/3积分计算误差：
 
 
由于3点积分（辛普森1/3积分）与4点积分（辛普森3/8积分）具有相同的误差阶，因此一般优先使用复合辛普森1/3积分，复合辛普森3/8积分应用极少，因此本文不做介绍。感兴趣的同学可以下去自己推一下复合辛普森3/8积分计算公式，应该也很简单。
 
（穿插）积分改进思路2
 
可以看出，将区间划分为更小的区间可以有效的提高积分计算的精度。那么我们是不是可以将积分区间不断地细分，从而来得到更精确的积分结果？理论上是可以的！！
 方法1：（区间折半法）
 a.计算2个区间的积分，用复合辛普森1/3法则（当然，也可以用梯形法则或者布尔法则等都可以，但是推荐复合辛普森1/3法则）;
 b.将两个区间等分为4个区间，用复合辛普森1/3法则计算四个子区间的积分和；
 c.比较a和b两步计算的结果相对误差是否接近一个很小的数:(Ia - Ib)/Ib<e。如果相对误差满足要求，则结束，否则继续细分区间，用复合辛普森1/3法则计算积分，直到满足要求为止。
 
方法2：改进方法1
	a.对每个子区间一分为二后，比较两个区间的积分值与一个区间的积分值是否相近;
	b.如果相近则返回该值；否则对每个子区间执行a。（递归实现）
 
上面说，理论上是可以的！！理论上是可以的！！但是实际不行！！！
 因为计算机浮点数运算是存在误差的。当区间数量很大的时候，浮点数计算的舍入误差（不理解的可以自行去学习了解一下）变得很大，会限制积分的提高。下面就是一个很好的例子。
 
 
1.3常用的高精度积分
 
上述区间二分的思想是很有用的。在龙贝格积分和自适应积分中都会找到它的影子。
 
1.3.1 龙贝格积分
 
龙贝格积分的思想，及推导过程：
 
 
 
 
 
 类推：联合两个误差为O(h^4)的积分，改进后可以得到一个误差为O(h ^ 6)的积分：
 
 由此产生龙贝格积分：
 
 
 
（穿插）积分改进思路3
 
龙贝格积分是用梯形公式推导得到的，那么我们是不是可以用辛普森公式做类似的推导，也能到达一套新的公式呢？答案是：理论上是可以的！！
 
由辛普森1/3公式推导出来的加入误差修正的积分计算式：
 
 那么问题就来了。
 我们的计算机表达的double值是有范围的，而对于辛普森外，采用递归，递归8次就是16^8 = 2 ^ 32。计算机递归计算的结果很容易就产生溢出了。除非采用多精度浮点数计算才行，而采用多精度浮点数可能存在一些计算效率的问题。。。更深入的问题我就不继续深入探讨了（有兴趣的可以试试）。。。好了，总之来说，用辛普森推导出来的方法不是很实用。还是梯形龙贝格积分经典好使。
 
1.3.2 自适应积分
 
龙贝格积分中有两点：
 a。每次都会将所有的步长折半；（那么我们能不能只对计算误差较大的那个区间进行折半成两个子区间再积分呢？）
 b。辛普森公式推导的方法，由于计算机double表示范围有限没用上。
 
下面都会在自适应积分中派上用场。自适应积分思想：（递归方法）
 step1.辛普森1/3法则计算区间[a,b]的积分；
 step2.区间步长折半，复合辛普森1/3法则计算区间[a,b]的积分；
 如果两次误差小于给定的误差限，停止执行，返回计算结果
 ；
 如果两次误差大于给定的误差限，对两个子区间分别执行step1和step2。
 
伪代码：
 
 
1.3.3 高斯积分
 
高斯积分应用比较广泛的是高斯-勒让德公式，此外还有高斯-切比雪夫公式；以及区间【0，正无穷】的广义积分 高斯-拉盖尔公式；以及区间【负无穷，正无穷】的广义积分 高斯-埃尔米特公式。
 
高斯-勒让德公式
 
介绍高斯积分的文档就很多了，本文不细讲推导的过程，给出计算方法：
 
 多点高斯-勒让德公式：
 
 
 高斯勒让德公式的区间为【-1，1】.对于任意区间为【a，b】；需要做一下简单的变量代换，
 
 
高斯积分误差：
 
高斯勒让德公式的区间为【-1，1】，误差：
 
 对于任意区间为【a，b】，误差：
 
 
 由上式继续往下推导，可以得出，当区间长度【b-a】> 2(n+1)时，误差会发散。因此对于高斯点为n的高斯积分，积分区间长度不要大于2(n+1)，否则积分结果很可能会不准确。
 此外，高斯积分的误差还与被积函数的光顺性有关，被积函数光顺性越差，误差越大。详细可参考《现代数值积分》第5章 数值积分与数值微分，P139页）。
 
高斯积分的优点：高斯积分时给定节点数下代数精度最高的求积公式。
 
高斯积分的不足：高斯积分每次改变积分点的个数，所以的积分点的函数值都需要重新计算。
 
高斯-切比雪夫公式
 
高斯点为切比雪夫节点。
 
（未完待续。。。）
 
广义积分：高斯-拉盖尔公式
 
广义积分：高斯-埃尔米特公式
 
2 多重积分方法
 
二重积分的辛普森公式
 
二重积分的蒙特卡洛方法
 
多重（高维）积分的蒙特卡洛方法

梯形多步法和辛普森积分
 具体算法以及相关定理参考以下链接：
 
 
http://wenku.baidu.com/link?url=ZAScYogajXHrTTRa5xjpUPtS7OQQXZ_LfXaWNkczTtWf2MJgx0RZFUuca4iRGUcPHAVXOyGOIydwz3j-rriBiyLxRMnio7m4fNFcJOKJqui
 
 
#include <stdio.h>
 
#include <string.h>
 #include <stdlib.h>
 #include <math.h>
 
 #define EPS 1.0E-14        //计算精度
 #define DIV 1000        //分割区间，值越大，精确值越高
 #define ITERATE 20        //二分区间迭代次数，区间分割为2^n，初始化应该小一点，否则会溢出
 
 #define RECTANGLE 0        //矩形近似
 #define TRAPEZOID 1        //梯形近似
 #define TRAPEZOID_FORMULA 2    //递推梯形公式
 #define SIMPSON_FORMULA 3     //辛普森公式
 #define BOOL_FORMULA 4        //布尔公式 
 
 double function1(double);       
 double function2(double);
 double function3(double); 
 void Integral(int, double f(double), double, double);      //矩形, 梯形逼近求定积分公式 
 void Trapezoid_Formula(double f(double), double a, double b);    //递推梯形公式
 void Simpson_Formula(double f(double), double a, double b);        //递推辛普森公式
 void Bool_Formula(double f(double), double a, double b);        //递推布尔公式
 void Formula(int, double f(double), double, double);
 
 double function1(double x)
 {
     double y;
     y = cos(x);
     return y;
 }
 
 double function2(double x)
 {
     double y;
     y=1/x;
    // y = 2+sin(2 * sqrt(x));
     return y;
 }
 
 double function3(double x)
 {
     double y;
     y = exp(x);
     return y;
 }
 
 int main()
 {    
     printf("y=cos(x) , x=[0, 1]\n");
     printf("Rectangle : "); Integral(RECTANGLE, function1, 0, 1);
     printf("Trapezoid : "); Integral(TRAPEZOID, function1, 0, 1);
     Formula(TRAPEZOID_FORMULA, function1, 0, 1);
     Formula(SIMPSON_FORMULA, function1, 0, 1);
     Formula(BOOL_FORMULA, function1, 0, 1);
     printf("\n"); 
     
     printf("y=1/x , x=[1, 5]\n");
     printf("Rectangle : "); Integral(RECTANGLE, function2, 1, 5);
     printf("Trapezoid : "); Integral(TRAPEZOID, function2, 1, 5);
     Formula(TRAPEZOID_FORMULA, function2, 1, 5);
     Formula(SIMPSON_FORMULA, function2, 1, 5);
     Formula(BOOL_FORMULA, function2, 1, 5); 
     printf("\n");
     
     printf("y=exp(x) , x=[0, 1]\n");
     printf("Rectangle : "); Integral(RECTANGLE, function3, 0, 1);
     printf("Trapezoid : "); Integral(TRAPEZOID, function3, 0, 1);
     Formula(TRAPEZOID_FORMULA, function3, 0, 1);
     Formula(SIMPSON_FORMULA, function3, 0, 1);
     Formula(BOOL_FORMULA, function3, 0, 1);
     
     return 0;
 }
 //积分：分割-近似-求和-取极限
 //利用黎曼和求解定积分
 void Integral(int method, double f(double),double a,double b)   //f(double)，f(x), double a，float b，区间两点
 {
     int i;
     double x;            
     double sum = 0;      //矩形总面积 
     double h;            //矩形宽度
     double t = 0;        //单个矩形面积
     
     h = (b-a) / DIV;
      
     for(i=1; i <= DIV; i++)       
     {
         x = a + h * i;     
         switch(method)
         {
             case RECTANGLE :
                     t = f(x) * h;  
                 break;
             case TRAPEZOID :
                     t = (f(x) + f(x - h)) * h / 2;
                 break;
             default:
                 break;
         }
         sum = sum + t;   //各个矩形面积相加
     }
     
     printf("the value of function is %.10f\n", sum);
 }
 
 //递推梯形公式
 void Trapezoid_Formula(double f(double), double a, double b)//递推梯形公式
 {
     unsigned int i, j, M, J=1; 
     double h;
     double F_2k_1 = 0;
     double T[32];
     double res = 0;
     
     T[0] = (f(a) + f(b)) * (b-a)/2;
     for(i=0; i<ITERATE; i++)
     {
         F_2k_1 = 0;
         J = 1;
          
         for(j=0; j<=i; j++)    //区间划分
             J <<=  1;    //2^J
         h = (b - a) /J;    //步长
         //M = J/2;            //2M表示将积分区域划分的块数
         for(j=1; j <= J; j += 2)    //
         {
             F_2k_1 += f(a + j*h);    //f(2k-1)累加和
         }     
         T[i+1] = T[i]/2 + h * F_2k_1;        //递推公式
         res =  T[i+1]; 
         
         if(fabs(T[i+1] - T[i]) < EPS) 
             if(i < 3) continue;
             else break;
     }
     
     printf("Trapezoid_Formula : the value of function is %.10f\n", res);
     //return res;
 }
 //辛普森公式
 void Simpson_Formula(double f(double), double a, double b)        //辛普森公式
 {
     unsigned int i, j, M, J=1; 
     double h;
     double F_2k_1 = 0;
     double T[32], S[32];
     double res_T = 0, res_S = 0;
     
     T[0] = (f(a) + f(b)) * (b-a)/2;
     for(i=0; i<ITERATE; i++)
     {
         F_2k_1 = 0;
         J = 1;
          
         for(j=0; j<=i; j++)    //区间划分
             J <<=  1;    //2^J
         h = (b - a) /J;    //步长
         //M = J/2;            //2M表示将积分区域划分的块数
         for(j=1; j <= J; j += 2)    //
         {
             F_2k_1 += f(a + j*h);    //f(2k-1)累加和
         }     
         T[i+1] = T[i]/2 + h * F_2k_1;        //递推梯形公式
         res_T =  T[i+1]; 
         
         if(fabs(T[i+1] - T[i]) < EPS) 
             if(i < 3) continue;
             else break;
     }
     
     for(i=1; i<ITERATE; i++)
     {
         S[i] = (4 * T[i] - T[i-1]) / 3;        //递推辛普森公式 
         res_S = S[i];
         if(fabs(S[i] - S[i-1]) < EPS) 
             if(i < 3) continue;
             else break;
     }
     
     printf("Simpson_Formula : the value of function is %.10f\n", res_S);
     //return res_S;
 }
 
 //布尔公式 
 void Bool_Formula(double f(double), double a, double b)        //布尔公式 
 {
     unsigned int i, j, M, J=1; 
     double h;
     double F_2k_1 = 0;
     double T[32], S[32], B[32];
     double res_T = 0, res_S = 0, res_B;
     
     T[0] = (f(a) + f(b)) * (b-a)/2;
     for(i=0; i<ITERATE; i++)
     {
         F_2k_1 = 0;
         J = 1;
          
         for(j=0; j<=i; j++)    //区间划分
             J <<=  1;    //2^J
         h = (b - a) /J;    //步长
         //M = J/2;            //2M表示将积分区域划分的块数
         for(j=1; j <= J; j += 2)    //
         {
             F_2k_1 += f(a + j*h);    //f(2k-1)累加和
         }     
         T[i+1] = T[i]/2 + h * F_2k_1;        //递推公式
         res_T =  T[i+1]; 
         
         if(fabs(T[i+1] - T[i]) < EPS) 
             if(i < 3) continue;
             else break;
     }
     
     for(i=1; i<ITERATE; i++)
     {
         S[i] = (4 * T[i] - T[i-1]) / 3;                //递推辛普森公式 
         res_S = S[i];
         if(fabs(S[i] - S[i-1]) < EPS) 
             if(i < 3) continue;
             else break;
     }
     
     for(i=1; i<ITERATE; i++)
     {
         B[i] = (16 * S[i] - S[i-1]) / 15;            //递推布尔公式
         res_B = B[i];
         if(fabs(B[i] - B[i-1]) < EPS) 
             if(i < 3) continue;
             else break;
     }
     
     printf("Bool_Formula : the value of function is %.10f\n", res_B);
     //return res_B;
 } 
 
 //采用二分区间的方法迭代，实际运行速度很慢
 void Formula(int method, double f(double), double a, double b)//递推梯形公式
 {
     unsigned int i, j, M, J=1; 
     double h;
     double F_2k_1 = 0;
     double T[32], S[32], B[32];
     double res_B = 0, res_S = 0, res_T = 0;
     
     T[0] = (f(a) + f(b)) * (b-a)/2;
     for(i=0; i<ITERATE; i++)
     {
         F_2k_1 = 0;
         J = 1;
          
         for(j=0; j<=i; j++)    //区间划分
             J <<=  1;    //2^J
         h = (b - a) /J;    //步长
         //M = J/2;            //2M表示将积分区域划分的块数
         for(j=1; j <= J; j += 2)    //
         {
             F_2k_1 += f(a + j*h);    //f(2k-1)累加和
         }     
         T[i+1] = T[i]/2 + h * F_2k_1;        //递推公式
         res_T =  T[i+1]; 
         
         if(fabs(T[i+1] - T[i]) < EPS) 
             if(i < 3) continue;
             else break;
     }
     
     switch(method)
     {
         default : 
         case TRAPEZOID_FORMULA :
                 printf("Trapezoid_Formula : the value of function is %.10f\n", res_T);
                 //return res_T;
                 break;
         case SIMPSON_FORMULA :
                 for(i=1; i<ITERATE; i++)
                     {
                         S[i] = (4 * T[i] - T[i-1]) / 3; 
                         res_S = S[i];
                         if(fabs(S[i] - S[i-1]) < EPS) 
                             if(i < 3) continue;
                             else break;
                     }
                 printf("Simpson_Formula : the value of function is %.10f\n", res_S);
                 //return res_S;
                 break;
         case BOOL_FORMULA : 
                 for(i=1; i<ITERATE; i++)
                 {
                     S[i] = (4 * T[i] - T[i-1]) / 3; 
                     res_S = S[i];
                     if(fabs(S[i] - S[i-1]) < EPS) 
                         if(i < 3) continue;
                         else break;
                 }
                 for(i=1; i<ITERATE; i++)
                 {
                     B[i] = (16 * S[i] - S[i-1]) / 15;
                     res_B = B[i];
                     if(fabs(B[i] - B[i-1]) < EPS) 
                         if(i < 3) continue;
                         else break;
                 }
                 
                 printf("Bool_Formula : the value of function is %.10f\n", res_B);
                 //return res_B;
             break;
     }
     
 }