数值积分方法

落日映苍穹つ 2022-08-01 03:15 1082阅读 0赞

数值积分是工程师和科学家经常使用的基本工具，用来计算无法解析求解的定积分的近似解。 
如：Φ(x)=∫x0t3et−1dt不存在Φ(x)的解析解，要求Φ(5)。 
那么我们就要通过数值积分的方法来计算，数值积分的目的是，通过在有限个采样点上计算f(x)的值来逼近 f(x)在区间\[a,b\]上的定积分. 
设a=x0<x1<…<xM=b. 称形如

![这里写图片描述][20150619101924570]  
且具有性质  ∫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时，称为闭型牛顿－科特斯公式。

![这里写图片描述][20150619105317172]

### 左/中/右矩形公式、梯形公式 ###

左/中/右矩形公式

![这里写图片描述][20150619101433460]  
梯形公式  
![这里写图片描述][20150619101549260]  
图形如下  
![这里写图片描述][20150619101711326]

--------------------

### 辛普森公式 ###

#### 推导过程 ####

若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\] 上的近似积分。

![这里写图片描述][20150619103203729]  
由假设我们有  
![这里写图片描述][20150619103119394]  
![这里写图片描述][20150619103330564]  
所以可得

∫baf(x)dx=∫x2x0\+∫x4x2\+⋯\+∫xnxn−2f(x)dx

![这里写图片描述][20150619103717965]

#### 误差估计 ####

若令Sn=Δx3\[f(x0)\+4f(x1)\+2f(x2)\+4f(x3)⋯\+2f(xn−2)\+4f(xn−1)\+f(xn)\]且f(4)(x)在\[a,b\]上连续，则我们可以估计出辛普森公式的误差值为

![这里写图片描述][20150619104727878]

#### 例题1 ####

试用辛普森公式估计∫10e−x2dx，取n=6  
解：  
令f(x)=e−x2,Δx=16则

![这里写图片描述][20150619104932965]

--------------------

# 拉格朗日插值 #

### 分段线性插值 ###

例如：函数f(x)=11\+x2如果在区间\[−5，5\]上取11个等距节点:xk=−5\+k(k=0,1,2,...,10),由Lagrange插值公式可得到f(x)的10次L10(x)。如图所示：

![这里写图片描述][5_302.gif]  
 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 .

![这里写图片描述][5_309.gif]

#### 误差估计 ####

根据拉格朗日一次插值函数的余项,可以得到分段线性插值函数的插值误差估计: 对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)近似值。

![这里写图片描述][5_312.gif]  
解:  
在每个分段区间  \[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) 它们的函数值满足下表：

![这里写图片描述][5_108.gif]  
因为  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 ####

已知：

利用此三值的二次插值多项式求lg(12)的近似值。  
解：设x0=10,x1=15,x2=20,则:

![这里写图片描述][20150619121019048]  
故:  
  
![这里写图片描述][20150619121107143]  
![这里写图片描述][20150619121115919]  
  
所以  
  
![这里写图片描述][20150619121127683]

#### 误差估计 ####

我们在\[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  
估计误差：

![这里写图片描述][5_135.gif]

--------------------

# 三次样条插值法 #

### 三次样条曲线原理 ###

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. 左右两端点处特性（自然边界，固定边界，非节点边界）
> 
> 根据定点，求出每段样条曲线方程中的系数，即可得到每段曲线的具体表达式。
> 
> 插值和连续性：  
> ![这里写图片描述][26192553-43497c728faa4a2e9d3d78c451133034.png], 其中 i=0,1,…,n−1  
> 微分连续性：  
> ![这里写图片描述][26192555-7dbb8641155e423b92302e508975ea16.png], 其中i=0,1,…,n−2  
> 样条曲线的微分式：  
> ![image][]  
> 将步长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  
> 则要求解的方程组可写为：  
> 
> 
> ![这里写图片描述][26192741-c629f65d48af4a07912bb6b2fbdad544.png]  
> ![这里写图片描述][26192742-67d70a20fe4c433181dcfebe39a06a93.png]  
> 
> 
> ##### b. 固定边界(Clamped) #####
> 
> 首尾两端点的微分值是被指定的，这里分别定为A和B。则可以推出  
> 
> 
>   
> ![这里写图片描述][26192746-828dbca3271941088a9df631cfbfdbc9.png]  
> ![这里写图片描述][26192749-3d8d4b6bb3e7460181bcfe7cd61525c9.png]  
>   
> 将上述两个公式带入方程组，新的方程组左侧为  
> ![这里写图片描述][26192806-f0180ec05e234bdbbfbb1ed5032e02a6.png]
> 
> ##### 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)
> 
>   
> 新的方程组系数矩阵可写为：  
> 
> 
> ![这里写图片描述][26192903-d7db9e3fffb74f4b86afa926827e549b.png]
> 
> 右下图可以看出不同的端点边界对样条曲线的影响：  
> 
> 
>   
> ![这里写图片描述][26193007-3c4400f6621b4f98bf26a085e36d3130.png]
> 
> #### 算法总结 ####
> 
> 假定有n+1个数据节点  
> 
> 
> (x0,y0),(x1,y1),(x2,y2),…,(xn,yn),
> 
>   
> a. 计算步长  hi=xi\+1−xi(i=0,1,…,n−1)  
> b. 将数据节点和指定的首位端点条件带入矩阵方程  
> c. 解矩阵方程，求得二次微分值  mi 。  
> d. 计算样条曲线的系数：  
> 
> 
> ![这里写图片描述][26192936-2b0349526ad54de296f83cf1001938ef.png]  
> 其中  i=0,1,…,n−1  
> e. 在每个子区间  xi≤x≤xi\+1 中，创建方程  
> 
> 
> gi(x)=ai\+bi(x−xi)\+ci(x−xi)2\+di(x−xi)3

[20150619101924570]: /images/20220731/c8fc89a49726443f8e04025efaa761d6.png
[20150619105317172]: /images/20220731/70c8029b7b7c40eabf967b053b83d24b.png
[20150619101433460]: /images/20220731/6d5ab28f60a049e7a18425acc7a5282f.png
[20150619101549260]: /images/20220731/1d73ac41af454aa4b060c0a5b2845d3c.png
[20150619101711326]: /images/20220731/51f7f5a7c192424b88062df95a9670b7.png
[20150619103203729]: /images/20220731/9126e88a702a4660b78289bb84669e69.png
[20150619103119394]: /images/20220731/f3ec4413e82c411aac4b8400e2249d71.png
[20150619103330564]: /images/20220731/14f5a9919bc840d88dfae0b1f660581a.png
[20150619103717965]: /images/20220731/fbcc02b28b914173a627b9f2f495fdf0.png
[20150619104727878]: /images/20220731/361bbb8e300b46dda76a7ca590faef85.png
[20150619104932965]: /images/20220731/97d35c80ab1c4dfb80d51d10c380c582.png
[5_302.gif]: /images/20220731/c13fe73141494a6188b10cc6e65d73ce.png
[5_309.gif]: /images/20220731/00e999201d2a42e4a73c34f8509800b0.png
[5_312.gif]: /images/20220731/8d33701104964d05998fe41574076ebf.png
[5_108.gif]: /images/20220731/5a15903b629a40eaa2257d48b8edb28e.png
[20150619121019048]: /images/20220731/082095fbae3b434fab3df054145ac708.png
[20150619121107143]: /images/20220731/c2bdf53ed8a5474a822607ea9e7e0eb5.png
[20150619121115919]: /images/20220731/dfb1c8d9667a4cd8819f845eb27246e8.png
[20150619121127683]: /images/20220731/10b9d509ce8a469bbcf8b2796f86af96.png
[5_135.gif]: /images/20220731/2c5cae78232c42f9adece5fa76cedb8e.png
[26192553-43497c728faa4a2e9d3d78c451133034.png]: /images/20220731/c43ea8a92bfb4068858e0667bbbfe5dc.png
[26192555-7dbb8641155e423b92302e508975ea16.png]: /images/20220731/9b3b79b72d6d48ea967453ab28298934.png
[image]: /images/20220731/1685a28df71f421180e4e83c72f7a66a.png
[26192741-c629f65d48af4a07912bb6b2fbdad544.png]: /images/20220731/18e465c54f194aa59233f10452552a88.png
[26192742-67d70a20fe4c433181dcfebe39a06a93.png]: /images/20220731/631af2e9204744c896b0d5d91bf21e70.png
[26192746-828dbca3271941088a9df631cfbfdbc9.png]: /images/20220731/5887517e3c8c4dc69dd5613a52d48baa.png
[26192749-3d8d4b6bb3e7460181bcfe7cd61525c9.png]: /images/20220731/78ddca08bfa84b22a3b75ed49240265f.png
[26192806-f0180ec05e234bdbbfbb1ed5032e02a6.png]: /images/20220731/d1465ad38be64639b1ad09871bba55cd.png
[26192903-d7db9e3fffb74f4b86afa926827e549b.png]: /images/20220731/787f78c50d37403490e030ffd25e835f.png
[26193007-3c4400f6621b4f98bf26a085e36d3130.png]: /images/20220731/ffe9f8f8b73643a99b9b922bb47937d3.png
[26192936-2b0349526ad54de296f83cf1001938ef.png]: /images/20220731/0446f3db6cfe4e17a9e98a77e1e2aa4e.png