數值積分是工程師和科學家經常使用的基本工具，用來計算無法解析求解的定積分的近似解。
如：Φ(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
