2.7三次样条插值

三次样条插值的方法和思路摘要：1.三次样条插值的基本概念2.三次样条插值的数学原理3.三次样条插值的实现步骤4.三次样条插值的优缺点5.三次样条插值在实际应用中的案例正文：在日常的科学研究和工程应用中，我们经常会遇到需要对一组数据进行插值的问题。

插值方法有很多，其中三次样条插值是一种常见且有效的方法。

本文将从基本概念、数学原理、实现步骤、优缺点以及实际应用案例等方面，全面介绍三次样条插值的方法和思路。

一、三次样条插值的基本概念三次样条插值（Cubic Spline Interpolation）是一种基于分段多项式的插值方法。

它通过在各个节点上构建一条三次多项式曲线，使得这条曲线在节点之间满足插值条件，从而达到拟合数据的目的。

二、三次样条插值的数学原理三次样条插值的数学原理可以分为两个部分：一是分段三次多项式的构建，二是插值条件的满足。

1.分段三次多项式的构建假设有一组数据点序列为（x0，y0），（x1，y1），（x2，y2），（x3，y3），我们可以将这些数据点连接起来，构建一条分段三次多项式曲线。

分段三次多项式在每个子区间上都是一个三次多项式，它们之间通过节点值进行连接。

2.插值条件的满足为了使分段三次多项式在节点之间满足插值条件，我们需要在每个子区间上满足以下四个条件：（1）端点条件：三次多项式在区间的端点上分别等于节点值；（2）二阶导数条件：三次多项式在区间内的二阶导数等于节点间的斜率；（3）三阶导数条件：三次多项式在区间内的三阶导数等于节点间的曲率；（4）内部点条件：三次多项式在区间内部满足插值函数的连续性。

通过求解这四个条件，我们可以得到分段三次多项式的系数，从而实现插值。

三、三次样条插值的实现步骤1.确定插值节点：根据数据点的位置，选取合适的节点；2.构建分段三次多项式：根据节点值和插值条件，求解分段三次多项式的系数；3.计算插值结果：将待插值点的横坐标代入分段三次多项式，得到插值结果。

三次样条插值多项式实验的目的及意义：为了取得理想结果：在不增加更多的插值条件下，能够求得一个插值多项式，既有良好的逼近效果，又有好的光滑性，引进三次样条插值 多项式。

如果已知函数y=f(x)在节点a=x0<x1<…<xn=b 处的函数值和导数值：()i i x f y =，i=0,1,2,…,n如果S(x)满足条件：1. S(x)是一个分段的三次多项式且()i i y x S =；2. S(x)在[a,b]具有二阶连续导数。

则称S(x)是三次样条插值函数。

S(x)的具体形式为：()()()()⎪⎪⎩⎪⎪⎨⎧∈∈∈=-]12,121,01,[,...............][,][,n n n x x x x s x x x x s x x x x s x s其中()x S i 在[]i i x x ,1-上是三次多项式()iiiiid x c x b x a x S +++=23由插值条件()ii y x S =，i=0,1,2,…,n ，得n+1个条件。

边界条件一：()()nn y x S y x S '',''00== 边界条件二：()()nn y x S y x S '''',''''00==数学公式：()()2211133[2]()[2()]()i i i i i i i i i i ih x x x x h x x x x H x y y h h ---+-----=++2211122()()()()i i i i i i i ix x x x x x x x m m h h -------+ 算法描述：Step1:输入未知数X 及(xi,yi),i=0,1,…,n ； Step2:计算步长H[i]; Step3:计算[][]()⎪⎪⎪⎩⎪⎪⎪⎨⎧+=-=+=-+++i i i i i i i ii i i i i x x f x x f u g u hh h ,,311111λλλStep4:根据边界条件，求解相应的方程得到m0,m1,…, mn Step5:判断X 属于[]i i x x ,1-，i=1,2,…,n 中的哪一个Step6:计算()x s y i i ≈Step7:输出y. 程序原代码如下： #include "stdio.h" #define N 5 void main() { int i,k; float X,s,y0,yn;float a[N][N+1],h[N],u[N],v[N],g[N],m[N],p[N],q[N],w[N];printf("please input X:"); //X 为未知数的大小scanf("%f",&X);printf("please input x:"); //输入x的大小for(i=0;i<N;i++)scanf("%f",&a[i][0]);printf("please input y:"); //输入y的大小for(i=0;i<N;i++)scanf("%f",&a[i][1]);for(i=1;i<N;i++)h[i]=a[i][0]-a[i-1][0]; //计算步长for(i=1;i<N;i++){v[i]=h[i+1]/(h[i]+h[i+1]);u[i]=1-v[i];g[i]=3*u[i]*(a[i+1][1]-a[i][1])/h[i+1]+3*v[i]*(a[i][1]-a[i-1][1])/h[i]; }printf("\t(1)已知边界条件1\n");printf("\t(2)已知边界条件2\n");printf("请选择边界条件序号:");scanf("%d",&k);if(k==1){printf("请输入y0和yn的一阶导:"); //输入边界条件一scanf("%f%f",&m[0],&m[N-1]);p[0]=0; //用追赶法求解m[N]q[0]=0;g[1]=g[1]-v[1]*m[0];g[N-2]=g[N-2]-u[N-2]*m[N-1];for(i=1;i<N;i++){w[i]=2-u[i]*p[i-1];p[i]=v[i]/w[i];q[i]=(g[i]-u[i]*q[i-1])/w[i];}m[N-2]=q[N-2];for(i=N-3;i>0;i--)m[i]=q[i]-p[i]*m[i+1];printf("输出m[i]的值:\n");for(i=0;i<N;i++)printf("%f\n",m[i]);for(i=1;i<N;i++) //计算最终结果if(X>a[i-1][0]&&X<a[i][0])s=(h[i]+2*(X-a[i-1][0]))*(X-a[i][0])*(X-a[i][0])*a[i-1][1]/(h[i]*h[i]*h[i]) +(h[i]-2*(X-a[i][0]))*(X-a[i-1][0])*(X-a[i-1][0])*a[i][1]/(h[i]*h[i]*h[i])+ (X-a[i-1][0])*(X-a[i][0])*(X-a[i][0])*m[i-1]/(h[i]*h[i])+(X-a[i-1][0])*(X-a[i-1][0])*(X-a[i][0])*m[i]/(h[i]*h[i]);printf("s(%f)=%f\n",X,s);}if(k==2){printf("请输入y0和yn的二阶导:"); //输入边界条件二scanf("%f%f",&y0,&yn);g[0]=3*(a[1][1]-a[0][1])/h[1]-h[1]*y0/2;g[N-1]=3*(a[N-1][1]-a[N-2][1])/h[N-1]+h[N-1]*yn/2;q[0]=g[0];u[0]=1;v[N-1]=1;w[0]=2;for(i=1;i<N;i++){w[i]=2-v[i]*u[i-1]/w[i-1];q[i]=g[i]-v[i]*q[i-1]/w[i-1];}m[N-1]=q[N-1]/w[N-1];for(i=N-2;i>=0;i--)m[i]=(q[i]-u[i]*m[i+1])/w[i];printf("输出m[i]的值:\n");for(i=0;i<N;i++)printf("%f\n",m[i]);for(i=1;i<N;i++)if(X>=a[i-1][0]&&X<=a[i][0])s=(h[i]+2*(X-a[i-1][0]))*(X-a[i][0])*(X-a[i][0])*a[i-1][1]/(h[i]*h[i]*h[i]) +(h[i]-2*(X-a[i][0]))*(X-a[i-1][0])*(X-a[i-1][0])*a[i][1]/(h[i]*h[i]*h[i])+ (X-a[i-1][0])*(X-a[i][0])*(X-a[i][0])*m[i-1]/(h[i]*h[i])+(X-a[i-1][0])*(X-a[i-1][0])*(X-a[i][0])*m[i]/(h[i]*h[i]);printf("s(%f)=%f\n",X,s);}}数值计算：已知y=f(x)的如下数值求三次样条插值函数S（x），满足条件1.s’(0)=0,s’(4)=482.s’’(0)=0,s’’(4)=24Please input X:2.5Please input x:0 1 2 3 4 Please input y:-8 -7 0 19 56(1)已知边界条件1(2)已知边界条件2请选择边界条件的序号:1请输入y0和yn的一阶导:0 48 0.0000003.00000012.00000027.00000048.000000s(2.500000)=7.625000press any key tocontinue请选择边界条件的序号:2请输入y0和yn的二阶导:0 24 -0.0000003.00000011.99999927.00000248.0000007.625000press any key tocontinue s(2.500000)=7.625000 对计算结果进行评价分析：()()443845h M x S x f ≤-三次样条插值函数与三次Hermite 插值函数相比，不仅光滑度有提高，而且要求求解时还不需要增加内节点处的导数值，因此比较实用。

三次样条插值法根据李庆阳的《数值分析》这本教材中的讲解编写的程序，使用的是第一边界条件，用追赶法求解了M矩阵。

为了调用方便，我将整个函数所有的信息构造成一个结构体，输入插值点的坐标和数量，定义边界条件后，将这个结构体的指针作为参数传给Spline3（）函数，就完成了函数计算，计算结果也存储在该结构体中。

程序如下：/*=================================================================== ====*///=====================三次样条插值的函数S(x)实现=============================// 创建人：汪雅楠// 北京交通大学 QQ312818820// 说明：根据研究生教材《数值分析》（李庆杨）第51页~第56页编写/* 初始条件： 1. 已知两端的一阶导数值2. 已知两端的二阶导数值3. 周期样条函数### 此函数选择1条件 ###函数建立： 1. 设定样条函数S(x)的一阶导数为变量ki,用分段三次Hermitte差值2. 设定样条函数S(x)的二阶导数为变量Ki,用分段积分### 此函数选择2方法 ###矩阵求解：追赶法求解严格三对角占优矩阵{M},根据教材第195页编写*//*=================================================================== ====*/#include <stdio.h>///////////////////////////////////////////////////////////////////// ///////////#define MAXNUM 50 //定义样条数据区间个数最多为50个typedef struct SPLINE //定义样条结构体，用于存储一条样条的所有信息{ //初始化数据输入float x[MAXNUM+1]; //存储样条上的点的x坐标，最多51个点float y[MAXNUM+1]; //存储样条上的点的y坐标，最多51个点unsigned int point_num; //存储样条上的实际的点的个数float begin_k1; //开始点的一阶导数信息float end_k1; //终止点的一阶导数信息//float begin_k2; //开始点的二阶导数信息//float end_k2; //终止点的二阶导数信息//计算所得的样条函数S(x)float k1[MAXNUM+1]; //所有点的一阶导数信息float k2[MAXNUM+1]; //所有点的二阶导数信息//51个点之间有50个段，func[]存储每段的函数系数float a3[MAXNUM],a1[MAXNUM];float b3[MAXNUM],b1[MAXNUM];//分段函数的形式为Si(x) = a3[i] * {x(i+1) - x}^3 + a1[i] * {x(i+1) - x} +// b3[i] * {x - x(i)}^3 + b1[i] * {x - x(i)}//xi为x[i]的值，xi_1为x[i+1]的值}SPLINE,*pSPLINE;typedef int RESULT; //返回函数执行的结果状态，下面为具体的返回选项#ifndef TRUE#define TRUE 1#endif#ifndef FALSE#define FALSE -1#endif#ifndef NULL#define NULL 0#endif#ifndef ERR#define ERR -2#endif///////////////////////////////////////////////////////////////////// //////////////*=================================================================== ============*** 函数名称：Spline3()*** 功能说明：完成三次样条差值，其中使用追赶法求解M矩阵*** 入口参数：(pSPLINE)pLine 样条结构体指针pLine中的x[],y[],num,begin_k1,end_k1*** 出口参数：(pSPLINE)pLine 样条结构体指针pLine中的函数参数*** 返回参数：返回程序执行结果的状态TRUE or FALSE===================================================================== ===========*/RESULT Spline3(pSPLINE pLine){float H[MAXNUM] = {0}; //小区间的步长float Fi[MAXNUM] = {0}; //中间量float U[MAXNUM+1] = {0}; //中间量float A[MAXNUM+1] = {0}; //中间量float D[MAXNUM+1] = {0}; //中间量float M[MAXNUM+1] = {0}; //M矩阵float B[MAXNUM+1] = {0}; //追赶法中间量float Y[MAXNUM+1] = {0}; //追赶法中间变量int i = 0;////////////////////////////////////////计算中间参数if((pLine->point_num < 3) || (pLine->point_num > MAXNUM + 1)){return ERR; //输入数据点个数太少或太多}for(i = 0;i <= pLine->point_num - 2;i++){ //求H[i]H[i] = pLine->x[i+1] - pLine->x[i];Fi[i] = (pLine->y[i+1] - pLine->y[i]) / H[i]; //求F[x(i),x(i+1)] }for(i = 1;i <= pLine->point_num - 2;i++){ //求U[i]和A[i]和D[i]U[i] = H[i-1] / (H[i-1] + H[i]);A[i] = H[i] / (H[i-1] + H[i]);D[i] = 6 * (Fi[i] - Fi[i-1]) / (H[i-1] + H[i]);}//若边界条件为1号条件，则U[i] = 1;A[0] = 1;D[0] = 6 * (Fi[0] - pLine->begin_k1) / H[0];D[i] = 6 * (pLine->end_k1 - Fi[i-1]) / H[i-1];//若边界条件为2号条件，则//U[i] = 0;//A[0] = 0;//D[0] = 2 * begin_k2;//D[i] = 2 * end_k2;/////////////////////////////////////////追赶法求解M矩阵B[0] = A[0] / 2;for(i = 1;i <= pLine->point_num - 2;i++){B[i] = A[i] / (2 - U[i] * B[i-1]);}Y[0] = D[0] / 2;for(i = 1;i <= pLine->point_num - 1;i++){Y[i] = (D[i] - U[i] * Y[i-1]) / (2 - U[i] * B[i-1]);}M[pLine->point_num - 1] = Y[pLine->point_num - 1];for(i = pLine->point_num - 1;i > 0;i--){M[i-1] = Y[i-1] - B[i-1] * M[i];}//////////////////////////////////////////计算方程组最终结果for(i = 0;i <= pLine->point_num - 2;i++){pLine->a3[i] = M[i] / (6 * H[i]);pLine->a1[i] = (pLine->y[i] - M[i] * H[i] * H[i] / 6) / H[i];pLine->b3[i] = M[i+1] / (6 * H[i]);pLine->b1[i] = (pLine->y[i+1] - M[i+1] * H[i] * H[i] / 6) /H[i]; }return TRUE;}///////////////////////////////////////////////////////////////////// /////////////SPLINE line1;pSPLINE pLine1 = &line1;///////////////////////////////////////////////////////////////////// /////////////main(){line1.x[0] = 27.7;line1.x[1] = 28;line1.x[2] = 29;line1.x[3] = 30;line1.y[0] = 4.1;line1.y[1] = 4.3;line1.y[2] = 4.1;line1.y[3] = 3.0;line1.point_num = 4;line1.begin_k1 = 3.0;line1.end_k1 = -4.0;Spline3(pLine1);return 0;}///////////////////////////////////////////////////////////////////// /////////////。