C语言实现三次样条插值

《数值分析》实验习题1

已知函数值如下表：

试用三次样条插值求f(5.472)的近似值。

/*****************************************************************《数值分析》实验习题1机械工程 现代制造技术教育部重点实验室*****************************************************************/#include<stdio.h>#include<math.h>double * cal_matrix(int , double *, double *, double *);int cpm_x(double , double *, int);int main(){ int n = 10, i;//printf("请输入数据个数：\n");//scanf("%d", &n);double h[n], u[n-1], d[n+2], *p, test = 5.472, result, M[n+1];double x[n] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};double fx[n+2] = {0, 0.69314718, 1.0986123, 1.3862944, 1.6094378, 1.7917595, 1.9459101, 2.079445, 2.1972246, 2.3025851, 1, 0.1};//printf("请输入x的值：\n");//for(i=0; i<n; i++)//scanf("%lf", &x[i]);//printf("请输入对应的fx的值：\n");//for(i=0; i<n; i++)//scanf("%lf", &fx[i]);//printf("补充边界条件：\n");////将边界条件的数值一并存入x数组 //scanf("%lf %lf", &fx[n], &fx[n+1]);//printf("请输入待计算的X0的值:");//scanf("%lf", &test); for(i=1; i<n; i++)h[i] =  x[i] - x[i-1];for(i=1; i<n-1; i++)u[i] = h[i] / (h[i] + h[i+1]);for(i=1; i<n-1; i++)d[i] = 6 / (h[i] + h[i+1]) * ((fx[i+1] - fx[i]) / (x[i+1] - x[i]) - (fx[i] - fx[i-1]) / (x[i] - x[i-1]));d[0] = 6 / h[1] * ((fx[1] - fx[0]) / (x[1] - x[0]) - fx[n]);d[n-1] = 6 / h[n-1] * (fx[n+1] - (fx[n-1] - fx[n-2]) / (x[n-1] - x[n-2]));p = cal_matrix(n, u, d, M); //求三弯矩方程中的M i = cpm_x(test, x, n); //取x0所在的区间 result = pow((x[i+1]-test), 3) * p[i+1] / 6 /h[i+1] +pow((test -x[i]), 3) * M[i+2] / 6 / h[i+1] + (fx[i] - M[i+1] * pow(h[i+1], 2) / 6) * (x[i+1] - test) / h[i+1] +(fx[i+1] - M[i+2] * pow(h[i+1], 2) / 6) * (test - x[i]) / h[i+1]; printf("result: %.8lf", result) ;return 0;}double * cal_matrix(int n, double *temp_a, double *temp_b, double *x){//用追赶法解三对角方程组 int i;double a[n+1], b[n+1], c[n], l[n+1], r[n+1], z[n+1];//处理下标不一致的问题 for(i=1; i<n-1; i++)a[i+1] = temp_a[i];a[n] = 1;c[1] = 1;for(i=1; i<n-1; i++)c[i+1] = 1 - temp_a[i];for(i=0; i<n; i++)b[i+1] = temp_b[i];//LU分解r[1] = 2;for(i=2; i<n+1; i++){l[i] = a[i] / r[i-1];r[i] = 2 - l[i] * c[i-1];};//解Lz=bz[1] = b[1];for(i=2; i<n+1; i++)z[i] = b[i] - l[i] * z[i-1];//解Ux=zx[n] = z[n] / r[n];for(i=n-1; i>0; i--)x[i] = (z[i] - c[i] * x[i-1]) / r[i];return x; }int cpm_x(double x0, double *x, int n){int i;for(i=0; i<n; i++){if (x0 < x[i])return i-1;else continue;}}