《新编MATLAB&Simulink自学一本通》课件第20章 多项式回归与数据插值.ppt

1、2023-5-5主要内容主要内容 多项式回归多项式回归 插值问题的数学描述插值问题的数学描述 一维插值一维插值 二维插值二维插值 高维插值高维插值2023-5-5第一节第一节 多项式回归多项式回归2023-5-5一、多项式回归模型一、多项式回归模型11212,(0,),1,2,.nniiininiiidiyp xp xp xpNim2023-5-5二、多项式回归的二、多项式回归的MATLAB实现实现 p,S,mu=polyfit(x,y,n)自变量观测值向量自变量观测值向量因变量观测值向量因变量观测值向量多项式阶数多项式阶数1.polyfit函数的用法函数的用法系数向量的估计值系数向量的估计值

2、用于误差估计的结构体变量用于误差估计的结构体变量均值和标准差均值和标准差2023-5-5 y,delta=polyval(p,x,S)用户指定的自变量取值向量用户指定的自变量取值向量因变量估计值向量因变量估计值向量2.polyval函数的用法函数的用法多项式系数向量多项式系数向量用于误差估计的结构体变量用于误差估计的结构体变量误差标准差的估计值向量误差标准差的估计值向量2023-5-5 r=poly2sym(p,v)多项式的符号表达式多项式的符号表达式3.poly2sym函数的用法函数的用法多项式系数向量多项式系数向量自变量符号或取值自变量符号或取值2023-5-5三、多项式回归案例三、多项式

3、回归案例【例例20.1-1】现有我国现有我国2007年年1月至月至2011年年11月的食品零售价格分月的食品零售价格分类指数数据，如表类指数数据，如表20.1-1所示。数据来源：中华人民共和国国家所示。数据来源：中华人民共和国国家统计局网站月度统计数据。试根据这统计局网站月度统计数据。试根据这59组统计数据研究全国食组统计数据研究全国食品零售价格分类指数品零售价格分类指数 y（上年同月（上年同月=100）和时间）和时间 x 之间的关系。之间的关系。序号序号统计月度统计月度上年同月上年同月=100上年同期上年同期=100全国全国城市城市农村农村全国全国城市城市农村农村12007年1月104.91

4、04.4105.9104.9104.4105.922007年2月105.8105.2106.9105.3104.8106.432007年3月107.7107.4108.3106.1105.7107572011年9月113.5113.4114112.6112.4113.1582011年10月111.9111.8112.2112.5112.3113592011年11月108.7108.8108.6112.2112112.62023-5-51.数据的散点图数据的散点图951001051101151201252007.12007.32007.52007.72007.92007.112008.12008

5、.32008.52008.72008.92008.112009.12009.32009.52009.72009.92009.112010.12010.32010.52010.72010.92010.112011.12011.32011.52011.72011.92011.11时间(x)食品零售价格分类指数(y)2023-5-52.四次多项式拟合四次多项式拟合 假设假设 y 关于关于 x 的理论回归方程为的理论回归方程为43212345 yp xp xp xp xp 调用调用polyfit函数求系数估计值函数求系数估计值 p4=polyfit(x,y,4)p4=-0.0001 0.0096 -0

6、.3985 5.5635 94.2769 r=poly2sym(p4);r=vpa(r,5)r=-0.000074268*x4+0.0096077*x3-0.39845*x2+5.5635*x+94.2774320.00010.00960.39855.563594.2769yxxxx 2023-5-53.更高次多项式拟合更高次多项式拟合 调用调用polyfit函数作更高次（大于函数作更高次（大于4次）多项式拟合，并把次）多项式拟合，并把多次拟合的残差的模加以对比，评价拟合的好坏。多次拟合的残差的模加以对比，评价拟合的好坏。p5,S5=polyfit(x,y,5);%5次多项式拟合次多项式拟合

7、S5.normr p6,S6=polyfit(x,y,6);%6次多项式拟合次多项式拟合 S6.normr p7,S7=polyfit(x,y,7);%7次多项式拟合次多项式拟合 S7.normr p8,S8=polyfit(x,y,8);%8次多项式拟合次多项式拟合 S8.normr p9,S9=polyfit(x,y,9);%9次多项式拟合次多项式拟合 S9.normr2023-5-54.拟合效果图拟合效果图 在以上拟合结果的基础上，可以调用在以上拟合结果的基础上，可以调用polyval函数计算给函数计算给定自变量定自变量 处的因变量处的因变量 的预测值，从而绘制拟合效果图，的预测值，从而

8、绘制拟合效果图，从拟合效果图上直观地看出拟合的准确性。从拟合效果图上直观地看出拟合的准确性。代码略。代码略。2023-5-5951001051101151201252007.12007.32007.52007.72007.92007.112008.12008.32008.52008.72008.92008.112009.12009.32009.52009.72009.92009.112010.12010.32010.52010.72010.92010.112011.12011.32011.52011.72011.92011.11时间(x)食品零售价格分类指数(y)原始散点4次多项式拟合6次多项

9、式拟合8次多项式拟合9次多项式拟合2023-5-5第二节第二节 插值问题的数学描述插值问题的数学描述2023-5-5一、什么是插值一、什么是插值 所谓插值就是在已知离散数据的基础上补插连续函数，所谓插值就是在已知离散数据的基础上补插连续函数，使得这条连续曲线（或曲面）通过全部已知的离散数据点，使得这条连续曲线（或曲面）通过全部已知的离散数据点，利用插值方法可通过函数在有限个点处的取值状况，估算利用插值方法可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。出函数在其他点处的近似值。一维插值问题的一维插值问题的数学描述数学描述为：为：已知某一函数已知某一函数 y=g(x)（g(x)的

10、解析表达式可能十分复杂的解析表达式可能十分复杂,也可以是未知的）在区间也可以是未知的）在区间a,b上上n+1个互异点个互异点 xj 处的函数值处的函数值 y j，j=0,1,n，还知道，还知道g(x)在在a,b上有若干阶导数，如何上有若干阶导数，如何求出求出g(x)在在a,b上任一点上任一点x的近似值。的近似值。二、一维插值问题的数学描述二、一维插值问题的数学描述数学描述数学描述 一维插值方法的一维插值方法的基本思想基本思想是：根据是：根据g(x)在区间在区间a,b上上n+1个互异点个互异点 x j（称为（称为节点节点）的函数值）的函数值 y j，j=0,1,n，求，求一个足够光滑、简单便于计

11、算的函数一个足够光滑、简单便于计算的函数f(x)（称为插值函数）（称为插值函数）作为作为g(x)的近似表达式，使得的近似表达式，使得f(x j)=y j，j=0，1，n。（1）然后计算然后计算f(x)在区间在区间a,b(称为插值区间称为插值区间)上点上点x（称为（称为插值插值点点）的值作为原函数）的值作为原函数g(x)（称为被插函数）在此点的近似（称为被插函数）在此点的近似值。求插值函数值。求插值函数f(x)的方法称为插值方法，（的方法称为插值方法，（1）式称为插值）式称为插值条件条件。代数多项式比较简单，常用多项式作为插值函数。代数多项式比较简单，常用多项式作为插值函数。2基本思想基本思想

12、常用的一维插值方法有：常用的一维插值方法有：分段线性插值分段线性插值、拉格朗日拉格朗日（Lagrange）多项式插值）多项式插值、牛顿（牛顿（Newton）插值）插值、Hermite插值插值、最近邻插值最近邻插值、三次样条插值三次样条插值和和B样条插值样条插值等。等。3常用一维插值方法常用一维插值方法 二二维插值问题的维插值问题的数学描述数学描述为：为：已知某已知某二元二元函数函数 z=G(x,y)（解析表达式可能十分复杂解析表达式可能十分复杂,也可以是未知的）在也可以是未知的）在平面区域平面区域D上上N个互异点个互异点(xi,yi)处的函数处的函数值值 zi，i=0,1,N，求，求一个足够光

13、滑、简单便于计算的插一个足够光滑、简单便于计算的插值函数值函数f(x,y)。由插值函数可以计算原函数在平面区域。由插值函数可以计算原函数在平面区域 上任上任意点处的近似值。意点处的近似值。常用的二维插值方法常用的二维插值方法有：分片线性插值、有：分片线性插值、双线性插值、最近邻插值、三次样条插值和双线性插值、最近邻插值、三次样条插值和B样条插值等。样条插值等。三、二维插值问题的数学描述三、二维插值问题的数学描述 是一个细的、可弯曲的木制或塑料条，在飞机或轮是一个细的、可弯曲的木制或塑料条，在飞机或轮船等的设计制造过程中为描绘出光滑的外形曲线船等的设计制造过程中为描绘出光滑的外形曲线(放样放样)

14、所用的工具所用的工具 从物理上讲，样条满足插值点的约束，同时使势能从物理上讲，样条满足插值点的约束，同时使势能达到最小达到最小 三次样条本质上是一段一段的三次多项式拼合而成三次样条本质上是一段一段的三次多项式拼合而成的曲线，在拼接处的曲线，在拼接处,不仅函数是连续的不仅函数是连续的,且一阶和二阶导且一阶和二阶导数也是连续的数也是连续的，19461946年年,Schoenberg,Schoenberg将样条引入数学将样条引入数学,即即所谓的所谓的样条函数样条函数1.什么是样条？什么是样条？四、三次样条插值的数学描述四、三次样条插值的数学描述()(),S xf xa b三三次次样样则则条条称称为为

15、在在上上的的插插值值函函数数(2)(),(),(),S xS xSxa b 都都在在区区间间上上连连续续1(1)(),iiS xxx 在在每每个个小小区区间间上上都都是是三三次次多多项项式式(),S xa b则则称称为为区区间间上上的的三三次次样样条条函函数数(3)()S x若若三三次次样样条条函函数数满满足足(),0,1,iif xy in(),0,1,iiS xy in01,naxxxba b为为区区间间的的一一个个分分割割(),:S xa b如如果果函函数数在在区区间间上上满满足足条条件件2.三次样条插值函数三次样条插值函数3.三次样条插值原理三次样条插值原理niyxSii,1,0,)(

16、在在n n个小区间构造个小区间构造S(x)S(x)，共有，共有n n个三次多项式，需确定个三次多项式，需确定4n4n个参数个参数1,1,)()(limniyxSxSiixxi1,1),()(limnixSxSixxi1,1,)()(lim niMxSxSiixxi在所有节点上在所有节点上n+1个方程个方程在除端点外的节点上在除端点外的节点上3(n-1)个方程个方程这样就得到这样就得到4n-2个方程，为保证待定参数的唯一性，还差两个个方程，为保证待定参数的唯一性，还差两个方程。为此，常用的方法是对边界节点除函数值外附加要求，方程。为此，常用的方法是对边界节点除函数值外附加要求，这就是所谓的边界条

17、件。根据实际问题的不同，三次样条插值这就是所谓的边界条件。根据实际问题的不同，三次样条插值常用到下列常用到下列三类边界条件三类边界条件。r 周期边界条件：周期边界条件：当当y=g(x)是以是以 b-a=x0 -xn 为周期的周期函数时，要求为周期的周期函数时，要求S(x)也是周期函数，故端点处要满足也是周期函数，故端点处要满足此条件称为周期条件。此条件称为周期条件。r m边界条件：边界条件：即给定端点处的一阶导数值。即给定端点处的一阶导数值。00(),()nnS xyS xyr M边界条件：边界条件：即给定端点处的二阶导数值。即给定端点处的二阶导数值。00(),()nnSxySxy(0)(0)

18、,(0)(0)S aS bS aS b2023-5-5第三节第三节 一维插值一维插值一、一、Lagrange（拉格朗日）（拉格朗日）插值插值1拉格朗日线性插值拉格朗日线性插值插值多项式：插值多项式：101()L xaa x011011kkkkyaa xyaa x直线的两点式表达式直线的两点式表达式1111kkkkkkkkxxxxyyxxxx1111(),()kkkkkkkkxxxxlxlxxxxx分别称为节点分别称为节点x0和和x1的一次的一次插值插值基函数基函数插值函数为基函数的线性组合，组合系数插值函数为基函数的线性组合，组合系数就是对应节点上的函数值就是对应节点上的函数值111()()k

19、kkkkkyyL xxxyxx在在 上，过点上，过点11(,),(,)kkkkxyxy1,kkxx2.拉格朗日二次插值拉格朗日二次插值二次插值多项式：二次插值多项式：21111()()()()kkk kkkL xylxy lxylx对于对于 ，可假定：，可假定：1,()()0,()iiixxl xxxx在节点上，但基函数满足：基函数满足：11111()()(),()()kkkkkkkxxxxlxxxxx基函数如何确定？基函数如何确定？1111()()(),()()kkkkkkkxxxxlxxxxx11111()()()()()kkkkkkkxxxxlxxxxx从而从而过点过点1111(,),(

20、,),(,)kkkkkkxyxyxyik11()()()kkklxA xxxx111()()kkkkAxxxxr 构造一组插值基函数（构造一组插值基函数（n次次多项式多项式）011011()()()()()()()()()iiniiiiiiinxxxxxxxxl xxxxxxxxx缺点：缺点：当当 n 比较大时，插值多项式比较大时，插值多项式Ln(x)的收敛性与稳定性变的收敛性与稳定性变 差，逼近效果不理想。差，逼近效果不理想。li(x)是是 n次次多项式多项式，且满足，且满足0()1ijijl xijr 由插值基函数由插值基函数li(x)构造构造 n次拉格朗日插值次拉格朗日插值多项式多项式如

21、下：如下：0()()nni iiL xyl x3.拉格朗日拉格朗日n 次插值次插值r 拉格朗日拉格朗日n次插值的误差估计次插值的误差估计拉格朗日插值产生的截断误差为拉格朗日插值产生的截断误差为Rn(x)，则，则()()()nnR xg xL x(1)1()()()()()(1)!nnnngR xg xL xxn10()()nniixxx1()()(1)!nnMR xxn(1)()ngM4.自编拉格朗日插值函数自编拉格朗日插值函数lagrangefunction y=lagrange(x0,y0,x)if numel(x0)=numel(y0)error(原始坐标点应等长原始坐标点应等长)end

22、x0=x0(:);if any(diff(sort(x0)=0)error(插值节点不能有重复值插值节点不能有重复值)endn=numel(x0);m=numel(x);y0=y0(:);x=x(:);y=zeros(n,m);for i=1:n y(i,:)=y0(i)*prod(repmat(x,n-1,1)-.repmat(x0(1:n=i),1,m)/prod(x0(i)-x0(1:n=i);endy=sum(y);2023-5-5【例例20.3-1】Runge（龙格）现象：用函数（龙格）现象：用函数 在在区间区间-1,1 上产生上产生11个等距节点，然后调用自编个等距节点，然后调用自

23、编lagrange函数作函数作拉格朗日插值。拉格朗日插值。2()1(125)f xx-1-0.500.51-0.500.511.52Xf(x)=11+25x2 插值节点原函数图像Lagrange插值 x0=linspace(-1,1,11);y0=1./(1+25*x0.2);x=linspace(-1,1,100);f=1./(1+25*x.2);y=lagrange(x0,y0,x);plot(x0,y0,ko);hold on;plot(x,f,k,linewidth,2);plot(x,y,k:,linewidth,2);xlabel(X);ylabel($f(x)=frac11+25

24、x2$,Interpreter,latex);legend(插值节点插值节点,原函数图像原函数图像,Lagrange插值插值)2023-5-5二、二、interp1函数函数 yi=interp1(x,Y,xi,method)自变量观测值向量自变量观测值向量因变量观测值向量因变量观测值向量用户指定的插值点横坐标用户指定的插值点横坐标1.interp1函数的用法函数的用法计算得到的近似函数值计算得到的近似函数值指定所用的插值方法指定所用的插值方法2023-5-5【例例20.3-1续续】针对例针对例20.3-1中的数据，调用中的数据，调用interp1函数作一维函数作一维插值。插值。ylin=int

25、erp1(x0,y0,x);yspl=interp1(x0,y0,x,spline);plot(x0,y0,ko);hold on;plot(x,f,k,linewidth,2);plot(x,ylin,:,linewidth,2);plot(x,yspl,r-.,linewidth,2);xlabel(X);ylabel($f(x)=frac11+25x2$,Interpreter,latex);legend(插值节点插值节点,原函数图像原函数图像,分段线性插值分段线性插值,三次样条插值三次样条插值)-1-0.500.5100.20.40.60.81Xf(x)=11+25x2 插值节点原函数

26、图像分段线性插值三次样条插值2023-5-5【例例20.3-2】在加工机翼的过程中，已有机翼断面轮廓线上的在加工机翼的过程中，已有机翼断面轮廓线上的20组坐标点数据，如表组坐标点数据，如表20.3-2所列，其中所列，其中(x,y1)和和(x,y2)分别对应轮分别对应轮廓线的上下线。假设需要得到廓线的上下线。假设需要得到 x坐标每改变坐标每改变0.1时的时的y 坐标，试通坐标，试通过插值方法计算加工所需的全部数据，并绘制机翼断面轮廓线，过插值方法计算加工所需的全部数据，并绘制机翼断面轮廓线，求加工断面的面积。求加工断面的面积。03579111213141501.82.22.73.03.12.92

27、.52.01.601.21.72.02.12.01.81.21.01.6x1y2y2023-5-5x0=0,3,5,7,9,11,12,13,14,15;y01=0,1.8,2.2,2.7,3.0,3.1,2.9,2.5,2.0,1.6;y02=0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6;x=0:0.1:15;ysp1=interp1(x0,y01,x,spline);pp=interp1(x0,y02,spline,pp);ysp2=ppval(pp,x);xx=x,fliplr(x);ysp=ysp1,fliplr(ysp2);plot(x0,x0,y01,

28、y02,o);hold on;plot(xx,ysp,r,linewidth,2);xlabel(X);ylabel(Y);legend(插值节点插值节点,三次样条插值三次样条插值,location,northwest);%截面面积：截面面积：S1=trapz(x,ysp1)-trapz(x,ysp2)S2=trapz(xx,ysp)05101500.511.522.533.5XY 插值节点三次样条插值2023-5-5三、三、spline函数函数 yy=spline(x,Y,xx)自变量观测值向量自变量观测值向量因变量观测值向量因变量观测值向量用户指定的插值点横坐标用户指定的插值点横坐标1.s

29、pline函数的用法函数的用法计算得到的近似函数值计算得到的近似函数值2023-5-5四、四、csape和和csapi函数函数 pp=csape(x,y,conds)自变量观测值向量自变量观测值向量因变量观测值向量因变量观测值向量边界条件参数边界条件参数1.csape函数的用法函数的用法分段多项式形式的插值结果分段多项式形式的插值结果conds参数取值参数取值说说 明明complete或clamped给定端点处的一阶导数值，默认为拉格朗日边界条件not-a-knot非纽结边界条件，csapi函数使用的就是这种边界条件periodic周期边界条件second给定端点处的二阶导数值。默认为0,0，

30、同variational情形variational设定端点处的二阶导数值为02023-5-5 values=csapi(x,y,xx)2.csapi函数的用法函数的用法自变量观测值向量自变量观测值向量因变量观测值向量因变量观测值向量用户指定的插值点横坐标用户指定的插值点横坐标计算得到的近似函数值计算得到的近似函数值2023-5-52023-5-5fun=(x)sin(pi*x/2).*(x=-1&x=1|x=-1&x=1|x=-1&x=1|x x0=linspace(0,2*pi,15);y0=sin(x0)+0.3*(rand(size(x0)-0.5);pp=csaps(x0,y0,0.9

31、);sp1=spaps(x0,y0,1e-3);sp2=spap2(3,4,x0,y0);plot(x0,y0,ko);hold on fnplt(pp,2,r:)fnplt(sp1,2,k-.)fnplt(sp2,2,-)xlabel(X);ylabel(Y=sin(x);legend(插值节点插值节点,三次光滑样条插值三次光滑样条插值,光滑光滑B样条插值样条插值,最小二乘最小二乘B样条近似样条近似);01234567-1.5-1-0.500.511.5XY=sin(x)插值节点三次光滑样条插值光滑B样条插值最小二乘B样条近似2023-5-5【例例20.3-6】MATLAB创始者创始者Cle

32、ve Moler博士写的博士写的MATLAB数值计算数值计算一书中有一个有趣的例子，就是利用样条插值方法一书中有一个有趣的例子，就是利用样条插值方法描绘手的轮廓线。思路就是首先运行描绘手的轮廓线。思路就是首先运行x,y=ginput命令（命令（ginput为取点函数），然后把手放在屏幕前，用鼠标沿着手的轮廓点为取点函数），然后把手放在屏幕前，用鼠标沿着手的轮廓点击取点，不用点太密，轮廓变化大的地方（例如指尖和两指相击取点，不用点太密，轮廓变化大的地方（例如指尖和两指相连处）可多取几个点。点击完毕后按连处）可多取几个点。点击完毕后按Enter键退出取点模式，键退出取点模式，MATLAB会记录下这

33、些点的坐标。对这些点调用会记录下这些点的坐标。对这些点调用cscvn函数进行函数进行顺序节点插值，最后就可画出手的轮廓。顺序节点插值，最后就可画出手的轮廓。2023-5-5 figure(position,get(0,screensize);axes(position,0 0 1 1);x,y=ginput;curve=cscvn(x;y);plot(x,y,ko);hold on;fnplt(curve);xlabel(X);ylabel(笔者的手笔者的手);0.20.30.40.50.60.70.80.9100.20.40.60.81X笔者的手2023-5-5第四节第四节 二维插值二维插值

34、2023-5-5一、网格节点插值一、网格节点插值1.网格节点插值图解网格节点插值图解2023-5-5 前面介绍的前面介绍的csape、csapi、spapi、csaps、spaps和和spap2函数均可以作二维和高维网格节点插值，除此之外，函数均可以作二维和高维网格节点插值，除此之外，MATLAB多项式函数工具箱中还提供了多项式函数工具箱中还提供了interp2函数，用函数，用来作二维网格节点插值。来作二维网格节点插值。2.网格节点插值的网格节点插值的MATLAB函数函数2023-5-5 csape、csapi、spapi、csaps、spaps和和spap2函数用法函数用法函数名函数名调用格

35、式调用格式功能及参数说明功能及参数说明csapepp=csape(x1,x2,Y)pp=csape(x1,x2,Y,conds1,conds2)使用指定边界条件作三次样条插值csapipp=csapi(x1,x2,Y)values=csapi(x1,x2,Y,xx1,xx2)使用非纽结边界条件作三次样条插值spapisp=spapi(knork1,knork2,x1,x2,Y)指定阶次的B样条插值csapspp=csaps(x1,x2,Y,p1,p2)values=csaps(x1,x2,Y,p1,p2,xx1,xx2)三次光滑样条插值spapssp=spaps(x1,x2,Y,tol1,to

36、l2)sp,values=spaps(x1,x2,Y,tol1,tol2)光滑B样条插值spap2sp=spap2(knork1,knork2,k,x1,x2,Y)最小二乘B样条近似2023-5-5（1）ZI=interp2(X,Y,Z,XI,YI)（2）ZI=interp2(Z,XI,YI)（3）ZI=interp2(Z,ntimes)（4）ZI=interp2(X,Y,Z,XI,YI,method)节点节点Y坐标坐标节点节点X坐标坐标插值点处近似函数值插值点处近似函数值 interp2函数用法函数用法节点节点Z坐标坐标插值点插值点X坐标坐标插值点插值点Y坐标坐标指定所用的插值方法指定所用的

37、插值方法2023-5-5【例例20.4-1】在一丘陵地带测量高程，在一丘陵地带测量高程，x和和y 方向每隔方向每隔100米测一米测一个点，得高程数据如表个点，得高程数据如表20.4-3所列，试用插值方法拟合出地形曲所列，试用插值方法拟合出地形曲面。面。3.网格节点插值举例网格节点插值举例高程高程100200300400500100450478624697636200420478630712698300400412598674680400310334552626662xy2023-5-5x=100:100:500;y=100:100:400;X,Y=meshgrid(x,y);Z=450 478

38、 624 697 636 420 478 630 712 698 400 412 598 674 680 310 334 552 626 662;xd=100:20:500;yd=100:20:400;Xd1,Yd1=meshgrid(xd,yd);Xd2,Yd2=ndgrid(xd,yd);figure;%新建图形窗口新建图形窗口%-调用调用interp2函数作三次样条插值函数作三次样条插值-Zd1=interp2(X,Y,Z,Xd1,Yd1,spline);subplot(2,3,1);surf(Xd1,Yd1,Zd1);xlabel(X);ylabel(Y);zlabel(Z);titl

39、e(interp2)05000200400200400600800Xinterp2YZ05000200400200400600800XcsapeYZ05000200400200400600800XcsapiYZ05000200400200400600800XspapiYZ05000200400200400600800XcsapsYZ05000200400200400600800XspapsYZ2023-5-5二、散乱节点插值二、散乱节点插值1.散乱节点插值图解散乱节点插值图解2023-5-5 MATLAB多项式函数工具箱中提供了多项式函数工具箱中提供了TriScatteredInterp和和g

40、riddata函数，其中前者用来作二维或三维散乱节点插值，后函数，其中前者用来作二维或三维散乱节点插值，后者只能用来作二维散乱节点插值，前者比后者的运行效率高。者只能用来作二维散乱节点插值，前者比后者的运行效率高。2.散乱节点插值的散乱节点插值的MATLAB函数函数函数名调用格式功能及参数说明griddataZI=griddata(x,y,z,XI,YI)XI,YI,ZI=griddata(x,y,z,XI,YI).=griddata(.,method)返回插值点XI,YI处的近似函数值矩阵ZI。输入参数x,y是节点坐标向量，z是相应的原函数值向量。XI，YI是用户指定的插值点坐标，可以是同型

41、矩阵，也可以是向量（XI为行向量，YI为列向量）。method参数用来指定所用的插值方法，其可用取值为：linear 线性插值（默认）cubic 立方插值nearest 最近邻插值v4 MATLAB 4 中用到的插值方法TriScatteredInterpF=TriScatteredInterp(X,V)返回TriScatteredInterp类变量F。输入参数X是m行n列的矩阵，其中m为节点个数，n为2或3，表示维数，V是节点处原函数值，是长度为m的列向量。此时由VI=F(XI)计算插值点XI处的近似函数值VI。F=TriScatteredInterp(X,Y,V)X,Y,V是等长的列向量，

42、X,Y为节点坐标，V为相应的原函数值。此时由VI=F(XI,YI)计算插值点XI,YI处的近似函数值VI。F=TriScatteredInterp(X,Y,Z,V)X,Y,Z,V是等长的列向量，X,Y,Z为节点坐标，V为相应的原函数值。此时由VI=F(XI,YI,ZI)计算插值点XI,YI,ZI处的近似函数值VI。F=TriScatteredInterp(.,method)用method参数指定所用的插值方法，其可用取值为：natural 自然近邻插值linear 线性插值（默认）nearest 最近邻插值2023-5-5【例例20.4-2】2011高教社杯全国大学生数学建模竞赛高教社杯全国大

43、学生数学建模竞赛A题中给出题中给出了某城市城区土壤地质环境调查数据，包括采样点的位置、海了某城市城区土壤地质环境调查数据，包括采样点的位置、海拔高度及其所属功能区等信息数据，以及拔高度及其所属功能区等信息数据，以及8种主要重金属元素在种主要重金属元素在采样点处的浓度、采样点处的浓度、8种主要重金属元素的背景值数据。试根据调种主要重金属元素的背景值数据。试根据调查数据中给出的采样点坐标和查数据中给出的采样点坐标和8种主要重金属元素的浓度数据绘种主要重金属元素的浓度数据绘制制Cd元素的空间分布图。元素的空间分布图。3.散乱节点插值举例散乱节点插值举例2023-5-5xyz=xlsread(cumc

44、m2011A.xls,1,B4:D322);Cd=xlsread(cumcm2011A.xls,2,C4:C322);x=xyz(:,1);y=xyz(:,2);z=xyz(:,3);xd=linspace(min(x),max(x),60);yd=linspace(min(y),max(y),60);Xd,Yd=meshgrid(xd,yd);%调用调用griddata函数作散乱节点插值函数作散乱节点插值Zd1=griddata(x,y,z,Xd,Yd);Cd1=griddata(x,y,Cd,Xd,Yd);figure;surf(Xd,Yd,Zd1,Cd1);shading interp;

45、xlabel(X);ylabel(Y);zlabel(Z（griddata）);colorbar;0123x 104012x 1040100200300 XY Z（griddata）200400600800100012002023-5-5第五节第五节 高维插值高维插值2023-5-5一、高维插值函数一、高维插值函数 高维插值同样分为网格节点插值和散乱节点插值，用到的高维插值同样分为网格节点插值和散乱节点插值，用到的MATLAB函数除了函数除了csape、csapi、spapi、csaps、spaps和和spap2函数外，还有函数外，还有MATLAB多项式函数工具箱中提供的多项式函数工具箱中提供

46、的interp3（三（三维网格节点插值）、维网格节点插值）、interpn（n 维网格节点插值）、维网格节点插值）、griddata3（三维散乱节点插值）和（三维散乱节点插值）和griddatan（n 维散乱节点插值）函数。维散乱节点插值）函数。2023-5-5二、高维插值举例二、高维插值举例【例例20.5-1】用函数用函数 产生定义区域上的离散网格数据，然后作三维网格节点插值，产生定义区域上的离散网格数据，然后作三维网格节点插值，给出拟合误差的最大值。为了观察函数值在空间中的分布情况，给出拟合误差的最大值。为了观察函数值在空间中的分布情况，绘制立体切片图和等值面图。绘制立体切片图和等值面图。

47、(,)coscoscos,V x y zxyzx y zpp 调用调用interp3函数作三维网格节点插值插值函数作三维网格节点插值插值 xyz0=linspace(-pi,pi,30);X0,Y0,Z0=meshgrid(xyz0);V0=cos(X0)+cos(Y0)+cos(Z0);xyz=linspace(-pi,pi,60);X,Y,Z=meshgrid(xyz);V=cos(X)+cos(Y)+cos(Z);V1=interp3(X0,Y0,Z0,V0,X,Y,Z);err=V1-V;max(err(:)2023-5-5 调用调用slice函数绘制立体切片图函数绘制立体切片图 sl

48、ice(X,Y,Z,V,-1,1,0,0);shading interp alpha(0.5);xlabel(X);ylabel(Y);zlabel(Z);set(gca,color,none);axis(-3.5,3.5,-3.5,3.5,-3.5,3.5);colorbar;2023-5-5 调用自编函数调用自编函数MyIsosurface绘制等值面图绘制等值面图figure;subplot(2,3,1);MyIsosurface(X,Y,Z,V,-1.2);title(V(x,y,z)=-1.2);subplot(2,3,2);MyIsosurface(X,Y,Z,V,-1);title(V(x,y,z)=-1);subplot(2,3,3);MyIsosurface(X,Y,Z,V,-0.9);title(V(x,y,z)=-0.9);subplot(2,3,4);MyIsosurface(X,Y,Z,V,0);title(V(x,y,z)=0);.