最小二乘法求多项式的拟合曲线

crv_fit.h

//多项式曲线拟合 f(x)=a0+a1x+a2x^2+a3x^3+...anx^n

class Crv_fit {
public :
	Crv_fit(void);
	void clear(void);
	//~Crv_fit(void);

public :
	void init_ary(int m); //m为样本点个数; 初始化数组ary_x, ary_y; 初始化ary_x_cp,之后和ary_x共同为计算ary_A提供数据
	void cal_ary (int n); //n+1矩阵A的阶数; 计算矩阵中元素,分别存放在ary_A, ary_b中
	void init_matA(void); //初始化系数矩阵(对称方阵)mat_A
	void cal_mat (); //求解非齐次线性方程组Ax=b
	void creat_fit(); //计算拟合曲线数值,并存入数组fit_y,同时将拟合曲线的点集输出到新文件中
	
private :
	int N, M; //M为样本个数, N+1为方阵(对角阵)的阶
	double *ary_x, *ary_y, *ary_x_cp; //ary_x存放样本x, ary_y存放样本y, ary_x_cp,之后和ary_x共同为计算ary_A提供数据
	double *fit_y; //拟合曲线上的y值
	double *ary_A, *ary_b; //ary_A存放方阵A中元素,ary_b存放矩阵b元素
	double **mat_A; //系数矩阵(对称方阵)
	double *result; //结果数组
};

 crv_fit_file.h

#include <iostream>
#include <math.h>
#include <malloc.h>
#include <stdlib.h>

#include "crv_fit.h"

char fileName[20]; //(x,y)样本点集文本文件
char nwFileName[20]; //拟合曲线上(x,y)点集文本文件
////////////////////////////////
Crv_fit::Crv_fit() {
	M=0;	N=0;
	ary_x=NULL; ary_x_cp=NULL;	ary_y=NULL;
	ary_A=NULL; ary_b=NULL;
	mat_A=NULL;
	result=NULL;
}

//////////////////////////////////
void Crv_fit::init_ary(int m) { //初始化数组ary_x, ary_y; 初始化ary_x_cp,之后和ary_x共同为计算ary_A提供数据
	M=m;
	ary_x=(double *)malloc(M*sizeof(double));
	ary_x_cp=(double *)malloc(M*sizeof(double));
	ary_y=(double *)malloc(M*sizeof(double));
 
	FILE *fp;
	char ch;
	if((fp=fopen(fileName,"rt"))==NULL) {                           
		printf("Cannot open file strike any key exit!");
		exit(0);
	}

	ch=fgetc(fp);
	std::string str=""; //这里str是类string的一个对象
	int k=0;
	while(k<M) {
		/*----开始读取x值了-----*/
		while(ch==' ' || ch=='p' || ch=='\n')
			ch=fgetc(fp);
		while(ch!=' ' && ch!='\n' && ch!='p') {
			str+=ch;
			ch=fgetc(fp);
		}
		double aa=atof(str.c_str()); //将字符串转换成double型
		*(ary_x+k) = aa; //文本中的一行为矩阵valAry中的一行
		str="";
		/*----读取x结束-----*/
		/*----开始读取y值了-----*/
		while(ch==' ' || ch=='p' || ch=='\n')
			ch=fgetc(fp);
		while(ch!=' ' && ch!='\n' && ch!='p') {
			str+=ch;
			ch=fgetc(fp);
		}
		aa=atof(str.c_str()); //将字符串转换成double型
		*(ary_y+k) = aa; //文本中的一行为矩阵valAry中的一行
		str="";
		/*----读取x结束-----*/
		*(ary_x_cp+k)=1.0;
		ch=fgetc(fp);
		k++;
	}
	std::cout << "read all the vectors" << '\n';
	fclose(fp);
}

/////////////////////////////////////////////////
//right!!!!!
void Crv_fit::cal_ary(int n) { //计算组成系数数组ary_A和结果矩阵ary_b的元素
	N=n;
	ary_A=(double *)malloc((2*N+1)*sizeof(double));
	ary_b=(double *)malloc((N+1)*sizeof(double));
	std::cout << "calculate the ary_A" << '\n';
	*(ary_A+0)=M;
	
	*(ary_b+0)=0;
	for(int i=0; i<M; i++)
		*(ary_b+0)+=*(ary_y+i);
	std::cout << "y的和为 " << *(ary_b+0) << '\n'; 
	for(i=1; i<=N; i++) {
		double sumA=0.0;
		double sumb=0.0;
		for(int j=0; j<M; j++) {
			*(ary_x_cp+j) *= *(ary_x+j);
			sumA+=*(ary_x_cp+j);

			sumb+=*(ary_x_cp+j)* (*(ary_y+j));
		}
		*(ary_A+i)=sumA;
		*(ary_b+i)=sumb;
	}
	std::cout << "x的和为 " << *(ary_A+1) << '\n';
	for(i=N+1; i<=2*N; i++) {
		double sumA=0.0;
		for(int j=0; j<M; j++) {
			*(ary_x_cp+j) *= *(ary_x+j);
			sumA+=*(ary_x_cp+j);
		}
		*(ary_A+i)=sumA;
	}
}

//////////////////////////////////////////////////////////
//right!!!!!!
void Crv_fit::init_matA() { //初始化系数矩阵(对称方阵)mat_A (N+1 阶哦!!)
	mat_A=(double **)malloc((N+1)*sizeof(double));
	for(int i=0; i<=N; i++)
		*(mat_A+i)=(double *)malloc((N+1)*sizeof(double));
	for(i=0; i<=N; i++) {
		int bg=i;
		for(int j=0; j<=i; j++) {
			*(*(mat_A+i)+j)=*(ary_A+bg);
			if(j!=i)
				*(*(mat_A+j)+i)=*(*(mat_A+i)+j);
			bg++;
		}
	}
} 

//////////////////////////////////////////////////////////
//求解非齐次线性方程组Ax=b
void Crv_fit::cal_mat() { //求解矩阵Ax=b (A为 N+1 阶对称方阵, b 为N+1 维列向量)

	double *temp; //当该行主元为0时需要与主元列元素不为0的行交换, 临时数组temp用于交换
	temp = (double *)malloc((N+1)*sizeof(double));

	/*----高斯消元游戏开始了----*/

	for(int i=0; i<N; i++) { //需要对前 N 行进行消元
		double keyData = *(*(mat_A+i)+i); //记录当前行的主元
		int keyRow = i; //记录当前主元行
		
		/*---------当前主元为0,需要找到新的主元行进行交换--------*/
		if(keyData == 0) { 
			for(int k=i+1; k<=N; k++)
				if(*(*(mat_A+k)+i) != 0) {
					keyRow = k;
					break; //找到当前主元列中当前元素不为0的行 k ,用于交换使其成成为真正的主元行
				}

			/*----交换两行-----------------------------*/
			for(int j=i; j<=N; j++) {
				*(temp+j) = *(*(mat_A+i)+j);
				*(*(mat_A+i)+j) = *(*(mat_A+keyRow)+j);
				*(*(mat_A+keyRow)+j) = *(temp+j);
			}
			/*----交换两行-----------------------------*/
			//交换后主元行元素更新, 故要重置keyRow,keyData
			keyRow = i; 
			keyData = *(*(mat_A+i)+i); 
		}
		/*---------当前主元为0,需要找到新的主元行进行交换--------*/

		/*------开始对当前行进行高斯消元----------*/
		//将主元列中主元以下的所有元素置0
		for(int k=i+1; k<=N; k++) { //i 为主元行行号, 主元以下消元, 故从 i+1 开始
			double elm = *(*(mat_A+k)+i)/keyData; //主元实际就是对角线元素, 故主元所在列等于行值, 故主元列为 i; elm是用于消元的除式因子
			/*----本行所有元素进行消元相关变化---*/
			for(int j=i+1; j<=N; j++) {
				*(*(mat_A+k)+j) -= *(*(mat_A+i)+j)*elm; //本行元素减去主元行元素乘以除式因子elm变为新的本行元素; i 是主元行的行号
			}
			*(ary_b+k) -= *(ary_b+i)*elm; //右值列向量矩阵也要发生相应变化 
			/*----本行所有元素进行消元相关变化---*/
		}
		/*------开始对当前行进行高斯消元----------*/
	}
	
	/*----高斯消元游戏结束了----*/

	/*------开始解方程喽----*/
	result = (double *)malloc((N+1)*sizeof(double));
	*(result+N) = (*(ary_b+N))/(*(*(mat_A+N)+N));
	int k=N-1; //标记mat_b的下标
	for(i=N-1; i>=0; i--) { //回代求解
		double sum = 0.0;
		for(int j=N; j>i; j--) 
			sum += *(*(mat_A+i)+j)* (*(result+j));
		*(result+i) = (*(ary_b+k)-sum)/(*(*(mat_A+i)+i));  
		k--;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            
	}
	std::cout << "out put the result\n";
	for(i=0; i<=N; i++)
		std::cout << "a" << i << "=" << *(result+i) << '\n';
	std::cout << '\n';
	/*------Ax=b解出来了----*/
}

/////////////////////////////////////////////////////////////
//计算拟合曲线上的纵坐标值
void Crv_fit::creat_fit() {
	fit_y = (double *)malloc(M*sizeof(double));
	/*----开始计算多项式的和----*/
	for(int i=0; i<M; i++) { //共M个样本点
		double sum=*(result+0);
		double t_x=1.0;
		for(int k=1; k<=N; k++) {//k为多项式的最高幂
			t_x *= *(ary_x+i);
			sum += *(result+k)*t_x; //sum+=an*x^n
		}
		*(fit_y+i) = sum;
	}
	/*----计算结束----*/

	/*----将拟合曲线的点集写入新文本文件中----*/
	FILE *fp;
	if((fp=fopen(nwFileName,"w"))==NULL) {
		printf("Cannot build file strike any key exit!");
		exit(0);
	}
	for(i=0; i<M; i++) {
		fprintf(fp,"%s"," v  ");
		fprintf(fp,"%f%s%f",*(ary_x+i),"  ",*(fit_y+i));
		fprintf(fp,"%s","\n");
	}	
	fclose(fp);
	std::cout << "write all the vectors" << '\n';
	/*----将拟合曲线的点集写入新文本文件中----*/
}

 crv_fit.cpp

#include <iostream>
#include "crv_fit_file.h"

///////////////////////////////////////////////////////////////////////////////
int main() {
	Crv_fit crv_fit;
	//crv_fit.Crv_fit();
	std::cout << "crvfilea.txt: sapNum=12; nNum=9\n" << "crvfileb.txt: sapNum=10; nNum=1\n";
	std::cout << "read the original vectors" << '\n';
	std::cout << "input the fileName \n";
	std::cin >> fileName;
	int sapNum;
	int nNum;
	std::cout << "sapNum= ";
	std::cin >> sapNum;
	std::cout << "nNum= ";
	std::cin >> nNum;
	crv_fit.init_ary(sapNum); 
	crv_fit.cal_ary(nNum); 
	crv_fit.init_matA();
	crv_fit.cal_mat();
	std::cout << "input the nwFileName \n";
	std::cin >> nwFileName;
	crv_fit.creat_fit();
	return 0;
}

///////////////////////////////////////////////////////////////////////////////

 验证

样本点集文件

拟合曲线上的点集文件

调试了好长时间也没能将曲线显示出来，先放放以后再说吧~