数学建模插值法

1、.,插值法,.,问题,一水库上游河段降暴雨，根据预报测算上游流入水库的流量为Q(t) (102立方米/秒) ： t (时) 8 12 16 24 30 44 48 56 60 Q（t） 36 54 78 92 101 35 25 16 13 通过这个预报值，估计14 和 20 时上游流入水库的流量.,在一天24小时内, 从零点开始每间隔2小时测得的环境温度为(摄氏度) 12,9,9,10,18,24,28,27,25,20,18,15,13 推测在每一分钟时的温度.,.,给定数据表：,计算ln(0.54)的值.,.,一、问题的提出,在生产和实验中，关于函数f(x)，经常存在两种情况： （1）其

2、表达式不便于计算； （2）无表达式， 而只有函数在给定点的函数值， 怎样计算其它点的函数值？,.,二、问题的解决,两种方法： （1）函数插值； （2）曲线拟合,.,三、插值法,定义：当精确函数 y = f(x) 非常复杂或未知时，在一系列节点 x0 xn 处测得函数值 y0 = f(x0), ，yn = f(xn)， 由此构造一个简单易算的近似函数 p(x) f(x)，满足条件 p(xi) = f(xi) (i = 0, n)- 插值条件 这里的 p(x) 称为f(x) 的插值函数； 构造插值函数的方法为插值法。,当 p(x)为多项式时，我们称p(x)为插值多项式；构造插值多项式的方法称为多项

3、式插值法.,.,多项式插值法包含：,拉格朗日插值 牛顿插值 三次埃尔米特插值法 分段线性插值 分段三次埃尔米特插值法 三次样条插值,.,1、 拉格朗日插值公式,（）定义 对给定的n+1个节点x0 , x1,x2,xn及对应的函数值y0 , y1,y2,yn， 构造一个n次插值多项式： 即为拉格朗日插值公式，其中,插值基函数,.,（2）拉格朗日插值的matlab实现,function y=lagrange(x0,y0,x) % x0插值节点， y0插值节点处的函数值，x要计算函数值的点； n=length(x0); %计算x0的长度 m=length(x); %计算x的长度 for i=1:m

4、s=0; z=x(i);,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; %计算在x(i)处的函数值（拉格朗日） end,.,（3） 函数的调用,给定数据表： 用lagrange插值公式计算 clear x0=0.4:0.1:0.8； y0=log(x0)； x=0.54； lagrange(x0,y0,x),ln(0.54)的函数值 ans = -0.6161,.,（4）练习,用lagrange插值公式计算x=2.101,4.234处的

5、函数值。,.,2、牛顿插值法,牛顿插值公式： Nn(x)=f(x0)+fx0,x1(x-x0)+fx0,x1,x2,xn(x-x0)(x-x1)(x-xn) 其中： fx0,x1 一阶差商 fx0,x1,x2,xn n阶差商 注：牛顿插值法与拉格朗日插值法，同一个多项式，不同的表达方式，但是计算量不一样，牛顿插值法的计算量小。 这两种方法在MATLAB 中没有现成的命令和函数，因为它们存在龙格现象。,.,牛顿插值法的matlab实现,function Y=newpoly(x,y,X) %牛顿插值 n=length(x); D=zeros(n); D(:,1)=y; for j=2:n %计算差

6、商矩阵 for k=j:n D(k,j)=(D(k,j-1)-D(k-1,j-1)/(x(k)-x(k-j+1); end end C=D(n,n); for k=(n-1):-1:1 %构造插值多项式 C=conv(C,poly(x(k); m=length(C); C(m)=C(m)+D(k,k); end Y=polyval(C,X);,.,龙格现象,Runge在上个世纪初发现： 在-5,5上用n+1个等距节点作n次插值多项式Pn(x)， 当在n时，插值多项式Pn(x)在区间中部趋于 f(x)=1/(1+x2) , 但对于3.63x5的x，Pn(x)严重发散。 用图形分析问题。,.,fo

7、r n=10:2:20 %从10等份到20等份 x0=-5:10/n:5; %插值节点 y0=1./(1+x0.2); %插值节点处的精确函数值 x=-5:0.1:5; %要进行计算函数值的点 y=lagrange(x0,y0,x); %调用函数计算x点的函数值 plot(x0,y0,*,x,1./(1+x.2),r,x,y) %绘制图形 pause %等待，按任意键 end,.,for n=4:2:10 %从4等份到10等份 x0=-5:10/n:5; %插值节点 y0=1./(1+x0.2); %插值节点处的精确函数值 x=-5:0.1:5; %要进行计算函数值的点 y=lagrange(

8、x0,y0,x); %调用函数计算x点的函数值 plot(x0,y0,b-*,x,1./(1+x.2),r,x,y,k) %绘制图形 pause %等待，按任意键 end,.,3、分段低次插值法,（1）分段线性插值 定义： 已知n+1个不同节点x0,x1,xn ,构造分段多项式I(x)，使之满足 lI(x)在a,b上连续； l I(xk)=yk； lI(x)在xi,xi+1上是一次多项式； I(x)=,.,（2）分段三次埃尔米特插值法 定义： 已知n+1个不同节点x0,x1,xn ,构造分段多项式I(x)，使之满足： l I(xk)=yk， I(xk)=yk， ； lI(x)在xi,xi+1上

9、是三次次多项式。,.,4、三次样条插值法,对于给定n+1个不同节点x0,x1,xn及函数值y0,y1,yn， 其中a=x0x1xn=b，构造三次样条插值函数S(x)。 S(x)称为三次样条函数时需满足： l S(x)在a,b上二阶导数连续； l S(xk)=yk (k=0,1,n)； l每个子区间xk,xk+1上S(x)是三次多项式(k=0,1,n)。,.,四、插值法的matlab实现,命令：interp1(x0,y0,x，method) 其中：x0：插值节点； y0：插值节点处的函数值； x：要计算函数值的点； method： l i n e a r ：分段线性插值； c u b i c ：

10、分段三次埃尔米特插值； s p l i n e ：三次样条插值。,.,五、插值法的应用,1、一水库上游河段降暴雨，根据预报测算上游流入水库的流量为Q(t) (102立方米/秒) ： t (时) 8 12 16 24 30 44 48 56 60 Q（t） 36 54 78 92 101 35 25 16 13 通过这个预报值，分别用不同的数值方法插值法来估计14 和20时上游流入水库的流量，并利用不同的插值方法描绘水流量曲线。,.,程序设计： clear x0=8,12,16,24,30,44,48,56,60; y0=36,54,78,92,101,35,25,16,13; x=14,20;

11、 xx=8:0.1:60; yl=interp1(x0,y0,x,linear) %分段线性插值法 yc=interp1(x0,y0,x,cubic) %分段三次埃尔米特插值法 ys=interp1(x0,y0,x,spline) %三次样条插值法 xyl=interp1(x0,y0,xx,linear); xyc=interp1(x0,y0,xx,cubic); xys=interp1(x0,y0,xx,spline); plot(x0,y0,*,xx,xyl,r,xx,xyc,k,xx,xys,g),.,2、在一天24小时内,从零点开始每间隔2小时测得的环境温度为(摄氏度) 12,9,9,

12、10,18,24,28,27,25,20,18,15,13推测在每一分钟时的温度.并利用不同的插值方法描绘温度曲线。,.,3、 下表给出的x、y数据位于机翼端面的轮廓线上，Y1和Y2分别对应轮廓的上下线。假设需要得到x坐标每改变0.1时的y坐标，试完成加工所需数据，画出曲线.,.,二维插值,前面的命令interp1()是一维插值，现实中有很多问题需要二维插值. 二维插值是基于与一维插值同样的基本思想.,.,二维插值的MATLAB实现,在MATLAB中，二维插值命令常用的有两个， 1、一个是网格节点插值： z=interp2(x0,y0,z0,x,y,method) 其中， z：被插值点处的函数

13、值； x0,y0,z0：插值节点， x0,y0为向量，z0是矩阵，其列数等于x0的长度，行数等于y0的长度； x,y：被插值点；,interp1(x0,y0,x，method),.,山体地貌,要在某山区方圆大约27平方公里范围内修建一条公路，从山脚出发经过一个居民区，再到达一个矿区。横向纵向分别每隔400米测量一次，得到一些地点的高程：,试做出该山区的地貌图，并对几种插值法进行比较.,.,程序设计： clear x0=1200:400:4000; y0=1200:400:3600; z0=1130 1250 1280 1230 1040 900 500 700; 1320 1450 1420

14、1400 1300 700 900 850; 1390 1500 1500 1400 900 1100 1060 950; 1500 1200 1100 1350 1450 1200 1150 1010; 1500 1200 1100 1550 1600 1550 1180 1070; 1500 1550 1600 1550 1600 1600 1600 1550; 1480 1500 1550 1510 1430 1300 1200 980; xi=1200:1:4000; %加密数据点 yi=1200:1:3600; zil=interp2(x0,y0,z0,xi,yi,linear);

15、%线性插值 zic=interp2(x0,y0,z0,xi,yi,cubic); %三次插值 zis=interp2(x0,y0,z0,xi,yi,spline); %样条插值 subplot(2,2,1) mesh(x0,y0,z0) subplot(2,2,2) mesh(xi,yi,zil) subplot(2,2,3) mesh(xi,yi,zic) subplot(2,2,4) mesh(xi,yi,zis),.,.,2、另一个是离散数据节点的插值命令： z=griddata(x0,y0,z0,x,y,method) 其中， z：被插值点处的函数值； x0,y0,z0：插值节点， x

16、0,y0,z0均为向量； x,y：被插值点； method：插值方法，包括： linear线性插值； cubic 三次插值；,.,船在该海域会搁浅吗？,在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）*（-50，150）里的哪些地方船要避免进入.,.,解决问题的步骤： 1.作出测量点的分布图； 2.求出矩形区域（75，200）*（-50，150）的细分网格节点之横、纵坐标向量； 3.利用MATLAB中的散点插值函数求网格节点的水深； 4.作出海底曲面图形和等高线图； 5.作出水深小于5的海域范围.,.,程序,clear x0=129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5; y0=7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5; z0=-4 -8 -6 -8 -6 -8 -8 -9 -9 -8 -8 -9 -4 -9; subplot(2,2,1)