范文3XXX毕业论文.doc

本科生毕业论文最小二乘法原理的应用及其Matlab实现姓名XXX学号20084070650专业信息与计算科学指导教师XXX 副教授2012年5月20日II摘 要最小二乘法（又称最小平方法）是一种数学优化技术.它通过最小化误差的平方来寻找数据的最佳函数匹配，在参数估计以及预测预报等诸多领域得到了极为广泛的应用.利用最小二乘法可以简便的求得未知的数据，并使得这些求得的数据和实际数据之间误差的平方和最小.本文详细给出了最小二乘法的理论推导过程，包括一般最小二乘法原理和正交最小二乘法原理，以及用最小二乘进行回归分析的参数估计.并以实例的形式说明了如何使用Matlab求解不同的拟合问题并给出相应的Matlab程序.主要介绍了最小二乘法在以下三方面的应用，处理物理实验数据上的应用，在给药方案上的应用以及在估计水塔流量上的应用.用Matlab实现最小二乘方法的数据处理，避免了传统数据处理方法的弊端，增加了数据处理的准确性及快捷性，Matlab的可视化功能也可以更加直观地反映数据结果.关键词：最小二乘法；拟合；最小二乘法的Matlab实现AbstractMethod of the least square is a mathematical optimization technology.It looks for the best data matching function by the way of minimizing error square.It is applied in the parameter estimation and prediction,and many other areas.Using the least square method,we can detemine unknown data easily and make these data and actual data obtained the least squares of the error.This paper presents theoretical deducing process of the least square method and estimates the parameter of the regression analysis with the least squares.In addition,this paper illustrates how to use Matlab to solve different fitting problems through some actual examples and gives the corresponding Matlab program.Matlab is efficient in data processing.The visual function of Matlab can also reflect the data results more directly.Key Words: least square method, fitting, the matlab programming of the least square目 录1 引言12 曲线拟合的最小二乘法12.1 一般最小二乘问题12.2 用正交函数作最小二乘拟合42.3 模型线性化52.3.1 指数模型52.3.2 幂律模型 72.4 用Gauss-Newton方法解非线性最小二乘72.4.1 Gauss-Newton方法简介72.4.2 带非线性系数的模型92.5 求解最小二乘数据拟合问题步骤93 用Matlab解拟合问题103.1 用Matlab作线性最小二乘拟合103.2 用Matlab作非线性最小二乘拟合113.2.1 lsqcurvefit用法简介113.2.2 lsqnonlin用法简介113.3 用Matlab进行回归分析123.3.1 回归分析的相关理论123.3.2 用Matlab进行回归分析相关函数调用133.4 具体实例应用153.4.1 测力计示值与测量时对应温度的关系拟合问题实例1153.3.2 给药方案拟合问题实例2163.3.3 估计水塔的流量拟合问题实例319参考文献22致谢23附录24 301 引 言对工程实践中测得的数据进行理论分析，用恰当的函数进行模拟数据原型是一类十分重要的问题，我们常用的数据处理方法有列表法、作图法、逐差法等，这些方法虽有许多优点，但带有一定的主观任意性，容易造成人为误差.最常用的逼近原则是让实测数据和估计数据之间的平方和最小，这即是最小二乘法.最小二乘概念是从19世纪早期Gauss和Legendre的开拓性工作开始的.这种回归和参数估计得关键技术在科学和工程中已经成为基本工具.最小二乘法是在诸多学科领域中获得广泛应用的数据处理方法，用数理统计的方法来处理实验数据，具有较高的可信度，但由于最小二乘法计算较复杂，在数据处理中常常不被采用.Matlab软件是一种集中了数值计算、符号运算及强大的图形处理、程序设计语言等功能于一体的科学计算语言.用Matlab处理实验数据非常简便，克服了最小二乘法计算量大的缺点，将繁琐、枯燥的数值计算过程简单化、直观化、可视化，并能较准确的标记处出实验数据点和绘制出拟合曲线.本文在详细给出最小二乘法的理论推导过程38的基础上针对不同类型的拟合问题，包括线性问题和一般的非线性问题，用具体实例的形式给出了如何使用Matlab去求解及程序实现15.2曲线拟合的最小二乘法2.1 一般最小二乘问题设给定如下数据：其中,称为权系数，它可以表示此点的重要程度，也可以表示此点的重复次数.今要求建立,之间的函数关系，可以用插值法1实现.众所周知，反映自然规律的数量关系的函数通常有三种表示方法，解析法、图像法和表格法.大量实际问题中的函数关系是用表格的形式给出的.从提供的部分离散的函数值去进行理论的分析和设计都是极不方便甚至是不可能的，因此需要设法寻找与已知函数相符而形式简单的插值函数.设函数在区间上有定义，且已知在点上的函数值.若存在一个简单的函数使成立，就称为的插值函数，点成为插值节点，成为插值区间，求的方法称为插值法.但由于这些数据往往是由实验得到的，会带有一定的误差4，插值法要求在插值节点的函数值相等，这就将误差带入函数关系中.另外，这样的数据往往较多，会使所要求的插值多项式的次数较高，给实际应用带来不便，甚至影响逼近效果.因此，要求所求的函数关系不经过点，只要求的给定点上的误差按某种度量标准最小.若记，常用(称为2范数)1来度量误差的大小，故通常称为最小二乘逼近.问题的一般提法是：对给定的一组数据及权系数在函数类中找一个函数，使误差平方最小化，即 , （2.1） 问题等价于求.使 , （2.2） 若令 , （2.3） 即相当于求多元函数的极小值问题.由多元函数求极小值的必要条件知： , （2.4） 其中，改写可得 ,即, , （2.5） 令,则（2.5）式变为 , （2.6） 由于线性无关，所以式（2.6）的系数行列式.上述行列式不等于0，从而保证了方程组有唯一解.得到的即为所求.此处由如下定理保证：在上线性无关的充要条件为.应该注意的是，若，即用多项式作最小二乘法时，次数不宜太大，否则方程组往往是病态的.2.2 用正交函数作最小二乘拟合上面已经提到当用高次（一般大于等于7时）多项式作最小二乘拟合时，往往会使得法方程是病态方程组，这会使得方程组的解有较大的误差.但如果是关于点，带权正交时，即,式（2.6）的解为.则所求解为 . （2.7）这就解决了病态方程组的问题.现我们根据给定节点及权函数，构造出带权正交多项式.这里，用递推公式表示，即 （2.8） 这里是首项系数为1的次多项式.由的正交性，得系数，和如下计算.利用归纳法可以证明这样构造的是正交的.由（2.7）式即可得到用正交多项式所做的最小二乘曲线拟合.2.3 模型线性化有些模型显然不是线性的，不能直接通过最小二乘进行拟合，所以就需要对模型进行线性化.下面主要介绍两种非线性模型的线性化.2.3.1指数模型指数模型即中不是线性的出现在模型方程中，所以不能直接通过最小二乘进行拟合.一旦把数据点带入模型，困难就明显：要求解系数的方程组不是线性的而且不能表示为线性方程组.处理非线性系数的问题有两种方法，比较难的方法是直接极小化最小二乘误差，即解非线性最小二乘问题.Gauss-Newton方法2的一种重要应用是对系数为非线性的模型进行拟合.设是数据点，是要拟合的函数，这里是一组参数，选择它们来极小化残差的平方和.如果参数以线性方式进入模型，那么这是一组关于的线性方程，因此正规方程2的解给出了参数的最优选择.如果参数在模型中是非线性的，那么进行同样处理会得到关于是非线性的方程组的结果.例如，对数据点拟合模型，产生非线性方程组，即,因此非线性的进入模型.比较简单的方法改变这个问题.我们可以解一个不同的问题，它通过“线性化”模型而与原问题相关，而不是解原来的最小二乘问题.在指数模型这种情形下，通过应用以下的自然算法把模型线性化：, 注意，对于指数模型,的图形是线性的.初看起来似乎我们仅仅把一个问题换成另外一个.现在系数在模型中是线性的，但是还不是.然而通过重新命名，我们就可以写成 , (2.9) 现在两个系数及在模型中都是线性的.对最佳的及求解正规方程之后，可得到相应的.我们提出的原来的最小二乘问题是对拟合数据，也就是求使得方程成立的残差的平方和 极小化的,.而现在，我们求解的问题是在“对数空间”极小化最小二乘误差，也就是通过求使得方程成立的残差平方和 (2.10) 极小化的,.这是两种不同的极小化，并有不同的解，这意味着它们一般情况下会得到系数,的不同值.前者是我们已经定义的最小二乘后者并不是.然而，依据数据的前后关系，两者中任意一个都是更自然的选取.使用者需要确定哪一种误差对极小化更重要，是原来意义下的误差还是在“对数空间”下的误差.事实上，对数模型是线性的.因此通过对数变换把数据变成线性关系之后评估模型的拟合才会比较自然.故一般情况下采用此种做法.2.3.2 幂律模型 幂律模型即也可以通过在两边取对数进行线性化而得到化简 , （2.11） 把数据带入模型将给出 , （2.12）结果得到矩阵形式： , （2.13） 正规方程得以确定以及,以及.2.4 用Gauss-Newton方法解非线性最小二乘2.4.1 Gauss-Newton方法我们可以用Gauss-Newton方法对数据进行具有非线性系数的模型拟合.下面介绍该方法的一些基本理论.考虑含个未知量和个方程的方程组用下面的函数（2.14）表示误差的平方和，这里.在定义中包括常数以简化后面的公式.为了极小化，令梯度为零：.（2.15）注意此处我们已经用了梯度的内积法则.设是个变量的标量值函数.的梯度是向量值函数，这里下标表示关于那个变量的偏导数.设是个变量的向量值函数.的Jacobi矩阵是.下述为矩阵代数中两种典型积的乘积法则.设和是向量值函数，而且令是矩阵函数.点积是标量值函数.下面第一个公式表示如何求它的梯度.矩阵和向量的乘积是一个向量，它的Jacobi矩阵表示由第二条法则给出.向量点积法则：.矩阵/向量乘积法则：，这里表示的第列.为计算公式（2.15），需要向量值函数的Jacobi矩阵.把行向量写成列向量.可以用矩阵/向量乘积法则得到，这里是的第列.注意（的二阶偏导数矩阵），.Gauss-Newton方法求解误差平方的梯度的根，尽管在极小值处梯度必须是零，但反之不一定正确.所以这种方法可能收敛于最大值或中性点.说明这种算法结果时必须小心.2.4.2 带非线性系数的模型Gauss-Newton方法的一种重要应用是对系数为非线性的模型进行拟合.设是数据点，是要拟合的函数，这里是一组参数，选择它们来极小化残差的平方和.如果参数是非线性的进入模型，为了解原来的最小二乘问题，用Gauss-Newton方法.它用于对作为参数的向量函数的误差函数(2.14)进行极小化.矩阵是误差关于参数的偏导数的矩阵，它们是.有了这些信息，就可以实现Gauss-Newton方法了.2.5 求解最小二乘数据拟合问题步骤给定一组数据点.首先选取模型.即拟合函数中函数的选取.对表征的关系，式中为参数.因自变量已知时已测得，故类型决定后，由误差方程，就可用最小二乘法决定经验公式中参数.下面就介绍经验公式的选取.对于经验公式类型的选取主要靠专业知识来决定，从数学方面去决定两个量的经验公式主要有三种方法：观察法、近似法和严格计算法.这三种方法中，观察法简单、直观，但相对来说近似法和严格计算法比较精确，但比较繁琐.对于大多数情况，都可以采用观察法来确定经验公式.观察法就是将数据作图，与下图进行比较，看所作图形与哪个类似，就取该类型为经验公式类型. 图2-2 多项式拟合1 图2-3 多项式拟合2图2-1 直线拟合图2-6 幂律模型2图2-5 幂律模型1图2-4 反比例函数模型 选取模型后就用模型拟合数据.具体做法是把数据点带入模型.每个数据点产生一个将未知量作为参数的方程，譬如在上述图1中直线模型中的与.结果产生方程组，在这里未知量表示未知知参数.最后参数的最小二乘解将作为正规方程组2的解而求解.如果是非线性模型，则用上述介绍的Gauss-Newton方法解决.3用Matlab解拟合问题3.1 用Matlab作线性最小二乘拟合若作多项式拟合，可利用已有程序，此处为输出拟合多项式的系数即，为输入同长度的数组，为拟合多项式的次数.拟合出的多项式在处的值可用以下命令计算：.3.2 用Matlab作非线性最小二乘拟合Matlab57提供了两个求非线性最小二乘拟合的函数：lsqcurvefit和lsqnonlin.两个命令都要先建立M文件，在其中定义函数，但两者定义的方式不同.详见如下.3.2.1 lsqcurvefit用法简介，为已知数据点， lsqcurvefit用以求含参量(向量)的向量值函数中的参变量（向量），使得最小.输入格式为,是事先建立的定义函数的文件，自变量为,为迭代初值,和为已知数据点.3.2.2 lsqnonlin用法简介,为已知数据点，lsqnonlin用以求含参量（向量）的向量值函数中的参量，使得最小，其中.调用格式为，其中是事先建立的定义函数的文件，自变量为,为迭代初值.3.3 用Matlab进行回归分析3.3.1 回归分析的相关理论一般地，称由确定的模型为一元线性回归模型12，记为,固定的未知参数、称为回归系数，自变量x也称为回归变量.，称为对的回归直线方程.一元线性回归分析的主要任务是用试验值（样本值）对、和作点估计；对回归系数、作假设检验；在处对作预测，对作区间估计.模型的参数估计可采用最小二乘估计回归系数，有组独立观测值，设 记最小二乘法就是选择和的估计,使.其中,经验回归方程为.拟合出回归方程后，需要对回归方程进行显著性检验，对回归方程的显著性检验，归结为对假设进行检验.假设被拒绝，则回归显著，认为与存在线性关系，所求的线性回归方程有意义，否则回归不显著，与的关系不能用一元线性回归模型来描述，所得的回归方程也无意义.常见的检验法12有以下几种：(1) 检验法 当成立时，,其中 （回归平方和）,故，拒绝，否则就接受.(2) t检验法当成立时，t（n-2），故，拒绝，否则就接受.其中.(3) 检验法记，当|r| r1-时，拒绝H0；否则就接受H0.其中.回归系数和的置信水平为的置信区间分别为：.3.3.2 用Matlab进行回归分析相关函数调用(1) 确定回归系数的点估计值时用命令,其中，对一元线性回归，取即可.(2) 求回归系数的点估计和区间估计并检验回归模型：，其中bint为回归系数的区间估计，为残差，rint为回归系数的置信区间，用于检验回归模型的统计量中有三个数值即相关系数,值,与值对应的概率.相关系数越接近于1，说明该回归方程越显著，时拒绝，越大，说明回归方程越显著，与对应的概率时拒绝，回归模型成立.画出残差及其置信区间调用.(3) 一元多项式回归先进行回归，则确定多项式系数的命令为,其中,是多项式的系数，S是一个矩阵，用来估计预测误差.一元多项式回归命令为.而进行预测和预测误差估计的命令为求所得的回归多项式在处的预测值，如下命令.为求所得的回归多项式在处的预测值及预测值的显著性为的置信区间；缺省时为0.5.(4) 非线性回归先进行回归分析，确定回归系数的命令为nlinfit，其中为估计出的回归系数，为残差，为Jacobian矩阵，输入数据,分别为矩阵和维列向量，对一元非线性回归，为维列向量.model为是事先用M文件定义的非线性函数, beta0为回归系数的初值.预测和预测误差估计的命令为，即求nlinfit所得的回归函数在处的预测值及预测值的显著性为的置信区间,缺省时为0.5.3.4 具体实例应用3.4.1 测力计示值与测量时对应温度的关系9拟合问题实例1用实验的方法研究测力计示值与测量时对应温度的关系.现测得值对应如下： 表3-1 测力计示值与测量时对应温度的关系t/15 18 21 24 27 30F/N43.6143.6343.68 43.7143.7443.78由表3-1数据进行一次拟合，由Matlab编程建立测力计示值和测量时对应温度的一次模型和残差分析图如下图3-1所示.具体Matlab编程见附录一.回归分析得到的一次曲线为.用Matlab编程运行得到的检验回归模型的统计量结果为stats=0.9920 ,499.1250 ,0.0000.（注 stats(1)为相关系数,stats(2)为检验的值,stats(3)为与对应的概率.相关系数越接近1，说明回归方程越显著；显然0.00000)；假设3：血液容积v，t=0注射剂量为d,血药浓度立即为d/v.，在此设d=300,t及c(t)在某点处的值见表3-2，下面拟合求出参数,.图3-2 血药浓度变化规律 方法一：用线性最小二乘拟合，先做如下替换进行转换先作线性拟合求出中的和，再间接求出参数,，即.具体的Matlab编程见附录二. 运行后如下图3-3： 图3-3 线性最小二乘拟合图示计算结果得= -0.2347,=2.9943 ,则可得.用Matlab编程得到的检验回归模型的统计stats=1.0e+003*0.0010 ,1.6113, 0.0000，显然0.00000.05(显著性水平)，所以可以认为该模型成立.第八个数据点可视为异常点，从残差图可以看出所有的数据的残差里零点均较近，且残差的置信区间均包含零点，这说明回归模型能较好的符合原始数据.方法二：直接用非线性最小二乘拟合，调用lsqcurvefit，具体的Matlab编程见附录三.运行结果如下图3-4： 图3-4 非线性最小二乘拟合图示 图3-4 非线性最小二乘拟合图示用Matlab编程得到的检验回归模型的统计量stats = 0.9393 ,108.3216 ,0.0000，显然0.0000 拟合水位-时间函数从测量记录看，一天有两个供水时段（以下称第1供水时段和第2供水时段），和3个水泵不工作时段（以下称第1时段t=0到t=8.97，第2时段t=10.95到t=20.84和第3时段t=23以后）.对第1,2时段的测量数据直接分别作多项式拟合，得到水位函数.所谓多项式拟合，假如数据满足的函数原型为次多项式，则函数利用上述最小二乘原理可以求出，其格式命令为,其中输入参数和分别是已知数据的横坐标和纵坐标数组，指明采用的拟合多项式的次数，输出参数是一个的行向量，具体来说，是的系数，是的系数，,是常数项.为使拟合曲线比较光滑，多项式次数不要太高，次数一般在3-6.由于第3时段只有3个测量记录，无法对这一时段的水位作出较好的拟合.2 确定流量-时间函数 对于第1,2时段只需将水位函数求导数即可，对于两个供水时段的流量，则用供水时段前后（水泵不工作时段）的流量拟合得到，并且将拟合得到的第2供水时段流量外推，将第3时段流量包含在第2供水时段内.(1)第一时段的水位时间函数，并导出流量,具体的编程见附录四.得到的流量函数为.如下图3-5所示 图3-5 第一时段水位时间函数图(2)拟合第二时段的水位时间函数并导出流量,具体的Matlab编程见附录五.得到的流量函数为.如下图3-6所示 图3-6 第二时段水位时间函数图(3)拟合供水时段的流量在第1供水时段（）之前（即第1时段）和之后（即第2时段）各取几点，其流量已经得到，用它们拟合第1供水时段的流量.为使流量函数在和连续，我们简单地只取4个点，拟合3次多项式，具体Matlab实现编程见附录六.拟合的流量函数为 . 图3-7 第一供水时段流量时间函数图在第2供水时段之前取20,20.8两点的流水量，在该时刻之后(第3时段)仅有3个水位记录，我们用差分得到流量，然后用这4个数值拟合第2供水时段的流量，具体编程见附录七.拟合的流量函数为：.如下图3-8所示： 图3-8 第二供水时段流量时间函数图3 一天总用水量的估计总用水量等于两个水泵不工作时段和两个供水时段用水量之和，它们都可以由流量对时间的积分得到.第1、2时段和第1、2供水时段流量的积分之和，就是一天总用水量.虽然诸时段的流量已表为多项式函数，积分可以解析地算出，这里用数值积分计算如下：是第1时段用水量，0.1为积分步长，是第2时段的用水量，是第1供水时段用水量，是第2供水时段用水量，是一天总用水量.具体的Matlab编程见附录八.计算结果：.4 流量及用水总量的检验计算出的各时刻的流量可用水位记录的数值微分来检验.用水量可用第1时段水位测量记录中下降高度968-822=146来检验，类似地，用1082-822=260检验.根据上述计算结果，显然吻合.供水时段流量的一种检验方法如下：供水时段的用水量加上水位上升值260是该时段泵入的水量，除以时段长度得到水泵的功率（单位时间泵入的水量），而两个供水时段水泵的功率应大致相等.第1、2时段水泵的功率可计算如下：是第1供水时段水泵的功率（水量仍以高度计），设；输出第2供水时段各时刻的流量，是第2供水时段水泵的功率（水量仍以高度计）.具体Matlab编程见附录九.计算结果：.综上所述，流量函数为： 流量曲线见图3-9： 图3-9 流量曲线图参考文献1黄铎,陈兰平,王凤.数值分析M.北京:科学出版社,2004.2Timothy Sauer. Numerical AnalysisM.北京:人民邮电出版社,2010. 3刘志平,石林英.最小二乘法原理应用及其Matlab实现J.中国西部科技,2008,17(7):33- 34. 4费业泰.误差理论与数据处理M.北京:机械工业出版社,2000.5程卫国.MATLAB5.3精要、编程及高级应用M.北京:机械工业出版社,2000.6陆健.最小二乘法及其应用J.中国西部科技,2007,12:19-21.7郑阿奇.MATLAB实用教程M.北京:电子工业出版,2005.8王福昌,曹慧荣,朱红霞.经典最小二乘与全最小二乘法及其参数估计J.统计决策,2009, 36(4):16-17.9徐敏,柳青,马延琴.基于Matlab的最小二乘拟合在大学物理实验中的应用J.西南民族大学学报,2010,36(9):1-3. 10袁占斌,聂玉峰,欧阳洁.基于泰勒基函数的移动最小二乘法及误差分析J.数值计算与计算机应用,2012,33(1):25-31.11唐晓虹.最小二乘原理在桥梁检测中的应用J.中国西部科技,2011,10(18):40-42.12戴朝寿.数理统计简明教程M.北京:高等教育出版社,2009.13胡二邦,王寒.应用最小二乘法原理对我国天气取样候选方案的筛选比较J.辐射防护通讯,2004,24(6):1-8.14李庆扬,王能超,易大义.数值分析M.武汉:华中科技大学出版社,1986.15黄其明,廖鸿志.数值计算方法及其程序设计M.昆明:云南大学出版社,1995.致 谢值此本科学位论文完成之际，首先要感谢我的导师XXX老师.XXX老师从一开始的论文方向的选定，到最后的整篇文论的完成，都非常耐心的对我进行指导.给我提供了大量数据资料和建议，告诉我应该注意的细节问题，细心的给我指出错误.在此，谨向XXX老师致以崇高的敬意和衷心的感谢！论文的完成代表着我大学四年的完结，心里有诸多不舍，因为有你们，使我的人生自信而不张扬，丰富而不纷乱，最后不禁再次向可敬的老师，亲爱的同学，体贴的父母表达我的感激之心，谢谢你们！附 录附录一：t=15 18 21 24 27 30; %输入温度t的实验数据F=43.61 43.63 43.68 43.71 43.74 43.78; %输入对应的测力计示值subplot(1,2,1)plot(t,F,r*) %画出初始数据点t和F的散点图，用*标示B=polyfit(t,F,1) %用最小二乘法对t,F进行一元线性回归y=polyval(B,t); %拟合出的曲线t对应的y值hold onplot(t,y,c) %画出拟合曲线图xlabel(t(t/);ylabel(F/N)x=ones(6,1) t;b, bint,r,rint,stats=regress(F,x) %回归分析检验,bint回归系数的区间估计，r为残 差,rint为置信区间，stats用于检验回归模型的统计量subplot(1,2,2)rcoplot(r,rint) %画出残差及其置信区间附录二：d=300; %t=0注射剂量t=0.25 0.5 1 1.5 2 3 4 6 8; %时间变化 c=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01;%浓度变化y=log(c); %为线性化对浓度取对数a=polyfit(t,y,1) %线性拟合k=-a(1)v=d/exp(a(2)y1=polyval(a,t); %拟合所得曲线对应的浓度值subplot(1,2,1)plot(t,y,*,t,y1) %作出线性化后的曲线图x=ones(9,1) t;b, bint,r,rint,stats=regress(y,x) %回归分析检验,bint回归系数的区间估计，r为残 差,rint为置信区间，stats用于检验回归模型的统计量subplot(1,2,2)rcoplot(r,rint) %画出残差及其置信区间附录三：用M-文件fun2.m定义函数，function f=fun2(x,t)d=300;f=(x(1)/d)*exp(-x(2)*t) ; % x(1)=v; x(2)=k主程序如下：t=0.25 0.5 1 1.5 2 3 4 6 8;c=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01; x0=10,0.5; %设置初始值x=lsqcurvefit(fun2,x0,t,c);f=fun2(x,t) ; subplot(1,2,1)plot(t,c,*,t,f)xlabel(t(h)ylabel(c(ug/mg)v=ones(9,1) t;b, bint,r,rint,stats=regress(f,v) %回归分析检验,bint回归系数的区间估计，r为残 差,rint为置信区间，stats用于检验回归模型的统计量subplot(1,2,2)rcoplot(r,rint) %画出残差及其置信区间附录四：t=0 0.92 1.84 2.95 3.87 4.98 5.9 7.01 7.93 8.97 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.9 20.84 23.88 24.99 25.91;h=968 948 931 913 898 881 869 852 839 822 1082 1050 1021 994 965 941 918 892 866 843 822 1059 1035 1018;c1=polyfit(t(1:10),h(1:10),3) %用3次多项式拟合第1时段水位，c1输出3次多项式的系数a1=polyder(c1) %a1输出多项式（系数为c1）导数的系数tp1=0:0.1:9;x1=-polyval(a1,tp1); %x1输出多项式在tp1时刻的流量（取负为正值）plot(tp1,x1,c.)hold onaxis(0 25 12 34)xlabel(hour)ylabel(cm/houor)附录五：t=0 0.92 1.84 2.95 3.87 4.98 5.9 7.01 7.93 8.97 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.9 20.84 23.88 24.99 25.91;h=968 948 931 913 898 881 869 852 839 822 1082 1050 1021 994 965 941 918 892 866 843 822 1059 1035 1018;c2=polyfit(t(10.9:21),h(10.9:21),4); %用3次多项式拟合第2时段水 位，c2输出4次多项式的系数a2=polyder(c2) % a2输出多项式（系数为c2）导数的系数tp2=10.9:0.1:21;x2=-polyval(a2,tp2); % x2输出多项式在tp2时刻的流量（取负后为正）plot(tp2,x2,c.)hold onaxis(0 25 12 34)xlabel(hour)ylabel(cm/houor)附录六：t=0 0.92 1.84 2.95 3.87 4.98 5.9 7.01 7.93 8.97 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.9 20.84 23.88 24.99 25.91;h=968 948 931 913 898 881 869 852 839 822 1082 1050 1021 994 965 941 918 892 866 843 822 1059 1035 1018;c1=polyfit(t(1:10),h(1:10),3)a1=polyder(c1);c2=polyfit(t(10.9:21),h(10.9:21),4)a2=polyder(c2);xx1=-polyval(a1,8 9); %取第1时段在t=8,9的流量xx2=-polyval(a2,11 12); %取第2时段在t=11,12的流量xx12=xx1 xx2;c12=polyfit(8 9 11 12,xx12,3) %拟合3次多项式a12=polyder(c12)tp12=9:0.1:11;x12=polyval(c12,tp12); % x12输出第1供水时段各时刻的流plot(tp12,x12,c)hold onaxis(0 25 12 34)xlabel(hour)ylabel(cm/houor)附录七：t=0 0.92 1.84 2.95 3.87 4.98 5.9 7.01 7.93 8.97 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.9 20.84 23.88 24.99 25.91;h=968 948 931 913 898 881 869 852 839 822 1082 1050 1021 994 965 941 918 892 866 843 822 1059 1035 1018;c2=polyfit(t(11:21),h(11:21),4);a2=polyder(c2);dt3=diff(t(22:24); %最后3个时刻的两两之差dh3=diff(h(22:24); %最后3个水位的两两之差dht3=-dh3./dt3; %t(22)和t(23)的流量t3=20 20.8 t(22) t(23)xx3=-polyval(a2,t3(1:2),dht3; %取t3各时刻的流量c3=-polyfit(t3,xx3,3); %拟合3次多项式a3=polyder(c3)tp3=20.8:0.1:24;x3=-polyval(c3,tp3); % x3输出第2供水时段（外推至t=24）各时刻的流量plot(tp3,x3,c.)hold on axis(0 25 12 34)xlabel(hour)ylabel(cm/houor)附录八：t=0 0.92 1.84 2.95 3.87 4.98 5.9 7.01 7.93 8.97 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.9 20.84 23.88 24.99 25.91;h=968 948 931 913 898 881 869 852 839 822 1082 1050 1021 994 965 941 918 892 866 843 822 1059 1035 1018;c