【插值】插值方法原理详解
插值问题详解
注明出处:http://www.cnblogs.com/duye/p/8671820.html
1.
我在具体的应用(如数学建模竞赛)中,常常需要根据已知的函数点进行数据、模型的处理和分析,而通常情况下现有的数据是极少的,不足以支撑分析的进行,这时就需要使用一些数学的方法,“模拟产生”一些新的但又比较靠谱的值来满足需求。一般来说,我可以去调用MATLAB或者Python的一些库函数来实现,这个功能就是“插值”。然而这有个非常让我苦恼的问题,我可以从手册上知道这个函数实现“三次多项式插值”,那个函数实现“样条插值”.......但究竟在什么情况下使用何种插值方法呢?若不对插值方法做深入的学习,这个疑团恐难以解开。
于是,在这个原因驱动之下,我决定对常见、常用的插值方法比较深入的学习一下。我希望读者也是基于这个原因来读这篇文章,希望我的总结能对你有所帮助。
2. 插值
简单讲,插值就是根据已知数据点(条件),来预测未知数据点值得方法。具体来说,假如你有n个已知条件,就可以求一个n-1次的插值函数P(x),使得P(x)接近未知原函数f(x),并由插值函数预测出你需要的未知点值。而又n个条件求n-1次P(x)的过程,实际上就是求n元一次线性方程组。
代数插值
代数插值就是多项式插值,即所求插值函数为多项式函数:
显然,系数a0.....an为所求。如果已知n+1个条件,需要n+1个方程组如下:
这时,便可以用待定系数求解。
一、泰勒插值
首先需要回顾泰勒多项式:
余项:
满足n阶可导这个条件实在是太苛刻,导致实际上泰勒插值并不常用,下面介绍拉格朗日插值与牛顿插值,这两种方法在本质上是相同的。
二、拉格朗日插值
上面引论中提到,一般来说多项式插值就是求n-1个线性方程的解,拉格朗日插值即是基于此思想。拉格朗日创造性的避开的方程组求解的复杂性,引入“基函数”这一概念,使得快速手工求解成为可能。
DEF:求作<=n 次多项式 pn(x),使满足条件pn(xi)= yi,i = 0,1,…,n.这就是所谓拉格朗日( Lagrange)插值
先以一次(线性)为例,介绍基函数方法求解,再推广到任意次多项式:
已知x0,x1;y0,y1,求P(x)= a0 + a1x,使得P(x)过这两点。则显然:
下面介绍基函数的一般形式:
对于要求的插值函数P(x),可以证明,均可以化简为以下形式:
要使得:
则要求:
从而基函数表为:
从而,给出基函数的一般形式为:
简而记之:
此即著名的拉格朗日插值公式。
余项:
可以由罗尔定理(高数讲过的定理,忘记请自行google)等证明,拉格朗日插值的余项为:
值得注意的是,拉格朗日插值的方法,在插值区间内插值的精度远远大于区间外的精度,故一般说,区间外,拉格朗日插值是不准确的。
三、牛顿插值
牛顿插值本质上和朗格朗日插值无异,但为什么牛顿也要提出这么一种插值方法呢?这是因为,拉格朗日插值每增加一个新节点,都要重新计算,换言之,它不具备承袭性。牛顿经过严密的推导,总结了下列具有承袭性的插值方法。为理解牛顿插值公式,需先了解下面的概念:
1) 差商是什么
DEF:设有函数f (x)以及自变量的一系列互不相等的x0, x1,…, xn(即在i ¹ j时,x i¹ xj)的值 f(xi),称
为f (x)在点xi , xi处的一阶差商,并记作f [xi , xj]。又称
为f (x)在点xi, xj, xk处的二阶差商;称:
为f (x)在点x0, x1,…, xn处的n阶差商。
根据定义得到差商表如下:
差商具有许多优美的性质,如对称性等,这里不做更多说明。
2) 差商形式的牛顿插值公式
根据上述差商定义,自然得到下面的差商公式表:(即把差商定义式展开)
从而,把后一项不断的代入前一项,就得到:
把最后一项去掉,作为余项(因其含有未知的x),就得到牛顿插上公式:
可以证明,这是关于x的n次多项式。
对于余项上面已经提到。
3) 牛顿插值与拉格朗日插值的比较
设拉格朗日插值函数为P(x),牛顿插值函数为N(x),显然二者均满足:
P(xi) = N(xi) = f(xi);
由代数多项式插值的唯一性(本质都在解最开始提到的方程),显然有:
P(x) = N(x);
因而,两个插值方法的余项也是相等的。这很有意思,两个余项风格迥异,形式完全不同,却证明了其相等。
当增加一个节点时,对于拉格朗日插值,必须摒弃前面的所有计算去重新计算,而牛顿插值公式却告诉我们,增加的节点只需要在其后再加一项。这种承袭性使得牛顿插值再某些情境下会比拉格朗日插值更加灵活易用。另外还需说明一点,计算余项时,牛顿插值公式余项由于不需要导数,故f(x)是由离散点或者导数不存在时仍然适用,这是拉格朗日余项计算所不能比拟的。
4) 差分与牛顿插值(简介)
差分概念请自行google。
当给出的节点等距时,使用差分能达到更优的效果,其计算更加简便。由于原理与上述基本相同,这里不做赘述。
四、埃尔米特插值(Hermite)
有时候,我们不仅要求插值函数在给定节点上函数值重合,而且要求若干阶导数也重合;即:要求插值函数φ(x)满足:
下面介绍两点三次Hermite插值。
已知:x0,x1;f(x0)= y0,f(x1) = y1;f’(x0) = y0’,f’(x1) = y1’;
求解:
模仿拉格朗日多项式的思想:
设
将插值条件带入可得:
以α0为例:
同理可得α1;
对于β0:
同理得β1;
因此,得到插值公式如下:
同理可得余项为:
NOTICE: 这里也有一个小小的技巧,导数的次数以及阶乘元就是你得到的条件的个数(这里是4个),而后面的连乘就是(x-xi),至于次数,每个xi用到了几个条件,自然是几次。这里不做证明,但事实上,这是快速写出余项的好方法。
一般情况下的埃尔米特插值公式:
根据上述计算思想,这里给出一般的公式:
n+1个节点,唯一确定2n+1次埃尔米特插值多项式:
余项为:
这就是埃尔米特插值问题。然而这有个问题,我们实际操作上并不可能得到每一个节点处的导数值,换言之,这种方法中看不中用!我将在后面的样条插值做更详细的阐述。
分段差值
什么是分段差值?
简而言之,分段差值就是对每一个分段区间(xi,xi+1)分别进行插值,则最后所得插值函数为一分段函数。
为什么用分段差值?
看一经典的例子:
在[-5,5]上,对函数f(x) = 1/(1+x^2)进行多项式插值,随着次数的增高,可以看到虽然与原函数重合的点越来越多,可以在端点附近抖动越大,这称为Runge现象。如下图:
在这种原因下,诞生了更加常用的使用分段插值。
一、分段线性插值
什么是分段线性插值?
在每个区间[xi,xi+1]上,用1阶多项式 (直线) 逼近 f (x):
根据前面的拉格朗日插值,在每个区间[xi,xi+1]上得:
图示如下:
余项:
R(x) = f''(x)/2! * (x-xi)(x-xi+1) <= (xi+1-xi)^2/8 *max(|f''(x)|);
二、分段三次(埃尔米特)插值
分段差值简单易行,又克服了Runge现象,但它却导致一阶导数不连续,有时候这可是不能容忍的!为了克服线性插值一阶导数不连续的缺点,可以采取下面的分段Hermite插值。
对于每一个区间,如果不仅仅知道端点处的函数值f(xi-1)、f(xi),还知道f‘(xi-1)、f’(xi),那么我们就有四个条件,自然可以进行三次插值。
求解方法:
1. 可以基于承袭求解。即基于两点求拉格朗日插值,由此构造,并根据导数条件求出未知系数。
2.埃尔米特方法。前面已经提到,这里不做赘述。
三、样条插值(主讲常用的三次样条插值)
前面的分段插值,本质上还是多项式插值。如果采用分段多项式插值, 则由于插值基函数只是局部活跃(它们的支集是局部紧致的), 结点上的误差可以被控制在小的范围内, 因而也带来了内在的高度稳定性. 这是分段插值的一大优势。
许多实际问题希望插值函数具有较高阶的整体光滑性. 此时, 高次Hermite插值或分段高次Hermite插值可以利用(注意:分段高次Lagrange插值和Newton插值等是做不到的,在插值结点上它们只能保证插值函数连续)。
但是,对于埃尔米特插值,必须知道每一个点的导数,在实际应用中这显然是不现实的,实际上这种求解也是中看不中用的。样条插值应运而生,解决了这个问题,并在各个领域大展身手,可以说是插值方法的里程碑。
什么是样条插值?
简而言之,就是依然对每一个小区间进行插值,但是我们不需要依赖于导数的已知;由于我们要做的就是使得端点出满足某种条件的光滑(一般来说,三次样条插值就是要满足二阶导数连续),根据这个要求,我们就可以在未知导数的情况下推导出样条函数。具体方法看下面。
样条函数的确定?
N个区间上,每个区间求三次样条插值函数,显然一共需要4n个条件。
1. 显式条件:
则n个区间各有2个条件,共2n个。
2. 隐式条件:
三次样条插值,一次导数连续且二次导数连续,即满足:
这样,n个区间可以找出n-1对导数相等条件。至此,可以找出2n+2(n-1) = 4n-2个条件。
3. 边界条件:
一般来说,在整个插值区间[a,b]上,会对边界[a,b]端点有状态的要求。这就是边界条件,有了这两个边界条件,就得到了需要的4n个条件。
边界条件一般有三种:
1) 一阶条件
即,给出端点处的一阶导数值。
2) 二阶条件
即,给出端点处的二阶导数值。
3) 循环(周期)条件
即,f(x)是一个周期函数时,端点a,b处的一阶导数值、二阶导数值分别相等。
确定了样条函数后,下面讲下具体的求解。
样条函数的计算?(二次的方法)
设分段样条函数s(x)在0到n个节点处的二阶导数值,为:
积分,得:每一区间上的样条函数 :
对于引入的两个常数c1,c2,使用两个边界条件,自然可以求解得到:
对所有区间,整理后得到关于Mj-1,Mj,Mj+1的方程:
其中:
这称作“三弯矩方程”。
最终,由着n-1个方程,再加两个边界条件(这里我使用一阶条件),直接带入,联立n+1个方程,得到:
至此,样条插值变为求解此矩阵形式的方程组。
可以看到,系数矩阵严格对角占优,矩阵可逆,方程组存在唯一解,从而可以解出M0…Mn。而此方程组的解法并不困难,计算非常快捷。
总结
这就是常见、常用的插值方法及其原理,如果有必要,我会对每个方法给出源码实现(不借助第三方库)并告诉你在MATLAB或者Python中如何快速的使用它们。这篇文章就写到这了,希望通过这篇文章你能对插值有更加深刻的认识并对你的学习和应用有所帮助。