1、1 11/26/2022 数学建模教程数学建模教程拟拟 合与合与 插插 值值 2 11/26/2022 q在大量的应用领域中,人们经常面临这样的问题:在大量的应用领域中,人们经常面临这样的问题:给定一批数据点,需确定满足特定要求的曲线或给定一批数据点,需确定满足特定要求的曲线或曲面。对这个问题有两种方法。曲面。对这个问题有两种方法。q 一种是一种是插值法插值法,数据假定是正确的,要求以某种方法描述数,数据假定是正确的,要求以某种方法描述数据点之间所发生的情况。据点之间所发生的情况。q 另一种方法是另一种方法是曲线拟合或回归曲线拟合或回归。人们设法找出某条光滑曲线,。人们设法找出某条光滑曲线,它
2、最佳地拟合数据,但不必要经过任何数据点。它最佳地拟合数据,但不必要经过任何数据点。q 本专题的目的之一是:了解插值和拟合的基本内容及本专题的目的之一是:了解插值和拟合的基本内容及方法;方法;函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似,函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似,由于近似的要求不同,二者的数学方法上是完全不同的。由于近似的要求不同,二者的数学方法上是完全不同的。3 11/26/2022 假设已获得某函数关系的假设已获得某函数关系的成批离散实验数据成批离散实验数据或或观观测数据测数据,拟合问题就是为这样的大量离散数据建立对,拟合问题就是为这样的大量离散数
3、据建立对应的、近似的连续模型的一种应用基础问题。所建立应的、近似的连续模型的一种应用基础问题。所建立的模型的基本形式是一条曲线(一元曲线),称为的模型的基本形式是一条曲线(一元曲线),称为拟拟合曲线或经验公式。合曲线或经验公式。通常采用通常采用“误差的平方和最小误差的平方和最小”的原则,即的原则,即最小最小二乘拟合问题。二乘拟合问题。它它不要求不要求目标模型(即拟合曲线)目标模型(即拟合曲线)精确地精确地过已知过已知的各离散点,只要求目标模型的各离散点,只要求目标模型符合符合已知离散点分布的已知离散点分布的总体轮廓总体轮廓,并与已知离散点的误差,并与已知离散点的误差按某种意义按某种意义尽量地尽
4、量地小小。一、拟合问题一、拟合问题4 11/26/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 11/26/2022 拟拟 合合 问问 题题 引引 例例 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.
5、32 7.45 5.24 3.01已知一室模型快速静脉注射下的血药浓度数据已知一室模型快速静脉注射下的血药浓度数据(t=0注射注射300mg)求血药浓度随时间的变化规律求血药浓度随时间的变化规律c(t).作半对数坐标系作半对数坐标系(semilogy)下的图形下的图形00(),ktc tc eck为待定系数024681001011026 11/26/2022 iixy(,)im 1,2,L已知一组观测数据:已知一组观测数据:x()要求在某特定函数类要求在某特定函数类 中寻找一个函数中寻找一个函数 作为作为x()的近似函数,使得二者在节点产生的残差的近似函数,使得二者在节点产生的残差yf x()
6、iiixf x()(),im 1,2,L按按某种度量标准为最小某种度量标准为最小。mmiiiiixy222200()常用原则常用原则:残差平方和最小:残差平方和最小min曲曲 线线 拟拟 合合 问问 题题 的的 提提 法法7 11/26/2022 线性最小二乘拟合函数线性最小二乘拟合函数的选取的选取 +=a1+a2x+=a1+a2x+a3x2+=a1+a2x+a3x2=a1+a2/x=aebx=ae-bxnnxaxaxax0011()()()()L1.1.通过机理分析建立数学模型来确定通过机理分析建立数学模型来确定 ;()x 2.2.将数据将数据 (xi,yi)i=1,n 作图,通过直观判断确
7、定作图,通过直观判断确定 :()x 8 11/26/2022 曲线拟合问题最常用的解法曲线拟合问题最常用的解法线性最小二乘法的基本思路线性最小二乘法的基本思路 第二步:确定确定a1,a2,an 的的准则(最小二乘准则):准则(最小二乘准则):使使n个点个点(xi,yi)与与曲线曲线 y=(x)的的距离距离 i 的平方和最小的平方和最小。记记 221211(,)()mmniiiiiJ aaaxy 问题归结为,求问题归结为,求 a1,a2,an 使使 J(a1,a2,an)最小。最小。第一步:先选定一组函数先选定一组函数 使使其中其中 a1,a2,an 为待定系数。为待定系数。nxxxxx01()
8、()(),(),()Lnnxaxaxax0011()()()()L9 11/26/2022 方程组没有通常意义下的解,这类方程组称为方程组没有通常意义下的解,这类方程组称为超定方超定方当线性方程组的方程个数多于未知数的个数时,当线性方程组的方程个数多于未知数的个数时,设线性方程组为设线性方程组为程组程组或或矛盾方程组矛盾方程组。11112211211222221122nnnnmmmnnma xa xa xba xa xaxbaxaxaxb LLLm nAxb 最小二乘法的求解:预备知识最小二乘法的求解:预备知识10 11/26/2022 若能找到一组向量若能找到一组向量*12(,)Tnxxxx
9、 L令令最小,其中最小,其中使得使得20mii 1niijjija xb 则称则称 为该为该超定方程组的最小二乘解超定方程组的最小二乘解。*x1,2,im L2212111(,)mmnniijjiiijJ xxxa xb L L由多元函数取极值的必要条件有由多元函数取极值的必要条件有求其最小值。求其最小值。0kJx 1,2,kn L11 11/26/2022 21miiJ 211112211()nna xa xa xbL221122222()nna xa xaxbL21122()mmmnnmaxaxaxbLL L LkJx 1111122112()knnaa xa xa xbL22112222
10、22()knnaa xa xaxbL11222()0mkmmmnnmaaxaxaxbLL L L即即1111122111()knnkaa xa xa xa bL2211222222()knnkaa xa xaxa bL1122()0mkmmmnnmkmaaxaxaxabLL L L12 11/26/2022 1111122111()knnkaa xa xa xa bL2211222222()knnkaa xa xaxa bL1122()0mkmmmnnmkmaaxaxaxabLL L L12,kkmkaaaL12mbbbM12,kkmkaaaL11112122122212nnnmmmnxaaa
11、xaaaxaaa LLMMMLML 0kJx 1,2,kn L故得故得13 11/26/2022 12mbbbM11112122122212nnnmmmnxaaaxaaaxaaa LLMMMLML112111222212mmnnmnaaaaaaaaaLLMMLML112111222212mmnnmnaaaaaaaaaLLMMLML 即即TTA AxA b 称为正则方程组。称为正则方程组。该方程组的解即为超定方程组的最小二乘解。该方程组的解即为超定方程组的最小二乘解。14 11/26/2022 (2)求解正则方程组得最小二乘解。求解正则方程组得最小二乘解。用最小二乘法解超定方程组的步骤:用最小二
12、乘法解超定方程组的步骤:(1)计算计算 和和 ,得正则方程组,得正则方程组 。TA bTA ATTA AxA b 15 11/26/2022 1.02,4abixiyyabx 例例解得最小二乘解为解得最小二乘解为得方程组得方程组解:解:已知试验数据,用最小二乘法求拟合直线已知试验数据,用最小二乘法求拟合直线0.00.20.40.60.80.91.92.83.34.210.00.910.21.910.42.810.63.310.84.2ab 故得拟合直线故得拟合直线1.024yx 16 11/26/2022 可线性化模型的最小二乘拟合可线性化模型的最小二乘拟合很多实际问题中,变量间并非线性关系,
13、但拟合很多实际问题中,变量间并非线性关系,但拟合曲线可视为曲线可视为 的形式,的形式,()()f yabg x 指数函数指数函数如双曲线如双曲线即将非线性化问题转化为线性问题。即将非线性化问题转化为线性问题。yabx$令令()xg x$则有则有()yf y$11abyx bxyae lnlnyabx 17 11/26/2022 例例给定如下观测数据,试用指数曲线给定如下观测数据,试用指数曲线 进行拟合。进行拟合。bxyae ixiy1.01.251.51.752.05.15.796.537.458.46解解:bxyae lnlnyabx 令令ln,lnyy aa$yabx$则有则有ixiy$1
14、.01.251.51.752.01.6291.7561.8762.0082.135故故1.6291.7561.252.1352.0ababab$L L$解此超定方程组得解此超定方程组得1.122,0.505ab$3.071,a 则拟合曲线为则拟合曲线为0.5053.071xye 18 11/26/2022 多变量数据拟合多变量数据拟合12(,)nyf xxx L有时变量间关系为多元函数关系,有时变量间关系为多元函数关系,有如下观测数据有如下观测数据观测次数12m1xnx2xy11x12x1mxMMMMMM21x22x2mx1nx2nxnmx1x2xmx假定变量假定变量y与与n个变量个变量xi间
15、为线性关系,间为线性关系,可设拟合方程为可设拟合方程为01122()nnxaa xa xa x L19 11/26/2022 第第i组观测数据对应的残差为组观测数据对应的残差为下面考虑用最小二乘原理确定拟合方程的系数下面考虑用最小二乘原理确定拟合方程的系数ai。()iiixy1mii 按照最小二乘原理,应使按照最小二乘原理,应使 最小。最小。01122iinniiaa xa xa xyL令令011111nnaa xa xyL011222nnaa xa xyLL L L011mnnmmaa xa xyL 求解该超定方程组的求解该超定方程组的最小二乘解即可。最小二乘解即可。01122()nnxaa
16、 xa xa x L20 11/26/2022 用用MATLAB解拟合问题解拟合问题1 1、线性最小二乘拟合、线性最小二乘拟合2 2、非线性最小二乘拟合、非线性最小二乘拟合21 11/26/2022 用用MATLAB作线性最小二乘拟合作线性最小二乘拟合1.1.作作多多项式项式f(x)=a1xm+amx+am+1拟合拟合,可利用已有命令可利用已有命令:a=polyfit(x,y,m)3.3.对超定方程组对超定方程组)(11nmyaRnmmn可得最小二乘意义下的解。可得最小二乘意义下的解。,用,用yRa2.2.多项式在多项式在x x处的值处的值y y的计算命令的计算命令:y=y=polyvalpo
17、lyval(a a,x x)输出拟合多项式系数输出拟合多项式系数a=a1,am,am+1(数组)数组)输入同长度输入同长度数组数组X,Y拟合多项式拟合多项式次数次数22 11/26/2022 即要求即要求 出二次多项式出二次多项式:3221)(axaxaxf中中 的的123(,)Aa a a 使得使得:1121()iiif xy最小例例 对下面一组数据作二次多项式拟合对下面一组数据作二次多项式拟合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 21111222223211
18、1111,1,1,1xxyayxxaayxx23 11/26/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此时解法解法1 1解超定方程的方法解超定方程的方法2)计算结果)计算结果:=-9.8108,20.1293,-0.03170317.01293.208108.9)(2xxxfRAy24 11/26/2022 25 11/26/2022 1)输入命令:)输入命令:x=0:0.1:1;
19、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,k+,x,z,r)%作出数据点和拟合曲线的图形作出数据点和拟合曲线的图形2)计算结果:)计算结果:=-9.8108,20.1293,-0.0317解法解法2用多项式拟合的命令用多项式拟合的命令00.20.40.60.81-20246810120317.01293.208108.9)(2xxxf26 11/26/2022 1.1.lsqcurvefitlsqcurvefit已知数据点数据点:xdat
20、axdata=(xdata1,xdata2,xdataxdatan n)ydataydata=(ydataydata1 1,ydataydata2 2,ydataydatan n)用用MATLAB作非线性最小二乘拟合作非线性最小二乘拟合两个求非线性最小二乘拟合的函数:两个求非线性最小二乘拟合的函数:lsqcurvefit、lsqnonlin。相同点和不同点:两相同点和不同点:两个命令都要先建立个命令都要先建立M-M-文件文件fun.mfun.m,定义函,定义函数数f(x)f(x),但定义但定义f(x)f(x)的方式不同的方式不同。211(,)2niiiF x xdataydata最小 lsqc
21、urvefitlsqcurvefit用以求含参量用以求含参量x x(向量)的向量值函数向量)的向量值函数F(x,xdataF(x,xdata)=)=(F F(x x,xdataxdata1 1),),F F(x x,xdataxdatan n)T T使得使得 27 11/26/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=lsq
22、curvefit(fun,x0,xdata,ydata,);(5)x,resnorm,residual=lsqcurvefit(fun,x0,xdata,ydata,);fun是一个事先建立的是一个事先建立的定义函数定义函数F(x,xdata)的的M-文件文件,自变量为自变量为x和和xdata说明:说明:x=lsqcurvefit(fun,x0,xdata,ydata,options);迭代初值迭代初值已知数据点已知数据点选项见无选项见无约束优化约束优化28 11/26/2022 lsqnonlin用以求含参量用以求含参量x x(向量)的向量值函数向量)的向量值函数 f(x)f(x)=(f=(
23、f1 1(x),f(x),f2 2(x),f(x),fn n(x)(x)T T ,使得,使得 最小。最小。其中其中 fi(x)=f(x,xdatai,ydatai)=F(x,xdatai)-ydatai22212()()()()()Tnfx f xf xfxfx2.lsqnonlin已知数据点:已知数据点:xdataxdata=(xdata1,xdata2,xdataxdatan n)ydataydata=(ydataydata1 1,ydataydata2 2,ydataydatan n)29 11/26/2022 输入格式:输入格式:1)x=lsqnonlin(fun,x0);2)x=ls
24、qnonlin(fun,x0,lb,ub);3)x=lsqnonlin(fun,x0,lb,ub,options);4)x,resnorm=lsqnonlin(fun,x0,);5)x,resnorm,residual=lsqnonlin(fun,x0,);说明:说明:x=lsqnonlinlsqnonlin(fun,x0,options););fun是一个事先建立的是一个事先建立的定义函数定义函数f(x)的的M-文件,文件,自变量为自变量为x迭代初值迭代初值选项见无选项见无约束优化约束优化30 11/26/2022 100200 30040050060070080090010004.54 4
25、.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59jt310jc100.022111min(,)22jktjjF a b kabec例例2 用下面一组数据拟合用下面一组数据拟合 中的参数中的参数a,b,k0.0.2()ktc tabe 该问题即解的最优化问题:该问题即解的最优化问题:31 11/26/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)
26、=k;2)输入命令输入命令tdatatdata=100:100:1000=100:100:1000cdatacdata=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=x=lsqcurvefitlsqcurvefit(curvefun1,x0,tdata,cdata)(curvefun1,x0,tdata,cdata)f=f=curvefun1(x,tdata)
27、F(x,tdata)=,x=(a,b,k)Tktktbeabea),(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.59jt310jc32 11/26/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
28、.0063 0.0062 0.0062 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 tee33 11/26/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,
29、5.90,6.10,6.26,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 ctadaab
30、ecabec x=x=(a a,b b,k k)34 11/26/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可以看出,两个命令的计算结果是相同的可以看出,两个命令的计
31、算结果是相同的。4)结论)结论:即拟合得即拟合得a=0.0063 b=-0.0034 k=0.25420.0063 b=-0.0034 k=0.254235 11/26/2022 设有一个年产设有一个年产240吨的食品加工厂吨的食品加工厂,需要统计其生产费用需要统计其生产费用,但由于该厂各项资料不全但由于该厂各项资料不全,无法统计。在这种情况下无法统计。在这种情况下,统计统计部门收集了设备、生产能力和该厂大致相同的五个食品加部门收集了设备、生产能力和该厂大致相同的五个食品加工厂的产量与生产费用资料如下表:工厂的产量与生产费用资料如下表:如何根据已知五个食品加工厂的资料较准确地估计该厂的如何根据
32、已知五个食品加工厂的资料较准确地估计该厂的生产费用呢?生产费用呢?引例引例插值问题插值问题36 11/26/2022 37 11/26/2022 (2)有的函数甚至没有表达式,只是一种表格函数,)有的函数甚至没有表达式,只是一种表格函数,(1)如果函数表达式本身比较复杂,且需要多次)如果函数表达式本身比较复杂,且需要多次实际问题中经常要涉及到函数值的计算问题:实际问题中经常要涉及到函数值的计算问题:重复计算时,计算量会很大;重复计算时,计算量会很大;方便且表达简单的函数来近似代替,这就是方便且表达简单的函数来近似代替,这就是插值插值问题。问题。实际对于这两种情况,我们都需要寻找一个计算实际对于
33、这两种情况,我们都需要寻找一个计算而我们需要的函数值可能不在该表格中。而我们需要的函数值可能不在该表格中。问题问题:插值与拟合的区别?:插值与拟合的区别?38 11/26/2022 一一 维维 插插 值值一、一、插值的定义插值的定义二、插值的方法二、插值的方法三、用三、用Matlab解插值问题解插值问题拉格朗日插值拉格朗日插值牛顿插值牛顿插值反插值反插值三次样条插值三次样条插值分段插值法分段插值法39 11/26/2022 二维插值二维插值一、一、二维插值定义二维插值定义二、二、网格节点插值法网格节点插值法三、三、用用MatlabMatlab解插值问题解插值问题最邻近插值最邻近插值分片线性插值
34、分片线性插值双线性插值双线性插值网格节点数据的插值网格节点数据的插值散点数据的插值散点数据的插值40 11/26/2022 一维插值的定义一维插值的定义已知已知 n+1个节点个节点,1,0(),(njyxjj其中其中jx互不相同,不妨设互不相同,不妨设),10bxxxan求任一插值点求任一插值点)(*jxx 处的插值处的插值.*y0 x1xnx0y1y*x*y41 11/26/2022 构造一个构造一个(相对简单的相对简单的)函数函数(),yx 通过全部节点通过全部节点,即即()(0,1,)jjxyjn 再用再用()x 计算插值,即计算插值,即*().yx 0 x1xnx0y1y*x*y42
35、11/26/2022 代数插值的唯一性代数插值的唯一性是惟一的。是惟一的。定理定理个不同节点,满足插值条件个不同节点,满足插值条件n 1ip xy()的的n次插值多项式次插值多项式in 1,2,Lnnnnp xa xaxa xa 1110()L当选择代数多项式作为插值函数时,称为代数多当选择代数多项式作为插值函数时,称为代数多项式插值。项式插值。代数插值代数插值43 11/26/2022 一、线性插值一、线性插值(n=1,一次插值),一次插值)求求解解 L1(x)=a1 x+a0ixiy0 x1x0y1y已知已知使得使得 L1(xi)=yi.(i=0,1)(1xLy )(xfy 0 x1xOx
36、y如果令如果令1010)(xxxxxl 0101)(xxxxxl 0)(,1)(1000 xlxl则称则称 l0(x),l1(x)为为x0,x1上的上的线性插值线性插值基函数基函数。这时。这时10100101yxxxxyxxxx 根据点斜式得到根据点斜式得到)()(0010101xxxxyyyxL 并称其为一次并称其为一次Lagrange插值多项式。插值多项式。1011()0,()1l xl xf(x)L1(x)=y0l0(x)+y1 l1(x)44 11/26/2022 二、抛物线插值二、抛物线插值(n=2,二次插值),二次插值)求求解解 L2(x)=a2x2+a1 x+a0使得使得 L2(
37、xi)=yi ,i=0,1,2.已知已知ixiy0 x1x0y1y2x2y45 11/26/2022 抛物插值抛物插值仿照线性函数的构造方法,构造仿照线性函数的构造方法,构造L xlx ylx ylx y001122()()()()且且lx 00()1 lx 01()0lx 10()0 lx 11()1lx 02()0lx 12()0lx 20()0 lx 21()0lx 22()1其中要求其中要求 均为二次多项式。均为二次多项式。il x i (),1,2,3设设lxA xxxx012()()()即即由由求求Alx 00()1A xxxx0102()()1Axxxx 01021()()xxx
38、xlxxxxx 1200102()()()()()46 11/26/2022 故故同理同理xxxxlxxxxx 1200102()()()()()xxxxlxxxxx 0211012()()()()()xxxxlxxxxx 0122021()()()()()L xlx ylx ylx y001122()()()()xxxxyxxxx 1200102()()()()xxxxyxxxx 0211012()()()()xxxxyxxxx 0122021()()()()47 11/26/2022 例题例题求线性插值函数。求线性插值函数。已知已知 的函数表的函数表 yf x ()13x12y2 1解:解
39、:L x ()12(1)xx(3)(2)xx2519822(13)(12)xx(1)(2)(31)(32)xx(1)(3)(21)(23)xx22.59.58fL (1.5)(1.5)0.62548 11/26/2022 多项式可使用类似方法。多项式可使用类似方法。分析分析由由 个不同节点个不同节点 构造构造 次次n 1niixy in(,)1,2,Lik nnL xlx ylx ylx y0011()()()()Lkknkkkkkkkknxxxxxxxxxxlxxxxxxxxxxx 12111211()()()()()()()()()()()LLLL其中其中且且kklx()1 kilx()0
40、上述多项式称为上述多项式称为拉格朗日插值多项式拉格朗日插值多项式。拉格朗日插值多项式拉格朗日插值多项式nkkklx y 0()49 11/26/2022 插值余项和误差估计插值余项和误差估计且依赖于且依赖于 。其中其中xab niixxx 0()(),nfR xf xL xxn (1)()()()()()(1)!设设f(x)在在a,b有有n+1阶导数,则阶导数,则n次插值多次插值多定理定理L x()项式项式 对任意对任意 ,有插值余项,有插值余项xa b ,50 11/26/2022 拉格朗日多项式插值的这种振荡现象叫 Runge现象现象55,11)(2xxxg 采用拉格朗日多项式插值:选取不
41、同插值节点个数n+1,其中n为插值多项式的次数,当n分别取2,4,6,8,10时,绘出插值结果图形.例例51 11/26/2022 52 11/26/2022 分段线性插值分段线性插值0111111()(),(),0,nnjjjjjjjjjjjjjjLxy lxxxxxxxxxxlxxxxxx 其 它0lim()(),nnnL xf xxxx xjxj-1xj+1x0 xnxoy53 11/26/2022 66,11)(2xxxg例例用分段线性插值法求插值用分段线性插值法求插值,并观察插值误差并观察插值误差.在在-6,6中平均选取中平均选取41个点作插值个点作插值,结果如图示结果如图示54 1
42、1/26/2022 比分段线性插值更光滑。比分段线性插值更光滑。xyxi-1 xiab 在数学上,光滑程度的定量描述是:函数在数学上,光滑程度的定量描述是:函数(曲线曲线)的的k阶阶导数存在且连续,则称导数存在且连续,则称该曲线具有该曲线具有k阶光滑性阶光滑性。光滑性的阶次越高,则越光滑。是否存在较低次的分光滑性的阶次越高,则越光滑。是否存在较低次的分段多项式达到较高阶光滑性的方法?三次样条插值就是一段多项式达到较高阶光滑性的方法?三次样条插值就是一个很好的例子。个很好的例子。三次样条插值三次样条插值55 11/26/2022 所谓样条函数,从数学上说,就是按照一定光滑性所谓样条函数,从数学上
43、说,就是按照一定光滑性它在每个内节点上具有直到它在每个内节点上具有直到m-1m-1阶连续导数。阶连续导数。要求要求“装配装配”起来的分段多项式。起来的分段多项式。它在每个子区间内是它在每个子区间内是m m次多项式,次多项式,三次样条插值三次样条插值三次样条函数满足:三次样条函数满足:在每个子区间在每个子区间 上是一个三次多项式。上是一个三次多项式。(1)S x()iix x 1,(2)S x()在每个内节点上具有直到二阶的连续导数。在每个内节点上具有直到二阶的连续导数。(3)iiS xy()in 0,1,L1012121nnnS xxx xS xxx xxS xxxx (),(),S()(),
44、M上可设上可设故在每一个小区间故在每一个小区间iixx 1,iiiiiS xa xb xc xd32()求解使用三弯矩法求解使用三弯矩法56 11/26/2022 例例66,11)(2xxxg用三次样条插值选取用三次样条插值选取11个基点计算插值个基点计算插值57 11/26/2022 用用MATLABMATLAB作插值计算作插值计算一维插值函数:一维插值函数:yi=interp1(x,y,xi,method)插值方法插值方法被插值点被插值点插值节点插值节点xixi处的插处的插值结果值结果nearest :最邻近插值最邻近插值linear :线性插值;线性插值;spline :三次样条插值;三
45、次样条插值;cubic :立方插值。立方插值。缺省时:缺省时:分段线性插值。分段线性插值。注意:所有的插值方法都要求注意:所有的插值方法都要求x是单调的,并且是单调的,并且xi不不能够超过能够超过x的范围。的范围。58 11/26/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
46、 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)59 11/26/2022 60 11/26/2022 xy机翼下轮廓线X035791 11 21 31 41 5Y01.21.72.02.12.01.81.21.01.6例例 已知飞机下轮廓线上数据如下,求已知飞机下轮廓线上数据如下,求x每改变每改变0.1时的时的y值
47、。值。61 11/26/2022 62 11/26/2022 二维插值的定义二维插值的定义 xyO O第一种(网格节点):第一种(网格节点):63 11/26/2022 已知已知 m n个节点个节点),2,1;,.,2,1(),(njmizyxijji 其中其中jiyx,互不相同,不妨设互不相同,不妨设bxxxam 21dyyycn 21 构造一个二元函数构造一个二元函数),(yxfz 通过全部已知节点通过全部已知节点,即即再用再用),(yxf计算插值,即计算插值,即).,(*yxfz ),1,0;,1,0(),(njmizyxfijji 64 11/26/2022 第二种(散乱节点):第二种
48、(散乱节点):yx0 065 11/26/2022 已知已知n个节点个节点),.,2,1(),(nizyxiii 其中其中),(iiyx互不相同,互不相同,构造一个二元函数构造一个二元函数),(yxfz 通过全部已知节点通过全部已知节点,即即),1,0(),(nizyxfiii 再用再用),(yxf计算插值,即计算插值,即).,(*yxfz 66 11/26/2022 注意:注意:最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。最邻近插值最邻近插值x y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O O 二维或高维情形的最邻近插值,与被插值点最邻近的节点的函数值即为所
49、求。67 11/26/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)=f468 11/26/2022 插值函数为:jii1ij1jy)xx(xxyyy)yy)(ff()xx)(ff(f)y,x(fj23i121第二片(上三角形区域):(x,y)满足iii1ij1jy)xx(xxyyy插值函数为:)xx)(ff()yy)(ff(f)y,x(fi43j141注意注意:(x,
50、y)当然应该是在插值节点所形成的矩形区域内。显然,分片线性插值函数是连续的;分两片的函数表达式如下:第一片(下三角形区域):(x,y)满足69 11/26/2022 双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下:(,)()()f x yaxb cyd 其中有四个待定系数,利用该函数在矩形的四个顶点(插值节点)的函数值,得到四个代数方程,正好确定四个系数。双线性插值双线性插值x y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O O70 11/26/2022 要求要求x0,y0 x0,y0单调;单调;x x,y y可取可取为矩阵,或为矩阵,或x x取取行向量,行向量