数值积分方法
最后更新于：2022-04-01 16:06:59
数值积分是工程师和科学家经常使用的基本工具，用来计算无法解析求解的定积分的近似解。
如：Φ(x)=∫x0t3et−1dt不存在Φ(x)的解析解，要求Φ(5)。
那么我们就要通过数值积分的方法来计算，数值积分的目的是，通过在有限个采样点上计算f(x)的值来逼近 f(x)在区间[a,b]上的定积分.
设a=x0<x1<…<xM=b. 称形如
且具有性质∫baf(x)dx=Q[f]+E[f]的公式为数值积分或面积公式。项E[f]称为积分的截断误差，值{xk}Mk=0称为面积节点， {wk}Mk=0称为权。
下面介绍几种常用的数值积分方法。
### 基于多项式差值的面积公式
通过M+1个等距点{(xk,f(xk))}Mk=0存在唯一的次数小于等于M的多项式PM(x)。当用该多项式来近似[a,b]上的f(x)时，PM(x)的积分就近似等于f(x)的积分，这类公式称为牛顿－科特斯公式。当使用采样点x0=a和xM=b时，称为闭型牛顿－科特斯公式。
### 左/中/右矩形公式、梯形公式
左/中/右矩形公式
梯形公式
图形如下
### 辛普森公式
#### 推导过程
若f(x)在[a,b]上有定义，将区间[a,b]分割成n等分（取n为偶数），则有a=x0<x1<…<xM=b，其中
xi=a+iΔx（∀i=0,1,…,n,Δx=b−an为步长）
这里我们想通过(x0,f(x0)),(x1,f(x1)),(x2,f(x2))三点抛物线g(x)=αx2+βx+γ来取代f(x)在[x0,x2]的定义，今儿求出它的近似积分值(如图)，最后用累加的方式求得f(x)在[a,b]上的近似积分。
由假设我们有
所以可得
∫baf(x)dx=∫x2x0+∫x4x2+⋯+∫xnxn−2f(x)dx
#### 误差估计
若令Sn=Δx3[f(x0)+4f(x1)+2f(x2)+4f(x3)⋯+2f(xn−2)+4f(xn−1)+f(xn)]且f(4)(x)在[a,b]上连续，则我们可以估计出辛普森公式的误差值为
#### 例题1
试用辛普森公式估计∫10e−x2dx，取n=6
解：
令f(x)=e−x2,Δx=16则
## 拉格朗日插值
### 分段线性插值
例如：函数f(x)=11+x2如果在区间[−5，5]上取11个等距节点:xk=−5+k(k=0,1,2,...,10),由Lagrange插值公式可得到f(x)的10次L10(x)。如图所示：
L10(x)仅在区间的中部能较好的逼近函数f(x),在其它部位差异较大，而且越接近端点，逼近效果越差。
可以证明，当节点无限加密时，Ln(x)也只能在很小的范围内收敛，这一现象称为Runge现象。它表明通过增加节点来提高逼近程度是不适宜的，因而不采用高次多项式插值。
#### 推导过程
已知函数f(x)在区间[xk,xk+1]的端点上的函数值yk=f(xk),yk+1=f(xk+1),求一个一次函数y=P1(x)使得yk=f(xk),yk+1=f(xk+1), 其几何意义是已知平面上两点(xk,yk),(xk+1,yk+1),求一条直线过该已知两点。
由直线的点斜式公式可知：
P1(x)=yk+yk+1−ykxk+1−xk(x−xk)
把此式按照yk和yk+1写成两项:
P1(x)=x−xk+1xk−xk+1yk+x−xkxk+1−xkyk+1
记
lk=x−xk+1xk−xk+1 lk+1=x−xkxk+1−xk
并称它们为一次插值基函数。该基函数的特点如下：
li(xk)={10k=ik!=i
对于i=0，
l0(x)=x−x1x0−x1,x∈[x0,x1]
其它点上，l0(x)=0;
对于i=1,2,…,n−1，
li(x)={x−xi−1xi−xi−1x−xi+1xi−xi+1x∈[xi−1,xi]x∈[xi,xi+1]
其它点上，li(x)=0; 对于i=n，
ln(x)=x−xn−1xn−xn−1,x∈[xn−1,xn]
其它点上，ln(x)=0. 于是，
P1(x)=∑k=0nyklk(x)
此表达式与前面的表达式是相同的,这是因为在区间[xk,xk+1]上, 只有lk(x),lk+1(x)是非零的，其它基函数均为零。
从而
P1(x)=yklk(x)+yk+1lk+1(x)
此形式称之为拉格朗日型插值多项式。其中, 插值基函数与yk、yk+1 无关，而由插值结点xk、xk+1 所决定。一次插值多项式是插值基函数的线性组合, 相应的组合系数是该点的函数值yk、yk+1 .
#### 误差估计
根据拉格朗日一次插值函数的余项,可以得到分段线性插值函数的插值误差估计: 对x∈[a,b],当x∈[xk,xk+1]时,
R(x)=12f′′(ξ)(x−xk)(x−xk+1)
则R(x)≤h28M，其中
h=max0≤k≤n−1|xk+1−xk|,M=maxx∈(a,b)f′′(x)
于是有下面的定理：
如果f(x)在[a,b]上二阶连续可微，则分段连续函数φ(x)的余项有以下误差估计：
|R(x)|=|f(x)−ϕ(x)|≤h28M
其中
h=max0≤k≤n−1|xk+1−xk|,M=maxx∈(a,b)f′′(x)
于是可以加密插值结点, 缩小插值区间, 使h减小, 从而减小插值误差。
#### 例题1
已知函数y=f(x)=11+x2在区间[0,5]上取等距插值节点(如下表)，
求区间上分段线性插值函数,并利用它求出f(4.5)近似值。
解:
在每个分段区间[k,k+1]上，
P1(x)=x−(k+1)k−(k+1)yk+x−(k)k+1−(k)yk+1=−yk(x−k−1)+yk+1(x−k)
于是
P1(4.5)=−0.05882×(4.5−5)+0.03846×(4.5−4)=0.04864
### 拉格朗日型二次插值多项式
#### 推导过程
已知函数y=f(x)在点xk−1,xk,xk+1上的函数值yk−1=f(xk−1),yk=f(xk),yk+1=f(xk+1), 求一个次数不超过二次的多项式P2(x), 使其满足,
P2(xk−1)=yk−1,P2(xk)=yk,P2(xk+1)=yk+1
其几何意义为:
已知平面上的三个点
(xk−1,yk−1),(xk,yk),(xk+1,yk+1),
求一个二次抛物线, 使得该抛物线经过这三点。
#### 插值基本多项式
有三个插值结点xk−1,xk,xk+1构造三个插值基本多项式，要求满足：
(1) 基本多项式为二次多项式；
(2) 它们的函数值满足下表：
因为lk−1(xk)=0,lk−1(xk+1)=0, 故有因子(x−xk)(x−xk+1), 而其已经是一个二次多项式, 仅相差一个常数倍, 可设
lk−1(x)=a(x−xk)(x−xk+1)
又因为
lk−1(xk−1)=1⇒a(xk−1−xk)(xk−1−xk+1)=1
得
a=1(xk−1−xk)(xk−1−xk+1)
从而
lk−1(x)=(x−xk)(x−xk+1)(xk−1−xk)(xk−1−xk+1)
同理得
lk(x)=(x−xk−1)(x−xk+1)(xk−xk−1)(xk−xk+1)
lk+1(x)=(x−xk−1)(x−xk)(xk+1−xk−1)(xk+1−xk)
#### 拉格朗日型二次插值多项式
由前述, 拉格朗日型二次插值多项式:
P2(x)=yk−1lk−1(x)+yklk(x)+yk+1lk+1(x)
P2(x)是三个二次插值多项式的线性组合,因而其是次数不超过二次的多项式,且满足:
P2(xi)=yi,(i=k−1,k,k+1)
#### 例题2
已知：
| xi | 10 | 15 | 20 |
|-----|-----|-----|-----|
| yi=lgxi | 1 | 1.1761 | 1.3010 |
利用此三值的二次插值多项式求lg(12)的近似值。
解：设x0=10,x1=15,x2=20,则:
故:
所以
#### 误差估计
我们在[a,b]上用多项式Pn(x) 来近似代替函数f(x), 其截断误差记作Rn(x)=f(x)−Pn(x)
当x在插值结点xi上时Rn(xi)=f(xi)−Pn(xi)=0,下面来估计截断误差：
定理1：
设函数y=f(x)的n阶导数y(n)=f(n)(x)在[a,b]上连续，
y(n+1)=f(n+1)(x)在(a,b)上存在；
插值结点为:a≤x0<x1<…<xn≤b,
Pn(x)是n次拉格朗日插值多项式；
则对任意x∈[a,b]有：
Rn(x)=1(n+1)!f(n+1)(ξ)ωn+1(x)
其中ξ∈(a,b), ξ依赖于x：ωn+1(x)=(x−x0)(x−x1)…(x−xn)
设maxa≤x≤b|fn+1(x)|≤Mn+1，则
|Rn(x)|≤1(n+1)!Mn+1ωn+1(x)
易知，线性插值的截断误差为:
R1(x)=12f′′(ξ)(x−x0)(x−x1)
二次插值的截断误差为:
R2(x)=16f(3)(ξ)(x−x0)(x−x1)(x−x2)
#### 例题3
在例2中，用lg10,lg15和lg20计算lg12.
P2(12)=1.0766,
e=|1.0792−1.0766|=0.0026
估计误差：
## 三次样条插值法
### 三次样条曲线原理
x:a=x0<x1<…<xM=by:y0 y1 … yn
#### 定义
样条曲线S(x)是一个分段定义的公式。给定n+1个数据点，共有n个区间，三次样条方程满足以下条件：
a. 在每个分段区间[xi,xi+1]（i=0,1,…,n−1，x递增），S(x)=Si(x) 都是一个三次多项式。
b. 满足S(xi)=yi（i=0,1,…,n）
c.S(x)，导数S′(x)，二阶导数S′′(x)在[a,b]区间都是连续的，即S(x)曲线是光滑的
所以n个三次多项式分段可以写作：
Si(x)=ai+bi(x−xi)+ci(x−xi)2+di(x−xi)3,i=0,1,…,n−1
其中ai,bi,ci,di代表4n个未知系数。
#### 求解
已知:
a. n+1个数据点[xi,yi],i=0,1,…,n
b. 每一分段都是三次多项式函数曲线进行逼近
c. 节点达到二阶连续
d. 左右两端点处特性（自然边界，固定边界，非节点边界）
根据定点，求出每段样条曲线方程中的系数，即可得到每段曲线的具体表达式。
插值和连续性：
, 其中 i=0,1,…,n−1
微分连续性：
, 其中i=0,1,…,n−2
样条曲线的微分式：
将步长hi=xi+1−xi带入样条曲线的条件：
a. 由Si(xi)=yi(i=0,1,…,n−1)推出
ai=yi
b. 由Si(xi+1)=yi+1(i=0,1,…,n−1)推出
ai+hibi+h2ici+h3idi=yi+1
c. 由S′i(xi+1)=S′i+1(xi+1)(i=0,1,…,n−2)推出
S′i(xi+1)=bi+2ci(xi+1−xi)+3di(xi+1−xi)2=bi+2cih+3dih2S′i+1(xi+1)=bi+1+2ci(xi+1−xi+1)+3di(xi+1−xi+1)2=bi+1
由此可得：
bi+2hici+3h2idi−bi+1=0
d. 由S′′i(xi+1)=S′′i+1(xi+1)(i=0,1,…,n−2)推出
2ci+6hidi−2ci+1=0
设mi=S′′i(xi)=2ci，则
a. 2ci+6hidi−ci+1=0可写为：
mi+6hidi−mi+1=0，推出
di=mi+1−mi6hi
b. 将ci,di带入 yi+hibi+h2ici+h3idi=yi+1可得：
bi=yi+1−yihi−hi2mi−hi6(mi+1−mi)
c. 将bi,ci,di带入bi+2hici+3hidi=bi+1(i=0,1,…,n−2)可得：
himi+2(hi+hi+1)mi+1+hi+1mi+2=6[yi+2−yi+1hi+1−yi+1−yihi]
#### 端点条件
由i的取值范围可知，共有n-1个公式， 但却有n+1个未知量m 。要想求解该方程组，还需另外两个式子。所以需要对两端点x0和xn的微分加些限制。 选择不是唯一的，3种比较常用的限制如下。
##### a. 自由边界(Natural)
首尾两端没有受到任何让它们弯曲的力，即S′′=0。具体表示为m0=0和 mn=0
则要求解的方程组可写为：
##### b. 固定边界(Clamped)
首尾两端点的微分值是被指定的，这里分别定为A和B。则可以推出
将上述两个公式带入方程组，新的方程组左侧为
##### c. 非节点边界(Not-A-Knot)
指定样条曲线的三次微分匹配，即
S′′′o(x1)=S′′′1(x1)S′′′n−2(xn−1)=S′′′n−1(xn−1)
根据S′′′i(x)=6di和di=mi+1−mi6hi，则上述条件变为
h1(m1−m0)=h0(m2−m1)hn−1(mn−1−mn−2)=hn−2(mn−mn−1)
新的方程组系数矩阵可写为：
右下图可以看出不同的端点边界对样条曲线的影响：
#### 算法总结
假定有n+1个数据节点
(x0,y0),(x1,y1),(x2,y2),…,(xn,yn),
a. 计算步长hi=xi+1−xi(i=0,1,…,n−1)
b. 将数据节点和指定的首位端点条件带入矩阵方程
c. 解矩阵方程，求得二次微分值mi。
d. 计算样条曲线的系数：
其中i=0,1,…,n−1
e. 在每个子区间xi≤x≤xi+1中，创建方程
gi(x)=ai+bi(x−xi)+ci(x−xi)2+di(x−xi)3