本文共 9403 字，大约阅读时间需要 31 分钟。

一，原理

二，根据原理得到的  python 源码。

import matplotlib.pyplot as plt  import math  import numpy  import random    fig = plt.figure()  ax = fig.add_subplot(111)    #阶数为9阶  order=9  #生成曲线上的各个点  x = numpy.arange(-1,1,0.02)  y = [((a*a-1)*(a*a-1)*(a*a-1)+0.5)*numpy.sin(a*2) for a in x]  #ax.plot(x,y,color='r',linestyle='-',marker='')  #,label="(a*a-1)*(a*a-1)*(a*a-1)+0.5"    #生成的曲线上的各个点偏移一下，并放入到xa,ya中去  i=0  xa=[]  ya=[]  for xx in x:      yy=y[i]      d=float(random.randint(70,130))/100      #ax.plot([xx*d],[yy*d],color='m',linestyle='',marker='.')      i+=1      xa.append(xx*d)      ya.append(yy*d)    '''''for i in range(0,5):     xx=float(random.randint(-100,100))/100     yy=float(random.randint(-60,60))/100     xa.append(xx)     ya.append(yy)'''  ax.plot(xa,ya,color='m',linestyle='',marker='.')  #plt.show()#print(ax)  #进行曲线拟合  matA=[]  for i in range(0,order+1):      matA1=[]      for j in range(0,order+1):          tx=0.0          for k in range(0,len(xa)):              dx=1.0              for l in range(0,j+i):                  dx=dx*xa[k]              tx+=dx          matA1.append(tx)      matA.append(matA1)    #print(len(xa))  #print(matA[0][0])  matA=numpy.array(matA)    matB=[]  for i in range(0,order+1):      ty=0.0      for k in range(0,len(xa)):          dy=1.0          for l in range(0,i):              dy=dy*xa[k]          ty+=ya[k]*dy      matB.append(ty)     matB=numpy.array(matB)    matAA=numpy.linalg.solve(matA,matB)  print (matAA)print (matA)print (matB)#画出拟合后的曲线  #print(matAA)  xxa= numpy.arange(-1,1.06,0.01)  yya=[]  for i in range(0,len(xxa)):      yy=0.0      for j in range(0,order+1):          dy=1.0          for k in range(0,j):              dy*=xxa[i]          dy*=matAA[j]          yy+=dy      yya.append(yy)  ax.plot(xxa,yya,color='g',linestyle='-',marker='')    ax.legend("aw")  plt.show()

python的结果

三，下面是c++实现。

c++第一部分：这里用了一个自己写的Matrix 类，作为头文件Matrix.h

#include
   
       #include 
    
           // std::ifstream#include 
     
            #include 
      
               using namespace std;class Matrix{private:	unsigned  row, col, size;	double *pmm;//数组指针    public:	Matrix(unsigned r, unsigned c) :row(r), col(c)//非方阵构造     	{		size = r*c;		if (size>0)		{			pmm = new double[size];			for (unsigned j = 0; j
       
        0)		{			pmm = new double[size];			for (unsigned j = 0; j
        
         >(istream&, Matrix&); friend ofstream &operator<<(ofstream &out, Matrix &obj); friend ostream &operator<<(ostream&, Matrix&); friend Matrix operator+(const Matrix&, const Matrix&); friend Matrix operator*(const Matrix&, const Matrix&); Matrix cov(_In_opt_ bool flag = true); //协方差阵 或者样本方差 double det(); //行列式 Matrix solveAb(Matrix &obj); // b是行向量或者列向量 Matrix diag(); //返回对角线元素 //Matrix asigndiag(); //对角线元素 Matrix T()const; //转置 void sort(bool);//true为从小到大 Matrix adjoint(); Matrix inverse(); void QR(_Out_ Matrix&, _Out_ Matrix&)const; Matrix eig_val(_In_opt_ unsigned _iters = 1000); Matrix eig_vect(_In_opt_ unsigned _iters = 1000); unsigned Row()const{ return row; } unsigned Col()const{ return col; } double norm1();//1范数 double norm2();//2范数 double*operator[](int i){ return pmm + i*col; }//注意this加括号， (*this)[i][j] void zeromean(_In_opt_ bool flag = true);//默认参数为true计算列 void normalize(_In_opt_ bool flag = true);//默认参数为true计算列 Matrix exponent(double x);//每个元素x次幂 Matrix eye();//对角阵 void maxlimit(double max,double set=0);//对角阵};//友元仅仅是指定了访问权限，不是一般意义的函数声明，最好在类外再进行一次函数声明。 //istream &operator>>(istream &is, Matrix &obj); //ostream &operator<<(ostream &is, Matrix &obj); //对称性运算符，如算术，相等，应该是普通非成员函数。 //Matrix operator*( const Matrix&,const Matrix& ); //Matrix operator+( const Matrix&,const Matrix& ); //dets递归调用 double dets(int n, double *&aa){ if (n == 1) return aa[0]; double *bb = new double[(n - 1)*(n - 1)];//创建n-1阶的代数余子式阵bb int mov = 0;//判断行是否移动 double sum = 0.0;//sum为行列式的值 for (int arow = 0; arow
         
           brow ? 0 : 1; //bb中小于arow的行，同行赋值，等于的错过，大于的加一 for (int j = 0; j
          
           max ? 0 : pmm[i*col + j]; } } }Matrix Matrix::eye()//对角阵{ for (unsigned i = 0; i< row; i++) { for (unsigned j = 0; j < col; j++) { if (i == j) { pmm[i*col + j] = 1.0; } } } return *this;}void Matrix::zeromean(_In_opt_ bool flag){ if (flag == true) //计算列均值 { double *mean = new double[col]; for (unsigned j = 0; j < col; j++) { mean[j] = 0.0; for (unsigned i = 0; i < row; i++) { mean[j] += pmm[i*col + j]; } mean[j] /= row; } for (unsigned j = 0; j < col; j++) { for (unsigned i = 0; i < row; i++) { pmm[i*col + j] -= mean[j]; } } delete[]mean; } else //计算行均值 { double *mean = new double[row]; for (unsigned i = 0; i< row; i++) { mean[i] = 0.0; for (unsigned j = 0; j < col; j++) { mean[i] += pmm[i*col + j]; } mean[i] /= col; } for (unsigned i = 0; i < row; i++) { for (unsigned j = 0; j < col; j++) { pmm[i*col + j] -= mean[i]; } } delete[]mean; }}void Matrix::normalize(_In_opt_ bool flag){ if (flag == true) //计算列均值 { double *mean = new double[col]; for (unsigned j = 0; j < col; j++) { mean[j] = 0.0; for (unsigned i = 0; i < row; i++) { mean[j] += pmm[i*col + j]; } mean[j] /= row; } for (unsigned j = 0; j < col; j++) { for (unsigned i = 0; i < row; i++) { pmm[i*col + j] -= mean[j]; } } ///计算标准差 for (unsigned j = 0; j < col; j++) { mean[j] = 0; for (unsigned i = 0; i < row; i++) { mean[j] += pow(pmm[i*col + j],2);//列平方和 } mean[j] = sqrt(mean[j] / row); // 开方 } for (unsigned j = 0; j < col; j++) { for (unsigned i = 0; i < row; i++) { pmm[i*col + j] /= mean[j];//列平方和 } } delete[]mean; } else //计算行均值 { double *mean = new double[row]; for (unsigned i = 0; i< row; i++) { mean[i] = 0.0; for (unsigned j = 0; j < col; j++) { mean[i] += pmm[i*col + j]; } mean[i] /= col; } for (unsigned i = 0; i < row; i++) { for (unsigned j = 0; j < col; j++) { pmm[i*col + j] -= mean[i]; } } ///计算标准差 for (unsigned i = 0; i< row; i++) { mean[i] = 0.0; for (unsigned j = 0; j < col; j++) { mean[i] += pow(pmm[i*col + j], 2);//列平方和 } mean[i] = sqrt(mean[i] / col); // 开方 } for (unsigned i = 0; i < row; i++) { for (unsigned j = 0; j < col; j++) { pmm[i*col + j] /= mean[i]; } } delete[]mean; }}double Matrix::det(){ if (col == row) return dets(row, pmm); else { cout << ("行列不相等无法计算") << endl; return 0; }}/ istream &operator>>(istream &is, Matrix &obj){ for (unsigned i = 0; i
           
            > obj.pmm[i]; } return is;}ostream &operator<<(ostream &out, Matrix &obj){ for (unsigned i = 0; i < obj.row; i++) //打印逆矩阵 { for (unsigned j = 0; j < obj.col; j++) { out << (obj[i][j]) << "\t"; } out << endl; } return out;}ofstream &operator<<(ofstream &out, Matrix &obj)//打印逆矩阵到文件 { for (unsigned i = 0; i < obj.row; i++) { for (unsigned j = 0; j < obj.col; j++) { out << (obj[i][j]) << "\t"; } out << endl; } return out;}Matrix operator+(const Matrix& lm, const Matrix& rm){ Matrix ret(lm.row, lm.col); for (unsigned i = 0; i
            
             pmm[j]) { tem = pmm[i]; pmm[i] = pmm[j]; pmm[j] = tem; } } else { if (pmm[i]
             
              << endl; return m; } Matrix m(row); for (unsigned i = 0; i
              
               = 0; --m) { sum = 0; for (unsigned j = m + 1; j < col; ++j) { sum += pmm[m * col + j] * ret.pmm[j * ret.col + count]; } sum = -sum / pmm[m * col + m]; midSum += sum * sum; ret.pmm[m * ret.col + count] = sum; } midSum = sqrt(midSum); for (unsigned i = 0; i < ret.row; ++i) { ret.pmm[i * ret.col + count] /= midSum; //每次求出一个列向量 } } *this = matcopy;//恢复原始矩阵； return ret;}Matrix Matrix::cov(bool flag){ //row 样本数，column 变量数 if (col == 0) { Matrix m(0); return m; } double *mean = new double[col]; //均值向量 for (unsigned j = 0; j
               
                bi) //ai行的代数余子式是：小于ai的行，aa与bb阵，同行赋值 pi = 0; else pi = 1; //大于等于ai的行，取aa阵的ai+1行赋值给阵bb的bi行 if (aj>bj) //ai行的代数余子式是：小于ai的行，aa与bb阵，同行赋值 pj = 0; else pj = 1; //大于等于ai的行，取aa阵的ai+1行赋值给阵bb的bi行 bb[bi*(n - 1) + bj] = pmm[(bi + pi)*n + bj + pj]; } } if ((ai + aj) % 2 == 0) q = 1;//因为列数为0，所以行数是偶数时候，代数余子式为-1. else q = (-1); ret.pmm[ai*n + aj] = q*dets(n - 1, bb); //加符号变为代数余子式 } } delete[]bb; return ret;}Matrix Matrix::inverse(){ double det_aa = det(); if (det_aa == 0) { cout << "行列式为0 ，不能计算逆矩阵。" << endl; Matrix rets(0); return rets; } Matrix adj = adjoint(); Matrix ret(row); for (unsigned i = 0; i
               
              
             
            
           
          
         
        
       
      
     
    
   

c++第二部分：主函数部分

#include "matrix.h"#include 
   
    #include
    
     #include
     
      const int N = 100;double fun(double x){	return (pow((x*x - 1), 3.0) + 0.5)*sin(x * 2);}int main(){	Matrix x(1, N);	Matrix y(1, N); 	ofstream ofs;	ofs.open("pydata.txt", ofstream::out);	double src = -1.0;	for (int i = 0; i < N; i++)	{		src += 0.02;		x[0][i] = src;		y[0][i] = fun(src);	}	Matrix sx(x);//加入随机因素前的x	for (int i = 0; i < N; i++)	{		double ran = rand() % 600 / 1000.0 + 0.7;  // 0.7~1.3 倍的随机		x[0][i] = x[0][i]*ran;		y[0][i] = y[0][i] * ran;	} 	//构建矩阵	int k = 8;  //次方数	Matrix xx(k + 1, k + 1); //累加矩阵k+1阶	for (size_t r = 0; r < k+1; r++)	{		for (size_t c = r; c < k+1; c++)		{			double sum = 0;			for (size_t i = 0; i < N; i++)			{				sum += pow(x[0][i], r+c);			}			xx[c][r]=xx[r][c] = sum;		}	}	cout << "xx:\n"<
     
    
   

c++第三部分： 根据c++的输出  pydata.txt 用python画图

import matplotlib.pyplot as plt  import numpy  fin = open("pydata.txt", "r")    line = fin.readline()   x=[]y=[]py=[]while line:        #print(line)    line=line.strip()        line=line.split()        x.append(float(line[0]))    y.append(float(line[1]))    py.append(float(line[2]))    line = fin.readline()    fin.close()    ax = plt.subplot(111)  xa= numpy.arange(-1,1,0.02)  ax.plot(xa,y,color='m',linestyle='',marker='.')  ax.plot(xa,py,color='g',linestyle='-',marker='')  plt.show()

画图的结果