1、第八章第八章 插值与拟合建模插值与拟合建模 重庆邮电大学重庆邮电大学 数理学院数理学院沈世云沈世云一一.插值插值主要内容主要内容2、掌握用数学软件包求解插值问题。、掌握用数学软件包求解插值问题。1、了解插值的基本内容。、了解插值的基本内容。11一维插值一维插值22二维插二维插值值33实验作业实验作业拉格朗日插值拉格朗日插值分段线性插值分段线性插值三次样条插值三次样条插值一一 维维 插插 值值一、一、插值的定义插值的定义二、插值的方法二、插值的方法三、用三、用Matlab解插值问题解插值问题返回返回返回返回二维插值二维插值一、一、二维插值定义二维插值定义二、网格节点插值法二、网格节点插值法三、三
2、、用用MatlabMatlab解插值问题解插值问题最邻近插值最邻近插值分片线性插值分片线性插值双线性插值双线性插值网格节点数据的插值网格节点数据的插值散点数据的插值散点数据的插值一维插值的定义一维插值的定义已知已知 n+1个节点个节点,1,0(),(njyxjj其中其中jx互不相同，不妨设互不相同，不妨设),10bxxxan求任一插值点求任一插值点)(*jxx 处的插值处的插值.*y0 x1xnx0y1y节点可视为由节点可视为由)(xgy 产生产生,，g表达式复杂表达式复杂,，或无封闭形式或无封闭形式,，或未知或未知.。*x*y 构造一个构造一个(相对简单的相对简单的)函数函数),(xfy 通
3、过全部节点通过全部节点,即即),1,0()(njyxfjj再用再用)(xf计算插值，即计算插值，即).(*xfy 0 x1xnx0y1y*x*y返回返回2.3 拉格朗日插值拉格朗日插值1、线性插值)(xf假设已知在区间,ba上的两个结点10,xx和它们的函数值 10 xxx )()()(10 xfxfxf求一个一次多项式xaaxP101)(，使得多项式)(1xP在结点上满足条件 1,0),()(1ixfxPii 这种插值方法称为线性插值方法（也称两点插值）。可以求出：)()()(101001011xfxxxxxfxxxxxP分段线性插值分段线性插值其它,0,)()()(1111110jjjjj
4、jjjjjjnjjjnxxxxxxxxxxxxxxxlxlyxL计算量与n无关;n越大，误差越小.nnnxxxxgxL0),()(limxjxj-1xj+1x0 xnxoy2、抛物插值)(xf已知在区间,ba上的三个结点210,xxx和它们的函数值 210 xxxx )()()()(210 xfxfxfxf求一个次数不超过2的多项式22102)(xaxaaxP，使得它在结点上满足条件2,1,0),()(2ixfxPii这种插值方法称为抛物线插值法，可求出：)()()()()()()()()()()(121012002010212xfxxxxxxxxxfxxxxxxxxxP)()()()()(2
5、021210 xfxxxxxxxx3、n次拉格朗日插值,ba假设取区间 上的1n个结点nxxx,10，并且已知函数)(xf在这此点的函数值 nxxxxx210 )()()()()(210nxfxfxfxfxf现在求一个次数不超过n的多项式nnnxaxaxaaxP2210)(，使得满足条件 nixfxPiin,1,0),()(这种插值方法称为n次多项式插值（或称代数插值），利用拉格朗日插值插值方法可得)()()()()()()()()(,2,21,10,0nnnnnnnxfxlxfxlxfxlxfxlxP0()0,1,njknjkjj kxxlxknxx其中上述多项式称为 n次拉格朗日（Lagr
6、ange）插值多项式，函数nkxlnk,1,0),(,称为拉格朗日插值基函数。当n=1,2时，n次拉格朗日（Lagrange）插值多项式即为线性插值多项式和抛物插值多项式。比分段线性插值更光滑。比分段线性插值更光滑。xyxi-1 xiab 在数学上，光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续，则称该曲线具有k阶光滑性。光滑性的阶次越高，则越光滑。是否存在较低次的分段多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例子。三次样条插值三次样条插值 三次样条插值,1,),()(1nixxxxsxSiii,)()3),1,0()()2),1()()10223niiiiiiixxCx
7、SniyxSnidxcxbxaxs)1,1()()(),()(),()(111 nixsxsxsxsxsxsiiiiiiiiiiii自然边界条件）(0)()()40 nxSxS)(,)4)3)2xSdcbaiiiig g(x x)为被插值函数为被插值函数。)()(limxgxSn4321x例1 已知函数发f(x)的函数表如下：31.28.15.2)(xf求其拉格朗日插值多项式，并求)5.2(f的近似值。解 由于给出了4个插值结点，所以可做出次数不超过3的拉格朗日插值多项式。43132361)41()31()21()4()3()2()(233,0 xxxxxxxl6219421)42()32()
8、22()4()3()2()(233,1xxxxxxxl472721)43()33()23()4()3()2()(233,2xxxxxxxl161161)44()34()24()4()3()2()(233,3xxxxxxxl 将上列4式代入n=3 的拉格朗日插值公式，可得所要求的插值多项式为6.4933.29.00667.0)(233xxxxP将x=2.5代入可得f(2.5)的近似值为1.8496。拉格朗日插值法适合节点较少的情况，当节点较多的大范围高次插值的逼近效果往往并不理想且当插值结点增加时，计算越来越繁。为了提高精度和减少计算，还有牛顿插值法下、三次样条插值等，具体可参阅有关书籍。用用M
9、ATLABMATLAB作插值计算作插值计算一维插值函数：一维插值函数：yi=interp1(x，y，xi，method)插值方法插值方法被插值点被插值点插值节点插值节点xixi处的插处的插值结果值结果nearest ：最邻近插值：最邻近插值linear ：线性插值；线性插值；spline ：三次样条插三次样条插值；值；cubic ：立方插值。立方插值。缺省时：缺省时：分段线性插值。分段线性插值。注意：所有的插值方法都要求注意：所有的插值方法都要求x x是单调的，并且是单调的，并且xi不不能够超过能够超过x的范围。的范围。例：在例：在1-121-12的的1111小时内，每隔小时内，每隔1 1小时
10、测量一次小时测量一次温度，测得的温度依次为：温度，测得的温度依次为：5 5，8 8，9 9，1515，2525，2929，3131，3030，2222，2525，2727，2424。试估计每隔。试估计每隔1/101/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),y
11、label(Degrees Celsius)例2 在用外接电源给电容器充电时，电容器两端的电压V将会随着充电时间t发生变化，已知在某一次实验时，通过测量得到下列观测值，分别用拉格朗日插值法、分段线性插值法、三次样条插值法画出V随着时间t变化的曲线图，分别计算当时间t=7s时，三种插值法各自算得电容器两端电压的近似值。1295.64321st 4.101.106.90.92.83.72.6vV解 由于MATLAB没有提供现成的拉格朗日插值命令，我们可以编写一个函数lglrcz.m来完成，其他两种插值法可用现成的命令。用MATLAB软件进行三种插值计算的程序为szczqx.m。程序程序lglrcz
12、.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*(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(
13、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,fdxx,g,t0,scyt,b)gtext(lglr)gtext(fdxx)gtext(scyt)执行结果是laglr=9.52988980716254fendxx=9.70000000000000sancyt=9.67118039327444图形如图1所示。图中曲线lglr表示拉格朗日插值曲线，scyt表示三次样条插值曲线，fdxx表示分段线性插值曲线。案例:水流量的估计
14、美国某州的用水管理机构要求各社区提供以每小时多少加仑计的用水量以及每天所用的总水量。许多社区没有测量流入或流出水塔水量的装置，只能代之以每小时测量水塔中的水位，其误差不超过5%。需要注意的是，当水塔中的水位下降到最低水位L时，水泵就自动向水塔输水直到最高水位H，此期间不能测量水泵的供水量，因此，当水泵正在输水时不容易建立水塔中水位和用水量之间的关系。水泵每天输水一次或两次，每次约2小时。2.1 问题问题 试估计任何时刻（包括水泵正在输水时间）从水塔流出的水流量f(t)，并估计一天的总用水量。已知该水塔是一个高为40ft（英尺），直径为57ft（英尺）的正圆柱，表5-1给出了某个小镇一天水塔水位
15、的真实数据，水位降至约27.00ft水泵开始工作，水位升到35.50ft停止工作。（注：1ft（英尺）=0.3048m（米）表1 某小镇某天水塔水位时间/s水位/0.01ft时间/s水位/0.01ft03175466363350331631104995332606635305453936316710619299457254308713937294760574301217921289264554292721240285068535284225223279571854276728543275275021269732284269779254水泵开动35932水泵开动82649水泵开动39332水泵开
16、动8596834753943535508995333974331834459327033402.2 问题分析 本实验所指流量可视为单位时间内流出水的体积。由于水塔是正圆柱形，横截面积是常数，所以在水泵不工作时段，流量很容易根据水位相对时间的变化算出。问题的难点在于如何估计水泵供水时段的流量。水泵供水时段的流量只能靠供水时段前后的流量经插值或拟合得到。作为用于插值或拟合的原始数据，我们希望水泵不工作时段的流量越准确越好。这此流量大体上可由两种方法计算，一是直接对表12-1中的水量用数值微分算出各时段的流量，用它们拟合其它时刻或连续时间的流量；二是先用表中数据拟合水位一时间函数，求导数即可得到连续
17、时间的流量。有了任何时刻的流量，就不难计算一天的总用水量。其实，水泵不工作时段的用水量可以由测量记录直接得到，由表1中下降水位乘以水搭的截面积就是这一时段的用水量。这个数值可以用来检验数据插值或拟合的结果。在具体给出本问题的解答之前，先介绍一个简单的数据插值方法。2.2 问题求解问题求解 为了表示方便，我们将2.1节问题中所给表1中的数据全部化为国际标准单位（表2），时间用小时（h），高度用米（m）：表12-2 一天内水塔水位记录时间/h水位/m时间/h水位/m09.6812.9510.210.929.4513.889.941.849.3114.989.652.959.1315.909.413
18、.878.9816.839.184.988.8117.938.925.908.6919.048.667.008.5219.968.437.938.3920.848.228.978.2222.02水泵开动9.98水泵开动22.96水泵开动10.93水泵开动23.8810.5910.9510.8224.9910.3512.0310.5025.9110.181模型假设（1）流量只取决于水位差，与水位本身无关，故由物理学中Torriceli定律：小孔流出的液体的流速正比于水面高度的平方根。题目给出水塔的最低和最高水位分别是8.1648m)3024.027(和10.7352m)3024.050.35(（
19、设出口的水位为零），因为sqrt1467.1)1648.87352.10(，约为1，所以可忽略水位对流速的影响。（2）将流量看作是时间的连续函数，为计算简单，不妨将流量定义成单位时间流出水的高度，即水位对时间变化率的绝对值（水位是下降的），水塔截面积为3475.2334*)3024.0*57(2S（m2），得到结果后乘以s即可。2流量估计方法 首先依照表2所给数据，用MATLAB作出时间水位散点图（图2）。下面来计算水箱流量与时间的关系。根据图1.，一种简单的处理方法为，将表2中的数据分为三段，然后对每一段的数据做如下处理：设某段数据),(,),(),(1100nnyxyxyx，相邻数据中点的
20、平均流速用下面的公式（流速=（右端点的水位右端点的水位）/区间长度）：iiiiiixxyyxxv111)2(每段数据首尾点的流速用下面的公式计算：)()43()(022100 xxyyyxv)()43()(221nnnnnnxxyyyxv用以上公式求得时间与流速之间的数据如表3。时间/h流速/cmh-1时间/h流速/cmh-1029.8912.4931.520.4621.7413.4229.031.3818.4814.4326.362.39516.2215.4426.093.4116.3016.3724.734.42515.3217.3823.645.4413.0418.4923.426.45
21、15.4519.5025.007.46513.9820.4023.868.4516.3520.8422.178.9719.2922.02水泵开动9.98水泵开动22.96水泵开动10.93水泵开动23.8827.0910.9533.5024.4321.6211.4929.6325.4518.4825.9113.30由表3作出时间流速散点图如图3。1）插值法 由表12-3，对水泵不工作时段1，2采取插值方法，可以得到任意时刻的流速，从而可以知道任意时刻的流量。我们分别采取拉格朗日插值法，分段线性插值法及三次样条插值法；对于水泵工作时段1应用前后时期的流速进行插值，由于最后一段水泵不工作时段数据太
22、少，我们将它与水泵工作时段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；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)；/*注：注：l
23、glrcz为一函数，程序同为一函数，程序同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)gtext(fdxx)gtext(scyt)运行结果lglrjf=145.6231；fdxxjf=147.1430；sancytjf=1412.6870 图4是对第1未供水段数据用三种不同方法得到的插
24、值函数图，图中曲线lglr、fdxx和scyt分别表示用拉格朗日插值法，分段线性插值法及三次样条插值法得到的曲线。由表2知，第1未供水时段的总用水高度为146(=968822)，可见上述三种插值方法计算的结果与实际值（146）相比都比较接近。考虑到三次样条插值方法具有更加良好的性质，建议采取该方法。其他三段的处理方法与第1未供水时段的处理方法类似，这里不再详细叙述，只给出数值结果和函数图像（图5图7），图中曲线标记同图4。图5 第一供水段时间流速示意图图6 第2未供水段时间流速示意图图7 混合时段时间流速示意图图8是用分段线性及三次样条插值方法得到的整个过程的时间流速函数示意图。表4 各时段及
25、一天的总用水量（用水高度）第1未供水段第2未供水段第3供水段混合时段全天拉格朗日插值法145.6231258.866454.068992.1337550.6921分段线性插值法147.1430258.969749.605176.4688532.1866三次样条插值法145.6870258.654753.333481.7699539.4450 表5是对一天中任取的4个时刻分别用3种方法得到的水塔水流量近似值。时间6.8810.8815.8822.8815.9826671234851433.7426009085346325.5662241818047734.7099621055169414.827
26、2413793103432.9976262626262725.4465591397849525.471578947368415.0527858196582033.7089553595325925.5490892055742329.41733175863551注：拉格朗日插值法分段线性插值法三次样条插值法2）拟合法（1）拟合水位时间函数 从表12-2中的测量记录看，一天有两次供水时段和三次未供水时段，分别对第1，2未供水时段的测量数据直接作多项式拟合，可得到水位函数（注意，根据多项式拟合的特点，此处拟合多项式的次数不宜过高，一般以36次为宜）。对第3未供水时段来说，数据过少不能得到很好的拟合。设
27、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,10.21,9.94,9.65,9.41,9.18,8.92,8
28、.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=-polyval(c2,tp2)；plot(tp2,x2)（2）确定流量时间函数对
29、于第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=-polyval(a1,7.93,8.97)；*/为下面
30、的程序准备数据为下面的程序准备数据/*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单位保持一致单位保持一致/*subpl
31、ot(1,2,2)plot(tp2,x23*100)*/与图与图12.10单位保持一致单位保持一致/*程序运行得到第1，2未供水时段的时间流量图如图10，可以看到与图8中用插值给出的曲线比较吻合。如果用5次多项式拟合则得图11的图形，显然较三次拟合的效果好。而第1供水时段的流量则用前后时期的流量进行拟合得到。为使流量函数在t=9和t=11连续，我们只取4个点，用三次多项式拟合得到第1供水时段的时间流量图形如图12，可以看到与图8中的相应部分比较吻合。图12.12图12.8dygsdsy=7.93,8.97,10.95,12.03；dygsdls=x14,x24；nhjg=polyfit(dyg
32、sdsj,dygsdls,3)；nhsj=7.93:0.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，
33、可以看到与图8中的相应部分也比较吻合。图12.13，图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）一天总用水量的估计 分别对供水的两个时段
34、和不供水的两个时段积分（流量对时间）并求和得到一天的总用水量约为526.8935（此数据是总用水高度，单位为cm）。表12-6列出各段用水量，与插值法算得的表12-4相比，二者较为吻合。表12-6时段第1未供水段第2未供水段第1供水段第2供水段全天用水用水高度145.67260.6646.6073.9635526.8935（4）流量及总用水量的检验 计算出各时刻的流量可用水位记录的数值微分来检验，各时段的用水高度可以用实际记录的水位下降高度来检验。例如，算得第1未供水段的用水量高度是145.67，而实际记录的水位下降高度为968822=146cm，两者是吻合的；算得第2未供水段的用水量高度是2
35、60.66cm，而实际记录的水位下降高度为1082822=260cm，两者也是吻合的。从算法设计和分析可知，计算结果与各时段所用的拟合多项式的次数有关。表12-7给出的是对第1，2未供水时段分别用五、六多项式拟合后得到的用水量结果。表12-7时段第1未供水段第2未供水段第1供水段第2供水段全天用水用水高度 146.5150257.760546.131776.3076526.71483）结果分析由表12-4可以看出，使用三次样条插值法得到的结果)6547.258,6870.145(260cm相差不大，说明插值结果与原始数据比较吻合。与表12-2中记录的下降高度 146cm，按三次样条插值法估计出
36、全天的用水量约为 42.1258781103475.2334450.539由表12-7可以得全天的用水量约为 82.1229075103475.2337148.5262.5 内容小结内容小结 本实验主要进行水塔水流量的估算。第一种估算方法为插值方法，我们用了三种不同的插值法进行估计，在求解的过程中，可以熟悉数据插值的理论和方法；第二种估算方法为数据拟合法，用多项式进行拟合，得到水塔水流量的估计。2二维插值命令二维插值命令interp2的具体使用格式的具体使用格式zz=interp2(x,y,z,xx,yy,method)该指令的意思是根据数据向量x，y，z按method指定的方法来做插值，然后
37、将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)它的具体含义跟前面的一、二维插值是相似的，在此不作解释，读者可在MATLAB工作空间中用help interp3命令获得。4样条插值命令样条插值命
38、令spline的具体使用格式的具体使用格式yy=spline(x,y,xx)它的意思等同于命令yy=interp1(x,y,xx,cubic)xy机翼下轮廓线X035791 11 21 31 41 5Y01.21.72.02.12.01.81.21.01.6例例 已知飞机下轮廓线上数据如下，求已知飞机下轮廓线上数据如下，求x每改变每改变0.1时的时的y值。值。To MATLAB(plane)返回返回二维插值的定义二维插值的定义 xyO O第一种（网格节点）：第一种（网格节点）：已知已知 m n个节点个节点),2,1;,.,2,1(),(njmizyxijji 其中其中jiyx,互不相同，不妨设
39、互不相同，不妨设bxxxam 21dyyycn 21 构造一个二元函数构造一个二元函数),(yxfz 通过全部已知节点通过全部已知节点,即即再用再用),(yxf计算插值，即计算插值，即).,(*yxfz ),1,0;,1,0(),(njmizyxfijji 第二种（散乱节点）：第二种（散乱节点）：yx0 0已知已知n个节点个节点),.,2,1(),(nizyxiii 其中其中),(iiyx互不相同，互不相同，构造一个二元函数构造一个二元函数),(yxfz 通过全部已知节点通过全部已知节点,即即),1,0(),(nizyxfiii 再用再用),(yxf计算插值，即计算插值，即).,(*yxfz
40、返回返回 注意：注意：最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。最邻近插值最邻近插值x y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O O 二维或高维情形的最邻近插值，与被插值点最邻近的节点的函数值即为所求。返回返回 将四个插值点（矩形的四个顶点）处的函数值依次简记为：分片线性插值分片线性插值xy (xi,yj)(xi,yj+1)(xi+1,yj)(xi+1,yj+1)O Of(xi,yj)=f1，f(xi+1,yj)=f2，f(xi+1,yj+1)=f3，f(xi,yj+1)=f4插值函数为：jii1ij1jy)xx(xxyyy)yy)(ff()xx)(ff
41、(f)y,x(fj23i121第二片(上三角形区域)：(x,y)满足iii1ij1jy)xx(xxyyy插值函数为：)xx)(ff()yy)(ff(f)y,x(fi43j141注意注意：(x,y)当然应该是在插值节点所形成的矩形区域内。显然，分片线性插值函数是连续的；分两片的函数表达式如下：第一片(下三角形区域)：(x,y)满足返回返回 双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下：)dcy)(bax()y,x(f其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。双线性插值双线性插值x y(x1,y1)(x1,y2)(x
42、2,y1)(x2,y2)O O返回返回 要求要求x0,y0 x0,y0单调；单调；x x，y y可取可取为矩阵，或为矩阵，或x x取取行向量，行向量，y y取为列向量，取为列向量，x,yx,y的值分别不能超出的值分别不能超出x0,y0 x0,y0的范围。的范围。z=interp2(x0,y0,z0,x,y,method)被插值点插值方法用用MATLAB作网格节点数据的插值作网格节点数据的插值插值节点被插值点的函数值nearest nearest 最邻近插值最邻近插值linear linear 双线性插值双线性插值cubic cubic 双三次插值双三次插值缺省时缺省时,双线性插值双线性插值例：
43、测得平板表面例：测得平板表面3 3*5 5网格点处的温度分别为：网格点处的温度分别为：82 81 80 82 84 82 81 80 82 84 79 63 61 65 81 79 63 61 65 81 84 8484 84 82 85 86 82 85 86 试作出平板表面的温度分布曲面试作出平板表面的温度分布曲面z=f(x,y)z=f(x,y)的图形。的图形。输入以下命令：x=1:5;y=1:3;temps=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86;mesh(x,y,temps)1.先在三维坐标画出原始数据，画出粗糙的温度分布曲图.2以平滑
44、数据,在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)例例 山区地貌：山区地貌：在某山区测得一些地点的高程如下表。平面区域为在某山区测得一些地点的高程如下表。平面区域为 1200=x=4000,1200=y x=1:5;n y=1:5;n t=100,100,100,100,100;105,120,122,125,122;n 110,130,155,157,130;115,133,1
45、57,160,140;n 113,132,149,154,128nt=n 100 100 100 100 100n 105 120 122 125 122n 110 130 155 157 130n 115 133 157 160 140n 113 132 149 154 128n mesh(x,y,t)数据原图数据原图n xx=1:0.1:5;n yy=1:0.1:5;n tt=interp2(x,y,t,xx,yy,cubic);n mesh(xx,yy,tt)二维插值得到的示意图二维插值得到的示意图 例例 在某海域测得一些点在某海域测得一些点(x,y)(x,y)处的水深处的水深z z由下
46、由下表给出，船的吃水深度为表给出，船的吃水深度为5 5英尺，在矩形区域（英尺，在矩形区域（7575，200200）*（-50-50，150150）里的哪些地方船要避免进入。）里的哪些地方船要避免进入。xyz129 140 103.5 88 185.5 195 1057.5 141.5 23 147 22.5 137.5 85.54 8 6 8 6 8 8xyz157.5 107.5 77 81 162 162 117.5-6.5 -81 3 56.5 -66.5 84 -33.59 9 8 8 9 4 9 )1(.150,50200,75.2hd三次插值法作二维插值在矩形区域.3 作海底曲面图
47、.1 输入插值基点数据To MATLAB hd1返回返回4.作出水深小于5的海域范围,即z=5的等高线.实验作业实验作业 山区地貌：山区地貌：在某山区测得一些地点的高程如下表：在某山区测得一些地点的高程如下表：(平平面区域面区域1200=x=4000,1200=y=3600)1200=x=4000,1200=y=3600)，试作出该山区的，试作出该山区的地貌图和等高线图，并对几种插值方法进行比较。地貌图和等高线图，并对几种插值方法进行比较。36003200280024002000160012001480 1500 1550 1510 1430 1300 1200 9801500 1550 16
48、00 1550 1600 1600 1600 15501500 1200 1100 1550 1600 1550 1380 10701500 1200 1100 1350 1450 1200 1150 10101390 1500 1500 1400 900 1100 1060 9501320 1450 1420 1400 1300 700 900 8501130 1250 1280 1230 1040 900 500 700Y/x1200 1600 2000 2400 2800 3200 3600 4000返回返回二二.拟合拟合主要内容主要内容2、掌握用数学软件求解拟合问题。、掌握用数学软件求
49、解拟合问题。1、直观了解拟合基本内容。、直观了解拟合基本内容。1 1、拟合问题引例及基本理论。拟合问题引例及基本理论。4 4、实验作业。实验作业。2、用数学软件求解拟合问题。用数学软件求解拟合问题。3、应用实例应用实例拟拟 合合2.2.拟合的基本原理拟合的基本原理1.拟合问题引例拟合问题引例拟拟 合合 问问 题题 引引 例例 1 1温度温度t(0C)20.5 32.7 51.0 73.0 95.7电阻电阻R()765 826 873 942 1032已知热敏电阻数据：已知热敏电阻数据：求求60600C时的电阻时的电阻R。2040608010070080090010001100 设设 R=at+
50、ba,b为待定系数为待定系数拟拟 合合 问问 题题 引引 例例 2 2 t(h)0.25 0.5 1 1.5 2 3 4 6 8c(g/ml)19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01已知一室模型快速静脉注射下的血药浓度数据已知一室模型快速静脉注射下的血药浓度数据(t=0注射注射300mg)求血药浓度随时间的变化规律求血药浓度随时间的变化规律c(t).作半对数坐标系作半对数坐标系(semilogy)下的图形下的图形为待定系数kcectckt,)(002468100101102MATLAB(aa1)曲曲 线线 拟拟 合合 问问 题题 的的
