三次样条插值

因为网上并没有三次样条插值的简单实现

所以写了一个最通俗易懂的实现方式

方法用途：

求出一条平滑的曲线 连接给出的一系列数据(横纵坐标)

方法说明：

• 每两个点之间构建一条曲线（ax³+bx²+cx+d）
• 因此有N个数据 就有4*（N-1）个未知数
• 两条曲线的端点处取值相等（连接）
• 每条曲线连接处的导数相等 (平滑）
• 最左端和最右端的二阶导数等于0（一般做法）

代码实现：

#include<stdio.h>
#include<math.h>

#define N 25//数据个数
#define M 4*(N-1)//未知数个数

double f1a(double x)
    return x*x*x;

double f1b(double x)
    return x*x;

double f1c(double x)
    return x;

double f2a(double x)
    return 3*x*x;

double f2b(double x)
    return 2*x;

double f3a(double x)
    return 6*x;

int main()
{
    double a[M][M+1]={0};//求解矩阵
    double data[2][N]={{0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24},
                    {5,5,4,4,4,3,3,3,5,5,9,11,12,13,13,14,13,11,10,9,8,8,7,6,5}};//数据矩阵
    int i;//行
    int j;//列
    int k;
    for(i=0;i<N-1;i++)//求每个左端点的函数值 (N-1)个方程
    {
        a[i][M]=data[1][i];//方程的右端＝y
        j=4*i;
        //求a,b,c的系数
        a[i][j]=f1a(data[0][i]);
        a[i][j+1]=f1b(data[0][i]);
        a[i][j+2]=f1c(data[0][i]);
        a[i][j+3]=1;//d的系数为1
    }

    for(k=0;k<N-1;k++,i++)//求每个右端点的函数值 (N-1)个方程
    {
        a[i][M]=data[1][k+1];//方程的右端＝y
        j=4*k;
        //求a,b,c的系数
        a[i][j]=f1a(data[0][k+1]);
        a[i][j+1]=f1b(data[0][k+1]);
        a[i][j+2]=f1c(data[0][k+1]);
        a[i][j+3]=1;//d的系数为1
    }

    for(k=0;k<N-2;k++,i++)//每段函数连接处的导数相等 (N-2)个方程
    {
        a[i][M]=0;//方程的右端＝0
        j=4*k;
        //求a,b,a1,b1的系数
        a[i][j]=f2a(data[0][k+1]);
        a[i][j+1]=f2b(data[0][k+1]);
        a[i][j+2]=1;//c的系数为1
        a[i][j+4]=-1*f2a(data[0][k+1]);
        a[i][j+5]=-1*f2b(data[0][k+1]);
        a[i][j+6]=-1;
    }

    for(k=0;k<N-2;k++,i++)//每段函数连接处的二阶导数相等 (N-2)个方程
    {
        a[i][M]=0;//方程的右端＝0
        j=4*k;
        //求a,a1的系数
        a[i][j]=f3a(data[0][k+1]);
        a[i][j+1]=2;//b的系数为2
        a[i][j+4]=-1*f3a(data[0][k+1]);
        a[i][j+5]=-2;
    }

    //最左端和最右端函数的二阶导数等于0 2个方程
    a[i][M]=0;
    a[i][0]=f3a(data[0][0]);
    a[i][1]=2;
    i++;
    a[i][M]=0;
    a[i][M-4]=f3a(data[0][N-1]);
    a[i][M-3]=2;
/*    //输出构造完的矩阵
    printf("构造的矩阵为：\n");
    for(j=0;j<M;j++)
    {
        for(i=0;i<M+1;i++)
            printf("%.3f\t",a[j][i]);
        printf("\n");
    }*/
    qz(a);

    return 0;
}

求解方程组：

int qz(double a[M][M+1])//列主元消去法
{
    int i,j,k,temp1;
    double temp,x,max;

    for(k=0;k<M-1;k++)
    {
        temp1=k;
        max=a[k][k];
        for(j=k+1;j<M;j++)//找出列主元
            if(max<fabs(a[j][k]))
            {
                max=fabs(a[j][k]);
                temp1=j;
            }
        for(j=k;j<M+1;j++)//第一行与列主元行交换
        {
            x=a[temp1][j];
            a[temp1][j]=a[k][j];
            a[k][j]=x;
        }
        for(j=k+1;j<M;j++)//使每行第一个数为0
        {
            temp=a[j][k]/a[k][k];
            for(i=k;i<=M+1;i++)
            {
                a[j][i]-=a[k][i]*temp;        
            }
        }
    }
    for(j=0;j<M;j++)//使对角线为1
    {
        temp=a[j][j];
        for(i=j;i<=M+1;i++)
            a[j][i]/=temp;
    }


    for(k=M-1;k>0;k--)//回代
    {
        for(j=k-1;j>=0;j--)
        {
            temp=a[j][k]/a[k][k];
            for(i=k;i<=M+1;i++)
            {
                a[j][i]-=a[k][i]*temp;        
            }
        }
    }

/*    printf("解得：\n");
    for(j=0;j<M;j++)
    {
        for(i=0;i<M+1;i++)
            printf("%.3f\t",a[j][i]);
        printf("\n");
    }*/
    for(i=0;i+3<M;i+=4)
    {
        printf("%fx^3+%fx^2+%fx+%f\n",a[i][M],a[i+1][M],a[i+2][M],a[i+3][M]);
    }
    return 1;
}