三次样条插值算法

1 总体说明

三次样条插值算法是一种计算量和效果都比较理想的插值算法。关于三次样条插值算法的原理这里不做过多的解释，下面的代码是我在网上收集了两种C++实现版本的基础上自己整合的一个版本。由于本人刚接触C++不久，水平有限。没有使用模板机制将代码做的更通用。关于算法实现有下面几点说明。

1. 所有有关的类都被包含到SplineSpace命名空间中。
2. SplineSpace中一个有三个类分别是异常类(SplineFailure),接口类(SplineInterface)和实现类(Spline)。有一个枚举类型说明边界条件(BoundaryCondition),取值为：GivenFirstOrder和GivenSecondOrder。分别对应I型边界条件和II型边界条件。
3. 接口类定义了Spline在实现的过程中必须要有的三个方法：单点插值、多点插值和自动生成插值序列。
4. 异常类是可能被实现类抛出的类，如果在实现类的运行过程中出现了已知数据过少构造失败、使用了外插值、设定输出点数过少等行为会抛出该类。因此应该将插值的过程用try...catch(SplineFailure sf)包裹起来。如：

double x0[2]={1,2};
double y0[2]={3,4};try
{SplineInterface* sp = new Spline(x0,y0,2);//...
}
catch(SplineFailure sf)
{cout<<sf.GetMessage()<<endl;
}

上面代码就会抛出异常并显示“构造失败,已知点数过少”。

2 插值方法调用

2.1单点插值

调用方法如下：

#include <iostream>
#include "Spline.h"using namespace std;
using namespace SplineSpace;int main(void)
{//单点插值测试double x0[5]={1,2,4,5,6};       //已知的数据点double y0[5]={1,3,4,2,5};try{//Spline sp(x0,y0,5,GivenSecondOrder,0,0);SplineInterface* sp = new Spline(x0,y0,5);  //使用接口，且使用默认边界条件double x=4.5;double y;sp->SinglePointInterp(x,y); //求x的插值结果ycout<<"x="<<x<<"时的插值结果为:"<<y<<endl;}catch(SplineFailure sf){cout<<sf.GetMessage()<<endl;}getchar();  //程序暂停
}

此时屏幕会输出"x=4.5时的插值结果为2.71107"。

1. 可以直接构造Spline对象进行使用，也可以将它转化成对应的接口使用。
2. 默认边界条件为II型边界条件(已知边界的二阶导数)，默认左右边界二阶导数为都0，即自然边界条件。
3. 已知的数据点数要至少为3（3个点对应一条曲线）。

2.2多点插值

调用方法如下（省略了和上面重复的代码部分）：

//多点插值测试double x0[5]={1,2,4,5,6};       //已知的数据点
double y0[5]={1,3,4,2,5};double x[4] = {1.5,2.5,3.5,4.5};    //插值点
double y[4];
double leftBound=0,RightBound=0;    //边界导数try
{Spline sp(x0,y0,5,GivenSecondOrder,leftBound,RightBound);sp.MultiPointInterp(x,4,y);           //求x的插值结果yfor(int i = 0;i < 4;i++){cout<<"x="<<x[i]<<"时的插值结果为:"<<y[i]<<endl;}
}
catch(SplineFailure sf)
{cout<<sf.GetMessage()<<endl;
}
getchar();  //程序暂停

显示结果如下：

x=1.5时的插值结果为:2.01383
x=2.5时的插值结果为:3.8978
x=3.5时的插值结果为:4.62372
x=4.5时的插值结果为:2.71107

2.3自动生成插值序列

自动生成插值序列方法会根据指定的插值点数，在插值自变量区间等间距的生成插值点。这个方法在绘图的时候非常好用。指定的插值点数越多，绘制出来的图形质量越好。调用方法如下：

//自动插值测试double x0[5]={1,2,4,5,6};       //已知的数据点
double y0[5]={1,3,4,2,5};double x[10];   //插值点
double y[10];try
{SplineInterface* sp = new Spline(x0,y0,5);    //使用接口，且使用默认边界条件sp->AutoInterp(10,x,y);           //求x的插值结果yfor(int i = 0;i < 10;i++)cout<<x[i]<<",";cout<<endl;for(int i = 0;i < 10;i++)cout<<y[i]<<",";
}
catch(SplineFailure sf)
{cout<<sf.GetMessage()<<endl;
}
getchar();  //程序暂停

显示结果如下：

1,1.55556,2.11111,2.66667,3.22222,3.77778,4.33333,4.88889,5.44444,6,
1,2.12528,3.21225,4.14572,4.639,4.38218,3.1524,2.02937,2.79183,5,

将上面结果复制粘贴到matlab中绘图，并且和matlab中的三次样条拟合结果作比较。matlab代码如下:

clear
clcx0=[1,2,4,5,6];         %已知的自变量
y0=[1,3,4,2,5];         %已知的因变量pp=csape(x0,y0,'second',[0,0]); %二阶边界条件,且导数都为0
x = linspace(1,6,200);
y = ppval(pp,x);%C++算法得出了的结果
xCpp = [1,1.55556,2.11111,2.66667,3.22222,3.77778,4.33333,4.88889,5.44444,6];
yCpp = [1,2.12528,3.21225,4.14572,4.639,4.38218,3.1524,2.02937,2.79183,5];h=figure;
plot(xCpp,yCpp,'*',x,y);
legend('c++插值结果','matlab插值结果')
title('插值结果')
h.Name = '插值结果';
h.NumberTitle = 'off';

结果如下：

可以看到C++中的算法得到的结果和matlab中得到的结果相同。
C++源程序代码：
有CSDN下载积分的童鞋可以使用右边的下载链接支持一下：三次样条插值C++代码
没有积分的童鞋也可以直接在本站下载源码：三次样条插值C++代码