三次插值多项式(三次样条插值原理)

三次样条(cubic spline)插值

1. 拟合与插值2. 分段三次样条插值3. 三对角线性系统的求解4. 固定边界三次样条的求解5. 附录5.1 龙格库塔现象

1. 拟合与插值

一个情形

有时我们获得了某条曲线上若干个离散点的信息,却不知道具体的曲线方程,又希望得到这些离散点之外的信息。

两种方法

拟合:不要求方程通过所有已知点,讲究神似,曲线整体趋势一致即可。

三次插值多项式(三次样条插值原理)

插值:要求方程通过所有已知点,讲究形似。

样条插值

在数值分析这个数学分支中,样条插值是使用一种名为样条的特殊分段多项式进行插值的形式。由于样条插值可以使用低阶多项式样条实现较小的插值误差,这样就避免了使用高阶多项式所出现的龙格现象(见附录1),所以样条插值得到了流行。常用的是三次样条插值。

2. 分段三次样条插值

所谓分段三次样条插值,顾名思义,就是将区间 分成 个区间

,共有 个点,其中端点值

。三次样条就是说每个分段区间上的曲线都是一个三次多项式方程 ,并满足以下条件。

1. 在每个分段区间

上,

都是一个三次多项式方程。

2. 在所有已知点上满足插值条件,即

3. 曲线光滑,即

平滑, 连续。

我们可将构造如下的三次多项式方程:

我们称这个方程为三次样条函数 。显然,每个分段区间上的三次样条函数均有 个未知系数,有 个分段区间,共有 个未知系数,需要 个方程来求解。

分析

我们需要寻找出 个方程去求解 个未知系数。

1. 所有已知点满足插值条件,即

。除端点外的 内部点

满足

,共有

,共有 个方程。总计 个方程。

2. 个内部点一阶导数连续,即

。总计 个方程。

3. 个内部点二阶导数连续,即

。总计 个方程。

这样,我们就找出了 个方程了,还差 个方程就可以解出所有未知系数了。这 个方程我们通过边界条件给出。

自然边界(natural spline):指定端点的二阶导数为 ,即

固定边界(clamped spline):指定端点一阶导数为特定值,即

非扭结边界(not-a-knot spline):使最前面两个插值点的三阶导数值相等、最后面两个插值点的三阶导数值相等,即

推导

为了方便推导,我们令:

可得:

并令

可得:

可得:

可得:

我们令

,则 ,由 可得

将 、 、

代入 中,可得:

至此,

均表示成二阶导的关系式,将其代入 可得:

这样,我们就能构造出一个以二阶导 为未知数的线性方程组了。

自然边界:

可以看出,左侧的系数矩阵为严格对角占优矩阵,故线性方程组有唯一解。再根据

关于二阶导的表达式就能求出每个分段区间上的具体方程了。

固定边界:

有: ,而

,则有:

有:

,将

的表达式带入可得:

则线性方程组为:

非扭结边界:

有:

,代入 的表达式可得:

则线性方程组为:

算法

输入: 个数据节点

[1] 计算步长

[2] 配置并求解线性方程组;

[3] 计算样条曲线的系数;

[4] 在每个分段区间

上创建三次多项式方程

3. 三对角线性系统的求解

对于系数矩阵为严格对角占优三对角矩阵的线性方程组,有唯一解。且其 分解、前代及回代过程,只需

次操作。

在编程方面,通常无需存储整个 系数矩阵,使用三个列向量来代替。

voidtridag(intn,doublea[],doubleb[],doublec[],doubler[],doubleu[]){doublebet;double*gma=vector(1,n-1,(double)(NULL));bet=b[0];u[0]=r[0]/bet;for(inti=1;i<n;i ){//分解及前代gma[i]=c[i-1]/bet;bet=b[i]-a[i]*gma[i];u[i]=(r[i]-a[i]*u[i-1])/bet;}for(inti=n-2;i>=0;i–){//回代u[i]-=gma[i 1]*u[i 1];}}4. 固定边界三次样条的求解

仿照前述三对角线性系统的求解方法,对于固定边界的三次样条,编写如下样条插值的功能块。

/***在处理列表函数时,程序spline仅调用一次。一旦调用过后,对任意的xx,插值函数的值通过splint获得。**//***给定数组x[0…n-1](严格单调递增)与y[0…n-1]构成一个列表函数,即y[i]= f(x[i])。* yp1与ypn表示插值函数分别在第1个与第n-1个点处的一阶导数。*函数修改y2[0…n-1],表示插值函数在表中点x_i处的二阶导数。**/voidspline(doublex[],doubley[],intn,doubleyp1,doubleypn,doubley2[]){doublebet;double*gma=vector(1,n-1,(double)(NULL));bet=2*(x[1]-x[0]);y2[0]=6*((y[1]-y[0])/(x[1]-x[0])-yp1)/bet;for(inti=1;i<n-1;i ){//分解及前代gma[i]=(x[i]-x[i-1])/bet;bet=2*(x[i 1]-x[i-1])-(x[i]-x[i-1])*gma[i];y2[i]=(6*((y[i 1]-y[i])/(x[i 1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]))-(x[i]-x[i-1])*y2[i-1])/bet;}gma[n-1]=(x[n-1]-x[n-2])/bet;bet=2*(x[n-1]-x[n-2])*gma[n-1];y2[n-1]=(6*(ypn-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]))-(x[n-1]-x[n-2])*y2[n-2])/bet;for(inti=n-2;i>=0;i–){//回代y2[i]-=gma[i 1]*y2[i 1];}}/***给定数组x[0…n-1]和y[0…n-1]作为列表函数,数组y2a[0…n-1]由spline程序输出给定。*给定要计算的xx值,程序通过yy返回其三次样条插值结果。**/voidsplint(doublex[],doubley[],doubley2[],intn,doublexx,double*yy){intk,kl=0,kh=n-1;doublea,b,c,d;//findthecorrespondingintervalwithbinarysearchingwhile(kh-kl>1){k=(kh kl)/2;if(x[k]>xx){kh=k;}else{kl=k;}}a=y[kl];b=(y[kh]-y[kl])/(x[kh]-x[kl])-(x[kh]-x[kl])*y2[kl]/2-(x[kh]-x[kl])*(y2[kh]-y2[kl])/6;c=y2[kl]/2;d=(y2[kh]-y2[kl])/(x[kh]-x[kl])/6;*yy=a b*(xx-x[kl]) c*(xx-x[kl])*(xx-x[kl]) d*(xx-x[kl])*(xx-x[kl])*(xx-x[kl]);}

测试

针对上述功能块

,进行如下测试:在

上选取若干个控制点,计算控制点之外的插值点。

#include<stdio.h>#include<stdlib.h>#include”nrutil.h”intmain(){FILE*fp;intn=21;doublexx;doubleyy;doubleyp1=-20;doubleypn=20;double*x=vector(0,n-1,(double)(NULL));double*y=vector(0,n-1,(double)(NULL));double*y2=vector(0,n-1,(double)(NULL));//selectncontrolpointsfromf(x)=x^2fp=fopen(“control.dat”,”w”);for(inti=0;i<n; i){x[i]=i-10.0;y[i]=(i-10.0)*(i-10.0);fprintf(fp,”%lf\t%lf\n”,x[i],y[i]);}fclose(fp);spline(x,y,n,yp1,ypn,y2);//calculatetheinterpolatingpointsfp=fopen(“interpolate.dat”,”w”);xx=-10.0;while(xx<=10.0){splint(x,y,y2,n,xx,&yy);fprintf(fp,”%lf\t%lf\n”,xx,yy);xx =0.1;}fclose(fp);return0;}

画图脚本:

setxlabel’x’setylabel’y’setmultiplotlayout1,1title’cubicspline’plot”control.dat”u1:2t’control’wpointspointtype7linecolor’red’linewidth1pointsize1,\”interpolate.dat”u1:2t’interpolatation’wpointspointtype6linecolor’blue’linewidth1pointsize1unsetmultiplot

效果图:

5. 附录5.1 龙格库塔现象

龙格现象是在一组等间距插值点上使用高次多项式插值时出现的区间边缘处的振荡问题,它表明使用高次多项式插值并不总是能提高准确性。随着多项式阶次的增加,以这种方式产生的多项式

实际上可能偏离 ,通常发生在靠近插值点的端点。

解决龙格现象的方法:

使用分段多项式样条可以避免龙格库塔现象。如果要减小插值误差,那么可以增加构成样条的多项式的数目,而不必是增加多项式的阶次。

发表评论

登录后才能评论