三次样条插值算法

  这是三次样条插值计算的简单算法代码，记录一下

  

#include<iostream>
#include<math.h>
using namespace std;
float Hermite(int n, float xx)
{
    float *x = new float[n+1],*ar = new float[n+1];
    float *y = new float[n+1],*br = new float[n+1];
    float *h = new float[n],m0,m1;
    float *a = new float[n+1], *b = new float[n+1];
    float *m = new float[n+2],fx;
    int i,choice=0;
    for(i=0; i<n+1; i++)//输入xi，yi
    {
        cout<<"input x"<<i<<":";
        cin>>x[i];
        cout<<"input y"<<i<<":";
        cin>>y[i];
    }
    for(i=0; i<n; i++)//计算hi
    {
        h[i] = x[i+1] - x[i];
    }
    for(i=1; i<n; i++)//i = 1,2,...n 计算α，β
    {
        ar[i] = h[i-1] / (h[i-1] + h[i]);
        br[i] = 3 * ((1-ar[i])/h[i-1]*(y[i]-y[i-1])+ar[i]/h[i]*(y[i+1]-y[i]));
    }
    CHOICE:
    cout<<"输入边界条件：
1.第一种边界条件.
2.第二种边界条件.
"<<endl;
    cin>>choice;//选择边界条件
    switch (choice)
    {
    case 1:
        cout<<"input M0:
";
        cin>>m0;
        cout<<"input Mn:
";
        cin>>m1;
        ar[0] = 0;        //第一种边界条件
        br[0] = 2 * m0;
        ar[n] = 1;
        br[n] = 2 * m1;
        break;
    case 2:
        ar[0] = 1;       //第二种边界条件
        br[0] = 3 / h[0] * (y[1] - y[0]);
        ar[n] = 0;
        br[n] = 3 / h[n-1] * (y[n] - y[n-1]);
        break;
    default:
        cout<<"input error!"<<endl;
        goto CHOICE;
        break;
    }
    //计算a,b
    a[0] = -ar[0]/2;
    b[0] = br[0]/2;
    for(i=1; i<n+1; i++)
    {
        a[i] = -ar[i]/(2+(1-ar[i])*a[i-1]);
        b[i] = (br[i]-(1-ar[i])*b[i-1])/(2+(1-ar[i])*a[i-1]);
    }
    //计算mi
    m[n+1] = 0.0;
    for(i=n; i>-1; i--)
    {
        m[i] = a[i] * m[i+1] + b[i];
    }
    if (xx < x[0] || xx > x[n])
    {
        cout<<"输入的数不在x取值范围内!"<<endl;
        return 0.0;
    }
      for(i=0; i<n+1; i++)//输出mi
    {
        cout<<"m"<<i<<"="<<m[i]<<" ";
        cout<<endl;
    }
    for(i=0; i<n; i++)
    {
        if(xx >= x[i] && xx < x[i+1])break;
    }
    //计算结果
    fx = (1+2*(xx-x[i])/(x[i+1]-x[i]))*pow(((xx-x[i+1])/(x[i]-x[i+1])),2)*y[i] +
         (1+2*(xx-x[i+1])/(x[i]-x[i+1]))*pow(((xx-x[i])/(x[i+1]-x[i])),2)*y[i+1] +
         (xx-x[i])*pow(((xx-x[i+1])/(x[i]-x[i+1])),2)*m[i] +
         (xx-x[i+1])*pow(((xx-x[i])/(x[i+1]-x[i])),2)*m[i+1];
    cout<<"The answer is:
"<<"s("<<xx<<")="<<fx<<endl;    //输出结果
    return fx;
}
int main()
{
    int n;
    float xx;
    cout<<"输入点的个数:
";
    cin>>n;
    cout<<"输入x的值:
";
    cin>>xx;
    Hermite(n-1,xx);
}