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=b
	z[1] = b[1];
	for(i=2; i<n+1; i++)
		z[i] = b[i] - l[i] * z[i-1];
	//解Ux=z
	x[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;
	}
}

/*****************************************************************
《数值分析》实验习题2
机械工程 
现代制造技术教育部重点实验室
*****************************************************************/

#include<stdio.h>
#include<malloc.h>
#include<math.h>
#define PI 3.14159265

typedef struct node  
{
	int num;
    double T;
    double S;
    double C;
    double R;
    struct node *next; 
}ParaList;

double fun(double );
double calT(ParaList* ,double ,double );
double calS(ParaList*);
double calC(ParaList*);
double calR(ParaList*);

int main(){
	double a = 0, b = 3, eps = 0.00001;
	
	//第二题参数 
//	double a = 1, b = 3, eps = 0.00001;

//	printf("请输入积分下限\n");
//	scanf("%lf", &a);
//	printf("请输入积分上限\n");
//	scanf("%lf", &b);
//	printf("请输入eps");
//	scanf("%lf", &eps);
	//算前四行的T，S，C的值并存入参数链表中 
	ParaList *phead, *p;
	int i;
	phead = (ParaList*)malloc(sizeof(ParaList));
	phead->num = 1;
	phead->T = (b - a) / 2 * (fun(a) + fun(b));
	p= phead;
	for(i=2; ;i++){
		p->next = (ParaList*)malloc(sizeof(ParaList));
		p->next->num = i;
		p->next->T = calT(p, a, b);
		p->next->S = calS(p);
		p->next->C = calC(p);
		p->next->R = calR(p);
		if(i>=5 && (fabs(p->next->R-p->R)<eps)){
			printf("计算结果为%.8lf，共进行了%d次迭代", p->next->R, p->next->num);
			break;
		}
		p = p->next;
	} 
}

double fun(double x){
	//第一题函数 
	double fresult;
	fresult = x * pow(1 + pow(x, 2), 0.5);
	return fresult;
}

//double fun(double x){
//	//第二题函数 
//	double fresult;
//	fresult = pow(3, x)*pow(x, 1.4)*(5*x+7)*sin(pow(x, 2));
//	return fresult;
//}

double calT(ParaList* s, double a, double b){
	//传入某一行的指针和积分上下限，用于计算下一行的T值 
	int i;
	double sum = 0, tresult;
	for(i=0; i<pow(2, s->num -1); i++)
		sum += fun(a + (i + 0.5) * (b - a) / pow(2, s->num -1));
	tresult = 0.5 * (s->T + (b - a) / pow(2, s->num -1) * sum);
	return tresult;
}

double calS(ParaList* t){
	//传入某一行的指针，用于计算下一行的S值 
	double sresult;
	sresult =(4 * t->next->T - t->T) / 3;
	return sresult; 
} 

double calC(ParaList* s){
	//传入某一行的指针，用于计算下一行的C值 
	double cresult;
	cresult = (16 * s->next->S - s->S) / 15;
	return cresult;
}

double calR(ParaList* c){
	//传入某一行的指针，用于计算下一行的R值 
	double rresult;
	rresult = (64 * c->next->C - c->C) / 63;
	return rresult;
}