将特殊段分割出去，对分割出的数据进行三次样条插值
（1），再用复化梯形公式
（2）求出特殊段s11,s12,s21,s22的值。
由于这四个特殊段面积的计算方法类同，所以这里就以s11的计算为例。
s11的计算：
A2=xlsread（'第一题数据','下边疆','A11:
B13'）;
x2=A2（:
1）;y2=A2（:
2）;
x21=max（x2）:
-0.001:
min（x2）;
y2i1=interp1（x2,y2,x21,'spline'）;
s11=-trapz（x21,y2i1）;%由于x21的值是降序的所以符号取反的结果才是面积
s11
计算结果为：
s11=2.7150e+003
s12=1.3585e+003
s21=3.8407e+003
s22=2.2610e+003
3．计算s1和s2
s1中包含s11和s12的图
s2中包含s21和s22的图
由于s1和s2的计算方法类同，所以这就以s1的计算为例。
对分段的数据进行三次样条插值（除特殊段以外，具体分割见附录1下边疆），每段都运用复化梯形公式并求出和记为s10。
观察图形可得出，s10与s1相比就只多算了两个特殊段的面积s11，s12，所以s1=s10-s11-s12。
计算s1的关键代码：
s10=trapz（x10,y1i）+trapz（x30,y3i）+trapz（x40,y4i）+trapz（x60,y6i）;
%用复化梯形公式求除特殊段以外的各段与X轴围成的面积并求和
计算结果为：
s1=7.4597e+004
s2=2.1257e+005
4．计算S
从以上的过程已经得出s=s2-s1，还要将s按比例转化所得就是国土面积S。
s=s2-s1；
S=s*100/9；
S
计算结果为：
S=1.5331e+006（平方公里）
5.计算相对误差t
相对误差的求解公式为：
计算结果为：
t=0.0213
6．误差分析
误差的产生主要来源于数据点的个数有点少以及数据点之间不够均匀，这样三次样条插值后的数据作图就会与实际地图和地形相比有较大的误差。
7．分析和总结
由于t的值较小，所以以上的计算结果S=1.5331e+006平方公里可以作为国土面积的近似值。
在以后的不规则图形面积的计算中，此方法可以考虑选择使用。
由于梯形公式只有一次代数精确度（3），所以会产生计算的结果不够精确。
但要很精确求出不规则图形面积，就必须测量出更多的数据点以及选用代数精确度更高的算法，才能够更加减少插值和计算上的误差。
至此，国土面积计算完成。
注释：
（1）三次样条插值：
参照参考文献[2]46~50页
（2）复化梯形公式：
参照参考文献[2]90~91页
（3）代数精确度：
参照参考文献[2]88页
参考文献：
[1]玉莉等，MATLAB函数速查手册，：
化学工业，2010
[2]袁东锦，计算方法—数值分析，：
师大学，2007
[3]蒲俊吉家锋伊良忠，Matlab6.0数学手册，浦东：
浦东电子，2002
附录：
1．对附件数据的分段
下边疆:
A03:
B10;
A11:
B13;%（特殊段s11）
A13:
B31;
A32:
B39;
A39:
B40;%（特殊段s12）
A40:
B56;
上边疆:
A03:
B11;
A11:
B14;%（特殊段s21）
A14:
B40;
A40:
B41;%（特殊段s22）
A41:
B49;
A50:
B52;
A53:
B64;
2．求国土面积及相对误差的完整代码
A1=xlsread（'第一题数据','下边疆','A03:
B10'）;
A2=xlsread（'第一题数据','下边疆','A11:
B13'）;
A3=xlsread（'第一题数据','下边疆','A13:
B31'）;
A4=xlsread（'第一题数据','下边疆','A32:
B39'）;
A5=xlsread（'第一题数据','下边疆','A39:
B40'）;
A6=xlsread（'第一题数据','下边疆','A40:
B56'）;
%导入下边疆的实验数据并分好计算的数据段
x1=A1（:
1）;x2=A2（:
1）;x3=A3（:
1）;
x4=A4（:
1）;x5=A5（:
1）;x6=A6（:
1）;
y1=A1（:
2）;y2=A2（:
2）;y3=A3（:
2）;
y4=A4（:
2）;y5=A5（:
2）;y6=A6（:
2）;%给相应的变量赋值
x10=min（x1）:
0.001:
max（x1）;%对每一段数据点按连结顺序进行点横坐标的加密处理
x20=max（x2）:
-0.001:
min（x2）;
x30=min（x3）:
0.001:
max（x3）;
x40=min（x4）:
0.001:
max（x4）;
x50=max（x5）:
-0.001:
min（x5）;
x60=min（x6）:
0.001:
max（x6）;
y1i=interp1（x1,y1,x10,'spline'）;%对数据进行三次样条插值
y2i=interp1（x2,y2,x20,'spline'）;
y3i=interp1（x3,y3,x30,'spline'）;
y4i=interp1（x4,y4,x40,'spline'）;
y5i=interp1（x5,y5,x50,'spline'）;
y6i=interp1（x6,y6,x60,'spline'）;
x=[x10x20x30x40x50x60];y=[y1iy2iy3iy4iy5iy6i];
%对三次样条插值后的数据按连结顺序合并
s10=trapz（x10,y1i）+trapz（x30,y3i）+trapz（x40,y4i）+trapz（x60,y6i）;
%用梯形公式求除特殊段以外的各段与X轴围成的面积并求和
s11=-trapz（x20,y2i）;%对特殊段面积s11的计算
s12=-trapz（x50,y5i）;%对特殊段面积s12的计算
s1=s10-s11-s12;%计算下疆界与X轴围成面积的精确值
plot（x,y）
holdon;
A1=xlsread（'第一题数据','上边疆','A03:
B11'）;
A2=xlsread（'第一题数据','上边疆','A11:
B14'）;
A3=xlsread（'第一题数据','上边疆','A14:
B40'）;
A4=xlsread（'第一题数据','上边疆','A40:
B41'）;
A5=xlsread（'第一题数据','上边疆','A41:
B49'）;
A6=xlsread（'第一题数据','上边疆','A50:
B52'）;
A7=xlsread（'第一题数据','上边疆','A53:
B64'）;
%导入上边疆的实验数据并分好计算的数据段
x1=A1（:
1）;x2=A2（:
1）;x3=A3（:
1）;x4=A4（:
1）;
x5=A5（:
1）;x6=A6（:
1）;x7=A7（:
1）;
y1=A1（:
2）;y2=A2（:
2）;y3=A3（:
2）;y4=A4（:
2）;
y5=A5（:
2）;y6=A6（:
2）;y7=A7（:
2）;%给相应的变量赋值
x10=min（x1）:
0.001:
max（x1）;%对每一段数据点按连结顺序进行点横坐标的加密处理
x20=max（x2）:
-0.001:
min（x2）;
x30=min（x3）:
0.001:
max（x3）;
x40=max（x4）:
-0.001:
min（x4）;
x50=min（x5）:
0.001:
max（x5）;
x60=min（x6）:
0.001:
max（x6）;
x70=min（x7）:
0.001:
max（x7）;
y1i=interp1（x1,y1,x10,'spline'）;%对数据进行三次样条插值
y2i=interp1（x2,y2,x20,'spline'）;
y3i=interp1（x3,y3,x30,'spline'）;
y4i=interp1（x4,y4,x40,'spline'）;
y5i=interp1（x5,y5,x50,'spline'）;
y6i=interp1（x6,y6,x60,'spline'）;
y7i=interp1（x7,y7,x70,'spline'）;
x=[x70x60x50x40x30x20x10];y=[y7iy6iy5iy4iy3iy2iy1i];
%对三次样条插值后的数据按连结顺序合并
s20=trapz（x10,y1i）+trapz（x30,y3i）+trapz（x50,y5i）+trapz（x60,y6i）+trapz（x70,y7i）;%用梯形公式求除特殊段以外的各段与X轴围成的面积并求和
s21=-trapz（x20,y2i）;%对特殊段面积s21的计算
s22=-trapz（x40,y4i）;%对特殊段面积s22的计算
s2=s20-s21-s22;%计算上疆界与X轴围成面积的精确值
s=s2-s1；
S=s*100/9；%计算实际国土面积
S
plot（x,y）
gridon;
t=abs（1.5665e+006-S）/1.5665e+006;%计算相对误差
t
升级会员 