最小二乘法求多項式的擬合曲線

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;
}

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

 驗證

樣本點集文件

擬合曲線上的點集文件

調試了好長時間也沒能將曲線顯示出來，先放放以后再說吧~