三次样条插值 SPLINE

.
//三次样条插值 SPLINE

#include <iostream>
#include <complex>//exp(const complex<Type>& _ComplexNum)函数
#include <math.h>//使用三角函数及log函数

void TSS(double a[], double b[], double c[], double d[], double x[], const int n);
//本算法解三对角线性方程组，系数存于数组a、b和c中，d是右端向量，计算结果存于x中，n为数组长度
void FINDK(double x[], int n, double x_aim, int& K);
//本算法找出x_aim所在区间(xk-1<=x<=xk),以xk的下标值K作为结果
void SPLINE(double x[], double y[], double M[], const int n, double r0, double d0, double un, double dn);
//本算法计算三次样条的参数Mi，存放于数组M中，参数r0,d0,un,dn根据不同边界条件给定;数组x,y,M的长度为n+1

 

int main()
{
 using namespace std;

 const int N = 6;
 double x
= {0};
 double y
= {0};
 double M
= {0};

 for(int i = 0; i < N; i++)
 {
  x[i] = -1.0 + 2.0 * i / ( N - 1);
  y[i] = 1.0 / ( 1.0 + 25.0 * x[i] * x[i]);
 }

 double x_input = -1 + 2.0 * 1 / 100;
 int K = 0;

 //确定边界条件
 double r0 = 0;
 double un = 0;
 double d0 = 0;
 double dn = 0;

/*
 const int N = 3;
 double x
= {2.2, 2.4, 2.6};
 double y
= {0.5207843, 0.5104147, 0.4813306};
 double M
= {0};

 double x_input = 2.5;
 int K = 0;

 //确定边界条件
 double r0 = 1;
 double un = 1;
 double d0 = 6 * ( 0.5104147 - 0.5207843 ) / ( 2.4 - 2.2 ) /  ( 2.4 - 2.2 );
 double dn = -6 * ( 0.4813306 - 0.5104147 ) / ( 2.6 - 2.4 ) / ( 2.6 - 2.4 );
*/

 SPLINE(x, y, M, N-1, r0, d0, un, dn);

 //输出结果
 for(int i = 0; i < N; i++)
  cout << "M[" << i << "] = " << M[i] << endl;

 FINDK(x, N, x_input, K);
 cout << "/nK = " << K << endl;

 double h_temp = x[K] - x[K-1];
 double x_r = x[K] - x_input;
 double x_l = x_input - x[K-1];
 double y_output = ( M[K-1] * pow(x_r, 3) / 6 + M[K] * pow(x_l, 3) / 6 + ( y[K-1] - M[K-1] * pow(h_temp, 2) / 6 ) * x_r + ( y[K] - M[K] * pow(h_temp, 2) / 6 ) * x_l ) / h_temp;

 cout << "插值结果为 ：" << y_output << endl;

 return 0;
}

void TSS(double a[], double b[], double c[], double d[], double x[], const int n)
{
 //申请动态数组
 double *u = new double
;
 double *l = new double
;
 double *y = new double
;

 //追赶法
 u[0] = b[0];
 y[0] = d[0];
 for(int k = 1; k < n; k++)
 {
  l[k] = a[k] / u[k-1];
  u[k] = b[k] - l[k] * c[k-1];
  y[k] = d[k] - l[k] * y[k-1];
 }
 x[n-1] = y[n-1] / u[n-1];
 for(int k = n-2; k >= 0; k--)
 {
  x[k] = (y[k] - c[k] * x[k+1]) / u[k];
 }
 
 //释放空间
 delete []u;
 delete []l;
 delete []y;
}

void FINDK(double x[], int n, double x_aim, int& K)
{
 K = 1;
 for(int i = 1; i < n; i++)
 {
  if(x_aim <= x[i])
  {
   K = i;
   break;
  }
  else
   K = i + 1;
 }

//插值点大于x[n-1]时，令K取n-1
 if(x_aim > x[n-1]) K = n-1;
}

void SPLINE(double x[], double y[], double M[], const int n, double r0, double d0, double un, double dn)
{

 //申请动态数组，实际上h,a,c只要用n个长度，b要用n+1个长度
 double *h = new double[n+1];
 double *a = new double[n+1];
 double *b = new double[n+1];
 double *c = new double[n+1];

 //1
 for(int i = 0; i <= n; i++)//此处要取到n，把y中所有值赋到M中
  M[i] = y[i];
 //2 求y的二阶差商，即二阶导数
 for(int k = 1; k <= 2; k++)
  for(int i = n; i >= k; i--)//计算二阶差商，存储到M2至Mn中
   M[i] = ( M[i] - M[i-1] ) / ( x[i] - x[i-k] );
 //3
 h[1] = x[1] - x[0];
 //4
 for(int i = 1; i < n; i++)
 {
  h[i+1] = x[i+1] - x[i];
  c[i] = h[i+1] / ( h[i] + h[i+1] );
  a[i] = 1 - c[i];
  b[i] = 2;
  M[i] = 6 * M[i+1];
 }
 //5
 M[0] = d0;
 M
= dn;
 c[0] = r0;
 b[0] = 2;
 a
= un;
 b
= 2;
 //6调用TSS(a, b, c, M, M, n+1)计算M
 TSS(a, b, c, M, M, n+1);

 //释放申请的内存
 delete []h;
 delete []a;
 delete []b;
 delete []c;

 //test
 std::cout <<"TEST: Can SPLINE function be executed?/n/n";
}

 

 

算法原理见《计算方法教程（第二版）》凌永祥、陈明奎编，西安交通大学出版社，第110页