样条（Spline）函数是由舍恩伯格于1946年提出的。样条是富有弹性的细木条或有机玻璃条，它的作用相当于“万能”曲线板。早期船舶、汽车、飞机放样时用铅压铁压住样条，使其通过一系列型值点，调整压铁达到设计要求后绘制其曲线，称为样条曲线。这样设计曲线的方法在20世纪六七十年代得到了广泛应用。
通常单一的曲线段或曲面片难以表达复杂的形状，必须将一些曲线段拼接成组合曲线，或将一些曲面片拼接成组合曲面，才能表达复杂的形状。为了保证在结合点处光滑过渡，需要满足连续性条件。连续性条件有两种：参数连续性（Parametric Continuity）与几何连续性（Geometric Continuity）。
零阶参数连续性，记为C0，指相邻两段曲线在结合点处具有相同的坐标。
一阶参数连续性，记为C1，指相邻两段曲线在结合点处具有相同的一阶导数
二阶参数连续性，记为C2，指相邻两段曲线在结合点处具有相同的一阶导数和二阶导数。
与参数连续性不同的是，几何连续性只要求参数成比例而非相等。
已知n个型值点Pi(xi, yi), i = 1, 2,…, n且a = x1 < x2 <…< xn = b。若 y=s(x)满足下列条件：
（1）型值点Pi在函数 y=s(x) 上。
（2） s(x) 在整个区间[a, b]上二阶连续可导。
（3）在每个子区间[xi, xi+1], i = 1, 2,…, n-1上，分段函数 都是参数x的三次多项式。则称函数 是过型值点的三次样条函数，由三次样条函数构成的曲线称为三次样条曲线。
第i段的 si(x) 表示为:
式中，ai，bi，ci，di为待定系数，i=1，2 … ，n-1。
第i段曲线的首端通过Pi（xi，yi），末端通过Pi+1（xi+1，yi+1）
Pi处的二阶导数为Mi。
其中各型值点的弯矩Mi的力学解释如下：
子区间示意图如下：
用型值点处的二阶导数Mi表示的三次B样条函数为:
见个人资源：三次样条插值函数求解过程.pdf
（1）读入n个型值点且满足其x坐标递增。
（2）根据实际情况确定三次样条曲线的边界条件。
（3）计算曲线的系数，将其表达为型值点二阶导数的函数。
（4）用追赶法求解三弯矩方程。
（5）将求解出的参数代入三次样条函数的表达式中，构造三次样条曲线。
（6）循环访问每个节点。在每个子区间内，按照精度要求，使用直线段连接各段内的若干等分点，即可绘制出三次样条曲线。
首先定义数据结构p2.cpp与头文件(主要如下)：
CP2::CP2(double x, double y)
{
this->x = x;
this->y = y;
}
之后在主函数中定义读取函数：
// CGeometricfiguretestViewmessage handlers
void CGeometricfiguretestView::ReadPoint()
{
P[1].x = -340, P[1].y = -200;
P[2].x = -150, P[2].y = 0;
P[3].x = 0, P[3].y = -50;
P[4].x = 100, P[4].y = -100;
P[5].x = 250, P[5].y = -100;
P[6].x = 350, P[6].y = -50;
}
void CGeometricfiguretestView::DrawDataPoint(CDC* pDC)//绘制型值点
{
CBrush NewBrush, *pOldBrush;
NewBrush.CreateSolidBrush(RGB(0, 0, 0));
pOldBrush = pDC->SelectObject(&NewBrush);
for( int i = 1; i < 7; i++)
pDC->Ellipse(ROUND(P[i].x - 5), ROUND(P[i].y - 5), ROUND(P[i].x + 5), ROUND(P[i].y + 5));
pDC->SelectObject(pOldBrush);
}
具体过程见注释：
void CGeometricfiguretestView::DrawCubicSpline(CDC* pDC)//三次样条曲线
{
int n = 6;
const int dim = 7;//二维数组维数
double b1 = 10, bn = 10;//边界条件:"夹持端",给出起点和终点的一阶导数
double h[dim], lambda[dim], mu[dim], D[dim];
double l[dim], m[dim], u[dim];
double M[dim], K[dim];
double a[dim], b[dim], c[dim], d[dim];
for(int i = 1; i < n; i++)//计算hi＝xi+1-xi
h[i] = P[i+1].x - P[i].x;
for(int i = 2; i < n; i++)
{
lambda[i] = h[i-1]/(h[i-1]+h[i]);//计算λ
mu[i] = h[i]/(h[i-1]+h[i]);//计算μ
D[i] = 6/(h[i-1]+h[i])*((P[i+1].y-P[i].y)/h[i]-(P[i].y-P[i-1].y)/h[i-1]);//计算D
}
D[1]=6*((P[2].y-P[1].y)/h[1]-b1)/h[1];//夹持端的D[1]
D[n]=6*(bn-(P[n].y-P[n-1].y)/h[n-1])/h[n-1];//夹持端的D[n]
mu[1]=1;//夹持端的μ[1],
lambda[n]=1;//夹持端的λ[n]
//追赶法求解三弯矩方程
l[1]=2;
u[1]=mu[1]/l[1];
for(int i=2; i <= n; i++)
{
m[i]=lambda[i];
l[i]=2-m[i]*u[i-1];
u[i]=mu[i]/l[i];
}
K[1] = D[1]/l[1];//解LK=D
for(int i = 2; i <= n;i++)
{
K[i]=(D[i]-m[i]*K[i-1])/l[i];
}
M[n] = K[n];//解UM=K
for(int i = n-1; i >= 1;i--)
{
M[i]=K[i]-u[i]*M[i+1];
}
//计算三次样条函数的系数
for(int i = 1; i < n; i++)
{
a[i] = P[i].y;
b[i] = (P[i+1].y-P[i].y)/h[i] - h[i]*(M[i]/3+M[i+1]/6);
c[i] = M[i]/2;
d[i] = (M[i+1]-M[i])/(6*h[i]);
}
pDC->MoveTo(ROUND(P[1].x), ROUND(P[1].y));
double xStep = 0.5;//x的步长
double x, y;//当前点
for(int i = 1; i < n; i++)//循环访问每个节点
{
for(x = P[i].x; x < P[i+1].x; x += xStep)//循环访问每个节点
{
y=a[i]+b[i]*(x-P[i].x)+c[i]*(x-P[i].x)*(x-P[i].x)+d[i]*(x-P[i].x)*(x-P[i].x)*(x-P[i].x);
pDC->LineTo(ROUND(x), ROUND(y));//绘制样条曲线
}
}
}
void CGeometricfiguretestView::OnDraw(CDC* pDC)
{
CTestDoc* pDoc = GetDocument();
ASSERT_VALID(pDoc);
if (!pDoc)
return;
// TODO: add draw code for native data here
CRect rect;//定义客户区矩形
GetClientRect(&rect);//获得客户区矩形的信息
pDC->SetMapMode(MM_ANISOTROPIC);//自定义二维坐标系
pDC->SetWindowExt(rect.Width(), rect.Height());//设置窗口范围
pDC->SetViewportExt(rect.Width(), -rect.Height());//设置视区范围，且x轴水平向右为正，y轴垂直向上为正
pDC->SetViewportOrg(rect.Width() / 2, rect.Height() / 2);//设置客户区中心为二维坐标系原点
rect.OffsetRect(-rect.Width() / 2, -rect.Height() / 2);//rect矩形与客户区重合
ReadPoint();
DrawDataPoint(pDC);
DrawCubicSpline(pDC);
}
编译运行，可见如下：
我们也可以改变其中一个边界约束条件：
double b1 = 10, bn = -5;//边界条件:"夹持端",给出起点和终点的一阶导数
运行可以看到如下：
