地球椭球体表面方程的不同近似刻画公式的优劣性比较.docx

地球椭球体表面方程的不同近似刻画公式的优劣性比较

一.实验目的

地球椭球体表面方程的不同近似刻画公式的优劣性比较。通过不同公式的计算结果，分析各个公式的计算误差及其优劣。

二.实验原理

1、计算并绘图

将角度等分为适合数量，分别计算出精确公式以及三个近似公式对应的半径值。利用画图软件画出“角度——半径”曲线，根据对称性，只需要0~π/2的半径即可。特别的，在计算椭球体的半径时候，可以选取yoz平面所在的截面进行计算。从所做图中粗略的比较异同。

三个近似公式为：

1r? a

1

1??2fcos2?

2r??a?1??fcos2??

2

? ? 2

1??

r3??R?1??f?cos

?? ?

3

? ? ??

?其中，a，f，R为常数。精确值计算方法：

?

地球椭球（单位：km）

2 2 2

x??yz

椭圆体方程 ? ??1，用yoz平面截取，得到

a2 b2

2y z2

2

?

a2 b2

??1；

又因为

tan??????y，

z

即可得到角度与半径（

2、误差比较分析

y2

y2?z2

计算出三个近似公式的绝对误差与相对误差，作出“角度——误差”曲线，更直观地比较三个近似公式的误差大小以及优劣性。

三.实验过程

1、计算半径并绘图

利用Matlab软件编写代码，作出四张“角度——半径”的曲线图。实验中，90°被平分为2500份。

代码如下：a=6378.135;b=6356.751;f=1/298.26;

x=0:pi/5000:pi/2;z=a./(sqrt(tan(x).*tan(x)+a*a/(b*b)));y=sqrt((1-z.*z/(b*b))*a*a);

subplot221plot(x,sqrt(y.*y+z.*z));

title(无近似);

xlabel(角度（弧度制）);ylabel(半径);

y1=a./(sqrt(1+2*f.*cos(x).*cos(x)));subplot222;

plot(x,y1,r);

title(近似1);xlabel(角度（弧度制）);ylabel(半径);

y2=a*(1-f*cos(x).*cos(x));subplot223;plot(x,y2,y);

title(近似2);

xlabel(角度（弧度制）);ylabel(半径);

y3=6371*(1-f*(cos(x).*cos(x)-1/3));subplot224;

plot(x,y3,g);

title(近似3);xlabel(角度（弧度制）);ylabel(半径);

2、计算绝对误差与相对误差并绘图

在此基础上，再改写代码，画出“角度——绝对误差”与“角度——相对误差”曲线，代码如下：

绝对误差代码：a=6378.135;b=6356.751;f=1/298.26;

x=0:pi/5000:pi/2;z=a./(sqrt(tan(x).*tan(x)+a*a/(b*b)));y=sqrt((1-z.*z/(b*b))*a*a);

y1=a./(sqrt(1+2*f.*cos(x).*cos(x)));y11=abs(sqrt(y.*y+z.*z)-y1);

subplot221;plot(x,y11,r);

title(第一种近似绝对误差);

xlabel(角度（弧度制）);ylabel(误差);

y2=a*(1-f*cos(x).*cos(x));y22=abs(sqrt(y.*y+z.*z)-y2);subplot222;

plot(x,y22,b);title(第二种近似绝对误差);xlabel(角度（弧度制）);ylabel(误差);

y3=6371*(1-f*(cos(x).*cos(x)-1/3));y33=abs(sqrt(y.*y+z.*z)-y3);

subplot223;plot(x,y33,g);

title(第三种近似绝对误差);

xlabel(角度（弧度制）);ylabel(误差);

相对误差代码：a=6378.135;b=6356.751;f=1/298.26;

x=0:pi/5000:pi/2;z=a./(sqrt(tan(x).*tan(x)+a*a/(b*b)));y=sqrt((1-z.*z/(b*b))*a*a);

y1=a./(sqrt(1+2*f.*cos(x).*cos(x)));y11=abs(sqrt(y.*y+z.*z)-y1)./sqrt(y.*y+z.*z);subplot221;

plot(x,y11,r);title(第一种近似相对误差)