三次样条插值算法
1 总体说明
三次样条插值算法是一种计算量和效果都比较理想的插值算法。关于三次样条插值算法的原理这里不做过多的解释，下面的代码是我在网上收集了两种C++实现版本的基础上自己整合的一个版本。由于本人刚接触C++不久，水平有限。没有使用模板机制将代码做的更通用。关于算法实现有下面几点说明。
- 所有有关的类都被包含到SplineSpace命名空间中。
- SplineSpace中一个有三个类分别是异常类(SplineFailure),接口类(SplineInterface)和实现类(Spline)。有一个枚举类型说明边界条件(BoundaryCondition),取值为：GivenFirstOrder和GivenSecondOrder。分别对应I型边界条件和II型边界条件。
- 接口类定义了Spline在实现的过程中必须要有的三个方法：单点插值、多点插值和自动生成插值序列。
- 异常类是可能被实现类抛出的类，如果在实现类的运行过程中出现了已知数据过少构造失败、使用了外插值、设定输出点数过少等行为会抛出该类。因此应该将插值的过程用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的插值结果y
cout<<"x="<<x<<"时的插值结果为:"<<y<<endl;
}
catch(SplineFailure sf)
{
cout<<sf.GetMessage()<<endl;
}
getchar(); //程序暂停
}
此时屏幕会输出"x=4.5时的插值结果为2.71107"。
- 可以直接构造Spline对象进行使用，也可以将它转化成对应的接口使用。
- 默认边界条件为II型边界条件(已知边界的二阶导数)，默认左右边界二阶导数为都0，即自然边界条件。
- 已知的数据点数要至少为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的插值结果y
for(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的插值结果y
for(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
clc
x0=[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++源程序代码：
下载源码：三次样条插值C++代码
