c - acos() 在 C 中有两个接近点的坏结果

Question

我正在创建一个函数，该函数使用两点的纬度/经度(以度数而非弧度为单位)和余弦球面定律来计算两点的距离。我遇到的问题是，由于函数 acos() 中的舍入误差，当两点彼此非常接近时，我获得的结果远非如此。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef struct{
 double lat;
 double lon;
} point;

#define DEG_TO_RAD 0.017453292519943295769236907684886
#define EARTH_RADIUS_IN_METERS 6372797.560856

double distance(point a, point b) {
  double arg=sin(a.lat * DEG_TO_RAD) * sin(b.lat * DEG_TO_RAD) + cos(a.lat * DEG_TO_RAD) * cos(b.lat * DEG_TO_RAD) * cos((a.lon-b.lon) * DEG_TO_RAD);
  if(arg>1) arg=1;
  else if (arg<-1) arg=-1;
  printf("arg=%.12f acos(arg)=%.12f\n",arg, acos(arg));    //to see the problem
  return acos(arg) * EARTH_RADIUS_IN_METERS;
}

int main(){
  point p1,p2;

  p1.lat=63.0;
  p1.lon=27.0;
  p2.lat=p1.lat;
  p2.lon=p1.lon;

  printf("dist=%.8f\n",distance(p1,p2));

  return 0;
}

输出是

arg=1.000000000000 acos(arg)=0.000000014901
dist=0.09496208

如您所见，当它计算 acos() 时，它应该给出零，但它确实给出了一些误差，这些误差在乘以地球半径后被极大地放大了。当两点不相等但非常接近时也会发生这种情况。如果它有任何用处，我的经纬度数据最多有 7 位小数。

Answer 1

你从 acos 得到的结果和它得到的一样好:问题是 arg 的计算总是会有一个小错误并返回一个值是稍微关闭。当两点相等或非常接近时，结果小于一，例如 1-10^-16。如果你看the graph of acos(x)您会看到它在 x=1 时几乎是垂直的，这意味着即使 arg 中最轻微的错误也会对相对错误产生巨大影响。换句话说，该算法在数值上是不稳定的。

您可以使用 haversine formula以获得更好的结果。