例:机床加工
待加工零件的外形根据工艺要求由一组数据(x,y)给出(平面情况),用程控铣床加工时每一刀只能沿x和y方向走一小步,这就需要从已知数据得到加工所要求的步长很小的(x,y)坐标。
下表给出x,y数据位于机翼断面的下轮廓线上,假设需要得到x坐标每改变0.1的y坐标,试完成加工所需数据线,画出曲线,并求出x=0处的曲线斜率和13x
15范围内的y的最小值,要求用分段线性和三次样条两种插值方法计算
x | 0 | 3 | 5 | 7 | 9 | 11 | 12 | 13 | 14 | 15 |
y | 0 | 1.2 | 1.7 | 2.0 | 2.1 | 2.0 | 1.8 | 1.2 | 1.0 | 1.6 |
clc,clear;
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;
y1 = interp1(x0,y0,x); %默认一维线性插值
y2 = interp1(x0,y0,x,'spline'); %立方样条插值
pp1 = csape(x0,y0); %三次样条插值返回值pp形式
y3 = fnval(pp1,x); %三次样条插值要求x处的函数值必须使用函数fnval
pp2 = csape(x0,y0,'second'); %边界为二阶导数,默认值[0,0]
y4 = fnval(pp2,x);
[x',y1',y2',y3',y4'] %转置 行变列输出
subplot(1,3,1);
plot(x0,y0,'+',x,y1)
title('分段线性插值');
subplot(1,3,2);
plot(x0,y0,'+',x,y2);
title('立方样条插值');
subplot(1,3,3);
plot(x0,y0,'+',x,y3);
title('三次样条插值');
dx = diff(x);
dy = diff(y3);
dy_dx = dy./dx;
dy_dx0 = dy_dx(1);
ytemp = y3(131:151);
ymin = min(ytemp);
index = find(y3 == ymin);
xmin = x(index);
[xmin,ymin]
运行结果
可以看出,分段线性插值的光滑性较差(特别是在x=14附近弯曲处),建议选用三次样条插值的结果)
ans =
0 0 0 0 0
0.1000 0.0400 0.0499 0.0497 0.0441
0.2000 0.0800 0.0990 0.0987 0.0881
0.3000 0.1200 0.1474 0.1470 0.1321
0.4000 0.1600 0.1951 0.1946 0.1760
0.5000 0.2000 0.2421 0.2415 0.2198
0.6000 0.2400 0.2884 0.2878 0.2635
0.7000 0.2800 0.3340 0.3333 0.3070
0.8000 0.3200 0.3788 0.3781 0.3503
0.9000 0.3600 0.4230 0.4223 0.3934
1.0000 0.4000 0.4665 0.4658 0.4362
1.1000 0.4400 0.5094 0.5086 0.4788
1.2000 0.4800 0.5515 0.5507 0.5211
1.3000 0.5200 0.5930 0.5922 0.5631
1.4000 0.5600 0.6338 0.6331 0.6046
1.5000 0.6000 0.6739 0.6732 0.6459
1.6000 0.6400 0.7134 0.7128 0.6867
1.7000 0.6800 0.7523 0.7516 0.7271
1.8000 0.7200 0.7904 0.7899 0.7670
1.9000 0.7600 0.8280 0.8275 0.8064
2.0000 0.8000 0.8649 0.8644 0.8453
2.1000 0.8400 0.9012 0.9007 0.8837
2.2000 0.8800 0.9368 0.9364 0.9215
2.3000 0.9200 0.9719 0.9715 0.9587
2.4000 0.9600 1.0063 1.0060 0.9952
2.5000 1.0000 1.0401 1.0398 1.0311
2.6000 1.0400 1.0732 1.0731 1.0664
2.7000 1.0800 1.1058 1.1057 1.1009
2.8000 1.1200 1.1378 1.1377 1.1347
2.9000 1.1600 1.1692 1.1692 1.1678
3.0000 1.2000 1.2000 1.2000 1.2000
3.1000 1.2250 1.2302 1.2302 1.2314
3.2000 1.2500 1.2599 1.2599 1.2621
3.3000 1.2750 1.2889 1.2890 1.2920
3.4000 1.3000 1.3174 1.3175 1.3211
3.5000 1.3250 1.3454 1.3455 1.3494
3.6000 1.3500 1.3727 1.3728 1.3771
3.7000 1.3750 1.3995 1.3997 1.4041
3.8000 1.4000 1.4258 1.4259 1.4303
3.9000 1.4250 1.4515 1.4517 1.4560
4.0000 1.4500 1.4767 1.4768 1.4810
4.1000 1.4750 1.5014 1.5015 1.5053
4.2000 1.5000 1.5255 1.5256 1.5291
4.3000 1.5250 1.5491 1.5492 1.5523
4.4000 1.5500 1.5722 1.5723 1.5749
4.5000 1.5750 1.5947 1.5948 1.5970
4.6000 1.6000 1.6168 1.6168 1.6186
4.7000 1.6250 1.6383 1.6384 1.6396
4.8000 1.6500 1.6594 1.6594 1.6602
4.9000 1.6750 1.6799 1.6799 1.6803
5.0000 1.7000 1.7000 1.7000 1.7000
5.1000 1.7150 1.7196 1.7196 1.7193
5.2000 1.7300 1.7387 1.7386 1.7381
5.3000 1.7450 1.7573 1.7572 1.7565
5.4000 1.7600 1.7754 1.7754 1.7744
5.5000 1.7750 1.7930 1.7930 1.7920
5.6000 1.7900 1.8102 1.8101 1.8090
5.7000 1.8050 1.8269 1.8268 1.8257
5.8000 1.8200 1.8430 1.8430 1.8419
5.9000 1.8350 1.8588 1.8587 1.8576
6.0000 1.8500 1.8740 1.8740 1.8729
6.1000 1.8650 1.8887 1.8887 1.8877
6.2000 1.8800 1.9030 1.9030 1.9021
6.3000 1.8950 1.9168 1.9168 1.9160
6.4000 1.9100 1.9301 1.9301 1.9294
6.5000 1.9250 1.9430 1.9429 1.9424
6.6000 1.9400 1.9553 1.9553 1.9549
6.7000 1.9550 1.9672 1.9672 1.9669
6.8000 1.9700 1.9786 1.9786 1.9784
6.9000 1.9850 1.9895 1.9895 1.9895
7.0000 2.0000 2.0000 2.0000 2.0000
7.1000 2.0050 2.0100 2.0100 2.0101
7.2000 2.0100 2.0195 2.0195 2.0196
7.3000 2.0150 2.0285 2.0285 2.0287
7.4000 2.0200 2.0370 2.0370 2.0372
7.5000 2.0250 2.0450 2.0451 2.0452
7.6000 2.0300 2.0525 2.0526 2.0527
7.7000 2.0350 2.0595 2.0596 2.0597
7.8000 2.0400 2.0660 2.0660 2.0662
7.9000 2.0450 2.0719 2.0720 2.0721
8.0000 2.0500 2.0773 2.0774 2.0774
8.1000 2.0550 2.0822 2.0822 2.0823
8.2000 2.0600 2.0865 2.0865 2.0865
8.3000 2.0650 2.0902 2.0902 2.0902
8.4000 2.0700 2.0933 2.0934 2.0934
8.5000 2.0750 2.0959 2.0960 2.0959
8.6000 2.0800 2.0979 2.0980 2.0979
8.7000 2.0850 2.0994 2.0994 2.0993
8.8000 2.0900 2.1002 2.1002 2.1001
8.9000 2.0950 2.1004 2.1004 2.1004
9.0000 2.1000 2.1000 2.1000 2.1000
9.1000 2.0950 2.0990 2.0990 2.0990
9.2000 2.0900 2.0974 2.0974 2.0975
9.3000 2.0850 2.0952 2.0952 2.0954
9.4000 2.0800 2.0925 2.0925 2.0927
9.5000 2.0750 2.0893 2.0893 2.0896
9.6000 2.0700 2.0857 2.0856 2.0860
9.7000 2.0650 2.0815 2.0814 2.0819
9.8000 2.0600 2.0770 2.0769 2.0775
9.9000 2.0550 2.0721 2.0719 2.0726
10.0000 2.0500 2.0668 2.0666 2.0673
10.1000 2.0450 2.0611 2.0610 2.0617
10.2000 2.0400 2.0552 2.0550 2.0558
10.3000 2.0350 2.0490 2.0488 2.0496
10.4000 2.0300 2.0425 2.0423 2.0431
10.5000 2.0250 2.0358 2.0357 2.0364
10.6000 2.0200 2.0289 2.0288 2.0294
10.7000 2.0150 2.0219 2.0218 2.0223
10.8000 2.0100 2.0147 2.0146 2.0150
10.9000 2.0050 2.0074 2.0073 2.0076
11.0000 2.0000 2.0000 2.0000 2.0000
11.1000 1.9800 1.9924 1.9925 1.9922
11.2000 1.9600 1.9841 1.9842 1.9837
11.3000 1.9400 1.9742 1.9744 1.9736
11.4000 1.9200 1.9621 1.9622 1.9612
11.5000 1.9000 1.9469 1.9471 1.9459
11.6000 1.8800 1.9280 1.9282 1.9270
11.7000 1.8600 1.9046 1.9048 1.9036
11.8000 1.8400 1.8759 1.8761 1.8751
11.9000 1.8200 1.8413 1.8414 1.8408
12.0000 1.8000 1.8000 1.8000 1.8000
12.1000 1.7400 1.7516 1.7514 1.7522
12.2000 1.6800 1.6970 1.6967 1.6984
12.3000 1.6200 1.6377 1.6372 1.6398
12.4000 1.5600 1.5749 1.5743 1.5777
12.5000 1.5000 1.5099 1.5092 1.5133
12.6000 1.4400 1.4442 1.4435 1.4478
12.7000 1.3800 1.3790 1.3783 1.3825
12.8000 1.3200 1.3157 1.3151 1.3186
12.9000 1.2600 1.2556 1.2552 1.2573
13.0000 1.2000 1.2000 1.2000 1.2000
13.1000 1.1800 1.1501 1.1506 1.1477
13.2000 1.1600 1.1063 1.1074 1.1011
13.3000 1.1400 1.0687 1.0704 1.0607
13.4000 1.1200 1.0377 1.0399 1.0271
13.5000 1.1000 1.0134 1.0160 1.0009
13.6000 1.0800 0.9960 0.9988 0.9827
13.7000 1.0600 0.9857 0.9885 0.9729
13.8000 1.0400 0.9828 0.9851 0.9721
13.9000 1.0200 0.9875 0.9889 0.9810
14.0000 1.0000 1.0000 1.0000 1.0000
14.1000 1.0600 1.0205 1.0186 1.0295
14.2000 1.1200 1.0492 1.0450 1.0686
14.3000 1.1800 1.0863 1.0799 1.1162
14.4000 1.2400 1.1320 1.1237 1.1714
14.5000 1.3000 1.1866 1.1768 1.2330
14.6000 1.3600 1.2503 1.2398 1.3000
14.7000 1.4200 1.3233 1.3131 1.3712
14.8000 1.4800 1.4057 1.3973 1.4457
14.9000 1.5400 1.4979 1.4927 1.5223
15.0000 1.6000 1.6000 1.6000 1.6000
ans =
13.8000 0.9851
例2:
已知速度曲线v(t)上的四个数据点如表所列
速度的4个观测值
t | 0.15 | 0.16 | 0.17 | 0.18 |
v(t) | 3.5 | 1.5 | 2.5 | 2.8 |
用三次样条插值求位移
clc,clear;
x0 = [0.15 0.16 0.17 0.18];
y0 = [3.5 1.5 2.5 2.8];
%三次样条插值
pp = csape(x0,y0);
format long g
xishu = pp.coefs %显示每个区间上三次多项式的系数
s = quadl(@(t)ppval(pp,t),0.15,0.18) %求积分
format % 恢复短小数的显示格式
运行结果
xishu =-616666.666666666 33500 -473.333333333333 3.5
-616666.666666666 15000 11.6666666666667 1.5
-616666.666666666 -3500.00000000001 126.666666666666 2.5
s =0.068625
二维插值
例:在一丘陵地带测量高度,x和y方向每隔100m测一个点,试插值一曲面,确定合适的模型,并由此找出最高点和该点的x
x y | 100 | 200 | 300 | 400 | 500 |
100 | 636 | 697 | 624 | 478 | 450 |
200 | 698 | 712 | 630 | 478 | 420 |
300 | 680 | 674 | 598 | 412 | 400 |
400 | 662 | 626 | 552 | 334 | 310 |
程序如下
clc,clear;
x = 100:100:500;
y = 100:100:400;
z = [636 697 624 478 450
698 712 630 478 420
680 674 598 412 400
662 626 552 334 310];
pp = csape({x,y},z')
xi = 100:10:500;
yi = 100:10:400;
cz = fnval(pp,{xi,yi});
[i,j] = find(cz == max(max(cz))); %寻找最高点坐标
x = xi(i),y = yi(j)
zmax = cz(i,j) %寻找最高点高度
运行结果:
x =
170
y =180
zmax =720.6252
即在坐标(170,180)最高,高720.6252
例四
在某海域测得一些点(x,y)的水深z由表给出,在适当的矩形区域内画出海底曲面的图形
x | 129 | 140 | 103.5 | 88 | 185.5 | 195 | 105 | 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 |
clc,clear;
x=[129,140,103.5,88,185.5,195,105,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];
xmm = minmax(x) %求x的最小值和最大值
ymm = minmax(y) %求y的最小值和最大值
xi = xmm(1):xmm(2);
yi = ymm(1):ymm(2);
zi1 = griddata(x,y,z,xi,yi','cubic'); %立方插值
zi2 = griddata(x,y,z,xi,yi','nearest'); %最近邻插值
zi = zi1; %立方插值和最近点插值的混合插值的初始值
zi(isnan(zi1))=zi2(isnan(zi1)) %把立方插值中的不确定值换成最近点插值的结果
subplot(1,2,1);
plot(x,y,'*')
subplot(1,2,2)
mesh(xi,yi,zi)
matla插值时外插值是不确定的,这里使用了混合插值,把不确定的插值换成了最近点插值的结果
运行结果