.
//三次样条插值 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[N] = {0};
double y[N] = {0};
double M[N] = {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[N] = {2.2, 2.4, 2.6};
double y[N] = {0.5207843, 0.5104147, 0.4813306};
double M[N] = {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[n];
double *l = new double[n];
double *y = new double[n];
//追赶法
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[n] = dn;
c[0] = r0;
b[0] = 2;
a[n] = un;
b[n] = 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页