欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

C语言实现三次样条插值

程序员文章站 2023-12-27 21:58:27
...

《数值分析》实验习题1

已知函数值如下表:


C语言实现三次样条插值

试用三次样条插值求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;
}


上一篇:

下一篇: