数学建模——拟合与插值

1、1 3/5/2022 数学建模教程数学建模教程拟拟 合与合与 插插 值值 2 3/5/2022 q在大量的应用领域中，人们经常面临这样的问题：在大量的应用领域中，人们经常面临这样的问题：给定一批数据点，需确定满足特定要求的曲线或给定一批数据点，需确定满足特定要求的曲线或曲面。对这个问题有两种方法。曲面。对这个问题有两种方法。q 一种是一种是插值法插值法，数据假定是正确的，要求以某种方法描述数，数据假定是正确的，要求以某种方法描述数据点之间所发生的情况。据点之间所发生的情况。q 另一种方法是另一种方法是曲线拟合或回归曲线拟合或回归。人们设法找出某条光滑曲线，。人们设法找出某条光滑曲线，它最佳地拟

2、合数据，但不必要经过任何数据点。它最佳地拟合数据，但不必要经过任何数据点。 q 本专题的主要目的是：了解插值和拟合的基本内容；本专题的主要目的是：了解插值和拟合的基本内容；掌握用掌握用MatlabMatlab求解插值与拟合问题的基本命令。求解插值与拟合问题的基本命令。 函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。由于近似的要求不同，二者的数学方法上是完全不同的。3 3/5/2022 内容提纲内容提纲1.拟合问题引例及基本理论2.Matlab求解拟合问题3.应用实例4.插值问

3、题引例及基本理论5.Maltab求解插值问题6.应用实例4 3/5/2022 拟拟 合合 问问 题题 引引 例例 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+ba,b为待定系数为待定系数一、拟合问题一、拟合问题5 3/5/2022 拟拟 合合 问问 题题 引引 例例 2 2 t (h) 0.25 0.5 1 1.5 2 3 4 6 8c ( g/ml) 19.21

4、 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01已知一室模型快速静脉注射下的血药浓度数据已知一室模型快速静脉注射下的血药浓度数据(t=0注射注射300mg)求血药浓度随时间的变化规律求血药浓度随时间的变化规律c(t).作半对数坐标系作半对数坐标系(semilogy)下的图形下的图形00( ),ktc tc eck为待定系数024681001011026 3/5/2022 曲曲 线线 拟拟 合合 问问 题题 的的 提提 法法已知一组（二维）数据，即平面上已知一组（二维）数据，即平面上 n个点个点（xi,yi) i=1,n, 寻求一个函数（曲线）寻求一个函数

5、（曲线）y=f(x), 使使 f(x) 在某种准则下与所有在某种准则下与所有数据点最为接近，即曲线拟合得最好。数据点最为接近，即曲线拟合得最好。 +xyy=f(x)(xi,yi)i i 为点为点（xi,yi) 与与曲线曲线 y=f(x) 的距离的距离7 3/5/2022 线性最小二乘拟合线性最小二乘拟合 f(x)=a1r1(x)+ +amrm(x)中中函数函数rr1 1(x), r(x), rm m(x)(x)的选取的选取 1. 1. 通过机理分析建立数学模型来确定通过机理分析建立数学模型来确定 f(x)f(x)；+f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+

6、a2/xf=aebxf=ae-bx2. 2. 将数据将数据 (xi,yi) i=1, n 作图，通过直观判断确定作图，通过直观判断确定 f(x)：8 3/5/2022 曲线拟合问题最常用的解法曲线拟合问题最常用的解法线性最小二乘法的基本思路线性最小二乘法的基本思路第一步: :先选定一组函数先选定一组函数 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)

7、 与与曲线曲线 y=f(x) 的距离的距离 i 的平方和最小的平方和最小 。记记 )2()()(),(211211221iiknimkkininiiimyxrayxfaaaJ 问题归结为，求问题归结为，求 a1,a2, am 使使 J(a1,a2, am) 最小。最小。9 3/5/2022 线性最小二乘法的求解：预备知识线性最小二乘法的求解：预备知识超定方程组超定方程组：方程个数大于未知量个数的方程组：方程个数大于未知量个数的方程组)( 221111212111mnyarararyarararnmnmnnmm即即 Ra=ynmnmnnmyyyaaarrrrrrR112111211,其中其中超定

8、方程一般是不存在解的矛盾方程组。超定方程一般是不存在解的矛盾方程组。 如果有向量如果有向量a使得使得 达到最小，达到最小，则称则称a为上述为上述超定方程的最小二乘解超定方程的最小二乘解。 212211)(imniimiiyararar10 3/5/2022 线性最小二乘法的求解线性最小二乘法的求解定理：定理：当当R RT TR R可逆时，超定方程组（可逆时，超定方程组（3 3）存在最小二乘解，）存在最小二乘解， 且即为方程组且即为方程组 R RT TRa=RRa=RT Ty -y -正则（正规）方程组正则（正规）方程组的解：的解：a=(Ra=(RT TR)R)-1-1R RT Ty y所以，曲

9、线拟合的最小二乘法要解决的问题，实际上就是求以所以，曲线拟合的最小二乘法要解决的问题，实际上就是求以下超定方程组的最小二乘解的问题。下超定方程组的最小二乘解的问题。nmnmnmyyyaaaxrxrxrxrR111111,)()()()(其中其中Ra=y （3）11 3/5/2022 用用MATLAB解拟合问题解拟合问题1 1、线性最小二乘拟合、线性最小二乘拟合2 2、非线性最小二乘拟合、非线性最小二乘拟合12 3/5/2022 用用MATLAB作线性最小二乘拟合作线性最小二乘拟合1. 1. 作多项式作多项式f(x)=a1xm+ +amx+am+1拟合拟合, ,可利用已有命令可利用已有命令:a=

10、polyfit(x,y,m)3.3.对超定方程组对超定方程组)(11nmyaRnmmn可得最小二乘意义下的解。可得最小二乘意义下的解。，用，用yRa2.2.多项式在多项式在x x处的值处的值y y的计算命令：的计算命令：y=polyvaly=polyval（a a，x x）输出拟合多项式系数输出拟合多项式系数a=a1,am,am+1 （数组）数组）输入同长度输入同长度数组数组X，Y拟合多项式拟合多项式次数次数13 3/5/2022 即要求即要求 出二次多项式出二次多项式:3221)(axaxaxf中中 的的123( ,)Aa a a 使得使得:1121() iiif xy最小例例 对下面一组数

11、据作二次多项式拟合对下面一组数据作二次多项式拟合xi 0.1 0.2 0.4 0.5 0.6 0.7 0.8 0.9 1 yi 1.978 3.28 6.16 7.34 7.66 9.58 9.48 9.30 11.2 211112222232111111,1,1,1xxyayxxaayxx14 3/5/2022 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=Ry11 11211121xxxxR此时MATLAB(zxec1

12、)解法解法1 1解超定方程的方法解超定方程的方法2）计算结果）计算结果： = -9.8108, 20.1293, -0.03170317.01293.208108.9)(2xxxfRAy15 3/5/2022 16 3/5/2022 2）计算结果：）计算结果： = -9.8108, 20.1293, -0.0317解法解法2用多项式拟合的命令用多项式拟合的命令MATLAB(zxec2)00.20.40.60.81-20246810120317.01293.208108.9)(2xxxf1）输入命令：）输入命令：x=0:0.1:1; y=-0.447,1.978,3.28,6.16,7.08,7

13、.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) %作出数据点和拟合曲线的图形作出数据点和拟合曲线的图形17 3/5/2022 1. 1. lsqcurvefitlsqcurvefit已知数据点数据点：xdata=xdata=（xdata1，xdata2，xdataxdatan n） ydata=ydata=（ydataydata1 1，ydataydata2 2，ydataydatan n）用用MATLAB作非线性最小二乘拟合作非线性最小二乘拟合两个求非线性最小二乘拟合的函数：两个求非

14、线性最小二乘拟合的函数：lsqcurvefitlsqcurvefit、lsqnonlinlsqnonlin。相同点和不同点：两个命令都要先建立相同点和不同点：两个命令都要先建立M-M-文件文件fun.mfun.m，定义函，定义函数数f(x)f(x)，但定义，但定义f(x)f(x)的方式不同的方式不同。211( ,) 2niiiF x xdataydata最小 lsqcurvefitlsqcurvefit用以求含参量用以求含参量x x（向量）的向量值函数（向量）的向量值函数F(x,xdata)=F(x,xdata)=（F F（x x，xdataxdata1 1），），F F（x x，xdatax

15、datan n）T T使得使得 18 3/5/2022 输入格式输入格式: : (1) x = lsqcurvefit (fun,x0,xdata,ydata); (2) x =lsqcurvefit (fun,x0,xdata,ydata,lb, ub); (3) x =lsqcurvefit (fun,x0,xdata,ydata, lb, ub, options); (4) x, resnorm = lsqcurvefit (fun,x0,xdata,ydata,); (5) x, resnorm, residual = lsqcurvefit (fun,x0,xdata,ydata,)

16、; fun是一个事先建立的是一个事先建立的定义函数定义函数F(x,xdata) 的的M-文件文件, 自变量为自变量为x和和xdata说明：说明：x = lsqcurvefit (fun,x0,xdata,ydata,options);迭代初值迭代初值已知数据点已知数据点选项见无选项见无约束优化约束优化19 3/5/2022 lsqnonlin用以求含参量用以求含参量x x（向量）的向量值函数（向量）的向量值函数 f(x)f(x)=(f=(f1 1(x),f(x),f2 2(x),f(x),fn n(x)(x)T T ，使得，使得 最小。最小。 其中其中 fi（x）=f（x，xdatai，yda

17、tai） =F(x,xdatai)-ydatai22212( ) ( )( )( )( )Tnfx f xf xfxfx2. lsqnonlin已知数据点：已知数据点： xdata=xdata=（xdata1，xdata2，xdataxdatan n） ydata=ydata=（ydataydata1 1，ydataydata2 2，ydataydatan n）20 3/5/2022 输入格式：输入格式： 1） x=lsqnonlin（fun，x0）； 2） x= lsqnonlin （fun，x0，lb，ub）； 3） x= lsqnonlin （fun，x0， ，lb，ub，options

18、）； 4） x， resnorm= lsqnonlin （fun，x0，）； 5） x， resnorm ， residual= lsqnonlin （fun，x0，）；说明：说明：x= lsqnonlinlsqnonlin （fun，x0，options）；）；fun是一个事先建立的是一个事先建立的定义函数定义函数f(x)的的M-文件，文件，自变量为自变量为x迭代初值迭代初值选项见无选项见无约束优化约束优化21 3/5/2022 100200 30040050060070080090010004.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59j

19、t310jc100.022111min( , , )22jktjjF a b kabec例例2 用下面一组数据拟合用下面一组数据拟合 中的参数中的参数a，b，k0.0.2( )ktc tabe 该问题即解的最优化问题：该问题即解的最优化问题：22 3/5/2022 1 1）编写）编写M-M-文件文件 curvefun1.mcurvefun1.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:1000tdata

20、=100:100:1000cdata=cdata=1e-03* *4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;6.50,6.59; x0=0.2,0.05,0.05; x0=0.2,0.05,0.05; x=lsqcurvefit (curvefun1,x0,tdata,cdata)x=lsqcurvefit (curvefun1,x0,tdata,cdata) f= f= curvefun1(x,tdata) F(x，tdata)= ，x=(a，b，k)Tk

21、tktbeabea),(10102. 002. 0解法解法1 1. 用命令用命令lsqcurvefitlsqcurvefit100200 30040050060070080090010004.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59jt310jc23 3/5/2022 3 3）运算结果）运算结果：f =0.0043 0.0051 0.0056 0.0059 0.0061 f =0.0043 0.0051 0.0056 0.0059 0.0061 0.0062 0.0062 0.0063 0.0063 0.0063 0.0062 0.0062

22、 0.0063 0.0063 0.0063 x =0.0063 -0.0034 0.2542 x =0.0063 -0.0034 0.25424）结论）结论：a=0.0063, b=-0.0034, k=0.2542 0.02 0.25420.0050840.00630.00340.00630.0034ttc tee24 3/5/2022 1）编写编写M-M-文件文件 curvefun2.mcurvefun2.m function f=curvefun2(x) tdata=100:100:1000; cdata=1e-03*4.54,4.99,5.35,5.65,5.90, 6.10,6.26

23、,6.39,6.50,6.59; f=x(1)+x(2)*exp(-0.02*x(3)*tdata)- cdata2）输入命令）输入命令: x0=0.2,0.05,0.05;x=lsqnonlin(curvefun2,x0)f= curvefun2(x)函数函数curvefun2的自变量是的自变量是x，cdata和和tdata是已知参数，故应是已知参数，故应将将cdata tdata的值写在的值写在curvefun2.m中中解法解法 2 2 用命令用命令lsqnonlinlsqnonlin 1010.020.02110,(,)ktktTf xF x tdata ctadaabecabec x=

24、 x=（a a，b b，k k）25 3/5/2022 3 3）运算结果为）运算结果为 f =1.0e-003 f =1.0e-003 * *（0.2322 -0.1243 -0.2495 -0.2413 0.2322 -0.1243 -0.2495 -0.2413 -0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792-0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792） x =0.0063 -0.0034 0.2542x =0.0063 -0.0034 0.2542可以看出，两个命令的计算结果是相同的可以看出，两个命令的计算

25、结果是相同的。4）结论）结论：即拟合得即拟合得a=0.0063 b=-0.0034 k=0.25420.0063 b=-0.0034 k=0.254226 3/5/2022 插值问题插值问题27 3/5/2022 拟合与插值的关系拟合与插值的关系说明：说明：函数插值与曲线拟合都是要根据一组数据构造一个函数作为近函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上完全不同。似，由于近似的要求不同，二者的数学方法上完全不同。 实例：实例：下面数据是某次实验所得，希望得到x和 f之间的关系？x124791 21 31 51 7f1 .53 .96 .611

26、.71 5 .61 8 .81 9 .62 0 .62 1 .1MATLAB(cn)问题：问题：给定一批数据点，需确定满足特定要求的曲线或曲面解决方案：解决方案：若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，就是数据拟合数据拟合，又称曲线拟合或曲面拟合。若要求所求曲线（面）通过所给所有数据点，就是插值问题插值问题；28 3/5/2022 最临近插值、线性插值、样条插值与曲线拟合结果：最临近插值、线性插值、样条插值与曲线拟合结果：0246810121416180510152025已知数据点spline三次多项式插值0246810121416180510152025已知数据点

27、linest三次多项式插值0246810121416180510152025已知数据点nearest三次多项式插值29 3/5/2022 一一 维维 插插 值值一、插值的定义一、插值的定义二、插值的方法二、插值的方法三、用三、用Matlab解插值问题解插值问题30 3/5/2022 返回返回二维插值二维插值一、二维插值定义一、二维插值定义二、网格节点插值法二、网格节点插值法三、用三、用MatlabMatlab解插值问题解插值问题最邻近插值最邻近插值分片线性插值分片线性插值双线性插值双线性插值网格节点数据的插值网格节点数据的插值散点数据的插值散点数据的插值31 3/5/2022 一维插值的定义一

28、维插值的定义已知已知 n+1个节点个节点, 1 , 0(),(njyxjj其中其中jx互不相同，不妨设互不相同，不妨设),10bxxxan求任一插值点求任一插值点)(*jxx 处的插值处的插值.*y0 x1xnx0y1y*x*y32 3/5/2022 构造一个构造一个(相对简单的相对简单的)函数函数),(xfy 通过全部节点通过全部节点, 即即), 1 ,0()(njyxfjj再用再用)(xf计算插值，即计算插值，即).(*xfy 0 x1xnx0y1y*x*y返回返回33 3/5/2022 的函数值的函数值 nyyy,10已知已知 y=f(x) 在在n+1 个点个点01na xxxbnkyx

29、fxpkkkn, 1 , 0,)()( 构造构造n次多项式次多项式 pn(x) ,使得使得从而得到从而得到 f(x) 的近似计算式的近似计算式 ,),()(baxxpxfn 34 3/5/2022 求解求解 L1(x)=a1 x+a0ixiy0 x1x0y1y已知已知使得使得 L1(xi) = yi . (i=0,1)(1xLy )(xfy 0 x1xOxy如果令如果令1010)(xxxxxl 0101)(xxxxxl 0)(, 1)(1000 xlxl则称则称 l0(x) , l1(x)为为x0, x1上的上的线性插值线性插值基函数基函数。这时。这时10100101yxxxxyxxxx 根据

30、点斜式得到根据点斜式得到)()(0010101xxxxyyyxL 并称其为一次并称其为一次Lagrange插值多项式。插值多项式。1011()0, ( )1l xl xf(x) L1(x)=y0l0(x)+y1 l1(x)35 3/5/2022 求解求解 L2(x)=a2x2+a1 x+a0使得使得 L2(xi)=yi , i=0,1,2.关于关于二次多项式二次多项式的构造采用如下方法：令的构造采用如下方法：令已知已知ixiy0 x1x0y1y2x2y并由插值条件并由插值条件得到得到)(20100 xxxxyA )(21011xxxxyB )(12022xxxxyC L2(x)=A(x-x1)

31、(x-x2)+B(x-x0)(x-x2)+C(x-x0)(x-x1)L2(x0)=y0 , L2(x1)=y1 , L2(x)=y236 3/5/2022 于是得到于是得到2120210121012002010212)()()()()()()(yxxxxxxxxyxxxxxxxxyxxxxxxxxxL 这时：这时：f(x) L2(x)=y0l0(x)+y1 l1(x)+ y2 l2(x)如果令如果令则有则有 jijixlijij, 0, 1)( )()()(2010210 xxxxxxxxxl )()()(2101201xxxxxxxxxl )()()(1202102xxxxxxxxxl 并称

32、其为并称其为二次二次Lagrange插值多项式。插值多项式。 基函数表示基函数表示 x0, x1, x2上的二次插值基函数37 3/5/2022 称为拉格朗日插值基函数拉格朗日插值基函数。n0iiiny)x(L)x(P 已知函数f(x)在n+1个点x0,x1,xn处的函数值为 y0,y1,yn 。求一n次多项式函数Pn(x)，使其满足： Pn(xi)=yi,i=0,1,n. 解决此问题的拉格朗日插值多项式公式如下其中Li(x) 为n次多项式：)xx()xx)(xx()xx)(xx()xx()xx)(xx()xx)(xx()x(Lni1ii1ii1i0in1i1i10i拉格朗日拉格朗日(Lagr

33、ange)插值插值38 3/5/2022 拉格朗日拉格朗日(Lagrange)插值插值特别地特别地:两点一次两点一次(线性线性)插值多项式插值多项式: 101001011yxxxxyxxxxxL三点二次三点二次(抛物抛物)插值多项式插值多项式: 2120210121012002010212yxxxxxxxxyxxxxxxxxyxxxxxxxxxL .,满足插值条件直接验证可知xLn39 3/5/2022 拉格朗日多项式插值的这种振荡现象叫 Runge现象现象55,11)(2xxxg 采用拉格朗日多项式插值：选取不同插值节点个数n+1，其中n为插值多项式的次数，当n分别取2,4,6,8,10时，

34、绘出插值结果图形.例例40 3/5/2022 41 3/5/2022 分段线性插值分段线性插值其它,0,)()()(1111110jjjjjjjjjjjnjjjnxxxxxxxxxxxxxxxlxlyxL计算量与n无关;n越大，误差越小.nnnxxxxgxL0),()(limxjxj-1xj+1x0 xnxoy42 3/5/2022 66,11)(2xxxg例例用分段线性插值法求插值用分段线性插值法求插值,并观察插值误差并观察插值误差.在在-6,6中平均选取中平均选取41个点作插值个点作插值,结果如图示结果如图示43 3/5/2022 比分段线性插值更光滑。比分段线性插值更光滑。xyxi-1

35、xiab 在数学上，光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续，则称该曲线具有k阶光滑性。 光滑性的阶次越高，则越光滑。是否存在较低次的分段多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例子。三次样条插值三次样条插值44 3/5/2022 , 1,),()(1nixxxxsxSiii,)()3), 1 ,0()()2), 1()()10223niiiiiiixxCxSniyxSnidxcxbxaxs) 1, 1()()(),()(),()(111 nixsxsxsxsxsxsiiiiiiiiiiii自然边界条件）(0)()()40 nxSxS)(,)4)3)2xSdcba

36、iiiig g( (x x) )为被插值函数为被插值函数。lim( )( )nS xg x 三次样条插值三次样条插值45 3/5/2022 例例66,11)(2xxxg用三次样条插值选取用三次样条插值选取11个基点计算插值个基点计算插值46 3/5/2022 用用MATLABMATLAB作插值计算作插值计算一维插值函数：一维插值函数：yi=interp1(x，y，xi，method)插值方法插值方法被插值点被插值点插值节点插值节点xixi处的插处的插值结果值结果nearest ：最邻近插值：最邻近插值linear ： 线性插值；线性插值；spline ： 三次样条插值；三次样条插值；cubic

37、 ： 立方插值。立方插值。缺省时：缺省时： 分段线性插值。分段线性插值。 注意：所有的插值方法都要求注意：所有的插值方法都要求x x是单调的，并且是单调的，并且xi不不能够超过能够超过x的范围。的范围。47 3/5/2022 例：在例：在1-121-12的的1111小时内，每隔小时内，每隔1 1小时测量一次温小时测量一次温度，测得的温度依次为：度，测得的温度依次为：5 5，8 8，9 9，1515，2525，2929，3131，3030，2222，2525，2727，2424。试估计每隔。试估计每隔1/101/10小时的温度小时的温度值。值。hours=1:12;temps=5 8 9 15

38、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)48 3/5/2022 49 3/5/2022 xy机翼下轮廓线X035791 11 21 31 41 5Y01 . 21 . 72 . 02 . 12 . 01 . 81 . 21 . 01 . 6例例 已知飞机下轮廓线上数据如下，求已知飞机下轮廓线上数据如下，求x每改变每改变0

39、.1时的时的y值。值。50 3/5/2022 51 3/5/2022 二维插值的定义二维插值的定义 xyO O第一种（网格节点）：第一种（网格节点）：52 3/5/2022 已知已知 m n个节点个节点 ),2 , 1;,.,2 , 1(),(njmizyxijji 其中其中jiyx ,互不相同，不妨设互不相同，不妨设bxxxam 21dyyycn 21 构造一个二元函数构造一个二元函数),(yxfz 通过全部已知节点通过全部已知节点,即即再用再用),(yxf计算插值，即计算插值，即).,(*yxfz ),1 ,0;,1 ,0(),(njmizyxfijji 53 3/5/2022 第二种（散

40、乱节点）：第二种（散乱节点）： yx0 054 3/5/2022 已知已知n个节点个节点),.,2 , 1(),(nizyxiii 其中其中),(iiyx互不相同，互不相同， 构造一个二元函数构造一个二元函数),(yxfz 通过全部已知节点通过全部已知节点,即即),1 ,0(),(nizyxfiii 再用再用),(yxf计算插值，即计算插值，即).,(*yxfz 55 3/5/2022 注意：注意：最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。最邻近插值最邻近插值x y(x1, y1)(x1, y2)(x2, y1)(x2, y2)O O 二维或高维情形的最邻近插值，与被插值点

41、最邻近的节点的函数值即为所求。56 3/5/2022 将四个插值点（矩形的四个顶点）处的函数值依次简记为： 分片线性插值分片线性插值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)=f457 3/5/2022 插值函数为：jii1ij1jy)xx(xxyyy)yy)(ff ()xx)(ff (f)y, x(fj23i121第二片(上三角形区域)：(x, y)满足iii1ij1jy)xx(xxyyy插值函数为：)xx)(ff ()yy

42、)(ff (f)y, x(fi43j141注意注意：(x, y)当然应该是在插值节点所形成的矩形区域内。显然，分片线性插值函数是连续的；分两片的函数表达式如下：第一片(下三角形区域)： (x, y)满足58 3/5/2022 双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下：)dcy)(bax()y, x(f其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。双线性插值双线性插值x y(x1, y1)(x1, y2)(x2, y1)(x2, y2)O O59 3/5/2022 要求要求x0,y0 x0,y0单调；单调；x x，

43、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 双三次插值双三次插值缺省时缺省时, , 双线性插值双线性插值60 3/5/2022 例：测得平板表面例：测得平板表面3 3* *5 5网格点处的温度

44、分别为：网格点处的温度分别为： 82 81 80 82 84 82 81 80 82 84 79 63 61 65 81 79 63 61 65 81 84 84 82 85 86 84 84 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.先在三维坐标画出原始数据，画出粗糙的温度分布曲图.61 3/5/2022 62 3/5/2022 再输入以下命

45、令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi)画出插值后的温度分布曲面图. 2以平滑数据,在x、y方向上每隔0.2个单位的地方进行插值.63 3/5/2022 64 3/5/2022 例例 山区地貌：山区地貌： 在某山区测得一些地点的高程如下表。平面区域为在某山区测得一些地点的高程如下表。平面区域为 1200=x=4000,1200=y0)k(0)模型假设模型假设1. 1. 机体看作一个房室，室内血药浓度均匀机体看作一个房室，室内血药浓度均匀一室模型一室模型模型建立模型建立d/c(0) 3得：由假设-kcdtdc 2得：由假设ktevdtc)( 在此，在此，d=300mg，t