1、 例：在例：在1-12的的11小时内，每隔小时内，每隔1小时测量一次温度，小时测量一次温度，测得的温度依次为：测得的温度依次为：5，8，9，15，25，29，31，30，22，25，27，24。试估计每隔。试估计每隔1/10小时的温度值。小时的温度值。0246810125101520253035HourDegrees Celsius第1页/共48页拉格朗日插值分段线性插值三次样条插值一、一、插值的定义插值的定义二、插值的方法二、插值的方法三、用三、用Matlab解插值问题解插值问题一维插值第2页/共48页一维插值的定义一维插值的定义已知已知 n+1个节点个节点,1,0(),(njyxjj其中其
2、中jx互不相同，不妨设互不相同，不妨设),10bxxxan求任一插值点求任一插值点)(*jxx 处的插值处的插值.*y0 x1xnx0y1y*x*y第3页/共48页 构造一个构造一个(相对简单的相对简单的)函数函数),(xfy 通过全部节点通过全部节点,即即),1,0()(njyxfjj再用再用)(xf计算插值，即计算插值，即).(*xfy 0 x1xnx0y1y*x*y第4页/共48页011110.,()().()(0,1.)nnnnnnnnjjaxxxb LxnLxa xaxa xaLxyjn设是 次多项式满足YXA次插值多项式个结点确定1式有唯一解，即0德蒙行列式，经列交换和转置后是范n
3、nXX第5页/共48页 称为拉格朗日插值基函数拉格朗日插值基函数。0()()nniiiP xL xy 已知函数f(x)在n+1个点x0,x1,xn处的函数值为 y0,y1,yn。求一n次多项式函数Pn(x)，使其满足：Pn(xi)=yi,i=0,1,n.解决此问题的拉格朗日插值多项式公式如下:其中Li(x)为n次多项式：01110111()()()()()()()()()()()iiniiiiiiiinxxxxxxxxxxL xxxxxxxxxxx拉格朗日拉格朗日(Lagrange)插值插值ijijxLji，当，当01)(第6页/共48页拉格朗日(Lagrange)插值特别地特别地:两点一次两
4、点一次(线性线性)插值多项式插值多项式:101001011yxxxxyxxxxxL三点二次三点二次(抛物抛物)插值多项式插值多项式:1202012012010210122021x xx xx xx xx xx xL xyyyxxxxxxxxxxxx第7页/共48页For examplex149161234取最接近x=5的点，x0=1，x1=4，x2=9为插值节点，运用插值公式L(5)=2.27.0 01 12 2L(x)=y l(x)+y()()l xy lx第8页/共48页For examplex757677782.76 2.83 2.90 2.97取最接近x=5的点，x0=1，x1=4，x
5、2=9为插值节点，运用插值公式L(5)=2.27.0 01 12 2L(x)=y l(x)+y()()l xy lx第9页/共48页21(),551g xxx 采用拉格朗日多项式插值：选取不同插值节点个数n+1，其中n为插值多项式的次数，当n分别取2,4,6,8,10时，绘出插值结果图形。例例第10页/共48页第11页/共48页分段线性插值分段线性插值其它,0,)()()(1111110jjjjjjjjjjjnjjjnxxxxxxxxxxxxxxxlxlyxL计算量与n无关;n越大，误差越小.nnnxxxxgxL0),()(limxjxj-1xj+1x0 xnxoy第12页/共48页xy机翼下
6、轮廓线X035791 11 21 31 41 5Y01.21.72.02.12.01.81.21.01.6例例 已知飞机下轮廓线上数据如下，求已知飞机下轮廓线上数据如下，求x每改变每改变0.1时的时的y值。值。第13页/共48页比分段线性插值更光滑。比分段线性插值更光滑。xyxi-1 xiab 在数学上，光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续，则称该曲线具有k阶光滑性。光滑性的阶次越高，则越光滑。是否存在较低次的分段多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例子。三次样条插值三次样条插值第14页/共48页 三次样条插值三次样条插值;.2,1,0,)()1(njjj
7、yxS上是三次多项式在每个小区间1.,2,1,01,)()2(njjjxxxS上有连续的二阶导数在,)()3(baxS第15页/共48页njjjyxS.,1,0,)(插值条件32)(xDxCxBAxSjjjjj1,.,2,1)()()3()()()2()()()1(njxSxSxSxSxSxSjjjjjj连接条件第16页/共48页00000(1)()(),()()(2)()(),()()()()0nnnnnmS xfxS xfxMSxfxSxfxSxSx边值条件：边值条件：边值条件：特别地：称为自然边界条件)(xS确定件、边界条件可以唯一利用插值条件、连接条第17页/共48页用用MATLAB作
8、插值计算作插值计算一维插值函数：一维插值函数：yi=interp1(x，y，xi，method)插值方法插值方法被插值点被插值点插值节点插值节点xi处的插处的插值结果值结果nearest ：最邻近插值最邻近插值linear ：线性插值线性插值；spline ：三次样条插值三次样条插值；cubic ：立方插值。立方插值。缺省时：缺省时：分段线性插值分段线性插值。注意注意：所有的插值方法都要求：所有的插值方法都要求x是单调的，并且是单调的，并且xi不不能够超过能够超过x的范围。的范围。第18页/共48页2.拟合的基本原理拟合的基本原理1.拟合问题引例拟合问题引例拟 合第19页/共48页拟拟 合合
9、问问 题题 引引 例例 1温度温度t(0C)20.5 32.7 51.0 73.0 95.7电阻电阻R()765 826 873 942 1032已知热敏电阻数据：已知热敏电阻数据：求求600C时的电阻时的电阻R。2040608010070080090010001100 设设 R=at+ba,b为待定系数为待定系数第20页/共48页 函数插值与曲线拟合都是要根据一组数据构造一个函数作函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同为近似，由于近似的要求不同，二者的数学方法上是完全不同的。的。实例：实例：下面数据是某次实验所得，希望得到X和
10、 f之间的关系？x124791 21 31 51 7f1.53.96.611.71 5.61 8.81 9.62 0.62 1.1问题问题：给定一批数据点，需确定满足特定要求的曲线或曲面解决方案解决方案：若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，这就是数据拟合数据拟合，又称曲线拟合或曲面拟合。若要求所求曲线（面）通过所给所有数据点，就是插值问题插值问题；拟合与插值的关系第21页/共48页线性插值、样条插值与曲线拟合结果：线性插值、样条插值与曲线拟合结果：第22页/共48页曲线拟合问题最常用的解法曲线拟合问题最常用的解法线性最小二乘法的基本思路线性最小二乘法的基本思路第
11、一步:先选定一组函数先选定一组函数 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 的平方和最小的平方和最小。记记)2()()(),(211211221iiknimkkininiiimyxrayxfaaaJ 问题归结为，求问题归结为，求 a1,a2,am 使使 J(a1,a2,am)最小。最小。第23页/共48页线性最小二乘法的
12、求解：预备知识线性最小二乘法的求解：预备知识超定方程组超定方程组：方程个数大于未知量个数的方程组：方程个数大于未知量个数的方程组)(221111212111mnyarararyarararnmnmnnmm即即 Ra=ynmnmnnmyyyaaarrrrrrR112111211,其中其中超定方程一般是不存在解的矛盾方程组。超定方程一般是不存在解的矛盾方程组。如果有向量如果有向量a使得使得 达到最小，达到最小，则称则称a为上述为上述超定方程的最小二乘解超定方程的最小二乘解。212211)(imniimiiyararar第24页/共48页 定理定理：当当RTR可逆时，超定方程组（可逆时，超定方程组（
13、3）存在最小二乘）存在最小二乘解，且即为方程组解，且即为方程组 RTRa=RTy的解：的解：a=(RTR)-1RTy 所以，曲线拟合的最小二乘法要解决的问题，实际上就是所以，曲线拟合的最小二乘法要解决的问题，实际上就是求以下超定方程组的最小二乘解的问题。求以下超定方程组的最小二乘解的问题。nmnmnmyyyaaaxrxrxrxrR111111,)()()()(其中其中Ra=y （3）线性最小二乘法的求解第25页/共48页线性最小二乘拟合线性最小二乘拟合 中函数中函数r1(x),rm(x)的选取的选取 1.通过机理分析建立数学模型来确定通过机理分析建立数学模型来确定 f(x)；+f=a1+a2x
14、f=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx 2.将数据将数据(xi,yi)i=1,n 作图，通过直观判断确定作图，通过直观判断确定 f(x)：第26页/共48页1.1.线性最小二乘拟合线性最小二乘拟合2.2.非线性最小二乘拟合非线性最小二乘拟合二、用MATLAB解拟合问题第27页/共48页1)作多项式作多项式f(x)=a1xm+amx+am+1拟合拟合,可利用已有程序可利用已有程序:a=polyfit(x,y,m)2)对超定方程组对超定方程组)(11nmyaRnmmn可得最小二乘意义下的解。可得最小二乘意义下的解。，用，用yRa3)多项式在
15、多项式在x处的值处的值y可用以下命令计算：可用以下命令计算：y=polyval（a，x）输出拟合多项式系数输出拟合多项式系数a=a1,am,am+1(数组数组)）输入同长度输入同长度的数组的数组X，Y拟合多项拟合多项式次数式次数1.用MATLAB作线性最小二乘拟合第28页/共48页即要求出二次多项式即要求出二次多项式:3221)(axaxaxf中中 的的),(321aaaA 使得使得:最小 )(1112iiiyxf例例 对下面一组数据作二次多项式拟合对下面一组数据作二次多项式拟合第29页/共48页1 1）输入以下命令输入以下命令：x=0:0.1:1;y=-0.447 1.978 3.28 6.
16、16 7.08 7.34 7.66 9.56 9.48 9.30 11.2;R=(x.2)x ones(11,1)；A=Ry11 11211121xxxxR此时解法解法1用解超定方程的方法用解超定方程的方法2 2）计算结果计算结果：=-9.8108 20.1293 -0.03170317.01293.208108.9)(2xxxf第30页/共48页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;A=polyfit(x,y,2)z=polyval(A,x);plot(x,y
17、,k+,x,z,r)%作出数据点和拟合曲线的图形作出数据点和拟合曲线的图形2 2）计算结果计算结果：=-9.8108 20.1293 -0.0317解法解法2用多项式拟合的命令用多项式拟合的命令00.20.40.60.81-20246810120317.01293.208108.9)(2xxxf第31页/共48页建模案例：黄河小浪底调水调沙问题 2004年6月至7月黄河进行了第三次调水调沙试验，特别是首次由小浪底、三门峡和万家寨三大水库联合调度，采用接力式防洪预泄放水，形成人造洪峰进行调沙试验获得成功。整个试验期为20多天，小浪底从6月19日开始预泄放水，至到7月13日恢复正常供水结束。小浪底
18、水利工程按设计拦沙量为75.5亿立方米，在这之前，小浪底共积泥沙达14.15亿吨。第32页/共48页建模案例：黄河小浪底调水调沙问题 这次调水调沙试验一个重要的目的就是由小浪底上游的三门峡和万家寨水库泄洪，在小浪底形成人造洪峰，冲刷小浪底库区沉积的泥沙，在小浪底水库开闸泄洪以后，从6月27日开始三门峡水库和万家寨水库陆续开闸放水，人造洪峰于29日先后到达小浪底，7月3日达到最大流量2700，使小浪底水库的排沙量也不断地增加。表7-1是由小浪底观测站从6月29日到7月10日检测到的试验数据。第33页/共48页表 7-1 试验观测数据（单位：水流为sm/3，含沙量为3kg/m）日期 6.29 6.
19、30 7.1 7.2 7.3 7.4 时间 8：00 20:00 8：00 20:00 8：00 20:00 8：00 20:00 8：00 20:00 8：00 20:00 水流量 1800 1900 2100 2200 2300 2400 2500 2600 2650 2700 2720 2650 含沙量 32 60 75 85 90 98 100 102 108 112 115 116 日期 7.5 7.6 7.7 7.8 7.9 7.10 时间 8：00 20:00 8：00 20:00 8：00 20:00 8：00 20:00 8：00 20:00 8：00 20:00 水流量 2
20、600 2500 2300 2200 2000 1850 1820 1800 1750 1500 1000 900 含沙量 118 120 118 105 80 60 50 30 26 20 8 5 注注：以上数据主要是根据媒体公开报道的结果整理而成，不一定与真实数据完全相符。现在，根据试验数据建立数学模型研究下面的问题：（1）给出估算任意时刻的排沙量及总排沙量的方法；（2）确定排沙量与水流量的变化关系。第34页/共48页1.模型假设1)水流量和排沙量都是连续的，不考虑上游泄洪所带来的含沙量和外界带来的含沙量。2)时间是连续变化的，所取时间点依次为1,2,3，,24,单位时间为12h。第35页
21、/共48页假设水流量和排沙量都是连续的，不考虑上游泄洪所带 的含沙量和外界(如雨水等)带 入的含沙量，时间是连续变化的，已知给定的观测时间间距是等间距的，以 6 月29 日 8 点位第一个观测节点，所取 t 依次为 1，2，3，2 4，单位时间为 12h。2.模型建立与求解问题一：给出估算任意时刻的排沙量及总排沙量的方法 根据试验数据，要计算任意时刻的排沙量，就要确定出排沙量随时间变化的规律，可以通过插值来实现，考虑到实际中排沙量应该是时间的连续函数，为了提高精度，我们采用三次样条函数来进行插值。第36页/共48页记第i次观测时水流量为3(/)iv ms，含沙量为3(/)ih kg m，则第i
22、次观测时的排排沙沙量量为 1）建立模型iiiyvh(/)kg s 已知给定的观测时刻是等间距的，以6月29 日零时刻开始计时，则各次观测时刻分别为：124,1,2,24itii 计时单位为小时。第一次观测时刻是第8小时，最后一次观测时刻是第284 小时。第37页/共48页由于时间间隔较大（12小时），直接计算排沙量误差很大，将排沙量的单位由(/)kg s，换算为吨/12小时（(/12)th），即 iiiyvh(/)kg s12 3.6iivh (/12)th 将数据进行插值，用三次样条函数对排沙量iy进行插值，设第i个观测节点到第1i个节点的三次样条插值函数为：321234()()(),1,1
23、,2,23iiiiiya t ia t ia t iati ii 则总总排排沙沙量量模模型型为：2313212341()()()iiiiiiiya tia tia tiadx 第38页/共48页下面利用matlab 进行插值，用csape 求出每两相邻观测节点间的三次样条函数，再每段积分求和得排沙量，程序如下：模型求解clc,clear i=1:24;%观测节点 v=12*3.6.*1800 1900,2100,2200,2300,2400,2500,2600,2650,2700,2720,.2650,2600,2500,2300,2200,2000,1850,1820,1800,1750,1
24、500,1000,900;h=32,60,75,85,90,98,100,102,108,112,115,116,118,120,118,105,80,60,50,30,26,20,8,5;%含沙量 y=v.*h;%计算排沙量 pp=csape(i,y);%进行三次样条插值 xishu=pp.coefs%求三次样条插值多项式的系数矩阵 z=0;syms x real;for i=1:23%节点i 到i+1 段的定积分来计算这段的排沙量 z=z+double(int(xishu(i,1)*(x-i)3+xishu(i,2)*(x-i)2+xishu(i,3)*(x-i)+xishu(i,4),i
25、,i+1);z end 第39页/共48页z=1.8440e+08 即总排沙量1.844 亿吨，此与媒体报道的排沙量几乎一样。计算结果第40页/共48页从试验数据可以看出，在排水排沙的过程中，随着时间的推移，水库内所剩泥沙越来越少，开始排沙量是随着水流量的增加而增加，而随后是随着水流量的减少而减少，显然变化规律并非是线性关系。因此，将问题分为两个部分，从开始水流量增加到最大值 2720m3/s 为第一阶段，从水流量最大值到结束为第二阶段，用拟合的方法分别来研究水流量与排沙量的关系。问题二：确定排沙量与水流量的变化关系。表 7-3 第一阶段试验观测数据 节节 点点 1 2 3 4 5 6 7 8
26、 9 10 11 水流量 1800 1900 2100 2200 2300 2400 2500 2600 2650 2700 2720 含沙量 32 60 75 85 90 98 100 102 108 112 115 表 7-4 第二阶段试验观测数据 节节 点点 12 13 14 15 16 17 18 19 20 21 22 23 24 水流量 2650 2600 2500 2300 2200 2000 1850 1820 1800 1750 1500 1000 900 含沙量 116 118 120 118 105 80 60 50 30 26 20 8 5 第41页/共48页先画出排沙
27、量与水流量的散点图：散点图图 7-18 第一、第二阶段排沙量与水流量散点图 第42页/共48页从散点图可以看出排沙量与水流量的大致趋势，对两个阶段我们都用多项式进行拟合，并计算它的最小误差S：误差计算方法设m次多项式拟合函数为 1miiiya x m次多项式拟合的最小误差为 1111()miijjiiSa xy （第一阶段）1311()miijjiiSa xy （第二阶段）最终采用的拟合多项式的次数m由最小误差S来确定，最小误差S越小，拟合效果越好。第43页/共48页从散点图看曲线接近于直线，故可利用 matlab 分别做 1 次、2 次、3 次、4 次多项式拟合，第一阶段多项式拟合结果如下：
28、第一阶段多项式拟合 第44页/共48页其最小误差为：S1=9.206181083635087e+008 S2=6.895433807257007e+008 S3=3.064729231715386e+008 S4=3.042796821308041e+008 从上图及最小误差来看，3 次多项式拟合与 4 次多项式效果差别不大，基本可以看出排沙量与水流量之间的关系。其 3 次多项式次多项式拟合与 4 次多项式次多项式分别为：43233.4683 102.41885814.84414590454.2078yvvv7433241.1 101.348 105.80910880.71887408344.7yvvvv 第45页/共48页从散点图看曲线接近于 3 次或 4 成多项式，故可利用 matlab只做 3 次、4 次多项式拟合，第二阶段多项式拟合结果 第二阶段多项式拟合 第46页/共48页其最小误差为：S3=3.612084193389563e+009 S4=1.748805174191605e+009 从拟合效果及最小误差来看，显然 4 次多项式拟合要好于 3次多项式拟合，4 次多项式次多项式分别为：7433242.8 101.811 104.0923891.04411322627.4967yvvvv 第47页/共48页谢谢您的观看！第48页/共48页
