《数模插值与拟合建模.ppt》由会员分享，可在线阅读，更多相关《数模插值与拟合建模.ppt（133页珍藏版）》请在三一办公上搜索。
1、第八章 插值与拟合建模,重庆邮电大学 数理学院沈世云,一.插值,主要内容,2、掌握用数学软件包求解插值问题。,1、了解插值的基本内容。,1一维插值,2二维插值,3实验作业,拉格朗日插值,分段线性插值,三次样条插值,一 维 插 值,一、插值的定义,二、插值的方法,三、用Matlab解插值问题,返回,返回,二维插值,一、二维插值定义,二、网格节点插值法,三、用Matlab解插值问题,最邻近插值,分片线性插值,双线性插值,网格节点数据的插值,散点数据的插值,一维插值的定义,返回,2.3 拉格朗日插值,1、线性插值,假设已知,在区间,上的两个结点,和它们的函数值,求一个一次多项式,，使得多项式,在结点
2、上满足条件,这种插值方法称为线性插值方法（也称两点插值）。,可以求出：,分段线性插值,计算量与n无关;n越大，误差越小.,2、抛物插值,已知,在区间,上的三个结点,和它们的函数值,求一个次数不超过2的多项式,，使得它在结点上满足条件,这种插值方法称为抛物线插值法，,可求出：,3、n次拉格朗日插值,假设取区间,上的,个结点,，并且已知函数,在这此点的函数值,现在求一个次数不超过,的多项式,，使得满足条件,这种插值方法称为,次多项式插值（或称代数插值），,利用拉格朗日插值插值方法可得,上述多项式称为 n次拉格朗日（Lagrange）插值多项式，,函数,称为拉格朗日插值基函数。,当n=1,2时，n次
3、拉格朗日（Lagrange）插值多项式即为线性插值多项式和抛物插值多项式。,比分段线性插值更光滑。,在数学上，光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续，则称该曲线具有k阶光滑性。光滑性的阶次越高，则越光滑。是否存在较低次的分段多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例子。,三次样条插值,三次样条插值,g(x)为被插值函数。,例1 已知函数发f(x)的函数表如下：,求其拉格朗日插值多项式，并求,的近似值。,解 由于给出了4个插值结点，所以可做出次数不超过3的拉格朗日插值多项式。,将上列4式代入n=3 的拉格朗日插值公式，可得所要求的插值多项式为,将x=2.5代入可得
4、f(2.5)的近似值为1.8496。,拉格朗日插值法适合节点较少的情况，当节点较多的大范围高次插值的逼近效果往往并不理想且当插值结点增加时，计算越来越繁。为了提高精度和减少计算，还有牛顿插值法下、三次样条插值等，具体可参阅有关书籍。,用MATLAB作插值计算,一维插值函数：,yi=interp1(x，y，xi，method),nearest：最邻近插值linear：线性插值；spline：三次样条插值；cubic：立方插值。缺省时：分段线性插值。,注意：所有的插值方法都要求x是单调的，并且xi不能够超过x的范围。,例：在1-12的11小时内，每隔1小时测量一次温度，测得的温度依次为：5，8，9
5、，15，25，29，31，30，22，25，27，24。试估计每隔1/10小时的温度值。,To MATLAB(temp),hours=1:12;temps=5 8 9 15 25 29 31 30 22 25 27 24;h=1:0.1:12;t=interp1(hours,temps,h,spline);(直接输出数据将是很多的)plot(hours,temps,+,h,t,hours,temps,r:)%作图xlabel(Hour),ylabel(Degrees Celsius),例2 在用外接电源给电容器充电时，电容器两端的电压V将会随着充电时间t发生变化，已知在某一次实验时，通过测量得
6、到下列观测值，分别用拉格朗日插值法、分段线性插值法、三次样条插值法画出V随着时间t变化的曲线图，分别计算当时间t=7s时，三种插值法各自算得电容器两端电压的近似值。,解 由于MATLAB没有提供现成的拉格朗日插值命令，我们可以编写一个函数lglrcz.m来完成，其他两种插值法可用现成的命令。,用MATLAB软件进行三种插值计算的程序为szczqx.m。,程序lglrcz.m：function y=lglrcz(x0,y0,x)n=length(x0)；m=length(x)；for i=1:m z=x(i)；s=0.0；for k=1:n p=1.0；for j=1:n if j=k p=p*
7、(z-x0(j)/(x0(k)-x0(j)；end end s=p*y0(k)+s；end y(i)=s；end,程序为azczf.mt=1,2,3,4,6.5,9,12；v=6.2,7.3,8.2,9.0,9.6,10.1,10.4；t0=0.2:0.1:12.5；lglr=lglrcz(t,v,t0)；laglr=lglrcz(t,v,7)；fdxx=interp1(t,v,t0)；fendxx=interp1(t,v,7)；scyt=interp1(t,v,t0,spline)；sancyt=interp1(t,v,7,spline)plot(t,v,*,t0,lglr,r,t0,fdx
8、x,g,t0,scyt,b)gtext(lglr)gtext(fdxx)gtext(scyt),执行结果是laglr=9.52988980716254fendxx=9.70000000000000sancyt=9.67118039327444,图形如图1所示。,图中曲线lglr表示拉格朗日插值曲线，scyt表示三次样条插值曲线，fdxx表示分段线性插值曲线。,案例:水流量的估计,美国某州的用水管理机构要求各社区提供以每小时多少加仑计的用水量以及每天所用的总水量。许多社区没有测量流入或流出水塔水量的装置，只能代之以每小时测量水塔中的水位，其误差不超过5%。需要注意的是，当水塔中的水位下降到最低水
9、位L时，水泵就自动向水塔输水直到最高水位H，此期间不能测量水泵的供水量，因此，当水泵正在输水时不容易建立水塔中水位和用水量之间的关系。水泵每天输水一次或两次，每次约2小时。,2.1 问题,试估计任何时刻（包括水泵正在输水时间）从水塔流出的水流量f(t)，并估计一天的总用水量。已知该水塔是一个高为40ft（英尺），直径为57ft（英尺）的正圆柱，表5-1给出了某个小镇一天水塔水位的真实数据，水位降至约27.00ft水泵开始工作，水位升到35.50ft停止工作。（注：1ft（英尺）=0.3048m（米）,表1 某小镇某天水塔水位,2.2 问题分析,本实验所指流量可视为单位时间内流出水的体积。由于水
10、塔是正圆柱形，横截面积是常数，所以在水泵不工作时段，流量很容易根据水位相对时间的变化算出。问题的难点在于如何估计水泵供水时段的流量。,水泵供水时段的流量只能靠供水时段前后的流量经插值或拟合得到。作为用于插值或拟合的原始数据，我们希望水泵不工作时段的流量越准确越好。这此流量大体上可由两种方法计算，一是直接对表12-1中的水量用数值微分算出各时段的流量，用它们拟合其它时刻或连续时间的流量；二是先用表中数据拟合水位一时间函数，求导数即可得到连续时间的流量。,有了任何时刻的流量，就不难计算一天的总用水量。其实，水泵不工作时段的用水量可以由测量记录直接得到，由表1中下降水位乘以水搭的截面积就是这一时段的
11、用水量。这个数值可以用来检验数据插值或拟合的结果。,在具体给出本问题的解答之前，先介绍一个简单的数据插值方法。,2.2 问题求解,为了表示方便，我们将2.1节问题中所给表1中的数据全部化为国际标准单位（表2），时间用小时（h），高度用米（m）：,表12-2 一天内水塔水位记录,1模型假设,（1）流量只取决于水位差，与水位本身无关，故由物理学中Torriceli定律：小孔流出的液体的流速正比于水面高度的平方根。题目给出水塔的最低和最高水位分别是8.1648m,和10.7352m,（设出口的水位为零），因为sqrt,，约为1，所以可忽略水位对流速的影响。,（2）将流量看作是时间的连续函数，为计算简
12、单，不妨将流量定义成单位时间流出水的高度，即水位对时间变化率的绝对值（水位是下降的），水塔截面积为,（m2），得到结果,后乘以s即可。,2流量估计方法,首先依照表2所给数据，用MATLAB作出时间水位散点图（图2）。,下面来计算水箱流量与时间的关系。,根据图1.，一种简单的处理方法为，将表2中的数据分为三段，然后对每一段的数据做如下处理：,设某段数据,，相邻数据中点的平均流速用下面的公式（流速=（右端点的水位右端点的水位）/区间长度）：,每段数据首尾点的流速用下面的公式计算：,用以上公式求得时间与流速之间的数据如表3。,由表3作出时间流速散点图如图3。,1）插值法,由表12-3，对水泵不工作时
13、段1，2采取插值方法，可以得到任意时刻的流速，从而可以知道任意时刻的流量。,我们分别采取拉格朗日插值法，分段线性插值法及三次样条插值法；对于水泵工作时段1应用前后时期的流速进行插值，由于最后一段水泵不工作时段数据太少，我们将它与水泵工作时段2合并一同进行插值处理（该段简称混合时段）。,我们总共需要对四段数据（第1，2未供水时段，第1供水时段，混合时段）进行插值处理，下面以第1未供水时段数据为例分别用三种方法算出流量函数和用水量（用水高度）。,下面是用MATLAB实现该过程的程序。,t=0,0.46,1.38,2.395,3.41,4.425,12.44,6.45,7.465,8.45,8.97
14、；v=29.89,21.74,18.48,16.22,16.30,15.32,13.04,15.45,13.98,16.35,19.29；t0=0:0.1:8.97；,lglr=lglrcz(t,v,t0)；/*注：lglrcz为一函数，程序同lglrcz.m*/lglrjf=0.1*trapz(lglr)fdxx=interpl(t,v,t0)；fdxxjf=0.1*trapz(fdxx)scyt=interpl(t,v,t0,spline)；sancytjf=0.1*trapz(scyt)plot(t,v,*,t0,lglr,r,t0,fdxx,g,t0,scyt,b)gtext(lglr
15、)gtext(fdxx)gtext(scyt),运行结果lglrjf=145.6231；fdxxjf=147.1430；sancytjf=1412.6870,图4是对第1未供水段数据用三种不同方法得到的插值函数图，,图中曲线lglr、fdxx和scyt分别表示用拉格朗日插值法，分段线性插值法及三次样条插值法得到的曲线。,由表2知，第1未供水时段的总用水高度为146(=968822)，可见上述三种插值方法计算的结果与实际值（146）相比都比较接近。考虑到三次样条插值方法具有更加良好的性质，建议采取该方法。,其他三段的处理方法与第1未供水时段的处理方法类似，这里不再详细叙述，只给出数值结果和函数图
16、像（图5图7），图中曲线标记同图4。,图5 第一供水段时间流速示意图,图6 第2未供水段时间流速示意图,图7 混合时段时间流速示意图,图8是用分段线性及三次样条插值方法得到的整个过程的时间流速函数示意图。,表4 各时段及一天的总用水量（用水高度）,表5是对一天中任取的4个时刻分别用3种方法得到的水塔水流量近似值。,注：拉格朗日插值法分段线性插值法三次样条插值法,2）拟合法,（1）拟合水位时间函数,从表12-2中的测量记录看，一天有两次供水时段和三次未供水时段，分别对第1，2未供水时段的测量数据直接作多项式拟合，可得到水位函数（注意，根据多项式拟合的特点，此处拟合多项式的次数不宜过高，一般以36
17、次为宜）。对第3未供水时段来说，数据过少不能得到很好的拟合。,设t，h分别为已输入的时刻和水位测量记录（由表12.2提供，水泵启动的4个时刻不输入），这样第1未供水时段各时刻的水位可由如下MATLAB程序完成：,t=0,0.92,1.84,2.95,3.87,4.98,5.90,7.00,7.93,8.97,10.95,12.03,12.95,13.88,14.98,15.90,16.83,17.93,19.04,19.96,20.84,23.8824.99,25.66h=9.68,9.48,9.31,9.13,8.98,8.81,8.69,8.52,8.39,8.22,10.82,10.50
18、,10.21,9.94,9.65,9.41,9.18,8.92,8.66,8.43,8.22,10.59,10.35,10.18；c1=polyfit(t(1:10),h(1:10),3)；tp1=0:0.1:8.9；x1=polyval(c1,tp1)；plot(tp1,x1),图12.9给出的是第1未供水时段的时间水位拟合函数图形。,变量x1中存放了以0.1为步长算出的各个时刻的水位高度。同样地，第2未供水时段时间水位图可由如下MATLAB程序完成，读者可自己上机运行查看。,c2=polyfit(t(11:21),h(11:21),3)；tp2=10.9:0.1:20.9；x2=-poly
19、val(c2,tp2)；plot(tp2,x2),（2）确定流量时间函数,对于第1，2未供水时段的流量可直接对水位函数求导，程序如下：,c1=polyfit(t(1:10),h(1:10),3)；c2=polyfit(t(11:21),h(11:21),3)；a1=polyder(c1)；a2=polyder(c2)；tp1=0:0.01:8.97；tp2=10.95:0.01:20.84；x13=-polyval(a1,tp1)；x113=-polyval(a1,0:0.01:8.97)；wgsysl1=100*trapz(tp1,x113)；*/计算第1未供水时段的总用水量/*x14=-p
20、olyval(a1,7.93,8.97)；*/为下面的程序准备数据/*x23=-polyval(a2,tp2)；x114=-polyval(a2,10.95:0.01:20.84)wgsysl2=100*trapz(tp2,x114)；*/计算第2未供水时段的总用水量/*x24=-polyval(a2,10.95,12.03)；*/为下面的程序准备数据/*x25=-polyval(a2,19.96,20.84)；*/为下面的程序准备数据/*subplot(1,2,1)plot(tp1,x13*100)*/与图12.10单位保持一致/*subplot(1,2,2)plot(tp2,x23*100
21、)*/与图12.10单位保持一致/*,程序运行得到第1，2未供水时段的时间流量图如图10，可以看到与图8中用插值给出的曲线比较吻合。,如果用5次多项式拟合则得图11的图形，显然较三次拟合的效果好。,而第1供水时段的流量则用前后时期的流量进行拟合得到。为使流量函数在t=9和t=11连续，我们只取4个点，用三次多项式拟合得到第1供水时段的时间流量图形如图12，可以看到与图8中的相应部分比较吻合。,图12.12,图12.8,dygsdsy=7.93,8.97,10.95,12.03；dygsdls=x14,x24；nhjg=polyfit(dygsdsj,dygsdls,3)；nhsj=7.93:0
22、.1:12.03；nhlsjg=polyval(nhjg,nhsj)；gssj1=8.97:0.01:10.95；gs1=polyval(nhjg,8.97:0.01:10.95)；gsysl1=100*trapz(gssj1,gs1)；*/该语句计算第1供水时段的总用水量/*plot(nhsj,100*nhlsjg),程序如下：,在第2供水时段之前取t=19.96,20.84两点的流量，用第3未供水时段的3个记录做差分得到两个流量数据21.62，18.48，然后用这4个数据做三次多项式拟合得到第2供水时段与第3未供水时段的时间流量图如图13，可以看到与图8中的相应部分也比较吻合。,图12.1
23、3，,图12.8,程序如下：,t3=19.96,20.84,t(22),t(23)；ls3=x25*100,21.62,18.48；nhhddxsxs=polyfit(t3,ls3,3)；tp3=19.96:0.01:25.91；xx3=polyval(nhhddxsxs,tp3)；gssj2=20.84:0.01:24；gs2=polyval(nhhddxsxs,20.84:0.01:24)；gsysl2=trapz(gssj2,gs2)；*/该语句计算第2供水时段的总用水量/*plot(tp3,xx3),（3）一天总用水量的估计,分别对供水的两个时段和不供水的两个时段积分（流量对时间）并求
24、和得到一天的总用水量约为526.8935（此数据是总用水高度，单位为cm）。表12-6列出各段用水量，与插值法算得的表12-4相比，二者较为吻合。,表12-6,（4）流量及总用水量的检验,计算出各时刻的流量可用水位记录的数值微分来检验，各时段的用水高度可以用实际记录的水位下降高度来检验。,例如，算得第1未供水段的用水量高度是145.67，而实际记录的水位下降高度为968822=146cm，两者是吻合的；,算得第2未供水段的用水量高度是260.66cm，而实际记录的水位下降高度为1082822=260cm，两者也是吻合的。,从算法设计和分析可知，计算结果与各时段所用的拟合多项式的次数有关。表12
25、-7给出的是对第1，2未供水时段分别用五、六多项式拟合后得到的用水量结果。,表12-7,3）结果分析,由表12-4可以看出，使用三次样条插值法得到的结果,260cm相差不大，说明插值结果与原始数据比较吻合。,与表12-2中记录的下降高度,146cm，,按三次样条插值法估计出全天的用水量约为,由表12-7可以得全天的用水量约为,2.5 内容小结,本实验主要进行水塔水流量的估算。第一种估算方法为插值方法，我们用了三种不同的插值法进行估计，在求解的过程中，可以熟悉数据插值的理论和方法；第二种估算方法为数据拟合法，用多项式进行拟合，得到水塔水流量的估计。,2二维插值命令interp2的具体使用格式,z
26、z=interp2(x,y,z,xx,yy,method),该指令的意思是根据数据向量x，y，z按method指定的方法来做插值，然后将xx，yy处插值函数的插值结点向量，如果xx，yy在插值范围之内，则返回值在zz中，否则返回值为空NaN。method是插值方法可选项，具体要求同一维插值的情况。,该命令还有以下几种省略格式：zz=interp2(z,xx,yy)zz=interp2(x,y,z,xx,yy)zz=interp2(z,ntimes),3三维插值命令interp3的具体使用格式,vi=interp3(x,y,z,v,xi,yi,zi,method),它的具体含义跟前面的一、二维插
27、值是相似的，在此不作解释，读者可在MATLAB工作空间中用help interp3命令获得。,4样条插值命令spline的具体使用格式,yy=spline(x,y,xx),它的意思等同于命令yy=interp1(x,y,xx,cubic),例 已知飞机下轮廓线上数据如下，求x每改变0.1时的y值。,To MATLAB(plane),返回,二维插值的定义,第一种（网格节点）：,已知 mn个节点,第二种（散乱节点）：,返回,注意：最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。,最邻近插值,二维或高维情形的最邻近插值，与被插值点最邻近的节点的函数值即为所求。,返回,将四个插值点（矩形
28、的四个顶点）处的函数值依次简记为：,分片线性插值,f(xi,yj)=f1，f(xi+1,yj)=f2，f(xi+1,yj+1)=f3，f(xi,yj+1)=f4,插值函数为：,第二片(上三角形区域)：(x,y)满足,插值函数为：,注意：(x,y)当然应该是在插值节点所形成的矩形区域内。显然，分片线性插值函数是连续的；,分两片的函数表达式如下：,第一片(下三角形区域)：(x,y)满足,返回,双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下：,其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。,双线性插值,返回,要求x0,y0单
29、调；x，y可取为矩阵，或x取行向量，y取为列向量，x,y的值分别不能超出x0,y0的范围。,z=interp2(x0,y0,z0,x,y,method),用MATLAB作网格节点数据的插值,nearest 最邻近插值linear 双线性插值cubic 双三次插值缺省时,双线性插值,例：测得平板表面3*5网格点处的温度分别为：82 81 80 82 84 79 63 61 65 81 84 84 82 85 86 试作出平板表面的温度分布曲面z=f(x,y)的图形。,输入以下命令：x=1:5;y=1:3;temps=82 81 80 82 84;79 63 61 65 81;84 84 82 8
30、5 86;mesh(x,y,temps),1.先在三维坐标画出原始数据，画出粗糙的温度分布曲图.,2以平滑数据,在x、y方向上每隔0.2个单位的地方进行插值.,再输入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi)画出插值后的温度分布曲面图.,To MATLAB(wendu),通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插值效果进行比较。,To MATLAB(moutain),返回,插值函数griddata格式为:,cz=griddata（x，y，z，cx，cy，method）,用
31、MATLAB作散点数据的插值计算,要求cx取行向量，cy取为列向量。,nearest 最邻近插值linear 双线性插值cubic 双三次插值v4-Matlab提供的插值方法缺省时,双线性插值,下面再给出一个二维插值应用的例子,x=1:5;y=1:5;t=100,100,100,100,100;105,120,122,125,122;110,130,155,157,130;115,133,157,160,140;113,132,149,154,128t=100 100 100 100 100 105 120 122 125 122 110 130 155 157 130 115 133 157
32、 160 140 113 132 149 154 128,mesh(x,y,t)数据原图,xx=1:0.1:5;yy=1:0.1:5;tt=interp2(x,y,t,xx,yy,cubic);mesh(xx,yy,tt)二维插值得到的示意图,例 在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）*（-50，150）里的哪些地方船要避免进入。,To MATLAB hd1,返回,4.作出水深小于5的海域范围,即z=5的等高线.,实验作业,山区地貌：在某山区测得一些地点的高程如下表：(平面区域1200=x=4000,1200=y=3600)，试作出该
33、山区的地貌图和等高线图，并对几种插值方法进行比较。,返回,二.拟合,主要内容,2、掌握用数学软件求解拟合问题。,1、直观了解拟合基本内容。,1、拟合问题引例及基本理论。,4、实验作业。,2、用数学软件求解拟合问题。,3、应用实例,拟 合,2.拟合的基本原理,1.拟合问题引例,拟 合 问 题 引 例 1,求600C时的电阻R。,设 R=at+ba,b为待定系数,拟 合 问 题 引 例 2,求血药浓度随时间的变化规律c(t).,作半对数坐标系(semilogy)下的图形,MATLAB(aa1),曲 线 拟 合 问 题 的 提 法,已知一组（二维）数据，即平面上 n个点（xi,yi)i=1,n,寻求
34、一个函数（曲线）y=f(x),使 f(x)在某种准则下与所有数据点最为接近，即曲线拟合得最好。,y=f(x),i 为点（xi,yi)与曲线 y=f(x)的距离,拟合与插值的关系,函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。,实例：下面数据是某次实验所得，希望得到X和 f之间的关系？,MATLAB(cn),问题：给定一批数据点，需确定满足特定要求的曲线或曲面,解决方案：,若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，这就是数据拟合，又称曲线拟合或曲面拟合。,若要求所求曲线（面）通过所给所有数据点，就是插值问题；
35、,最临近插值、线性插值、样条插值与曲线拟合结果：,曲线拟合问题最常用的解法线性最小二乘法的基本思路,第一步:先选定一组函数 r1(x),r2(x),rm(x),mn,令 f(x)=a1r1(x)+a2r2(x)+amrm(x)（1）其中 a1,a2,am 为待定系数。,第二步:确定a1,a2,am 的准则（最小二乘准则）：使n个点（xi,yi)与曲线 y=f(x)的距离i 的平方和最小。,记,问题归结为，求 a1,a2,am 使 J(a1,a2,am)最小。,线性最小二乘法的求解：预备知识,超定方程组：方程个数大于未知量个数的方程组,超定方程一般是不存在解的矛盾方程组。,如果有向量a使得 达到
36、最小，则称a为上述超定方程的最小二乘解。,线性最小二乘法的求解,定理：当RTR可逆时，超定方程组（3）存在最小二乘解，且即为方程组 RTRa=RTy的解：a=(RTR)-1RTy,所以，曲线拟合的最小二乘法要解决的问题，实际上就是求以下超定方程组的最小二乘解的问题。,线性最小二乘拟合 f(x)=a1r1(x)+amrm(x)中函数r1(x),rm(x)的选取,1.通过机理分析建立数学模型来确定 f(x)；,2.将数据(xi,yi)i=1,n 作图，通过直观判断确定 f(x)：,用MATLAB解拟合问题,1、线性最小二乘拟合,2、非线性最小二乘拟合,用MATLAB作线性最小二乘拟合,1.作多项式
37、f(x)=a1xm+amx+am+1拟合,可利用已有程序:,a=polyfit(x,y,m),2.对超定方程组,3.多项式在x处的值y可用以下命令计算：y=polyval（a，x）,例 对下面一组数据作二次多项式拟合,1）输入以下命令：x=0:0.1:1;y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2;R=(x.2)x ones(11,1)；A=Ry,MATLAB(zxec1),解法1用解超定方程的方法,2）计算结果：=-9.8108 20.1293-0.0317,1）输入以下命令：x=0:0.1:1;y=-0.447 1
38、.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2;A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,k+,x,z,r)%作出数据点和拟合曲线的图形,2）计算结果：=-9.8108 20.1293-0.0317,解法2用多项式拟合的命令,MATLAB(zxec2),1.lsqcurvefit已知数据点：xdata=（xdata1，xdata2，xdatan），ydata=（ydata1，ydata2，ydatan）,用MATLAB作非线性最小二乘拟合,Matlab的提供了两个求非线性最小二乘拟合的函数：lsqcurv
39、efit和lsqnonlin。两个命令都要先建立M-文件fun.m，在其中定义函数f(x)，但两者定义f(x)的方式是不同的,可参考例题.,lsqcurvefit用以求含参量x（向量）的向量值函数F(x,xdata)=（F（x，xdata1），F（x，xdatan）T中的参变量x(向量),使得,输入格式为:(1)x=lsqcurvefit(fun,x0,xdata,ydata);(2)x=lsqcurvefit(fun,x0,xdata,ydata,options);(3)x=lsqcurvefit(fun,x0,xdata,ydata,options,grad);(4)x,options=l
40、sqcurvefit(fun,x0,xdata,ydata,);(5)x,options,funval=lsqcurvefit(fun,x0,xdata,ydata,);(6)x,options,funval,Jacob=lsqcurvefit(fun,x0,xdata,ydata,);,说明：x=lsqcurvefit(fun,x0,xdata,ydata,options);,lsqnonlin用以求含参量x（向量）的向量值函数 f(x)=(f1(x),f2(x),fn(x)T 中的参量x，使得 最小。其中 fi（x）=f（x，xdatai，ydatai）=F(x,xdatai)-ydata
41、i,2.lsqnonlin,已知数据点：xdata=（xdata1，xdata2，xdatan）ydata=（ydata1，ydata2，ydatan）,输入格式为：1）x=lsqnonlin（fun，x0）；2）x=lsqnonlin（fun，x0，options）；3）x=lsqnonlin（fun，x0，options，grad）；4）x，options=lsqnonlin（fun，x0，）；5）x，options，funval=lsqnonlin（fun，x0，）；,说明：x=lsqnonlin（fun，x0，options）；,例2 用下面一组数据拟合 中的参数a，b，k,该问题即解
42、最优化问题：,MATLAB(fzxec1),1）编写M-文件 curvefun1.m function f=curvefun1(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata)%其中 x(1)=a;x(2)=b；x(3)=k;,2）输入命令tdata=100:100:1000cdata=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;x0=0.2,0.05,0.05;x=lsqcurvefit(curvefun1,x0,tdata,cdata)f=curvefun1(x,tdata),F(x，td
43、ata)=，x=(a，b，k),解法1.用命令lsqcurvefit,3）运算结果为：f=0.0043 0.0051 0.0056 0.0059 0.0061 0.0062 0.0062 0.0063 0.0063 0.0063 x=0.0063-0.0034 0.2542,4）结论：a=0.0063,b=-0.0034,k=0.2542,MATLAB(fzxec2),解法 2 用命令lsqnonlin f(x)=F(x,tdata,ctada)=x=（a，b，k）,1）编写M-文件 curvefun2.m function f=curvefun2(x)tdata=100:100:1000;c
44、data=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;f=x(1)+x(2)*exp(-0.02*x(3)*tdata)-cdata,2）输入命令:x0=0.2,0.05,0.05;x=lsqnonlin(curvefun2,x0)f=curvefun2(x),函数curvefun2的自变量是x，cdata和tdata是已知参数，故应将cdata tdata的值写在curvefun2.m中,3）运算结果为 f=1.0e-003*(0.2322-0.1243-0.2495-0.2413-0.1668-0.0724 0.0241
45、0.1159 0.2030 0.2792 x=0.0063-0.0034 0.2542,可以看出,两个命令的计算结果是相同的.,4）结论：即拟合得a=0.0063 b=-0.0034 k=0.2542,MATLAB解应用问题实例,1、电阻问题,2、给药方案问题,*3、水塔流量估计问题,MATLAB(dianzu1),电阻问题,得到 a1=3.3940,a2=702.4918,方法2.直接用,结果相同。,MATLAB(dianzu2),一室模型：将整个机体看作一个房室，称中心室，室内血药浓度是均匀的。快速静脉注射后，浓度立即上升；然后迅速下降。当浓度太低时，达不到预期的治疗效果；当浓度太高，又可
46、能导致药物中毒或副作用太强。临床上，每种药物有一个最小有效浓度c1和一个最大有效浓度c2。设计给药方案时，要使血药浓度 保持在c1c2之间。本题设c1=10，c2=25(ug/ml).,一种新药用于临床之前，必须设计给药方案.,药物进入机体后血液输送到全身，在这个过程中不断地被吸收、分布、代谢，最终排出体外，药物在血液中的浓度，即单位体积血液中的药物含量，称为血药浓度。,要设计给药方案,必须知道给药后血药浓度随时间变化的规律。从实验和理论两方面着手：,给药方案,1.在快速静脉注射的给药方式下，研究血药浓度（单位体积血液中的药物含量）的变化规律。,t,问题,2.给定药物的最小有效浓度和最大治疗浓
47、度，设计给药方案：每次注射剂量多大；间隔时间多长。,分析,理论：用一室模型研究血药浓度变化规律,实验：对血药浓度数据作拟合，符合负指数变化规律,3.血液容积v,t=0注射剂量d,血药浓度立即为d/v.,2.药物排除速率与血药浓度成正比，比例系数 k(0),模型假设,1.机体看作一个房室，室内血药浓度均匀一室模型,模型建立,在此，d=300mg，t及c（t）在某些点处的值见前表，需经拟合求出参数k、v,用线性最小二乘拟合c(t),MATLAB(lihe1),计算结果：,用非线性最小二乘拟合c(t),给药方案 设计,设每次注射剂量D,间隔时间,血药浓度c(t)应c1 c(t)c2,初次剂量D0 应
48、加大,给药方案记为：,2、,1、,计算结果：,给药方案：,c1=10,c2=25k=0.2347v=15.02,故可制定给药方案：,即:首次注射375mg，其余每次注射225mg，注射的间隔时间为4小时。,估计水塔的流量,2、解题思路,3、算法设计与编程,1、问题,某居民区有一供居民用水的园柱形水塔，一般可以通过测量其水位来估计水的流量，但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间无法测量水塔的水位和水泵的供水量通常水泵每天供水一两次，每次约两小时.水塔是一个高12.2米，直径17.4米的正园柱按照设计，水塔水位降至约8.2米时
49、，水泵自动启动，水位升到约10.8米时水泵停止工作表1 是某一天的水位测量记录，试估计任何时刻（包括水泵正供水时）从水塔流出的水流量，及一天的总用水量,流量估计的解题思路,拟合水位时间函数,确定流量时间函数,估计一天总用水量,拟合水位时间函数 测量记录看，一天有两个供水时段（以下称第1供水时段和第2供水时段），和3个水泵不工作时段（以下称第1时段t=0到t=8.97，第2次时段t=10.95到t=20.84和第3时段t=23以后）对第1、2时段的测量数据直接分别作多项式拟合，得到水位函数为使拟合曲线比较光滑，多项式次数不要太高，一般在36由于第3时段只有3个测量记录，无法对这一时段的水位作出较
50、好的拟合,2、确定流量时间函数 对于第1、2时段只需将水位函数求导数即可，对于两个供水时段的流量，则用供水时段前后（水泵不工作时段）的流量拟合得到，并且将拟合得到的第2供水时段流量外推，将第3时段流量包含在第2供水时段内,3、一天总用水量的估计 总用水量等于两个水泵不工作时段和两个供水时段用水量之和，它们都可以由流量对时间的积分得到。,算法设计与编程,1、拟合第1、2时段的水位，并导出流量,2、拟合供水时段的流量,3、估计一天总用水量,4、流量及总用水量的检验,1、拟合第1时段的水位，并导出流量 设t，h为已输入的时刻和水位测量记录（水泵启动的4个时刻不输入），第1时段各时刻的流量可如下得：1
