中南林业科技大学计算方法实验-插值和曲线拟合

拉格朗日差值

#include<iostream>

int main(){
	double fun(double xnum,double x[],double y[],int size);
	double x[4]={-2,0,4,5};
	double y[4]={5,1,-3,1};
	double ans1=fun(-2,x,y,4);
	double ans2=fun(0,x,y,4);
	printf("%lf\t%lf",ans1,ans2);
	
	//分段三次插值 
	/*double x1[3]={0.8,0.95,1.05};
	double y1[3]={0.90,1.00,1.25382};
	double ans1=fun(0.569,x,y,3);
	double ans2=fun(0.99,x1,y1,3);
	printf("%lf\t%lf",ans1,ans2);*/
	
	/*double x2[7]={1,2,3,4,5,6,7};
	double y2[7]={0.368,0.135,0.050,0.018,0.007,0.002,0.001};
	double ans1=fun(1.8,x,y,7);
	double ans2=fun(6.15,x,y,7);
	printf("%lf\t%lf",ans1,ans2);*/
	
}

double fun(double xnum,double x[],double y[],int size){
	double i=0;
	for(int m=0;m<=size-1;m++){
		double x1=1;
		double x2=1;
		for(int n=0;n<=size-1;n++){
			//printf("%d\t%d\n",x1,x2);
			if(n==m) continue;
			x1*=xnum-x[n];
			x2*=x[m]-x[n];
		}
		i+=y[m]*x1/x2;
	}
	return i; 
}

牛顿差值

#include<stdio.h>

int main(){
	double x[6]={0.4,0.55,0.65,0.8,0.95,1.05};
	double y[6]={0.41075,0.57815,0.69675,0.90,1.00,1.25382};
	
	/*double x[6]={1,2,3,4,5,6};
	double y[6]={0.368,0.135,0.050,0.018,0.007,0.002};*/

	double f[6][6];
	double xx=0.596;
	for(int i=0;i<=5;i++){
		f[i][0]=y[i];
	}
	for(int j=1;j<=5;j++){
		for(int i=j;i<=5;i++){
			f[i][j]=(f[i][j-1]-f[i-1][j-1])/(x[i]-x[i-j]);
		}
	}
	for(int i=0;i<=5;i++){
		printf("%lf\t",f[i][i]);
	}
	printf("\n");
	
	double ans=0;
	double ss=1;
	for(int i=0;i<=5;i++){
		for(int j=0;j<i;j++){
			ss*=xx-x[j];
		}
		ans+=f[i][i]*ss;
	}
	printf("%lf",ans);
}

曲线拟合

#include<stdio.h>
#include<math.h>

int main(){
	double fun(double a[],double b[],int size);
	double x[12]={0,5,10,15,20,25,30,35,40,45,50,55};
	double y[12]={0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.02,4.64};
	int size=11;
	int nhcs=4;		//拟合次数 

		double p=0,a[12],fai[12][12],f[12]; 
		a[0]=0;f[1]=0;
		for(int i=1;i<=size;i++){
			f[1]+=x[i];
			a[0]+=y[i];
			fai[0][i]=1;
		}
		f[1]=f[1]/size;
		a[0]=a[0]/size;
		
		for(int i=1;i<=size;i++){
			fai[1][i]=x[i]-f[1];
		}
		a[1]=fun(fai[1],y,size)/fun(fai[1],fai[1],size);
		
		double xfai[12],beida[12];
		for(int i=1;i<nhcs;i++){
			for(int k=1;k<=size;k++){
				xfai[k]=x[k]*fai[i][k];
			}
			f[i+1]=fun(fai[i],xfai,size)/fun(fai[i],fai[i],size);
			beida[i]=fun(fai[i],fai[i],size)/fun(fai[i-1],fai[i-1],size);
			for(int k=1;k<=size;k++){
				fai[i+1][k]=(x[k]-f[i+1])*fai[i][k]-beida[i]*fai[i-1][i];
			}
			a[i+1]=fun(fai[i+1],y,size)/fun(fai[i+1],fai[i+1],size);
		}
		
		printf("%d次拟合的结果为\n",nhcs);
		double sum=0;
		for(int i=0;i<=nhcs;i++){
			printf("a[%d]=%f",i,a[i]); 
			sum+=a[i]*pow(10,i);
			printf("\n%lf",sum);
		}
		
		for(int i=1;i<=nhcs;i++){
			printf("f[%d]=%f\n",i,f[i]);
		}
		for(int i=1;i<nhcs;i++){
			printf("beida[%d]=%f\n",i,beida[i]);
		}
		
		
	
} 


double fun(double a[],double b[],int size){
	int m=0;
	for(int i=1;i<size;i++){
		m+=a[i]*b[i];
	}
	return m;
}

最小二乘法

#include<stdio.h>
#include<math.h>

int main(){
//	void fun1(double a[],double b[],int num);
	double xx[6]={0,1,2,3,4,5};
	double yy[6]={5,2,1,1,2,3};
	double aa[3][3],bb[3];
	aa[0][0]=6;
	for(int i=0;i<3;i++){
		for(int j=0;j<3;j++){
			if(i==0&&j==0) continue;
			double sum=0;
			for(int n=0;n<6;n++){
				sum+=pow(xx[n],i+j); 
			}
			aa[i][j]=sum;
		}
	}
	
	for(int i=0;i<3;i++){
		double sum=0;
		for(int j=0;j<6;j++){
			sum+=pow(xx[j],i)*yy[j];
		}
		bb[i]=sum;
	}
	
//	fun1(a[3],b,3);
	
	double a[4][4],b[4];
	for(int i=0;i<=3;i++){
		b[i+1]=bb[i];
		for(int j=0;j<=3;j++){
			a[i+1][j+1]=aa[i][j];
		}
	}
	
	double x[4]={0},y[4]={0};
	double exp=0.0001;
	double s=1;
	while(s>exp){
		s=0;
		for(int i=1;i<=3;i++) y[i]=x[i];
		for(int i=1;i<=3;i++){
			double sum=0;
			for(int j=1;j<=i-1;j++){
				sum+=a[i][j]*x[j];
			}
			for(int j=i+1;j<=3;j++){
				sum+=a[i][j]*y[j];
			}
			x[i]=(b[i]-sum)/a[i][i];
			printf("%lf\t",x[i]);
		}
		printf("\n");
		for(int i=1;i<=3;i++)if(fabs(x[i]-y[i])>s) s=fabs(x[i]-y[i]);
	}
	
	printf("解：\n");
	for(int m=1;m<=3;m++){
		printf("%.5lf\t",x[m]);
	}
	
/*	for(int i=0;i<6;i++){
		double ss=0;
		for(int j=1;j<=3;j++){
			ss+=pow(x[i],j-1)*x[j];
		}
		printf("\n%lf",ss);
	}*/
}