在已知曲线函数的情况下，对已知数据使用最小二乘法进行参数求解，从而实现曲线拟合的方法最为简便高效。

数据如下，匹配函数为 y = a0+a1*x .

0	1	2	3	4	5
x	1.5	2.5	3.5	4.5	5.5
y	4	6	8	10	12

高斯消元法的声明文件头

#ifndef GAUSSEQUATION_H
#define GAUSSEQUATION_H
#include <vector>

class GaussEquation
{
    public:
        GaussEquation();
        GaussEquation(int M,double *matrix,double *b);
        virtual ~GaussEquation();
        bool Resolve(std::vector<double>& xValue);
    protected:

    private:
        double *m_matrixA;
        int m_DIM;
        double *m_bVal;
        int SelectPivotalElement(int col); //find the row number of max value in column i;
        bool IsPrecisionZero(double aij);
        void SimplePivotalRow(int i,int j);
        void SwapRow(int i,int j);
        void RowElimination(int i,int j,int k);
};

#endif // GUASSEQUATION_H

高斯消元法的类定义源文件

#include "GaussEquation.h"
#include <cassert>
#include <cmath>
GaussEquation::GaussEquation()
{
    //ctor
}

GaussEquation::~GaussEquation()
{
    //dtor
    delete[] m_matrixA;
    delete[] m_bVal;
}


GaussEquation::GaussEquation(int M,double *matrix,double *b)
{
    m_DIM = M;
    m_matrixA = new double[M*M];
    m_bVal = new double[M];
    for(int i = 0; i <M; i++)
    {   m_bVal[i] = b[i];
        for(int j = 0; j < M; j++)
        {
        m_matrixA[i*m_DIM+j] = matrix[i*m_DIM+j];
        }
    }
}


bool
GaussEquation::Resolve(std::vector<double>& xValue)
{

    assert(static_cast<int>(xValue.size()) == m_DIM);

    for(int i = 0; i<m_DIM;i++)
    {

        int pivotRow = SelectPivotalElement(i);
        if(pivotRow != i)
        {
            SwapRow(i,pivotRow);
        }
        if(IsPrecisionZero(m_matrixA[i*m_DIM+i]))
        {
            return false;
        }

        SimplePivotalRow(i,i);

        for(int j = i + 1; j < m_DIM; j++)
        {
            RowElimination(i,j,i);
        }
    }

    m_matrixA[(m_DIM-1)*m_DIM+m_DIM-1] = m_bVal[m_DIM-1];///m_matrixA[(m_DIM-1)*m_DIM+m_DIM-1];
    for(int i =m_DIM-2;i>=0;i--)//行 i
    {
        double totalCof = 0.0;
        for(int j=i+1;j<m_DIM;j++)//列 j 
        {
            totalCof += m_matrixA[i*m_DIM+j]*m_matrixA[j*m_DIM+j];
        }
        m_matrixA[i*m_DIM +i] = (m_bVal[i]-totalCof); // 行首已经归一化 

    }

    for(int i = 0; i < m_DIM; i++)
    {
        xValue[i] = m_matrixA[i*m_DIM+i];
    }

    return true;
}


bool
GaussEquation::IsPrecisionZero(double aij)
{
    return aij<0.001&&aij>-0.001;
}


//Elimination row j by row i and start this operation from column k
void
GaussEquation::RowElimination(int r_i,int r_j,int c_k)
{
    double scale;
    scale = m_matrixA[r_j*m_DIM+c_k]/m_matrixA[r_i*m_DIM+c_k];
    for(int n=c_k;n<m_DIM;n++)
    {
        m_matrixA[r_j*m_DIM+n] -= m_matrixA[r_i*m_DIM+n]*scale;
    }
    m_bVal[r_j] -= m_bVal[r_i]*scale;
}


//
int
GaussEquation::SelectPivotalElement(int col)
{
    int r=col;
    for(int i=col+1;i<m_DIM;i++)
    {
        if (fabs(m_matrixA[r*m_DIM+col])<fabs(m_matrixA[i*m_DIM+col]))
            r = i;
    }
    return r;
}


//==Transform Ar_i,c_j to 1;
void
GaussEquation::SimplePivotalRow(int r_i,int c_j)
{
    double scale = 1.0/m_matrixA[r_i*m_DIM+c_j];

    for(int i=c_j;i<m_DIM;i++)
    {
        m_matrixA[r_i*m_DIM+i]  *= scale;
    }
	m_bVal[r_i] *= scale;
}


//==Swap row r_i with row r_j
void
GaussEquation::SwapRow(int r_i,int r_j)
{
    double a_tmp,b_tmp;
    for(int i = 0; i < m_DIM; i++)
    {
        a_tmp = m_matrixA[r_i*m_DIM + i];
        m_matrixA[r_i*m_DIM + i] = m_matrixA[r_j*m_DIM + i];
        m_matrixA[r_j*m_DIM + i] = a_tmp;
    }

    b_tmp = m_bVal[r_i];
    m_bVal[r_i] = m_bVal[r_j];
    m_bVal[r_j] = b_tmp;

}

main程序

#include <iostream>
#include <vector>
#include "GaussEquation.h"
using namespace std;
double *InputArray(double* x, int size);
double *OutputArray(double *x, double *y, int size);

int main()
{
    int size;
    vector<double> xValue(2);
    vector<double>::iterator it = xValue.begin();
    bool finished;

    cout<<"Input the number of data set:";
    cin>>size;
    //scanf("%d",&size);

    double x[size],y[size];
    for(int i = 0; i < size; i++)
    {

        cin>>x[i]>>y[i];//scanf("%f,%f",x+i,y+i);
    }

    double *A = InputArray(x,size);
    double *b = OutputArray(x,y,size);

    GaussEquation equation(2,A,b);
    delete[] A;
    delete[] b;
    finished = equation.Resolve(xValue);
    cout<<"status:"<<finished<<endl;


    cout<<"a0: "<<(*it)<<",";
    it++;
    cout<<"a1: "<<(*it)<<endl;
}

// A = aT*a
double *InputArray(double* x, int size){
//    std::vector<double> A(2*2,0);
//    std::vector<double> a(size*2,1);
    double *A = new double[2*2];
    double  *a = new double[size*2]; // a = [(1,x),...]

    A[0] = 0;
    A[1] = 0;
    A[2] = 0;
    A[3] = 0;
    for(int i = 0; i < size; i++){
        a[i*2+1] = x[i];
        a[i*2+0] = 1;
    }

    for(int i = 0; i < 2; i++)
    {
        for(int j = 0; j < 2; j++)
        {
            for(int k = 0; k < size; k++)
            {
            A[i*2+j] += a[k*2+i]*a[k*2+j];  //aT[i,j] = a[j,i] --> aT[i,k] = a[k,i]
            }
        }
    }
    delete[] a;
    return A;

}

// b = aT*y
double *OutputArray(double *x, double *y, int size)
{
    double *b = new double[2];
    double  *a = new double[size*2];
    b[0] = 0;
    b[1] = 0;
    // a = [(1,x),...]
    // aT[i,j] =a [j,i]
    for(int i = 0; i < size; i++){
        a[i*2+1] = x[i];
        a[i*2+0] = 1;
    }

    // b = aT*y
    for(int i = 0; i < 2; i++)
    {
        for(int j = 0; j < size; j++)
        {
        b[i] += a[j*2+i]*y[j];
        }
    }
    delete[] a;
    return b;

}

程序执行结果：