差值与逼近（实验二）

爱被打了一巴掌 2022-12-29 14:24 25阅读 0赞

p.s 文末附试验测试工程下载链接

### 拉格朗日插值法（封装函数） ###

###### 函数接口 ######

double LagrangeInterPol(double arrX[],double arrY[], int n, double x)

其中arrX和arrY表示待插值数据点，n表示数据点个数，x表示待求点，return返回所求点的值。

###### 实现代码 ######

#include "LagrangeInterPol.h"
    #include "vetMtxBase.h"
    double LagrangeInterPol(double arrX[], double arrY[], int n, double x) 
    { 
    	double yResult = 0.0;
    	double *LValue = newVector(n);         //LValue[n]存放的是每次求解的插值基函数的在x处的值 
    	int k, m;	                           //循环上下标变量k,m
    	double temp1, temp2;    	           //插值基函数中的上下累乘temp1,temp2
    	//计算各基函数的值
    	for (k = 0; k < n; k++)
    	{ 
    		temp1=1.0;
    		temp2=1.0;
    		for(m=0;m<n;m++)
    		{ 
    			if(m==k)
    				continue;
    			temp1*=(x-arrX[m]);
    			temp2*=(arrX[k]-arrX[m]);
    		}
    		LValue[k]=temp1/temp2;
    	}
    	//基函数值与y各元素相乘再求和
    	for (int i = 0; i < n; i++)
    	{ 
     		yResult+=arrY[i]*LValue[i];
    	}
    	deleteVector(LValue); 
    	return yResult;
    }

### 牛顿插值法（封装函数） ###

###### 函数接口 ######

double NewtonInterPol(double arrX[], double arrY[], int n, double x)

其中arrX和arrY表示待插值数据点，n表示数据点个数，x表示待求点，return返回所求点的值。

###### 实现代码 ######

#include "NewtonInterPol.h"
    #include "vetMtxBase.h"
    double NewtonInterPol(double arrX[], double arrY[], int n, double x)
    { 
    	int i;
    	double  **table = newMtx(n+1,n+1);           //这里是n+1*n+1，存放Newton插值的差商表 
    	for ( i = 0; i <= n; i++) 
    	{  
    		table[0][i] = arrY[i];               //第一列放y值 
    	} 
    	//计算各阶差商
    	for ( i = 1; i <= n; i++) 
    	{ 
    		for(int j=i;j<=n;j++)
    		{ 
    			table[i][j]=(table[i-1][j]-table[i-1][j-1])/(arrX[j]-arrX[j-i]);
    		}
    	}
    	//根据Newton插值公式根据差商与各项的乘积，再求和，获得插值点的值
    	double res = 0.0;
    	for ( i = 0; i < n; i++)
    	{ 
    		double temp=table[i][i];
    		for(int j=0;j<i;j++)
    		{ 
    			temp*=x-arrX[j];
    		}
    		res+=temp;
    	}
    	deleteMtx(table,n);
    	return res;
    }

### 多项式拟合最小二乘法（封装函数） ###

###### 函数接口 ######

double Polyfit(double arrX[],double arrY[], int order, int n, double x, double Coeffi[])

其中arrX和arrY表示待插值数据点，order表示多项式阶次，n表示数据点个数，x表示待求点，Coeffi表示多项式系数，return返回所求点的值。

###### 实现代码 ######

#include "Polyfit.h"
    #include "vetMtxBase.h" 
    #include "SolverEqGauss.h"
    double Polyfit(double arrX[],double arrY[], int order, int n, double x, double Coeffi[])
    { 
    	double **A = newMtx(order+1,order+1);   //待求得线性方程组系数矩阵分配内存
    	double  *b = newVector(order + 1);      //待求得线性方程组常向量分配内存
    	int i,j,k; 
    	//求系数矩阵A
    	for ( i = 0; i <= order; i++)
    	{ 
    		for ( j = 0; j <= order; j++)
    		{ 
    			//核心代码
    			A[i][j]=0;
    			for(int k=0;k<n;k++)
    			{ 
    				A[i][j]=A[i][j]+pow(arrX[k],i+j);
    			}
    		}
    	}
    	//求方程右边的常数向量b
    	for ( i = 0; i <= order; i++)
    	{ 
    		b[i] = 0;
    		for ( k = 0; k < n; k++)
    		{ 
    			//学生需要写的核心代码
    			b[i]+=arrY[k]*pow(arrX[k],i);
    		}
    	}
    	SolverEqGauss(A, b, order + 1, 1.0e-10, Coeffi);     //调用高斯列主元消去法解法方程组(线性方程组)
    	double res = 0;                                      //用拟合后的多项式计算校验点数据 
    	for (int i=0;i<order+1;i++)
    	{ 
    		//核心代码
    		res+=Coeffi[i]*pow(x,i);
    	}
    	//释放内存
    	deleteMtx(A,order+1);
    	deleteVector(b);
    	return res;
    }

//附录实现上面功能的高斯列主元消去法函数(封装函数)
    //函数接口：void SolverEqCholesky(const double* const *A, const double *b, int n, double *x, double **L);
    #include "SolverEqGauss.h"
    #include "vetMtxBase.h"
    void SolverEqGauss(const double* const *A, const double *b, int n,double eps, double *x)
    { 
    	int i, j, k;
    	double **a    = newMtx(n,n+1);   //这里是n*n+1，因为要将b向量存入一起处理
    	double *btemp = newVector(n);    //临时的b向量存储，便于回代时处理
    	double d;                        //用于存放行乘子
    	double temp;                     //用于元素交换的临时变量
    	//赋值给内部矩阵与向量
    	for (i = 0; i < n; i++)
    	{ 
    		for (j = 0; j < n; j++)
    		{ 
    			a[i][j] = A[i][j];
    		}
    		a[i][n] = b[i];
    	}
    	if ( n <= 0)
    	{ 
    		cout << "输入的格式有误！\n";
    	}
    	for (k = 0; k < n - 1; k++) //找列主元最大值
    	{ 
    		double maxNum = 0;                    //存储主元素，即绝对值最大的元素
    		int row = k;                          //记录主元素所在的行
    		for (i = k; i < n; i++)
    		{ 
    			if (fabs(a[i][k]) > maxNum)
    			{ 
    				maxNum = fabs(a[i][k]);
    				row = i;
    			}
    		}
    		if (fabs(a[row][k]) < eps)
    		{ 
    			cout << "无法计算" << endl;
    			deleteMtx(a,n);
    			deleteVector(btemp);
    			return;
    		}
    		if (k != row) //换行
    		{ 
    			for (i = 0; i < n + 1; i++)
    			{ 
    				temp = a[k][i];
    				a[k][i] = a[row][i];
    				a[row][i] = temp;
    			}
    		}
    		for (i = k + 1; i < n; i++) //消元
    		{ 
    			d = a[i][k] / a[k][k];                   //计算行乘子
    			for (j = k; j < n + 1; j++)
    			{ 
    				a[i][j] = a[i][j] - d * a[k][j];    //更新i行的各元素
    			}
    		}
    	}
    	for (i = n - 1; i >= 0; i--)    //回代
    	{ 
    		d = 0;
    		for (j = i+1; j < n; j++)
    		{ 
    			d = d + a[i][j]*x[j];
    		}
    		x[i] = (a[i][n] - d) / a[i][i];
    	}
    	deleteMtx(a,n);
    	deleteVector(btemp);
    }

[实验二：插值与逼近测试用项目c++文件][c]

[c]: https://download.csdn.net/download/huangziguang/13763328