精华内容
参与话题
问答
  • 最小二乘法曲线拟合以及Matlab实现

    万次阅读 多人点赞 2017-12-28 17:31:23
    最小二乘法曲线拟合以及Matlab实现 在实际工程中,我们常会遇到这种问题:已知一组点的横纵坐标,需要绘制出一条尽可能逼近这些点的曲线(或直线),以进行进一步进行加工或者分析两个变量之间的相互关系。而获取这...

    最小二乘法曲线拟合以及Matlab实现

    在实际工程中,我们常会遇到这种问题:已知一组点的横纵坐标,需要绘制出一条尽可能逼近这些点的曲线(或直线),以进行进一步进行加工或者分析两个变量之间的相互关系。而获取这个曲线方程的过程就是曲线拟合。

    目录

    • 最小二乘法直线拟合原理
    • 曲线拟合
    • Matlab实现代码

    最小二乘法直线线拟合原理

    首先,我们从曲线拟合的最简单情况——直线拟合来引入问题。如果待拟合点集近似排列在一条直线上时,我们可以设直线 y=ax+by=ax+by=ax+b 为其拟合方程,系数 A=[a,b]A=[a,b]A=[a,b] 为待求解项,已知:
    这里写图片描述

    用矩阵形式表达为: Y=X0AY=X_{0}AY=X0A,其中:
    这里写图片描述
    要求解A,可在方程两边同时左乘 X0X_{0}X0 的逆矩阵,如果它是一个方阵且非奇异的话。

    但是,一般情况下 X0X_{0}X0 连方阵都不是,所以我们在此需要用 X0X_{0}X0 构造一个方阵,即方程两边同时左乘 X0X_{0}X0 的转置矩阵,得到方程: X0TY=X0TX0AX_{0}^{T}Y=X_{0}^{T}X_{0}AX0TY=X0TX0A

    此时,方程的系数矩阵 X0TX0X_{0}^{T}X_{0}X0TX0 为方阵,所以两边同时左乘新系数矩阵 X0TX0X_{0}^{T}X_{0}X0TX0 的逆矩阵,便可求得系数向量A ,即:(X0TX0)−1X0TY=A(X_{0}^{T}X_{0})^{-1}X_{0}^{T}Y=A(X0TX0)1X0TY=A

    方程A=(X0TX0)−1X0TYA=(X_{0}^{T}X_{0})^{-1}X_{0}^{T}YA=(X0TX0)1X0TY 右边各部分均已知,所以可直接求解得到拟合直线的方程系数向量A。

    曲线线拟合

    当样本点的分布不为直线时,我们可用多项式曲线拟合,即拟合曲线方程为n阶多项式

    y=∑i=0naixi=anxn+an−1xn−1+...+a1x+a0y=\sum_{i=0}^{n}a_ix^i=a_nx^n+a_{n-1}x^{n-1}+...+a_1x+a_0y=i=0naixi=anxn+an1xn1+...+a1x+a0

    用矩阵形式表示为: Y=X0AY=X_0AY=X0A ,其中:

    X0

    待求解项为系数向量A=[an,an−1,...,a2,a1,a0]TA=[a_n,a_{n-1},...,a_2,a_1,a_0]^TA=[an,an1,...,a2,a1,a0]T

    曲线拟合方程Y=X0AY=X_0AY=X0A 的求解方法与上面直线的求解方法一样,也是在方程Y=X0AY=X_0AY=X0A 两边同左乘X0X_0X0的转置矩阵得到: X0TY=X0TX0AX_{0}^{T}Y=X_{0}^{T}X_{0}AX0TY=X0TX0A

    再同时在新方程两边同时左乘X0TX0X_{0}^{T}X_{0}X0TX0 的逆矩阵,得到:(X0TX0)−1X0TY=A(X_{0}^{T}X_{0})^{-1}X_{0}^{T}Y=A(X0TX0)1X0TY=A

    上式左边各部分均已知,所以可直接求解得拟合曲线方程的系数向量A。

    Matlab实现代码

    %by hanlestudy@163.com
    clear
    clc
    x=[2,4,5,6,6.8,7.5,9,12,13.3,15];
    y=[-10,-6.9,-4.2,-2,0,2.1,3,5.2,6.4,4.5];
    [~,k]=size(x);
    for n=1:9
        X0=zeros(n+1,k);
        for k0=1:k           %构造矩阵X0
            for n0=1:n+1
                X0(n0,k0)=x(k0)^(n+1-n0);
            end
        end
        X=X0';
        ANSS=(X'*X)\X'*y';
        for i=1:n+1          %answer矩阵存储每次求得的方程系数,按列存储
           answer(i,n)=ANSS(i);
       end
        x0=0:0.01:17;
        y0=ANSS(1)*x0.^n    ;%根据求得的系数初始化并构造多项式方程
        for num=2:1:n+1     
            y0=y0+ANSS(num)*x0.^(n+1-num);
        end
        subplot(3,3,n)
        plot(x,y,'*')
        hold on
        plot(x0,y0)
    end
    suptitle('不同次数方程曲线拟合结果,从1到9阶')
    

    运行结果

    拟合曲线结果:

    可以看出看来,当多项式的阶数过小是,曲线并不能很好地反映出样本点的分布情况;但阶数过高时,会出现过拟合的情况。
    这里写图片描述

    系数矩阵answer:

    这里写图片描述

    Matlab自带函数——polyfit

    在matlab中,也有现成的曲线拟合函数polyfit,其也是基于最小二乘原理实现的,具体用法为:ans=polyfit(x,y,n). 其中x,y为待拟合点的坐标向量,n为多项式的阶数。

    下面代码是用polyfit函数来做曲线拟合:

    clear
    clc
    x=[2,4,5,6,6.8,7.5,9,12,13.3,15];
    [~,k]=size(x);
    y=[-10,-6.9,-4.2,-2,0,2.1,3,5.2,6.4,4.5];
    for n=1:9
        ANSS=polyfit(x,y,n);  %用polyfit拟合曲线
        for i=1:n+1           %answer矩阵存储每次求得的方程系数,按列存储
           answer(i,n)=ANSS(i);
       end
        x0=0:0.01:17;
        y0=ANSS(1)*x0.^n    ; %根据求得的系数初始化并构造多项式方程
        for num=2:1:n+1     
            y0=y0+ANSS(num)*x0.^(n+1-num);
        end
        subplot(3,3,n)
        plot(x,y,'*')
        hold on
        plot(x0,y0)
    end
    suptitle('不同次数方程曲线拟合结果,从1到9阶')
    

    运行结果:

    用polyfit拟合的结果与第一份代码运行的结果基本一样
    这里写图片描述

    申明

    本文为本人原创,转载请注明出处!

    展开全文
  • 算法系列之二十一:实验数据与曲线拟合

    万次阅读 多人点赞 2013-10-16 22:17:15
    曲线拟合(Curve Fitting)的数学定义是指用连续曲线近似地刻画或比拟平面上一组离散点所表示的坐标之间的函数关系,是一种用解析表达式逼近离散数据的方法。曲线拟合通俗的说法就是“拉曲线”,也就是将现有数据透过...

    12.1 曲线拟合

    12.1.1 曲线拟合的定义

            曲线拟合(Curve Fitting)的数学定义是指用连续曲线近似地刻画或比拟平面上一组离散点所表示的坐标之间的函数关系,是一种用解析表达式逼近离散数据的方法。曲线拟合通俗的说法就是“拉曲线”,也就是将现有数据透过数学方法来代入一条数学方程式的表示方法。科学和工程遇到的很多问题,往往只能通过诸如采样、实验等方法获得若干离散的数据,根据这些数据,如果能够找到一个连续的函数(也就是曲线)或者更加密集的离散方程,使得实验数据与方程的曲线能够在最大程度上近似吻合,就可以根据曲线方程对数据进行数学计算,对实验结果进行理论分析,甚至对某些不具备测量条件的位置的结果进行估算。

    12.1.2 简单线性数据拟合的例子

            回想一下中学物理课的“速度与加速度”实验:假设某物体正在做加速运动,加速度未知,某实验人员从时间t0 = 3秒时刻开始,以1秒时间间隔对这个物体连续进行了12次测速,得到一组速度和时间的离散数据,请根据实验结果推算该物体的加速度。

    时间 ()

    3

    4

    5

    6

    7

    8

    9

    10

    速度(米/秒)

    8.41

    9.94

    11.58

    13.02

    14.33

    15.92

    17.54

    19.22

    时间 ()

    11

    12

    13

    14

     

     

     

     

    速度(米/秒)

    20.49

    22.01

    23.53

    24.47

     

     

     

     

    12 – 1 物体速度和时间的测量关系表

            在选择了合适的坐标刻度之后,我们就可以在坐标纸上画出这些点。如图121所示,排除偏差明显偏大的测量值后,可以看出测量结果呈现典型的线性特征。沿着该线性特征画一条直线,使尽量多的测量点能够位于直线上,或与直线的偏差尽量小,这条直线就是我们根据测量结果拟合的速度与时间的函数关系。最后在坐标纸上测量出直线的斜率KK就是被测物体的加速度,经过测量,我们实验测到的物体加速度值是1.48米/2

     

    12 – 1 实验法测量加速度的过程

     

    12.2 最小二乘法曲线拟合

            使用数学分析进行曲线拟合有很多常用的方法,这一节我们先介绍一下最简单的最小二乘法,并给出使用最小二乘法解决上一节给出的速度与加速度实验问题。

     

    12.2.1 最小二乘法原理

            最小二乘法(又称最小平方法)通过最小化误差的平方和寻找数据的最佳函数匹配,利用最小二乘法可以简便地求得未知的数据,并使得这些求得的数据与实际数据之间误差的平方和为最小,当然,做为一种插值方法使用时,最小二乘法也可以用于曲线拟合。使用最小二乘法进行曲线拟合是曲线拟合种早期的一种常用方法,不过,最小二乘法理论简单,计算量小,即便是在使用三次样条曲线或RBF(Radial Basis Function)进行曲线拟合大行其道的今天,最小二乘法在多项式曲线或直线的拟合问题上,仍然得到广泛地应用。使用最小二乘法,选取的匹配函数的模式非常重要,如果离散数据呈现的是指数变化规律,则应该选择指数形式的匹配函数模式,如果是多项式变化规律,则应该选择多项式匹配模式,如果选择的模式不对,拟合的效果就会很差,这也是使用最小二乘法进行曲线拟合时需要特别注意的一个地方。

            下面以多项式模式为例,介绍一下使用最小二乘法进行曲线拟合的完整步骤。假设选择的拟合多项式模式是:

    这m个等式相当于m个方程,a0,a1,…amm个未知量,因此这m个方程组成的方程组是可解的,最小二乘法的第二步处理就是将其整理为针对a0,a1,…am的正规方程组。最终整理的方程组如下:

    最小二乘法的第三步处理就是求解这个多元一次方程组,得到多项式的系数a0,a1,…am,就可以得到曲线的拟合多项式函数。求解多元一次方程组的方法很多,高斯消元法是最常用的一种方法,下一节就简单介绍一下最小二乘算法实现所用的高斯消元法算法。

    12.2.2 高斯消元法求解方程组

            在数学上,高斯消元法是线性代数中的一个算法,可用来求解多元一次线性方程组,也可以用来求矩阵的秩,以及求可逆方阵的逆矩阵。高斯消元法虽然以数学家高斯的名字命名,但是最早出现在文献资料中应该是中国的《九章算术》。

            高斯消元法的主要思想是通过对系数矩阵进行行变换,将方程组的系数矩阵由对称矩阵变为三角矩阵,从而达到消元的目的,最后通过回代逐个获得方程组的解。在消元的过程中,如果某一行的对角线元素的值太小,在计算过程中就会出现很大的数除以很小的数的情况,有除法溢出的可能,因此在消元的过程中,通常都会增加一个主元选择的步骤,通过行交换操作,将当前列绝对值最大的行交换到当前行位置,避免了除法溢出问题,增加了算法的稳定性。

            高斯消元法算法实现简单,主要有两个步骤组成,第一个步骤就是通过选择主元,逐行消元,最终行程方程组系数矩阵的三角矩阵形式,第二个步骤就是逐步回代的过程,最终矩阵的对角线上的元素就是方程组的解。下面就给出高斯消元法的一个算法实现:

     76 /*带列主元的高斯消去法解方程组,最后的解在matrixA的对角线上*/

     77 bool GuassEquation::Resolve(std::vector<double>& xValue)

     78 {

     79     assert(xValue.size() == m_DIM);

     80 

     81     /*消元,得到上三角阵*/

     82     for(int i = 0; i < m_DIM - 1; i++)

     83     {

     84         /*按列选主元*/

     85         int pivotRow = SelectPivotalElement(i);

     86         if(pivotRow != i)/*如果有必要,交换行*/

     87         {

     88             SwapRow(i, pivotRow);

     89         }

     90         if(IsPrecisionZero(m_matrixA[i * m_DIM + i])) /*主元是0? 不存在唯一解*/

     91         {

     92             return false;

     93         }

     94         /*对系数归一化处理,使行第一个系数是1.0*/

     95         SimplePivotalRow(i, i);

     96         /*逐行进行消元*/

     97         for(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_DIM - 1) * m_DIM + m_DIM - 1];

    104     for(int i = m_DIM - 2; i >= 0; i--)

    105     {

    106         double totalCof = 0.0;

    107         for(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     }

    113 

    114     /*将对角线元素的解逐个存入解向量*/

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

    116     {

    117         xValue[i] = m_matrixA[i * m_DIM + i];

    118     }

    119 

    120     return true;

    121 }

     

            GuassEquation::Resolve()函数中m_matrixA是以一维数组形式存放的系数矩阵,m_DIM是矩阵的维数,SelectPivotalElement()函数从系数矩阵的第i列中选择绝对值最大的那个值所在的行,并返回行号,SwapRow()函数负责交换系数矩阵两个行的所有值,SimplePivotalRow()函数是归一化处理函数,通过除法操作将指定的行的对角线元素变换为1.0,以便简化随后的消元操作。

    12.2.3 最小二乘法解决“速度与加速度”实验

            根据12.2.1节对最小二乘法原理的分析,用程序实现最小二乘法曲线拟合的算法主要由两个步骤组成,第一个步骤就是根据给出的测量值生成关于拟合多项式系数的方程组,第二个步骤就是解这个方程组,求出拟合多项式的各个系数。根据对上文最终整理的正规方程组的分析,可以看出其系数有一定的关系,就是每一个方程式都比前一个方程式多乘了一个xi。因此,只需要完整计算出第一个方程式的系数,其他方程式的系数只是将前一个方程式的系数依次左移一位,然后单独计算出最后一个系数就可以了,此方法可以减少很多无谓的计算。求解多元一次方程组的方法就使用12.2.2节介绍的高斯消元法,其算法上一节已经给出。

            这里给出一个最小二乘算法的完整实现,以12.1.2节的数据为例,因为数据结果明显呈现线性方程的特征,因此选择拟合多项式为v = v0 + at,v0和a就是要求解的拟合多项式系数。

     99 bool LeastSquare(const std::vector<double>& x_value, const std::vector<double>& y_value,

    100                  int M, std::vector<double>& a_value)

    101 {

    102     assert(x_value.size() == y_value.size());

    103     assert(a_value.size() == M);

    104 

    105     double *matrix = new double[M * M];

    106     double *b= new double[M];

    107 

    108     std::vector<double> x_m(x_value.size(), 1.0);

    109     std::vector<double> y_i(y_value.size(), 0.0);

    110     for(int i = 0; i < M; i++)

    111     {

    112         matrix[ARR_INDEX(0, i, M)] = std::accumulate(x_m.begin(), x_m.end(), 0.0);

    113         for(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);

    118         for(int k = 0; k < static_cast<int>(x_m.size()); k++)

    119         {

    120             x_m[k] *= x_value[k];

    121         }

    122     }

    123     for(int row = 1; row < M; row++)

    124     {

    125         for(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);

    130         for(int k = 0; k < static_cast<int>(x_m.size()); k++)

    131         {

    132             x_m[k] *= x_value[k];

    133         }

    134     }

    135 

    136     GuassEquation equation(M, matrix, b);

    137     delete[] matrix;

    138     delete[] b;

    139 

    140     return equation.Resolve(a_value);

    141 }

            将表121的数据带入算法,计算得到v0 = 4.05545455,a = 1.48818182,比作图法得到的结果更精确。以上算法是根据最小二乘法的理论推导系数方程,并求解系数方程得到拟合多项式的系数的一种实现方法,除此之外,还可以利用预先计算好的最小二乘解析理论直接求得拟合多项式的系数,读者可自行学习相关的实现算法。

     

    12.3 三次样条曲线拟合

            曲线拟合基本上就是一个插值计算的过程,除了最小二乘法,其他插值方法也可以被用于曲线拟合。常用的曲线拟合方法还有基于RBFRadial Basis Function)的曲线拟合和三次样条曲线拟合。最小二乘法方法简单,便于实现,但是如果拟合模式选择不当,会产生较大的偏差,特别是对于复杂曲线的拟合,如果选错了模式,拟合的效果就很差。基于RBFRadial Basis Function)的曲线拟合方法需要高深的数学基础,涉及多维空间理论,将低维的模式输入数据转换到高维空间中,使得低维空间内的线性不可分问题在高维空间内变得线性可分,这种数学分析方法非常强大,但是这种方法不宜得到拟合函数,因此在需要求解拟合函数的情况下使用起来不是很方便。

            样条插值是一种工业设计中常用的、得到平滑曲线的一种插值方法,三次样条又是其中用的较为广泛的一种。使用三次样条曲线进行曲线拟合可以得到非常高精度的拟合结果,并且很容易得到拟合函数,本节的内容将重点介绍三次样条曲线拟合的原理和算法实现,并通过一个具体的例子将三次样条函数拟合的曲线与原始曲线对比显示,体会一下三次样条曲线拟合的惊人效果。

     

    12.3.1 插值函数

            前文提到过,曲线拟合的实质就是各种插值计算,因此,插值函数的选择决定了曲线拟合的效果。那么插值函数的数学定义是什么呢?若在[a, b]上给出n + 1个点a ≤ x0 < x1 < ⋯ < xn ≤ b,f(x) 是[a, b]上的实值函数, 要求一个具有n + 1个参量的函数s(x; a0,...,an) 使它满足

            s(xi; a0,..., an) = f(xi) , i = 0, 1, ⋯, n                          (3.5)

    则称s(x)为f(x) 在[a, b]上的插值函数. 若s(x) 关于参量a0, a1,...,an是线性关系, 即:

            s(x) = a0s0(x) + a1s1(x) + ⋯ + ansn(x)                                (3.6)

    s(x)就是多项式插值函数,如果si(x)是三角函数,则s(x)就是三角插值函数。

            比较常用的多项式插值函数是牛顿插值多项式和拉格朗日插值多项式,但是在多项式的次数比较高的情况下,插值点数n过多会导致多项式插值在收敛性和稳定性上失去保证,因此,当插值点数n较大的情况下,一般不使用多项式插值,而采用样条插值或次数较低的最小二乘法插值。

    12.3.2 样条函数的定义

            在所有能够保证收敛性和稳定性的插值函数中,最常用的,也是最重要的插值函数就是样条插值函数。采用样条函数计算出的插值曲线和曲面在飞机、轮船和汽车等精密机械设计中都得到了广泛的应用。样条插值函数的数学定义是这样的:

    设区间[a, b]上选取n - 1个节点(包括区间端点a和b共n + 1个节点),将其划分为n个子区间 a = x0 < x1 < ⋯ < xn = b, 如果存在函数s(x),使得s(x)满足以下两个条件:

    (1) s(x) 在整个区间[a, b]上具有m - 1阶连续导数;

    (2) s(x)在每个子区间[xi-1, xi], i= 1, 2, ⋯, n 上是m 次代数多项式(最高次数为m次);

    则称s(x)是区间[a, b]上的m 次样条函数。假如区间[a, b]上存在实值函数f(x),使得每个节点处的值f(xi)与s(xi)相等,即

            s(xi) = f(xi), i = 0, 1, …, n                                          (3.7)

    则称s(x)是实值函数f(x)的m 次样条插值函数。

            当m = 1时,样条插值函数就是分段线性插值, 此时虽然s(x)是属于区间[a, b]上的函数, 但它不光滑(连一阶连续导数性质都不具备),不能满足工程设计要求。工程设计通常使用较多的是m = 3时的三次样条插值函数,此时样条函数具有二阶连续导数性质。

            根据三次样条函数的定义,s(x)在每个子区间上的样条函数si(x)都是一个三次多项式,也就是说,三次样条函数s(x)由n个区间上的n个三次多项式组成,每个三次多项式可描述为以下形式:

            si(x) = aix3 + bix2 + cix + di      i = 1, 2, …, n                      (3.8)

    因此,要确定完整的样条函数s(x)需要确定ai、bi、ci和di公4n个系数。根据样条函数的定义,s(x)在区间内的n - 1个节点处都是连续的,并且其一阶导数si‘(x)和二阶导数si“(x)都是连续的,根据连续函数的性质(xi的左右导数相等),我们可以得到3(n - 1)个条件:

            si(xi - 0) = si+1(xi + 0)    i = 1, 2, …, n-1

            si‘(xi - 0) = si+1‘(xi + 0)    i = 1, 2, …, n-1

            si“(xi - 0) = si+1“(xi + 0)    i = 1, 2, …, n-1                          (3.9)

    再加上插值函数在包括区间端点a(就是x0),b(就是xn)在内的n + 1个节点处满足s(xi) = f(xi),又可以得到n + 1个条件,这样就具备了4n – 2个条件。

    12.3.3 边界条件

            为了解决4n个系数组成的方程组,最终确定的s(x),需要再补充两个边界条件使之满足4n个条件。常用的边界条件有以下几种:

            第一类边界条件,即满足s‘(x0) = f’(x0),s‘(xn) = f’(xn)两个条件,其中f(x)是实值函数。

            第二类边界条件,即满足s”(x0) = f”(x0),s”(xn) = f”(xn)两个条件,其中f(x)是实值函数。特别情况下,当f”(x0) = f”(xn) = 0的时候,也就是s”(x0) = s”(xn) = 0的情况下,第二类边界条件又被称为自然边界条件。

            当样条函数的实值函数f(x)是以[a, b]为周期的周期函数时,三次样条函数s(x)在两个端点处满足s‘(x0 - 0) = s‘(xn + 0)和s”(x0 - 0) = s”(xn + 0),这种情况又被成为第三类边界条件。

            工程技术中常用的是第一类边界条件和第二类边界条件,以及第二类边界条件的特殊情况自然边界条件。理想情况下,也就是实值函数已知的情况下,可以通过实值函数直接计算出边界条件的值,否则的话,就只能通过测量和计算得到边界条件的值,有时候甚至只能给出经验估计值,工程技术中通常根据实际情况灵活使用各类边界条件。

    12.3.4 推导三次样条函数

            求三次样条插值函数s(x)的方法很多,其基本原理都是首先求出由待定系数组成的s(x),以及其一阶导数s’(x)和二阶导数s”(x),然后将其带入到12.3.2和12.3.3节列举的4n个条件中,得到关于待定系数的方程组,最后求解方程组得到待定系数,并最终确定插值函数s(x)。

            求三次样条插值函数s(x)常用的方法是“三转角法”和“三弯矩法”。根据三次样条函数的性质,s(x)的一阶导数s’(x)是二次多项式, 二阶导数s”(x)是一次多项式(线性函数), “三转角法”和“三弯矩法”的主要区别是利用这两个特性推导插值函数s(x)、s’(x)和s”(x)的方式不同。“三转角法”利用s(x)的一阶导数s’(x)是二次多项式这个特性,对于子区间[xi, xi+1],利用抛物线插值公式获得一个通过xi和xi+1两个点的二次多项式做为s’(x),然后对s’(x)进行积分和微分(求导)运算,分别得到s(x),和s”(x),最后将它们带入4n个条件中求解系数方程组。“三弯矩法”则是利用s(x)的二阶导数s”(x)是一次多项式(线性函数)这个特性,对于子区间[xi, xi+1],首先假设一个通过xi和xi+1两个点的线性函数做为s”(x),然后对s”(x)进行连续两次积分运算得到s(x),再对s(x)进行求导得运算到s’(x),最后将它们带入4n个条件中求解系数方程组。这两种方法的本质是一样的,只是对s(x)的推导过程不同,接下来就介绍使用“三弯矩法”求解三次样条函数的方法。

            三次样条函数的求解过程就是系数方程组的推导过程,使用“三弯矩法”推导系数方程组,首先要确定插值函数的二阶导数s”(x)。根据三次样条函数的性质,在每个子区间[xi, xi+1]上,其二阶导数s”(x)是个线性方程,现在假设在xi和xi+1两个端点的二阶导数值分别是Mi和Mi+1,也就是s”(xi) = Mi,s”(xi+1) = Mi+1,则经过xi和xi+1的两点式直线方程是:

    从M0到Mn,有n+1个Mi的值需要求解,但是(3.20)只有n-1个等式,此时就需要用到两个边界条件了。

            如果是使用第二类边界条件,则直接可以得到以下两个条件等式:

            s”(x0) = M0 = f”(x0) = y0’                                                (3.21)

            s”(xn) = Mn = f”(xn) = yn’                                                (3.22)

    令d0 = 2y0’,dn = 2yn’,可以得到由第二类边界条件确定的两个方程:

            2M0 = d0                                                                    (3.23)

            2Mn = dn                                                                    (3.24)

            如果是使用第一类边界条件,即s‘(x0) = f’(x0),s‘(xn) = f’(xn)两个条件,则需要将这两个条件代入(3.16),通过计算得到两个条件等式。将s‘(x0) = y0’代入(3.16),得到:

    将第二类边界条件得到的(3.23)和(3.24)或第一类边界条件得到的(3.27)和(3.28)与(3.20)中的n – 1个等式组合在一起就得到一个关于Mi的方程组,求解此方程组可以得到Mi的值,代入到(3.15)即可得到三次样条函数方程。以第一类边界条件得到的(3.27)和(3.28)为例,与(3.20)连立得到以下方程组:

    这就是三弯矩方程组,其中Mi,i=0,1,…n就是三次样条函数s(x)的矩。根据(3.27)和(3.28),un = 1,v0 = 1,其余各系数可以通过(3.19)中的系数计算出来。这个方程组的系数矩阵是一个对角线矩阵,并且是一个严格对角占优的对角阵(ui和vi的值均小于主对角线的值,也就是ui和vi的值皆小于2),可以使用追赶法求解。下一节将介绍如何使用追赶法求解方程组,并给出求解的算法实现。

    12.3.5 追赶法求解方程组

            任意矩阵A,都可以通过克洛脱(Crout)分解得到两个三角矩阵:

    其中y1=d1/l1,其余各项的递推计算关系是:

    yi = (di – miyi-1)/li, i = 2,3,… ,n

    对于第二个方程,求解最终结果xi

    其中xn = yn,其余各项的递推求解关系是:

    xi = yi uixi+1, i = n 1, n 2, , 1

    递推计算yi和xi的过程分别被形象地形容为追的过程和赶的过程,这也是追赶法得名的原因,实际上这种方法在国际上叫做托马斯(Thomas)法。在这里需要强调一下,对三角矩阵的克洛脱分解需要满足几个条件,否则无法进行,这几个条件分别是:

            下面就给出一个追赶法求解方程组的通用算法实现,在使用之前需要判断系数矩阵是否是对三角矩阵,并且满足上述三个条件,相关的判断请读者自行添加:

     76 /*追赶法求对三角矩阵方程组的解*/

     77 bool ThomasEquation::Resolve(std::vector<double>& xValue)

     78 {

     79     assert(xValue.size() == m_DIM);

     80 

     81     std::vector<double> L(m_DIM);

     82     std::vector<double> M(m_DIM);

     83     std::vector<double> U(m_DIM);

     84     std::vector<double> Y(m_DIM);

     85 

     86     /*消元,追的过程*/

     87     L[0] = m_matrixA[ARR_INDEX(0, 0, m_DIM)];

     88     U[0] = m_matrixA[ARR_INDEX(0, 1, m_DIM)] / L[0];

     89     Y[0] = m_bVal[0] / L[0];

     90     for(int i = 1; i < m_DIM; i++)

     91     {

     92         if(IsPrecisionZero(m_matrixA[ARR_INDEX(i, i, m_DIM)]))

     93         {

     94             return false;

     95         }

     96         M[i] = m_matrixA[ARR_INDEX(i, i - 1, m_DIM)];

     97         L[i] = m_matrixA[ARR_INDEX(i, i, m_DIM)] - M[i] * U[i - 1];

     98         Y[i] = (m_bVal[i] - M[i] * Y[i - 1]) / L[i];

     99         if((i + 1) < m_DIM)

    100         {

    101             U[i] = m_matrixA[ARR_INDEX(i, i + 1, m_DIM)] / L[i];

    102         }

    103     }

    104     /*回代求解,赶的过程*/

    105     xValue[m_DIM - 1] = Y[m_DIM - 1];

    106     for(int i = m_DIM - 2; i >= 0; i--)

    107     {

    108         xValue[i] = Y[i] - U[i] * xValue[i + 1];

    109     }

    110 

    111     return true;

    112 }

     

    12.3.6 三次样条曲线拟合算法实现

            根据12.3.4节对三次样条函数的推导分析,三次样条曲线拟合算法的核心可分为三部分,第一部分是根据推导结果计算关于三次样条函数的“矩”的方程组的系数矩阵,第二部分就是用追赶法求解方程组,得到各个区间的三次样条函数,第三部分就是根据每个拟合点的输入的值xi,确定使用哪个区间的三次样条函数,并使用该区间的三次样条函数计算出三次样条插值yi,最后得到的一系列(xi, yi)组成的曲线就是三次样条拟合曲线。拟合算法也是按照上面的分析,分三个步骤计算插值:

            第一步是计算系数矩阵,其中u0、v0、d0和dn的值需要单独计算,其余的值可以通过(3.19)递推计算出来。

            第二步是将系数矩阵代入12.3.5节给出的追赶法通用算法,求出Mi的值。求解之前,先证明一下第一步得到系数矩阵是否满足追赶法的条件。首先,主对角线元素的值都是2,满足12.3.5节的条件(1)。其次,由ui和vi的计算条件可知,|ui| < 1,|vi| < 1,这也满足12.3.5节的条件(2)。最后,因为ai = 2,且ui和vi的和是1,所以12.3.5节的条件(3)也得到了满足。由上判断可知,求解三次样条函数的“矩”的系数矩阵满足使用追赶法求解的条件,可以使用追赶法求解。

            第三步是计算插值,需要将第二步计算得到的Mi代入(3.15),并选择合适的子区间样条函数计算出插值点的值。

            下面就给出采用三弯矩法实现的三次样条曲线拟合算法,CalcSpline()函数的参数Xi和Yi是n个插值点(包括起点和终点)的值,boundType是边界条件类型,b1和b2分别是对应的两个边界条件,这个算法支持第一类和第二类边界条件(包括自然边界条件)。内部的矩阵matrixA就是按照(3.29)构造的Mi方程组的系数矩阵,可用于直接用追赶法求解方程组。CalcSpline()函数的大部分代码都是在构造Mi方程组的系数矩阵,首先根据边界条件确定un 、v0 、d0和dn,其他系数则根据(3.19)的递推关系,在for(int i = 1; i < (m_valN - 1); i++)循环中依次计算出来,最后是利用12.3.5节给出的追赶法算法求出Mi。GetValue()函数负责计算给定区间内任意位置的插值,首先根据x的值确定使用哪个子区间的样条函数,然后根据(3.12)和(3.14)给出的关系计算插值。

     16 void SplineFitting::CalcSpline(double *Xi, double *Yi, int n, int boundType, double b1, double b2)

     17 {

     18     assert((boundType == 1) || (boundType == 2));

     19 

     20     double *matrixA = new double[n * n];

     21     if(matrixA == NULL)

     22     {

     23         return;

     24     }

     25     double *d = new double[n];

     26     if(d == NULL)

     27     {

     28         delete[] matrixA;

     29         return;

     30     }

     31 

     32     m_valN = n;

     33     m_valXi.assign(Xi, Xi + m_valN);

     34     m_valYi.assign(Yi, Yi + m_valN);

     35     m_valMi.resize(m_valN);

     36     memset(matrixA, 0, sizeof(double) * n * n);

     37 

     38     matrixA[ARR_INDEX(0, 0, m_valN)] = 2.0;

     39     matrixA[ARR_INDEX(m_valN - 1, m_valN - 1, m_valN)] = 2.0;

     40     if(boundType == 1) /*第一类边界条件*/

     41     {

     42         matrixA[ARR_INDEX(0, 1, m_valN)] = 1.0; //v0

     43         matrixA[ARR_INDEX(m_valN - 1, m_valN - 2, m_valN)] = 1.0; //un

     44         double h0 = Xi[1] - Xi[0];

     45         d[0] = 6 * ((Yi[1] - Yi[0]) / h0 - b1) / h0; //d0

     46         double hn_1 = Xi[m_valN - 1] - Xi[m_valN - 2];

     47         d[m_valN - 1] = 6 * (b2 - (Yi[m_valN - 1] - Yi[m_valN - 2]) / hn_1) / hn_1; //dn

     48     }

     49     else /*第二类边界条件*/

     50     {

     51         matrixA[ARR_INDEX(0, 1, m_valN)] = 0.0; //v0

     52         matrixA[ARR_INDEX(m_valN - 1, m_valN - 2, m_valN)] = 0.0; //un

     53         d[0] = 2 * b1; //d0

     54         d[m_valN - 1] = 2 * b2; //dn

     55     }

     56     /*计算ui,vi,di,i = 2,3,...,n-1*/

     57     for(int i = 1; i < (m_valN - 1); i++)

     58     {

     59         double hi_1 = Xi[i] - Xi[i - 1];

     60         double hi = Xi[i + 1] - Xi[i];

     61         matrixA[ARR_INDEX(i, i - 1, m_valN)] = hi_1 / (hi_1 + hi); //ui

     62         matrixA[ARR_INDEX(i, i, m_valN)] = 2.0;

     63         matrixA[ARR_INDEX(i, i + 1, m_valN)] = 1 - matrixA[ARR_INDEX(i, i - 1, m_valN)]; //vi = 1 - ui

     64         d[i] = 6 * ((Yi[i + 1] - Yi[i]) / hi - (Yi[i] - Yi[i - 1]) / hi_1) / (hi_1 + hi); //di

     65     }

     66 

     67     ThomasEquation equation(m_valN, matrixA, d);

     68     equation.Resolve(m_valMi);

     69     m_bCalcCompleted = true;

     70 

     71     delete[] matrixA;

     72     delete[] d;

     73 }

     74 

     75 double SplineFitting::GetValue(double x)

     76 {

     77     if(!m_bCalcCompleted)

     78     {

     79         return 0.0;

     80     }

     81     if((x < m_valXi[0]) || (x > m_valXi[m_valN - 1]))

     82     {

     83         return 0.0;

     84     }

     85     int i = 0;

     86     for(i = 0; i < (m_valN - 1); i++)

     87     {

     88         if((x >= m_valXi[i]) && (x < m_valXi[i + 1]))

     89             break;

     90     }

     91     double hi = m_valXi[i + 1] - m_valXi[i];

     92     double xi_1 = m_valXi[i + 1] - x;

     93     double xi = x - m_valXi[i];

     94 

     95     double y = xi_1 * xi_1 * xi_1 * m_valMi[i] / (6 * hi);

     96     y += (xi * xi * xi * m_valMi[i + 1] / (6 * hi));

     97 

     98     double Ai = (m_valYi[i + 1] - m_valYi[i]) / hi - (m_valMi[i + 1] - m_valMi[i]) * hi / 6.0;

     99     y += Ai * x;

    100     double Bi = m_valYi[i + 1] - m_valMi[i + 1] * hi * hi / 6.0 - Ai * m_valXi[i + 1];

    101     y += Bi;

    102     return y;

    103 }

    12.3.7 三次样条曲线拟合的效果

    本节将用定义一个原始函数,从原始函数的某个区间上抽取9个插值点,根据这9个插值点和原函数的边界条件,利用三次样条曲线插值进行曲线拟合,并将原始曲线和拟合曲线做对比,展示一下三次样条曲线拟合的效果。

    首先定义原始函数:

    f(x) = 3/(1 + x2)

    选择区间[0.0, 8.0]上的9个点做为插值点,计算各点的值为:

    x

    0.0

    1.0

    2.0

    3.0

    4.0

    5.0

    6.0

    7.0

    8.0

    y

    3.0

    1.5

    0.6

    0.3

    0.1765

    0.1154

    0.0811

    0.06

    0.0462

    12 – 2 原函数f(x)在各插值点的值

    求f(x)的导函数f’(x):

    f’(x) = -6x/(1 + x2)2

    根据f’(x)计算出在区间端点处的两个第一类边界条件f’(0.0) = 0.0, f’(8.0) = -0.01136。利用表12-2中的数据和这两个边界条件,计算出三次样条插值函数,并从0.0开始,以0.01为步长,连续求800个点的插值,将这些点连成曲线得到拟合曲线。为了作对比,同样从0.0开始,以0.01为步长,用f(x)函数连续计算800个点的原值,将这些点连成曲线得到原始曲线。分别用不同的颜色画出这两条曲线,如图12-2所示:

    图12-2 拟合曲线和原始曲线对比

            从图12-2可以看到,三次样条曲线拟合的效果非常好。同样在[0.0, 8.0]区间上,如果增加插值点的个数,将获得更好的拟合效果。比如以0.5为单位,将插值点增加到17个,则拟合的曲线与原始曲线几乎完全重合。

     

    展开全文
  • matlab 万能实用的非线性曲线拟合方法

    万次阅读 多人点赞 2018-08-13 09:43:41
    在科学计算和工程应用中,经常会遇到需要拟合一系列的离散数据,最近找了很多相关的文章方法,在这里进行总结一下其中最完整、几乎能解决所有离散参数非线性拟合的方法   第一步:得到散点数据 根据你的实际...

            在科学计算和工程应用中,经常会遇到需要拟合一系列的离散数据,最近找了很多相关的文章方法,在这里进行总结一下其中最完整、几乎能解决所有离散参数非线性拟合的方法


     

    第一步:得到散点数据

    根据你的实际问题得到一系列的散点

    例如:

    x=[3.2,3.6,3.8,4,4.2,4.8,5,5.4,6.2,6.4,6.6,6.9,7.1]';%加上一撇表示对矩阵的转置
    y=[0.38,0.66,1,0.77,0.5,0.66,0.83,1,0.71,0.71,1,0.87,0.83]';
    


    第二步:确定函数模型

    根据上述的实际散点确定应该使用什么样的曲线,或者说是想要模拟的曲线

    t=[3.2,3.6,3.8,4,4.2,4.8,5,5.4,6.2,6.4,6.6,6.9,7.1]';
    tt=[0.38,0.66,1,0.77,0.5,0.66,0.83,1,0.71,0.71,1,0.87,0.83]';
     
    plot(t,tt,'.');%得到散点图
    

    散点图如下所示:

    我们已知现存的几种典型的(也是绝大多数情况下的函数模型)



     

    选定一个与散点图像相匹配的函数模型,在此例中我们选择典型的S型曲线模型y= 1/(a+b*e^(-x)),其实此处的函数模型可以任意。


     

    第三步:确定选用函数模型中的未知参数

          首先了解一下matlab中的inline函数,inline是用来定义内联函数的
    比如说:

     

    y=inline('sin(x)','x') %第一个参数是表达式,第二个参数是函数变量
    y(0) %计算sin(0)的值
    y(pi) %计算sin(pi)的值
    q=quad(y,0,1); %计算sin(x) 在0到1上的积分
    

           之后,我们在代码中进行函数的定义

    x=[3.2,3.6,3.8,4,4.2,4.5,4.8,5,5.3,5.4,5.6,5.8,6,6.2,6.4,6.6,6.9,7.1]';
    y=[0.38,0.66,1,0.77,0.5,0.33,0.66,0.83,0.33,1,0.33,0.5,0.33,0.71,0.71,1,0.87,0.83]';
     
    myfunc = inline('1./(beta(1)+beta(2).*exp(-x))','beta','x');%三个参数分别为:函数模型(注意需要使用点除和点乘),待定系数,自变量
     
    beta0 = [0.2,0.2]';%待定系数的预估值
    beta = nlinfit(x,y,myfunc,beta0);%
    



     

    其中,beta返回了非线性拟合之后的待定系数,beta(1)和beta(2)表示待定系数,可以为任意数量的扩展beta(n),也就说明了选择函数模型的自由性,甚至可以有100个参数!

    beta0表示的是函数模型中待定系数的预估值,可以任意设定


     

            matlab 中的nlinfit(x,y,f,a)函数:用于拟合非线性表达式的函数
    f:符号函数句柄,如果是以m文件的形式调用的时候,别忘记加@.这里需要注意,f函数的返回值是和y匹对的,即拟合参数的标准是(f-y)^2取最小值,具体看下面的例子
    a:最开始预估的值(预拟合的未知参数的估计值)。如上面的问题如果我们预估A为1,B为2,则a=[1 2]
    x:我们已经获知的x的值
    y:我们已经获知的x对应的y的值(
    这部分不懂的在matlab中help命令进行了解


     

    求解出beta的大小:beta(1) = 1.1562  beta(2) = 15.1875;


     

    画图:使用plot()函数

    t=[3.2,3.6,3.8,4,4.2,4.8,5,5.4,6.2,6.4,6.6,6.9,7.1]';
    tt=[0.38,0.66,1,0.77,0.5,0.66,0.83,1,0.71,0.71,1,0.87,0.83]';
     
    plot(t,tt,'.');
    hold on%保证同时显示
     
    x = 0:0.01:8;
    y = 1./(1.1562+15.1875.*exp(-x));
     
    plot(x,y);
    



     

    结果:是不是很棒!~,另外可以自行加上对应的横纵坐标内容,这里就不多说了。


     

    总结一下matlab非线性拟合散点图的过程:得到散点数据=>确定函数模型=>求解函数模型的待定系数=>得到拟合函数的具体形式=>画出拟合图像

    =========================================================================================

    用过Matlab的拟合、优化和统计等工具箱的网友,会经常遇到下面几个名词:

    SSE(和方差、误差平方和):The sum of squares due to error
    MSE(均方差、方差):Mean squared error
    RMSE(均方根、标准差):Root mean squared error
    R-square(确定系数):Coefficient of determination
    Adjusted R-square:Degree-of-freedom adjusted coefficient of determination

    下面我对以上几个名词进行详细的解释下,相信能给大家带来一定的帮助!!

    一、SSE(和方差)
    该统计参数计算的是拟合数据和原始数据对应点的误差的平方和,计算公式如下

     

    SSE越接近于0,说明模型选择和拟合更好,数据预测也越成功。接下来的MSE和RMSE因为和SSE是同出一宗,所以效果一样

    二、MSE(均方差)
    该统计参数是预测数据和原始数据对应点误差的平方和的均值,也就是SSE/n,和SSE没有太大的区别,计算公式如下

     

    三、RMSE(均方根)
    该统计参数,也叫回归系统的拟合标准差,是MSE的平方根,就算公式如下

     

    在这之前,我们所有的误差参数都是基于预测值(y_hat)和原始值(y)之间的误差(即点对点)。从下面开始是所有的误差都是相对原始数据平均值(y_ba)而展开的(即点对全)!!!

    四、R-square(确定系数)
    在讲确定系数之前,我们需要介绍另外两个参数SSR和SST,因为确定系数就是由它们两个决定的
    (1)SSR:Sum
    of squares of the regression,即预测数据与原始数据均值之差的平方和,公式如下

    (2)SST:Total sum of squares,即原始数据和均值之差的平方和,公式如下

    细心的网友会发现,SST=SSE+SSR,呵呵只是一个有趣的问题。而我们的“确定系数”是定义为SSR和SST的比值,故

     

    其实“确定系数”是通过数据的变化来表征一个拟合的好坏。由上面的表达式可以知道“确定系数”的正常取值范围为[0
    1],越接近1,表明方程的变量对y的解释能力越强,这个模型对数据拟合的也较好。


    展开全文
  • 曲线拟合

    万次阅读 2016-08-29 10:26:14
    12.1 曲线拟合 12.1.1 曲线拟合的定义  曲线拟合(Curve Fitting)的数学定义是指用连续曲线近似地刻画或比拟平面上一组离散点所表示的坐标之间的函数关系,是一种用解析表达式逼近离散数据的方法。曲线拟合...

    12.1 曲线拟合

    12.1.1 曲线拟合的定义

            曲线拟合(Curve Fitting)的数学定义是指用连续曲线近似地刻画或比拟平面上一组离散点所表示的坐标之间的函数关系,是一种用解析表达式逼近离散数据的方法。曲线拟合通俗的说法就是“拉曲线”,也就是将现有数据透过数学方法来代入一条数学方程式的表示方法。科学和工程遇到的很多问题,往往只能通过诸如采样、实验等方法获得若干离散的数据,根据这些数据,如果能够找到一个连续的函数(也就是曲线)或者更加密集的离散方程,使得实验数据与方程的曲线能够在最大程度上近似吻合,就可以根据曲线方程对数据进行数学计算,对实验结果进行理论分析,甚至对某些不具备测量条件的位置的结果进行估算。

    12.1.2 简单线性数据拟合的例子

            回想一下中学物理课的“速度与加速度”实验:假设某物体正在做加速运动,加速度未知,某实验人员从时间t0 = 3秒时刻开始,以1秒时间间隔对这个物体连续进行了12次测速,得到一组速度和时间的离散数据,请根据实验结果推算该物体的加速度。

    时间 ()

    3

    4

    5

    6

    7

    8

    9

    10

    速度(米/秒)

    8.41

    9.94

    11.58

    13.02

    14.33

    15.92

    17.54

    19.22

    时间 ()

    11

    12

    13

    14

     

     

     

     

    速度(米/秒)

    20.49

    22.01

    23.53

    24.47

     

     

     

     

     12 – 1 物体速度和时间的测量关系表

            在选择了合适的坐标刻度之后,我们就可以在坐标纸上画出这些点。如图121所示,排除偏差明显偏大的测量值后,可以看出测量结果呈现典型的线性特征。沿着该线性特征画一条直线,使尽量多的测量点能够位于直线上,或与直线的偏差尽量小,这条直线就是我们根据测量结果拟合的速度与时间的函数关系。最后在坐标纸上测量出直线的斜率KK就是被测物体的加速度,经过测量,我们实验测到的物体加速度值是1.48米/2

     

     12 – 1 实验法测量加速度的过程

     

    12.2 最小二乘法曲线拟合

            使用数学分析进行曲线拟合有很多常用的方法,这一节我们先介绍一下最简单的最小二乘法,并给出使用最小二乘法解决上一节给出的速度与加速度实验问题。

     

    12.2.1 最小二乘法原理

            最小二乘法(又称最小平方法)通过最小化误差的平方和寻找数据的最佳函数匹配,利用最小二乘法可以简便地求得未知的数据,并使得这些求得的数据与实际数据之间误差的平方和为最小,当然,做为一种插值方法使用时,最小二乘法也可以用于曲线拟合。使用最小二乘法进行曲线拟合是曲线拟合种早期的一种常用方法,不过,最小二乘法理论简单,计算量小,即便是在使用三次样条曲线或RBF(Radial Basis Function)进行曲线拟合大行其道的今天,最小二乘法在多项式曲线或直线的拟合问题上,仍然得到广泛地应用。使用最小二乘法,选取的匹配函数的模式非常重要,如果离散数据呈现的是指数变化规律,则应该选择指数形式的匹配函数模式,如果是多项式变化规律,则应该选择多项式匹配模式,如果选择的模式不对,拟合的效果就会很差,这也是使用最小二乘法进行曲线拟合时需要特别注意的一个地方。

            下面以多项式模式为例,介绍一下使用最小二乘法进行曲线拟合的完整步骤。假设选择的拟合多项式模式是:

    这m个等式相当于m个方程,a0,a1,…amm个未知量,因此这m个方程组成的方程组是可解的,最小二乘法的第二步处理就是将其整理为针对a0,a1,…am的正规方程组。最终整理的方程组如下:

    最小二乘法的第三步处理就是求解这个多元一次方程组,得到多项式的系数a0,a1,…am,就可以得到曲线的拟合多项式函数。求解多元一次方程组的方法很多,高斯消元法是最常用的一种方法,下一节就简单介绍一下最小二乘算法实现所用的高斯消元法算法。

    12.2.2 高斯消元法求解方程组

            在数学上,高斯消元法是线性代数中的一个算法,可用来求解多元一次线性方程组,也可以用来求矩阵的秩,以及求可逆方阵的逆矩阵。高斯消元法虽然以数学家高斯的名字命名,但是最早出现在文献资料中应该是中国的《九章算术》。

            高斯消元法的主要思想是通过对系数矩阵进行行变换,将方程组的系数矩阵由对称矩阵变为三角矩阵,从而达到消元的目的,最后通过回代逐个获得方程组的解。在消元的过程中,如果某一行的对角线元素的值太小,在计算过程中就会出现很大的数除以很小的数的情况,有除法溢出的可能,因此在消元的过程中,通常都会增加一个主元选择的步骤,通过行交换操作,将当前列绝对值最大的行交换到当前行位置,避免了除法溢出问题,增加了算法的稳定性。

            高斯消元法算法实现简单,主要有两个步骤组成,第一个步骤就是通过选择主元,逐行消元,最终行程方程组系数矩阵的三角矩阵形式,第二个步骤就是逐步回代的过程,最终矩阵的对角线上的元素就是方程组的解。下面就给出高斯消元法的一个算法实现:

     76 /*带列主元的高斯消去法解方程组,最后的解在matrixA的对角线上*/

     77 bool GuassEquation::Resolve(std::vector<double>& xValue)

     78 {

     79     assert(xValue.size() == m_DIM);

     80 

     81     /*消元,得到上三角阵*/

     82     for(int i = 0; i < m_DIM - 1; i++)

     83     {

     84         /*按列选主元*/

     85         int pivotRow = SelectPivotalElement(i);

     86         if(pivotRow != i)/*如果有必要,交换行*/

     87         {

     88             SwapRow(i, pivotRow);

     89         }

     90         if(IsPrecisionZero(m_matrixA[* m_DIM + i])) /*主元是0? 不存在唯一解*/

     91         {

     92             return false;

     93         }

     94         /*对系数归一化处理,使行第一个系数是1.0*/

     95         SimplePivotalRow(i, i);

     96         /*逐行进行消元*/

     97         for(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_DIM - 1) * m_DIM + m_DIM - 1];

    104     for(int i = m_DIM - 2; i >= 0; i--)

    105     {

    106         double totalCof = 0.0;

    107         for(int j = i + 1; j < m_DIM; j++)

    108         {

    109             totalCof += m_matrixA[* m_DIM + j] * m_matrixA[* m_DIM + j];

    110         }

    111         m_matrixA[* m_DIM + i] = (m_bVal[i] - totalCof) / m_matrixA[* m_DIM+ i];

    112     }

    113 

    114     /*将对角线元素的解逐个存入解向量*/

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

    116     {

    117         xValue[i] = m_matrixA[* m_DIM + i];

    118     }

    119 

    120     return true;

    121 }

     

            GuassEquation::Resolve()函数中m_matrixA是以一维数组形式存放的系数矩阵,m_DIM是矩阵的维数,SelectPivotalElement()函数从系数矩阵的第i列中选择绝对值最大的那个值所在的行,并返回行号,SwapRow()函数负责交换系数矩阵两个行的所有值,SimplePivotalRow()函数是归一化处理函数,通过除法操作将指定的行的对角线元素变换为1.0,以便简化随后的消元操作。

    12.2.3 最小二乘法解决“速度与加速度”实验

            根据12.2.1节对最小二乘法原理的分析,用程序实现最小二乘法曲线拟合的算法主要由两个步骤组成,第一个步骤就是根据给出的测量值生成关于拟合多项式系数的方程组,第二个步骤就是解这个方程组,求出拟合多项式的各个系数。根据对上文最终整理的正规方程组的分析,可以看出其系数有一定的关系,就是每一个方程式都比前一个方程式多乘了一个xi。因此,只需要完整计算出第一个方程式的系数,其他方程式的系数只是将前一个方程式的系数依次左移一位,然后单独计算出最后一个系数就可以了,此方法可以减少很多无谓的计算。求解多元一次方程组的方法就使用12.2.2节介绍的高斯消元法,其算法上一节已经给出。

            这里给出一个最小二乘算法的完整实现,以12.1.2节的数据为例,因为数据结果明显呈现线性方程的特征,因此选择拟合多项式为v = v0 + at,v0和a就是要求解的拟合多项式系数。

     99 bool LeastSquare(const std::vector<double>& x_value, const std::vector<double>&y_value,

    100                  int M, std::vector<double>& a_value)

    101 {

    102     assert(x_value.size() == y_value.size());

    103     assert(a_value.size() == M);

    104 

    105     double *matrix = new double[* M];

    106     double *b= new double[M];

    107 

    108     std::vector<double> x_m(x_value.size(), 1.0);

    109     std::vector<double> y_i(y_value.size(), 0.0);

    110     for(int i = 0; i < M; i++)

    111     {

    112         matrix[ARR_INDEX(0, i, M)] = std::accumulate(x_m.begin(), x_m.end(),0.0);

    113         for(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);

    118         for(int k = 0; k < static_cast<int>(x_m.size()); k++)

    119         {

    120             x_m[k] *= x_value[k];

    121         }

    122     }

    123     for(int row = 1; row < M; row++)

    124     {

    125         for(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);

    130         for(int k = 0; k < static_cast<int>(x_m.size()); k++)

    131         {

    132             x_m[k] *= x_value[k];

    133         }

    134     }

    135 

    136     GuassEquation equation(M, matrix, b);

    137     delete[] matrix;

    138     delete[] b;

    139 

    140     return equation.Resolve(a_value);

    141 }

            将表121的数据带入算法,计算得到v= 4.05545455,a = 1.48818182,比作图法得到的结果更精确。以上算法是根据最小二乘法的理论推导系数方程,并求解系数方程得到拟合多项式的系数的一种实现方法,除此之外,还可以利用预先计算好的最小二乘解析理论直接求得拟合多项式的系数,读者可自行学习相关的实现算法。

     

    12.3 三次样条曲线拟合

            曲线拟合基本上就是一个插值计算的过程,除了最小二乘法,其他插值方法也可以被用于曲线拟合。常用的曲线拟合方法还有基于RBFRadial Basis Function)的曲线拟合和三次样条曲线拟合。最小二乘法方法简单,便于实现,但是如果拟合模式选择不当,会产生较大的偏差,特别是对于复杂曲线的拟合,如果选错了模式,拟合的效果就很差。基于RBFRadial Basis Function)的曲线拟合方法需要高深的数学基础,涉及多维空间理论,将低维的模式输入数据转换到高维空间中,使得低维空间内的线性不可分问题在高维空间内变得线性可分,这种数学分析方法非常强大,但是这种方法不宜得到拟合函数,因此在需要求解拟合函数的情况下使用起来不是很方便。

            样条插值是一种工业设计中常用的、得到平滑曲线的一种插值方法,三次样条又是其中用的较为广泛的一种。使用三次样条曲线进行曲线拟合可以得到非常高精度的拟合结果,并且很容易得到拟合函数,本节的内容将重点介绍三次样条曲线拟合的原理和算法实现,并通过一个具体的例子将三次样条函数拟合的曲线与原始曲线对比显示,体会一下三次样条曲线拟合的惊人效果。

     

    12.3.1 插值函数

            前文提到过,曲线拟合的实质就是各种插值计算,因此,插值函数的选择决定了曲线拟合的效果。那么插值函数的数学定义是什么呢?若在[a, b]上给出n + 1个点a ≤ x< x< ⋯ < xn ≤ b,f(x) 是[a, b]上的实值函数, 要求一个具有n + 1个参量的函数s(x; a0,...,an) 使它满足

            s(xi; a0,..., an) = f(xi) , i = 0, 1, ⋯, n                          (3.5)

    则称s(x)为f(x) 在[a, b]上的插值函数. 若s(x) 关于参量a0, a1,...,an是线性关系, 即:

            s(x) = a0s0(x) + a1s1(x) + ⋯ + ansn(x)                                (3.6)

    s(x)就是多项式插值函数,如果si(x)是三角函数,则s(x)就是三角插值函数。

            比较常用的多项式插值函数是牛顿插值多项式和拉格朗日插值多项式,但是在多项式的次数比较高的情况下,插值点数n过多会导致多项式插值在收敛性和稳定性上失去保证,因此,当插值点数n较大的情况下,一般不使用多项式插值,而采用样条插值或次数较低的最小二乘法插值。

    12.3.2 样条函数的定义

            在所有能够保证收敛性和稳定性的插值函数中,最常用的,也是最重要的插值函数就是样条插值函数。采用样条函数计算出的插值曲线和曲面在飞机、轮船和汽车等精密机械设计中都得到了广泛的应用。样条插值函数的数学定义是这样的:

    设区间[a, b]上选取n - 1个节点(包括区间端点a和b共n + 1个节点),将其划分为n个子区间 a = x0 < x1 < ⋯ < xn = b, 如果存在函数s(x),使得s(x)满足以下两个条件:

    (1) s(x) 在整个区间[a, b]上具有m - 1阶连续导数;

    (2) s(x)在每个子区间[xi-1, xi], i= 1, 2, ⋯, n 上是m 次代数多项式(最高次数为m次);

    则称s(x)是区间[a, b]上的m 次样条函数。假如区间[a, b]上存在实值函数f(x),使得每个节点处的值f(xi)与s(xi)相等,即

            s(xi) = f(xi), i = 0, 1, …, n                                          (3.7)

    则称s(x)是实值函数f(x)的m 次样条插值函数。

            当m = 1时,样条插值函数就是分段线性插值, 此时虽然s(x)是属于区间[a, b]上的函数, 但它不光滑(连一阶连续导数性质都不具备),不能满足工程设计要求。工程设计通常使用较多的是m = 3时的三次样条插值函数,此时样条函数具有二阶连续导数性质。

            根据三次样条函数的定义,s(x)在每个子区间上的样条函数si(x)都是一个三次多项式,也就是说,三次样条函数s(x)由n个区间上的n个三次多项式组成,每个三次多项式可描述为以下形式:

            si(x) = aix3 + bix2 + cix + di      i = 1, 2, …, n                      (3.8)

    因此,要确定完整的样条函数s(x)需要确定ai、bi、ci和di公4n个系数。根据样条函数的定义,s(x)在区间内的n - 1个节点处都是连续的,并且其一阶导数si‘(x)和二阶导数si“(x)都是连续的,根据连续函数的性质(xi的左右导数相等),我们可以得到3(n - 1)个条件:

            si(xi - 0) = si+1(xi + 0)    i = 1, 2, …, n-1

            si‘(xi - 0) = si+1‘(xi + 0)    i = 1, 2, …, n-1

            si“(xi - 0) = si+1“(xi + 0)    i = 1, 2, …, n-1                          (3.9)

    再加上插值函数在包括区间端点a(就是x0),b(就是xn)在内的n + 1个节点处满足s(xi) = f(xi),又可以得到n + 1个条件,这样就具备了4n – 2个条件。

    12.3.3 边界条件

            为了解决4n个系数组成的方程组,最终确定的s(x),需要再补充两个边界条件使之满足4n个条件。常用的边界条件有以下几种:

            第一类边界条件,即满足s‘(x0) = f’(x0),s‘(xn) = f’(xn)两个条件,其中f(x)是实值函数。

            第二类边界条件,即满足s”(x0) = f”(x0),s”(xn) = f”(xn)两个条件,其中f(x)是实值函数。特别情况下,当f”(x0) = f”(xn) = 0的时候,也就是s”(x0) = s”(xn) = 0的情况下,第二类边界条件又被称为自然边界条件。

            当样条函数的实值函数f(x)是以[a, b]为周期的周期函数时,三次样条函数s(x)在两个端点处满足s‘(x0 - 0) = s‘(xn + 0)和s”(x0 - 0) = s”(xn + 0),这种情况又被成为第三类边界条件。

            工程技术中常用的是第一类边界条件和第二类边界条件,以及第二类边界条件的特殊情况自然边界条件。理想情况下,也就是实值函数已知的情况下,可以通过实值函数直接计算出边界条件的值,否则的话,就只能通过测量和计算得到边界条件的值,有时候甚至只能给出经验估计值,工程技术中通常根据实际情况灵活使用各类边界条件。

    12.3.4 推导三次样条函数

            求三次样条插值函数s(x)的方法很多,其基本原理都是首先求出由待定系数组成的s(x),以及其一阶导数s’(x)和二阶导数s”(x),然后将其带入到12.3.2和12.3.3节列举的4n个条件中,得到关于待定系数的方程组,最后求解方程组得到待定系数,并最终确定插值函数s(x)。

            求三次样条插值函数s(x)常用的方法是“三转角法”和“三弯矩法”。根据三次样条函数的性质,s(x)的一阶导数s’(x)是二次多项式, 二阶导数s”(x)是一次多项式(线性函数), “三转角法”和“三弯矩法”的主要区别是利用这两个特性推导插值函数s(x)、s’(x)和s”(x)的方式不同。“三转角法”利用s(x)的一阶导数s’(x)是二次多项式这个特性,对于子区间[xi, xi+1],利用抛物线插值公式获得一个通过xi和xi+1两个点的二次多项式做为s’(x),然后对s’(x)进行积分和微分(求导)运算,分别得到s(x),和s”(x),最后将它们带入4n个条件中求解系数方程组。“三弯矩法”则是利用s(x)的二阶导数s”(x)是一次多项式(线性函数)这个特性,对于子区间[xi, xi+1],首先假设一个通过xi和xi+1两个点的线性函数做为s”(x),然后对s”(x)进行连续两次积分运算得到s(x),再对s(x)进行求导得运算到s’(x),最后将它们带入4n个条件中求解系数方程组。这两种方法的本质是一样的,只是对s(x)的推导过程不同,接下来就介绍使用“三弯矩法”求解三次样条函数的方法。

            三次样条函数的求解过程就是系数方程组的推导过程,使用“三弯矩法”推导系数方程组,首先要确定插值函数的二阶导数s”(x)。根据三次样条函数的性质,在每个子区间[xi, xi+1]上,其二阶导数s”(x)是个线性方程,现在假设在xi和xi+1两个端点的二阶导数值分别是Mi和Mi+1,也就是s”(xi) = Mi,s”(xi+1) = Mi+1,则经过xi和xi+1的两点式直线方程是:

    从M0到Mn,有n+1个Mi的值需要求解,但是(3.20)只有n-1个等式,此时就需要用到两个边界条件了。

            如果是使用第二类边界条件,则直接可以得到以下两个条件等式:

            s”(x0) = M0 = f”(x0) = y0’                                                (3.21)

            s”(xn) = Mn = f”(xn) = yn’                                                (3.22)

    令d0 = 2y0’,dn = 2yn’,可以得到由第二类边界条件确定的两个方程:

            2M0 = d                                                                   (3.23)

            2Mn = dn                                                                    (3.24)

            如果是使用第一类边界条件,即s‘(x0) = f’(x0),s‘(xn) = f’(xn)两个条件,则需要将这两个条件代入(3.16),通过计算得到两个条件等式。将s‘(x0) = y0’代入(3.16),得到:

    将第二类边界条件得到的(3.23)和(3.24)或第一类边界条件得到的(3.27)和(3.28)与(3.20)中的n – 1个等式组合在一起就得到一个关于Mi的方程组,求解此方程组可以得到Mi的值,代入到(3.15)即可得到三次样条函数方程。以第一类边界条件得到的(3.27)和(3.28)为例,与(3.20)连立得到以下方程组:

    这就是三弯矩方程组,其中Mi,i=0,1,…n就是三次样条函数s(x)的矩。根据(3.27)和(3.28),un = 1,v0 = 1,其余各系数可以通过(3.19)中的系数计算出来。这个方程组的系数矩阵是一个对角线矩阵,并且是一个严格对角占优的对角阵(ui和vi的值均小于主对角线的值,也就是ui和vi的值皆小于2),可以使用追赶法求解。下一节将介绍如何使用追赶法求解方程组,并给出求解的算法实现。

    12.3.5 追赶法求解方程组

            任意矩阵A,都可以通过克洛脱(Crout)分解得到两个三角矩阵:

    其中y1=d1/l1,其余各项的递推计算关系是:

    yi = (di – miyi-1)/li, i = 2,3,… ,n

    对于第二个方程,求解最终结果xi

    其中xn = yn,其余各项的递推求解关系是:

    xi = yi  uixi+1, i = n  1, n  2, , 1

    递推计算yi和xi的过程分别被形象地形容为追的过程和赶的过程,这也是追赶法得名的原因,实际上这种方法在国际上叫做托马斯(Thomas)法。在这里需要强调一下,对三角矩阵的克洛脱分解需要满足几个条件,否则无法进行,这几个条件分别是:

            下面就给出一个追赶法求解方程组的通用算法实现,在使用之前需要判断系数矩阵是否是对三角矩阵,并且满足上述三个条件,相关的判断请读者自行添加:

     76 /*追赶法求对三角矩阵方程组的解*/

     77 bool ThomasEquation::Resolve(std::vector<double>& xValue)

     78 {

     79     assert(xValue.size() == m_DIM);

     80 

     81     std::vector<double> L(m_DIM);

     82     std::vector<double> M(m_DIM);

     83     std::vector<double> U(m_DIM);

     84     std::vector<double> Y(m_DIM);

     85 

     86     /*消元,追的过程*/

     87     L[0] = m_matrixA[ARR_INDEX(0, 0, m_DIM)];

     88     U[0] = m_matrixA[ARR_INDEX(0, 1, m_DIM)] / L[0];

     89     Y[0] = m_bVal[0] / L[0];

     90     for(int i = 1; i < m_DIM; i++)

     91     {

     92         if(IsPrecisionZero(m_matrixA[ARR_INDEX(i, i, m_DIM)]))

     93         {

     94             return false;

     95         }

     96         M[i] = m_matrixA[ARR_INDEX(i, i - 1, m_DIM)];

     97         L[i] = m_matrixA[ARR_INDEX(i, i, m_DIM)] - M[i] * U[- 1];

     98         Y[i] = (m_bVal[i] - M[i] * Y[- 1]) / L[i];

     99         if((+ 1) < m_DIM)

    100         {

    101             U[i] = m_matrixA[ARR_INDEX(i, i + 1, m_DIM)] / L[i];

    102         }

    103     }

    104     /*回代求解,赶的过程*/

    105     xValue[m_DIM - 1] = Y[m_DIM - 1];

    106     for(int i = m_DIM - 2; i >= 0; i--)

    107     {

    108         xValue[i] = Y[i] - U[i] * xValue[+ 1];

    109     }

    110 

    111     return true;

    112 }

     

    12.3.6 三次样条曲线拟合算法实现

            根据12.3.4节对三次样条函数的推导分析,三次样条曲线拟合算法的核心可分为三部分,第一部分是根据推导结果计算关于三次样条函数的“矩”的方程组的系数矩阵,第二部分就是用追赶法求解方程组,得到各个区间的三次样条函数,第三部分就是根据每个拟合点的输入的值xi,确定使用哪个区间的三次样条函数,并使用该区间的三次样条函数计算出三次样条插值yi,最后得到的一系列(xi, yi)组成的曲线就是三次样条拟合曲线。拟合算法也是按照上面的分析,分三个步骤计算插值:

            第一步是计算系数矩阵,其中u0、v0、d0和dn的值需要单独计算,其余的值可以通过(3.19)递推计算出来。

            第二步是将系数矩阵代入12.3.5节给出的追赶法通用算法,求出Mi的值。求解之前,先证明一下第一步得到系数矩阵是否满足追赶法的条件。首先,主对角线元素的值都是2,满足12.3.5节的条件(1)。其次,由ui和vi的计算条件可知,|ui| < 1,|vi| < 1,这也满足12.3.5节的条件(2)。最后,因为ai = 2,且ui和vi的和是1,所以12.3.5节的条件(3)也得到了满足。由上判断可知,求解三次样条函数的“矩”的系数矩阵满足使用追赶法求解的条件,可以使用追赶法求解。

            第三步是计算插值,需要将第二步计算得到的Mi代入(3.15),并选择合适的子区间样条函数计算出插值点的值。

            下面就给出采用三弯矩法实现的三次样条曲线拟合算法,CalcSpline()函数的参数Xi和Yi是n个插值点(包括起点和终点)的值,boundType是边界条件类型,b1和b2分别是对应的两个边界条件,这个算法支持第一类和第二类边界条件(包括自然边界条件)。内部的矩阵matrixA就是按照(3.29)构造的Mi方程组的系数矩阵,可用于直接用追赶法求解方程组。CalcSpline()函数的大部分代码都是在构造Mi方程组的系数矩阵,首先根据边界条件确定un 、v0 、d0和dn,其他系数则根据(3.19)的递推关系,在for(int i = 1; i < (m_valN - 1); i++)循环中依次计算出来,最后是利用12.3.5节给出的追赶法算法求出Mi。GetValue()函数负责计算给定区间内任意位置的插值,首先根据x的值确定使用哪个子区间的样条函数,然后根据(3.12)和(3.14)给出的关系计算插值。

     16 void SplineFitting::CalcSpline(double *Xi, double *Yi, int n, int boundType,double b1, double b2)

     17 {

     18     assert((boundType == 1) || (boundType == 2));

     19 

     20     double *matrixA = new double[* n];

     21     if(matrixA == NULL)

     22     {

     23         return;

     24     }

     25     double *= new double[n];

     26     if(== NULL)

     27     {

     28         delete[] matrixA;

     29         return;

     30     }

     31 

     32     m_valN = n;

     33     m_valXi.assign(Xi, Xi + m_valN);

     34     m_valYi.assign(Yi, Yi + m_valN);

     35     m_valMi.resize(m_valN);

     36     memset(matrixA, 0, sizeof(double) * n * n);

     37 

     38     matrixA[ARR_INDEX(0, 0, m_valN)] = 2.0;

     39     matrixA[ARR_INDEX(m_valN - 1, m_valN - 1, m_valN)] = 2.0;

     40     if(boundType == 1) /*第一类边界条件*/

     41     {

     42         matrixA[ARR_INDEX(0, 1, m_valN)] = 1.0; //v0

     43         matrixA[ARR_INDEX(m_valN - 1, m_valN - 2, m_valN)] = 1.0; //un

     44         double h0 = Xi[1] - Xi[0];

     45         d[0] = 6 * ((Yi[1] - Yi[0]) / h0 - b1) / h0; //d0

     46         double hn_1 = Xi[m_valN - 1] - Xi[m_valN - 2];

     47         d[m_valN - 1] = 6 * (b2 - (Yi[m_valN - 1] - Yi[m_valN - 2]) / hn_1) /hn_1; //dn

     48     }

     49     else /*第二类边界条件*/

     50     {

     51         matrixA[ARR_INDEX(0, 1, m_valN)] = 0.0; //v0

     52         matrixA[ARR_INDEX(m_valN - 1, m_valN - 2, m_valN)] = 0.0; //un

     53         d[0] = 2 * b1; //d0

     54         d[m_valN - 1] = 2 * b2; //dn

     55     }

     56     /*计算ui,vi,di,i = 2,3,...,n-1*/

     57     for(int i = 1; i < (m_valN - 1); i++)

     58     {

     59         double hi_1 = Xi[i] - Xi[- 1];

     60         double hi = Xi[+ 1] - Xi[i];

     61         matrixA[ARR_INDEX(i, i - 1, m_valN)] = hi_1 / (hi_1 + hi); //ui

     62         matrixA[ARR_INDEX(i, i, m_valN)] = 2.0;

     63         matrixA[ARR_INDEX(i, i + 1, m_valN)] = 1 - matrixA[ARR_INDEX(i, i - 1,m_valN)]; //vi = 1 - ui

     64         d[i] = 6 * ((Yi[+ 1] - Yi[i]) / hi - (Yi[i] - Yi[- 1]) / hi_1) /(hi_1 + hi); //di

     65     }

     66 

     67     ThomasEquation equation(m_valN, matrixA, d);

     68     equation.Resolve(m_valMi);

     69     m_bCalcCompleted = true;

     70 

     71     delete[] matrixA;

     72     delete[] d;

     73 }

     74 

     75 double SplineFitting::GetValue(double x)

     76 {

     77     if(!m_bCalcCompleted)

     78     {

     79         return 0.0;

     80     }

     81     if((< m_valXi[0]) || (> m_valXi[m_valN - 1]))

     82     {

     83         return 0.0;

     84     }

     85     int i = 0;

     86     for(= 0; i < (m_valN - 1); i++)

     87     {

     88         if((>= m_valXi[i]) && (< m_valXi[+ 1]))

     89             break;

     90     }

     91     double hi = m_valXi[+ 1] - m_valXi[i];

     92     double xi_1 = m_valXi[+ 1] - x;

     93     double xi = x - m_valXi[i];

     94 

     95     double y = xi_1 * xi_1 * xi_1 * m_valMi[i] / (6 * hi);

     96     y += (xi * xi * xi * m_valMi[+ 1] / (6 * hi));

     97 

     98     double Ai = (m_valYi[+ 1] - m_valYi[i]) / hi - (m_valMi[+ 1] -m_valMi[i]) * hi / 6.0;

     99     y += Ai * x;

    100     double Bi = m_valYi[+ 1] - m_valMi[+ 1] * hi * hi / 6.0 - Ai * m_valXi[i+ 1];

    101     y += Bi;

    102     return y;

    103 }

    12.3.7 三次样条曲线拟合的效果

    本节将用定义一个原始函数,从原始函数的某个区间上抽取9个插值点,根据这9个插值点和原函数的边界条件,利用三次样条曲线插值进行曲线拟合,并将原始曲线和拟合曲线做对比,展示一下三次样条曲线拟合的效果。

    首先定义原始函数:

    f(x) = 3/(1 + x2)

    选择区间[0.0, 8.0]上的9个点做为插值点,计算各点的值为:

    x

    0.0

    1.0

    2.0

    3.0

    4.0

    5.0

    6.0

    7.0

    8.0

    y

    3.0

    1.5

    0.6

    0.3

    0.1765

    0.1154

    0.0811

    0.06

    0.0462

     12 – 2 原函数f(x)在各插值点的值

    求f(x)的导函数f’(x):

    f’(x) = -6x/(1 + x2)2

    根据f’(x)计算出在区间端点处的两个第一类边界条件f’(0.0) = 0.0, f’(8.0) = -0.01136。利用表12-2中的数据和这两个边界条件,计算出三次样条插值函数,并从0.0开始,以0.01为步长,连续求800个点的插值,将这些点连成曲线得到拟合曲线。为了作对比,同样从0.0开始,以0.01为步长,用f(x)函数连续计算800个点的原值,将这些点连成曲线得到原始曲线。分别用不同的颜色画出这两条曲线,如图12-2所示:

    图12-2 拟合曲线和原始曲线对比

            从图12-2可以看到,三次样条曲线拟合的效果非常好。同样在[0.0, 8.0]区间上,如果增加插值点的个数,将获得更好的拟合效果。比如以0.5为单位,将插值点增加到17个,则拟合的曲线与原始曲线几乎完全重合。


    网址:http://blog.csdn.net/kuikuijia/article/details/45405911

    展开全文
  • 网上搜集的最小二乘法曲线拟合演示程序,和对任意若干点进行曲线拟合,可选择拟合多项式次数 网上搜集的最小二乘法曲线拟合演示程序,和对任意若干点进行曲线拟合,可选择拟合多项式次数
  • 曲线拟合问题的提法是,已知一组(二维)数据,即平面上的n个点 , 互不相同,寻求一个函数(曲线),使 在某种准则下与所有数据点最为接近,即曲线拟合得最好。 最小二乘准则 系数 的确定 1.2 函数的选取 ...
  • 第一种是进行多项式拟合,数学上可以证明,任意函数都可以表示为多项式形式。具体示例如下。 ###拟合年龄 import numpy as np import matplotlib.pyplot as plt   #定义x、y散点坐标 x = [10,20,30,40,50,60,70,80...
  • 曲线拟合

    2019-10-27 14:15:45
    股票预测问题 plot函数讲解: ...第三个参数要加单引号 比如上面这个问题 >> x=[2,3,4,5,8,9,10,11,12,15,16,17,18,19,22,23,24,25,26,29,30]; >> y=[7.74,7.84,7.82,7.7...
  • 二维曲线拟合

    千次阅读 2018-11-20 16:49:03
    二维曲线拟合相关基础理论知识最小二乘法广义逆(伪逆矩阵)矩阵分解特征值分解(Eigen Value Decomposition,EVD)Householder变换法(QR Factorization,QR分解)奇异值分解(singular value decomposition,SVD)...
  • 最小二乘法的曲线拟合

    万次阅读 2018-08-26 14:08:50
    最小二乘法解决的问题:Ax=C 无解下的最优解 例子1: 一条过原点的直线OA,C是直线外一点,求C在OA上的投影点P 例子1   ...已知三个不在一条直线上的点A,B,C,求一条直线,使A,B,C到直线的距离和最小 ...
  • 最小二乘曲线拟合——C语言算法实现一

    万次阅读 多人点赞 2014-05-29 20:45:01
    给定一组数据,我们要对这组数据进行曲线拟合
  • 曲线拟合最小二乘法c语言

    热门讨论 2012-12-16 14:30:34
    曲线拟合最小二乘法c语言 计算方法 #include #include #define N 9 #define M 3 void main() { int i,j; float a[2][N],b[5][N],c[7]={0,0,0,0,0,0,0};; printf("请输入%d个点的X坐标\n",N); for(j=0;j;j++) ...
  • 最小二乘法曲线拟合(C++)

    热门讨论 2010-10-23 16:49:15
    用C++编写的程序,用最小二乘法实现对曲线拟合拟合的多项式达到六阶。
  • 最小二乘法入门(Matlab直线和曲线拟合

    万次阅读 多人点赞 2019-02-28 14:13:47
    参考博客:... 多的就不多说了,持续脱发中!!! 最小二乘法历史起源之类的:https://baike.baidu.com/item/%E6%9C%80%E5%B0%8F%E4%BA%8C%E4%B9%98%E6%B3%95/2522346?fr=aladdin ...
  • 最小二乘法多项式曲线拟合原理与实现

    万次阅读 多人点赞 2012-04-27 16:38:24
    最小二乘法多项式曲线拟合,根据给定的m个点,并不要求这条曲线精确地经过这些点,而是曲线y=f(x)的近似曲线y= φ(x)。 原理 [原理部分由个人根据互联网上的资料进行总结,希望对大家能有用]  给定数据点pi...
  • 多项式拟合最小二乘法曲线拟合C语言实现,有详细描述文档和代码
  • Python实现最小二乘法曲线拟合

    千次阅读 2018-12-18 18:48:18
    感谢 最小二乘法曲线拟合以及Matlab实现  作者 乐乐lelele 建议先看一下该作者的讲解,比较详细 最小二乘法拟合最终通过求解N次方程的系数来实现的 Y=X*A A=INV(X)*Y #coding=utf-8 #独立测试脚本 #最小...
  • 最小二乘法曲线拟合C语言可执行代码
  • 最小二乘法进行曲线拟合@TOC 先假设拟合曲线的阶次,再利用最小二乘方法估计待估计参数。 clc; clear all; z=[66.6,84.9,88.6,78.0,96.8,105.2,93.2,111.6,88.3,117.0,115.2]; plot(z,‘o–’); xlabel(‘第k年’); ...
  • 最小二乘法多项式曲线拟合及其python实现多项式曲线拟合问题描述最小二乘法针对overfitting,加入正则项python实现运行结果 多项式曲线拟合问题描述 问题描述:给定一些数据点,用一个多项式尽可能好的拟合出这些点...
  • 最小二乘法(又称最小平方法)是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。...最小二乘法还可用于曲线拟合。其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。
  • 对给定数据点{(Xi,Yi)}(i=0,1,…,m),在取定的函数类Φ 中,求p(x)∈Φ ,使误差的平方和E^2最小,E^2=∑[p(Xi)-Yi]^2。...函数p(x)称为拟合函数或最小二乘解,求拟合函数p(x)的方法称为曲线拟合最小二乘法
  • 一、实验内容已知一组实验数据如下表,求它的拟合曲线。 x(i) 1 2 3 4 5 f(i) 4 4.5 6 8 8.5 w(i) 2 1 3 1 1 二、程序清单与...% mypolyfit输出通过最小二乘法求得的拟合曲线并绘图验证 % 例如: % X=[1 2 3 4 5];...
  • 最小二乘法曲线拟合+C代码

    千次阅读 2015-07-30 21:54:54
    最小二乘法曲线的系数求解过程是解一个正规方程组的解的过程,下图是数值分析课本上,最小二乘法拟合的原理: 课本中的例子如下: c代码如下: #include "stdio.h" #include "stdlib.h" #...
  • 给出了最小二乘法在多元正交基函数拟合中的计算机实现方法,以常见的二次曲线拟合为例说明了程序编制的要点,在实验的数据处理中具有实用价值。
  • 曲线拟合最小二乘法在数据采集中的应用曲线拟合最小二乘法在数据采集中的应用
  • VB实现最小二乘法多次曲线拟合 VB实现最小二乘法多次曲线拟合 VB实现最小二乘法多次曲线拟合 VB实现最小二乘法多次曲线拟合
  • 目录 最小二乘法的原理 实例: 求解结果: c++程序源代码: 最小二乘法的原理 拟合函数: 式中:s(x)为拟合函数, 为拟合系数,为函数族 平方误差: ...
  • 最小二乘法进行最高3次曲线拟合

    千次阅读 2016-12-20 13:46:34
    最近在做跟踪时,需要预测被跟踪...下面为最小二乘法的核心代码,有需要可以参考: bool CNXMinSquare::Calc(std::vector &vtCoef, std::vector &vtPoint) { vtCoef.clear(); // 1. 根据函数值对 1 x x^2 x^3 施
  • 最小二乘法曲线拟合 C语言实现

    万次阅读 2014-04-26 11:55:00
    简单思路如下: 1,采用目标函数对多项式系数求偏导,得到最优值条件,组成一个方程组; 2,方程组的解法采用行列式变换(两次变换:普通行列式——三角行列式——对角行列式——求解),行列式的求解算法上优化过...

空空如也

1 2 3 4 5 ... 20
收藏数 69,554
精华内容 27,821
关键字:

曲线拟合