实验数据与曲线拟合1. 曲线拟合1. 曲线拟合的定义2. 简单线性数据拟合的例子2. 最小二乘法曲线拟合1. 最小二乘法原理2. 高斯消元法求解方程组3. 最小二乘法解决速度与加速度实验3. 三次样条曲线拟合1. 插值函数2. 样条函数的定义3. 边界条件4. 推导三次样条函数5. 追赶法求解方程组6. 三次样条曲线拟合算法实现7. 三次样条曲线拟合的效果4. 12.1 曲线拟合5. 12.1.1 曲线拟合的定义6. 曲线拟合(Curve Fitting)的数学定义是指用连续曲线近似地刻画或比拟平面上一组离散点所表示的坐标之间的函数关系，是一种用解析表达式逼近离散数据的方法。
曲线拟合通俗的说法就是“拉曲线”，也就是将现有数据透过数学方法来代入一条数学方程式的表示方法。
科学和工程遇到的很多问题，往往只能通过诸如采样、实验等方法获得若干离散的数据，根据这些数据，如果能够找到一个连续的函数（也就是曲线）或者更加密集的离散方程，使得实验数据与方程的曲线能够在最大程度上近似吻合，就可以根据曲线方程对数据进行数学计算，对实验结果进行理论分析，甚至对某些不具备测量条件的位置的结果进行估算。
7. 12.1.2 简单线性数据拟合的例子8. 回想一下中学物理课的“速度与加速度”实验：假设某物体正在做加速运动，加速度未知，某实验人员从时间t0 = 3秒时刻开始，以1秒时间间隔对这个物体连续进行了12次测速，得到一组速度和时间的离散数据，请根据实验结果推算该物体的加速度。
9. 表 12 – 1 物体速度和时间的测量关系表10. 在选择了合适的坐标刻度之后，我们就可以在坐标纸上画出这些点。
如图12–1所示，排除偏差明显偏大的测量值后，可以看出测量结果呈现典型的线性特征。
沿着该线性特征画一条直线，使尽量多的测量点能够位于直线上，或与直线的偏差尽量小，这条直线就是我们根据测量结果拟合的速度与时间的函数关系。
最后在坐标纸上测量出直线的斜率K，K就是被测物体的加速度，经过测量，我们实验测到的物体加速度值是1.48米/秒2。
11.12. 图 12 – 1 实验法测量加速度的过程13.14. 12.2 最小二乘法曲线拟合15. 使用数学分析进行曲线拟合有很多常用的方法，这一节我们先介绍一下最简单的最小二乘法，并给出使用最小二乘法解决上一节给出的速度与加速度实验问题。
16.17. 12.2.1 最小二乘法原理18. 最小二乘法（又称最小平方法）通过最小化误差的平方和寻找数据的最佳函数匹配，利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小，当然，做为一种插值方法使用时，最小二乘法也可以用于曲线拟合。
使用最小二乘法进行曲线拟合是曲线拟合种早期的一种常用方法，不过，最小二乘法理论简单，计算量小，即便是在使用三次样条曲线或RBF（Radial Basis Function）进行曲线拟合大行其道的今天，最小二乘法在多项式曲线或直线的拟合问题上，仍然得到广泛地应用。
使用最小二乘法，选取的匹配函数的模式非常重要，如果离散数据呈现的是指数变化规律，则应该选择指数形式的匹配函数模式，如果是多项式变化规律，则应该选择多项式匹配模式，如果选择的模式不对，拟合的效果就会很差，这也是使用最小二乘法进行曲线拟合时需要特别注意的一个地方。
19. 下面以多项式模式为例，介绍一下使用最小二乘法进行曲线拟合的完整步骤。
假设选择的拟合多项式模式是：20.21. 这m个等式相当于m个方程，a0,a1,…a m是m个未知量，因此这m个方程组成的方程组是可解的，最小二乘法的第二步处理就是将其整理为针对a0,a1,…a m的正规方程组。
最终整理的方程组如下：22.23. 最小二乘法的第三步处理就是求解这个多元一次方程组，得到多项式的系数a0,a1,…a m，，就可以得到曲线的拟合多项式函数。
求解多元一次方程组的方法很多，高斯消元法是最常用的一种方法，下一节就简单介绍一下最小二乘算法实现所用的高斯消元法算法。
24. 12.2.2 高斯消元法求解方程组25. 在数学上，高斯消元法是线性代数中的一个算法，可用来求解多元一次线性方程组，也可以用来求矩阵的秩，以及求可逆方阵的逆矩阵。
高斯消元法虽然以数学家高斯的名字命名，但是最早出现在文献资料中应该是中国的《九章算术》。
26. 高斯消元法的主要思想是通过对系数矩阵进行行变换，将方程组的系数矩阵由对称矩阵变为三角矩阵，从而达到消元的目的，最后通过回代逐个获得方程组的解。
在消元的过程中，如果某一行的对角线元素的值太小，在计算过程中就会出现很大的数除以很小的数的情况，有除法溢出的可能，因此在消元的过程中，通常都会增加一个主元选择的步骤，通过行交换操作，将当前列绝对值最大的行交换到当前行位置，避免了除法溢出问题，增加了算法的稳定性。
27. 高斯消元法算法实现简单，主要有两个步骤组成，第一个步骤就是通过选择主元，逐行消元，最终行程方程组系数矩阵的三角矩阵形式，第二个步骤就是逐步回代的过程，最终矩阵的对角线上的元素就是方程组的解。
下面就给出高斯消元法的一个算法实现：76/*带列主元的高斯消去法解方程组，最后的解在matrixA的对角线上*/77bool GuassEquation::Resolve(std::vector<double>& xValue)78{79 assert(xValue.size()== m_DIM);8081/*消元，得到上三角阵*/82for(int i =0; i < m_DIM -1; i++)83{84/*按列选主元*/85int pivotRow = SelectPivotalElement(i);86if(pivotRow != i)/*如果有必要，交换行*/87{88 SwapRow(i, pivotRow);89}90if(IsPrecisionZero(m_matrixA[i * m_DIM + i]))/*主元是0? 不存在唯一解*/91{92return false;93}94/*对系数归一化处理，使行第一个系数是1.0*/95 SimplePivotalRow(i, i);96/*逐行进行消元*/97for(int j = i +1; j < m_DIM; j++)98{99 RowElimination(i, j, i);100}101}102/*回代求解*/103 m_matrixA[(m_DIM -1)* m_DIM + m_DIM -1]= m_bVal[m_DIM -1]/m_matrixA[(m_D IM -1)* m_DIM + m_DIM -1];104for(int i = m_DIM -2; i >=0; i--)105{106double totalCof =0.0;107for(int j = i +1; j < m_DIM; j++)108{109 totalCof += m_matrixA[i * m_DIM + j]* m_matrixA[j * m_DIM + j];110}111 m_matrixA[i * m_DIM + i]=(m_bVal[i]- totalCof)/ m_matrixA[i * m_DIM+ i ];112}113114/*将对角线元素的解逐个存入解向量*/115for(int i =0; i < m_DIM; i++)116{117 xValue[i]= m_matrixA[i * m_DIM + i];118}119120return true;121}28.29. GuassEquation::Resolve()函数中m_matrixA是以一维数组形式存放的系数矩阵，m_DIM是矩阵的维数，SelectPivotalElement()函数从系数矩阵的第i列中选择绝对值最大的那个值所在的行，并返回行号，SwapRow()函数负责交换系数矩阵两个行的所有值，SimplePivotalRow()函数是归一化处理函数，通过除法操作将指定的行的对角线元素变换为1.0，以便简化随后的消元操作。
30. 12.2.3 最小二乘法解决“速度与加速度”实验31. 根据12.2.1节对最小二乘法原理的分析，用程序实现最小二乘法曲线拟合的算法主要由两个步骤组成，第一个步骤就是根据给出的测量值生成关于拟合多项式系数的方程组，第二个步骤就是解这个方程组，求出拟合多项式的各个系数。
根据对上文最终整理的正规方程组的分析，可以看出其系数有一定的关系，就是每一个方程式都比前一个方程式多乘了一个x i。
因此，只需要完整计算出第一个方程式的系数，其他方程式的系数只是将前一个方程式的系数依次左移一位，然后单独计算出最后一个系数就可以了，此方法可以减少很多无谓的计算。
求解多元一次方程组的方法就使用12.2.2节介绍的高斯消元法，其算法上一节已经给出。
32. 这里给出一个最小二乘算法的完整实现，以12.1.2节的数据为例，因为数据结果明显呈现线性方程的特征，因此选择拟合多项式为v = v0 + at，v0和a就是要求解的拟合多项式系数。
99bool LeastSquare(const std::vector<double>& x_value,const std::vector<double>&y_v alue,100int M, std::vector<double>& a_value)101{102 assert(x_value.size()== y_value.size());103 assert(a_value.size()== M);104105double*matrix =new double[M * M];106double*b=new double[M];107108 std::vector<double> x_m(x_value.size(),1.0);109 std::vector<double> y_i(y_value.size(),0.0);110for(int i =0; i < M; i++)111{112 matrix[ARR_INDEX(0, i, M)]= std::accumulate(x_m.begin(), x_m.end(),0.0); 113for(int j =0; j <static_cast<int>(y_value.size()); j++)114{115 y_i[j]= x_m[j]* y_value[j];116}117 b[i]= std::accumulate(y_i.begin(), y_i.end(),0.0);118for(int k =0; k <static_cast<int>(x_m.size()); k++)119{120 x_m[k]*= x_value[k];121}122}123for(int row =1; row < M; row++)124{125for(int i =0; i < M -1; i++)126{127 matrix[ARR_INDEX(row, i, M)]= matrix[ARR_INDEX(row -1, i +1,M)]; 128}129 matrix[ARR_INDEX(row, M -1, M)]= std::accumulate(x_m.begin(),x_m.end(), 0.0);130for(int k =0; k <static_cast<int>(x_m.size()); k++)131{132 x_m[k]*= x_value[k];133}134}135136 GuassEquation equation(M, matrix, b);137delete[] matrix;138delete[] b;139140return equation.Resolve(a_value);141}33. 将表12－1的数据带入算法，计算得到v0 = 4.05545455，a = 1.48818182，比作图法得到的结果更精确。