三次样条c++代码实现

￡神魔★判官ぃ 2022-05-23 00:38 414阅读 0赞

三次样条是一种效果很好的分段插值方法，不仅光滑性好，同时也避免了高次插值龙格振荡现象的发生。

首先先对三次样条插值思想做简单的介绍

![70][]

干货见下：

使用的例子是龙格函数的数据

#include<iostream>
    
    #include<math.h>
    
    using namespace std;
    
    
    
    int main()
    
    {
    
        double a[21][3];//第一列已知点的横坐标，2.3列是其二階均差
    
        double b[21];//已知点纵坐标
    
        double h[21];//存放步长
    
        double c1[21];//μ
    
        double c2[21];//λ
    
        double c3[21];//β
    
        double c4[21];//Ly=f,存放y的值
    
        double d[21];//d
    
        double M[21];//Ux=y,存放x的值，也就是原方程要求出的M
    
        double x[25]={-5,-4.75,-4.5,-4.25,-4,-3.75,-3.5,-3.25,-3,-2.75,-2.5,-2.25,                  
    
                              -2,-1.75,-1.5,-1.25,-1,-0.75,-0.5,-0.25,0};//原始函数横坐标
    
        double y[25];//三次样条所得计算值
    
        double z[25];//龙格函数横坐标对应的值
    
        double r[25];//三次样条计算值与原函数值相对误差
    
        int k;
    
        cout<<"请输入本次实验您使用的等距节点的个数：";
    
        cin>>k;
    
        cout<<"请输入已知点的坐标(x,y):"<<endl;
    
        for(int i=0; i<=k; i++)
    
        {
    
            cin>>b[i]>>a[i][0];
    
            if(i>=1)
    
            {
    
                h[i-1]=b[i]-b[i-1];
    
            }
    
        }
    
        for(int i=0; i<=k; i++)
    
        {
    
            if(i==0)
    
            {
    
                c2[0]=1;
    
            }
    
            else if(i==k)
    
            {
    
                c1[k]=1;
    
            }
    
            else
    
            {
    
                c1[i]=h[i-1]/(h[i-1]+h[i]);
    
                c2[i]=h[i]/(h[i-1]+h[i]);
    
            }
    
        }
    
        //2阶均差的计算
    
        for(int j=1; j<3; j++)
    
        {
    
            for(int i=j; i<=k; i++)
    
            {
    
               a[i][j]=(a[i][j-1]-a[i-1][j-1])/(b[i]-b[i-j]);
    
            }
    
        }
    
        for(int i=0;i<=k;i++)
    
        {
    
            if (i==0)
    
                d[i]=(6*(a[1][1]-5/338))/h[i];// 5/338是函数在x0处的导数值
    
            else if(i==k)
    
                 d[i]=(6*(5/338-a[k][1]))/h[i-1];
    
            else d[i]=6*a[2][i];
    
        }
    
        //追赶法解对角占优矩阵
    
        for(int i=0;i<=k;i++)
    
        {
    
            if(i==0)
    
            {
    
               c3[0]=c2[0]/2;     //2是对角占优矩阵对角线的第一个元素
    
               c4[0]=d[0]/2;
    
            }
    
            else if(i<k)
    
             {
    
                  c3[i]=c2[i]/(2-c1[i]*c3[i-1]);
    
                  c4[i]=(d[i]-c1[i]*c4[i-1])/(2-c1[i]*c3[i-1]);
    
             }
    
            else if(i==k)
    
            {
    
                 c4[i]=(d[i]-c1[i]*c4[i-1])/(2-c1[i]*c3[i-1]);
    
                 M[k]=c4[k];
    
            }
    
        }
    
        for(int i=k-1;i>=0;i--)
    
        {
    
            M[i]=c4[i]-c3[i]*M[i+1];
    
        }
    
         //利用所得写出三次样条表达式并计算函数值，并与原函数值进行比较
    
        for(int i=0;i<21;i++)
    
        {
    
            int j;
    
            for(int w=0;w<k;w++)
    
            {
    
              if((x[i]>=b[w])&&(x[i]<b[w+1]))
    
                j=w;
    
            }
    
     y[i]=M[j]*(b[j]-x[i])*(b[j]-x[i])*(b[j]-x[i])/(6*h[j]+M[j+1]*(x[i]-b[j])*(x[i]-b[j])*(x[i]-b[j])/(6*h[j]+(a[j][0]-(M[j]*h[j]*h[j]/6))*((b[j+1]-x[i])/h[j])
    
    +(a[j+1][0]-(M[j+1]*h[j]*h[j]/6))*((x[i]-b[j])/h[j]);
    
             z[i]=1/(x[i]*x[i]+1);
    
             r[i]=fabs((y[i]-z[i])/z[i]);
    
            cout<<"x的值："<<x[i]<<"   理论值："<<z[i]<<"     计算所得："<<y[i]<<"       误差为："<<r[i]<<endl;
    
        }
    
        return 0;
    
    }

[70]: /images/20220523/93d513cfb652460791dba489632fee66.png