实验四插值和拟合
插值和拟合
实验目的：了解数值分析建模的方法，掌握用Matlab进行曲线拟合的方法，理解用插值法建模的思想，运用Matlab一些命令及编程实现插值建模。
实验要求：理解曲线拟合和插值方法的思想，熟悉Matlab相关的命令，完成相应的练习，并将操作过程、程序及结果记录下来。 实验内容： 一、插值
1．插值的基本思想
·已知有n +1个节点（xj,yj），j = 0,1,…, n，其中xj互不相同，节点(xj, yj)可看成由某个函数 y= f（x）产生；
·构造一个相对简单的函数 y=P(x)；
·使P通过全部节点，即 P (xk) = yk，k=0,1,…, n ； ·用P (x)作为函数f ( x )的近似。 2．用MATLAB作一维插值计算 yi=interp1(x，y，xi，'method')
注：yi—xi处的插值结果；x，y—插值节点；xi—被插值点；method—插值方法（‘nearest’ ：最邻近插值；‘linear’ ：线性插值；‘spline’ ：三次样条插值；‘cubic’ ：立方插值；缺省时：线性插值）。 注意：所有的插值方法都要求x是单调的，并且xi不能够超过x的范围。 练习1：机床加工问题
机翼断面下的轮廓线上的数据如下表： x y 0 0 3 1.2 5 1.7 7 2.0 9 2.1 11 2.0 12 1.8 13 1.2 14 1.0 15 1.6 用程控铣床加工机翼断面的下轮廓线时 每一刀只能沿x方向和y方向走非常小的一步。 表3-1给出了下轮廓线上的部分数据
但工艺要求铣床沿x方向每次只能移动0.1单位. 这时需求出当x坐标每改变0.1单位时的y坐标。 试完成加工所需的数据，画出曲线. 步骤1：用x0，y0两向量表示插值节点；
步骤2：被插值点x=0:0.1:15; y=interp1(x0,y0,x,'spline'); 步骤3：plot(x0,y0,'k+',x,y,'r')
grid on
>> x0=[0 3 5 7 9 11 12 13 14 15 ];
>> y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6 ];
>> x=0:0.1:15;y=interp1(x0,y0,x,'spline');plot(x0,y0,'k+',x,y,'r') grid on
1
实验四插值和拟合
2.521.510.50051015
3．用MATLAB作网格节点数据的插值(二维) z=interp2(x0,y0,z0,x,y,’method’) 注：z—被插点值的函数值；x0，y0，z0—插值节点；x，y—被插值点；method—插值方法（‘nearest’ ：最邻近插值；‘linear’ ：双线性插值； ‘cubic’ ：双三次插值；缺省时：双线性插值）。 注意：要求x0,y0单调；x，y可取为矩阵，或x取行向量，y取为列向量，x,y的值分别不能超出x0,y0的范围。
4．用MATLAB作散点数据的插值计算 cz =griddata（x，y，z，cx，cy，‘method’）
注：cz—被插点值的函数值；x，y，z—插值节点；cx，cy—被插值点；method—插值方法（‘nearest’ ：最邻近插值；‘linear’ ：双线性插值； ‘cubic’ ：双三次插值；'v4‘：Matlab提供的插值方法；缺省时：双线性插值）。
练习2：航行区域的警示线
某海域上频繁地有各种吨位的船只经过。 为保证船只的航行安全，有关机构在低潮时对水深进行了测量，下表是他们提供的测量数据： 水道水深的测量数据
x 129.0 140.0 103.5 88.0 185.5 195.0 105.5 y 7.5 141.5 23.0 147.0 22.5 137.5 85.5 z 4 8 6 8 6 8 8
x 157.5 107.5 77.0 81.0 162.0 162.0 117.5 y -6.5 -81.0 3.0 56.5 -66.5 84.0 -33.5 z 9 9 8 8 9 4 9
其中（x, y）为测量点，z为（x, y）处的水深(英尺)，水深z是区域坐标（x, y）的函数z= z (x, y），
船的吨位可以用其吃水深度来反映，分为 4英尺、4.5英尺、5英尺和 5.5英尺 4 档。
2
实验四插值和拟合
航运部门要在矩形海域（75，200）×（－50，150）上为不同吨位的航船设置警示标记。 请根据测量的数据描述该海域的地貌，并绘制不同吨位的警示线，供航运部门使用。 x=[129 140 103.5 88 185.5 195 105.5 157.5 107.5 77 81 162 162 117.5]; y=[7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5]; z=[-4 -8 -6 -8 -6 -8 -8 -9 -9 -8 -8 -9 -4 -9];
cx=75:0.5:200; cy=-70:0.5:150;
cz=griddata(x,y,z,cx,cy','cubic');
meshz(cx,cy,cz),rotate3d
xlabel('X'),ylabel('Y'),zlabel('Z') %pause
figure(2),contour(cx,cy,cz,[-5 -5]);grid on, hold on plot(x,y,'+')
xlabel('X'),ylabel('Y')
140120100806040200-20-40-6080100120140X160180200Y
练习3：估计水塔的水流量
问题描述见教材P.91—93，请绘出三次样条插值曲线，并计算一天的总的用水量。
x0=[0.46 1.38 2.40 3.41 4.43 5.44 6.45 7.47 8.45 11.49 12.49 13.42 14.43 15.44 16.37 17.38 18.49 19.50 20.40 24.43 25.32];
y0=[11.20 9.7 8.6 8.1 9.3 7.2 7.9 7.4 8.4 15.6 16.4 15.5 13.4 13.8 12.9 12.2 12.2 12.9 12.6 12.2 3.5]; x=0:0.01:25.32;y=interp1(x0,y0,x,'spline');plot(x0,y0,'k+',x,y,'r');grid on
3
实验四插值和拟合
18161412108642051015202530
二、曲线拟合
已知一组（二维）数据，即平面上 n个点（xi,yi) i=1,…n, 寻求一个函数（曲线）y=f(x), 使 f(x) 在某种准则下与所有数据点最为接近，即曲线拟合得最好。
最常用的方法是线性最小二乘拟合 1．多项式拟合
? 对给定的数据（xj,yj），j = 0,1,…, n；
? 选取适当阶数的多项式，如二次多项式g(x)=ax^2+bx+c；
? 使g(x)尽可能逼近（拟合）这些数据，但是不要求经过给定的数据（xj,yj）； 2．多项式拟合指令
1）多项式f(x)=a1xm+ …+amx+am+1拟合指令:
a=polyfit(x,y,m)
a：输出多项式拟合系数a[a1,a2,…,am]；x,y：输出长度相同的数组；m：多项式的次数。 2）多项式在x处的值y的计算命令：
y=polyval（a,x）
练习4：对下面一组数据作二次多项式拟合 xi 0 0.1 0.2 0.3 0.4 6.16 0.5 7.34 0.6 7.66 0.7 9.58 0.8 9.48 0.9 9.30 1 11.2 yi -0.447 1.978 3.28 6.16 写出拟合命令： 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,'k+',x,z,'r')
4
实验四插值和拟合
作出数据点和拟合曲线：
121086420-200.10.20.30.40.50.60.70.80.91
3．可化为多项式的非线性拟和
曲线改直是工程中又一常用的判断曲线形式的方法,许多常见的函数都可以通过适当的变换转化为线性函数。
（1）幂函数 y?axb?c
lny?c?lna?blnx
（2）指数函数 y?ab?c
lny?c?lna?xlnb
（3）抛物函数 y?ax?bx?,(c 2xx?0 )y?c?ax?b x
三 实验内容
1 已知x=[1.2,1.8,2.1,2.4,2.6,3.0,3.3],y=[4.85,5.2,5.6,6.2,6.5,7.0,7.5],求对x和y进行6阶多项式拟合的系数.
x=[1.2,1.8,2.1,2.4,2.6,3.0,3.3]; >> y=[4.85,5.2,5.6,6.2,6.5,7.0,7.5]; >> polyfit(x,y,6)
2 分别用2,3,4,5阶多项式来逼近[0,3]上的正弦函数sin x,并做出拟合曲线及sin x函数曲线图,了解多项式的逼近程度和有效拟合区间随多项式的阶数有何变化. x=(0:0.1*pi:2*pi);
y=[0,0.3090,0.5878,0.8090,0.9511,1.0000,0.9511,0.8090,0.5878,0.3090,0.0000,-0.3090,-0.5878
5
实验四插值和拟合
,-0.8090,-0.9511,-1.0000,-0.9511,-0.8090,-0.5878,-0.3090,-0.0000]; y2=polyfit(x,y,2);z2=polyval(y2,x); y3=polyfit(x,y,3);z3=polyval(y3,x); y4=polyfit(x,y,4);z4=polyval(y4,x); y5=polyfit(x,y,5);z5=polyval(y5,x);
plot(x,y,'k+',x,z2,'r*',x,z3,'g:',x,z4,'b-',x,z5,'y.')
3 已知x=[0.1,0.8,1.3,1.9,2.5,3.1],y=[1.2,1.6,2.7,2.0,1.3,0.5],用不同的方法求x=2点的插值,并分析所得结果有何不同.
x0=[0.1,0.8,1.3,1.9,2.5,3.1]; y0=[1.2,1.6,2.7,2.0,1.3,0.5]; yi0=interp1(x0,y0,2,'nearest') yi1=interp1(x0,y0,2,'spline') yi2=interp1(x0,y0,2,'linear') 4 已知热敏电阻数据
温度t(0C) 20.5 32.7 51.0 73.0 95.7 电阻R(?) 765 826 873 942 1032 求60C时的电阻R。（用拟合方法） t=[20.5,32.7,51.0,73.0,95.7]; r=[765,826,873,942,1032]; A=polyfit(t,r,2); r=polyval(A,60) 5 对函数f(x)?1，x?[-5,5]，分别用分段线性插值和三次样条插值作插值（其中插值节点1?x2不少于20），并分别作出每种插值方法的误差曲线.
x0=[-5,-4.5,-4,-3.5,-3,-2.5,-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5];
y0=[0.0385,0.0471,0.0588,0.0755,0.1000,0.1379,0.2000,0.3077,0.5000,0.8000,1.0000,0.8000,
0.5000,0.3077,0.2000,0.1379,0.1000,0.0755,0.0588,0.0471,0.0385];
x=-5:0.5:5;y=interp1(x0,y0,x,'spline');subplot(1,2,1);plot(x0,y0,'m+',x,y,'g'); x=-5:0.5:5;y=interp1(x0,y0,x,'linear');subplot(1,2,2);plot(x0,y0,'k+',x,y,'r')
6