こんにちは。
回転楕円体面上（地球上）の二点間測地線距離 (geodesic distance) を計算する各方法の最大相対誤差 (maximum relative accuracy) の距離依存性を調べ比較しました（評価結果は下記のチャート）。

評価のための処理および各計算法は「測地線距離計算式・計算ライブラリの精度評価」（330K INFO）からそっくり拝借して使いました（ありがとうございます）。
なお今回は与えられた距離に対する最大（相対）誤差評価です（緯度、方位角依存性は調べません¹）。

計算方法選定ガイド

評価結果を参考にすると、簡素性優先かつ比較的精度が良い方法を求める方には、下記あたりが候補だろうと思いました。

短距離向け：Straight-tunnel-line、Gauss-mid-latitude、...
長距離向け：Andoyer、...

精度評価用入力データ

精度評価用入力データは楕円体面上に稠密に二点を発生させました（15万対）。その真距離を得るには GeographicLib（の GeodSolve）を利用しました。追加として"GeodTest.dat" (GeographicLib) も使いました（距離 10,000 km 以上）。

各計算法の紹介²

Spherical: Menelaus of Alexandria (about 100); Great-circle distance - Wikipedia; 完全な球面近似であり回転楕円体面を考慮しない形; 通称 haversine（haversine formula を用いて計算する場合）
Flat-surface: Newton, Cassini (1713); 起源探究は断念しました; Geographical distance - Wikipedia; 擬似短距離一次近似（平面近似）; 日本での通称は Hubeny（簡易式）（ただし不適切）
Hubeny (1954, 1959) 擬似短距離近似。Flat-surface よりも近似精度が改善されている。
Straight-tunnel-line: 「短距離近似式」- wikipedia
Gauss-mid-latitude: Gauss (1843, 1846); Lambert and Swick (1935); Rapp (1991) §6.4; Geographical distance - Wikipedia; 「回転楕円体面に沿う最短距離の式」- wikipedia
Bowring-short: Bowring (1981); Rapp (1991) §6.5; Geographical distance - Wikipedia
Andoyer: Forsyth (1895), Andoyer (1932); Thomas (1965)³
Andoyer-Lambert: Lambert (1942); Thomas (1965)³; Geographical distance - Wikipedia
Andoyer-Lambert-Thomas: Thomas (1970)³
Bowring-inverse: Bowring (1996); used by Geospatial Information Authority of Japan (GSI) based on this method, 距離と方位角の計算　計算式（国土地理院）
Vincenty: Vincentry (1975, 1976); 回転楕円体上の測地線及び航程線の算出について（海洋情報部研究報告）
Karney: Karney (2011); GeographicLib

なお今回私自身で書き下して用いたソースコード（下記）はこれらのうち spherical, flatsurface, straightline, gaussmidlatitude, andoyer のみです（どれも簡素な計算方法）。

短評

今回の評価結果に対する短評をまとめました。

spherical は誤差 5 ‰。距離依存性もほぼなく、対蹠点まで安定している。
flatsurface と hubeny は普遍的な短距離近似ではなく、緯度依存による計算精度変動が極めて大きいので、二点間測地線が高緯度・極付近を通過する条件では、短距離でも極めて大きい誤差を生じる⁴。
Straight-tunnel-line の誤差は、100 km で 10 ppm、1,000 km で 80 ppm、3,000 km で 1 %（spherical に逆転される）。
gaussmidlatitude の誤差は、100 km で 1 ppm、1,000 km で 8 ppm、10,000 km で 3 ‰（spherical に逆転される）。10,000 km 以上は対蹠点まで安定している。
andoyer⁵ の誤差は、10,000 km 以下では常に 15 ppm。ただし 10,000 km 以上で急激に増大する（対蹠点で 3 ‰）。
andoyerlambert の誤差は、10,000 km 以下では常に 1.5 ppm⁶。ただし 10,000 km 以上で急激に増大する（対蹠点で 3 ‰）。
andoyerlambertthomas の誤差は、10,000 km 以下では常に 3 ppb。ただし 10,000 km 以上で急激に増大する（対蹠点で 3 ‰）。
vencenty と bowring-inverse の誤差は安定して小さいが、対蹠点直近で誤差が急増する。反復型であり計算速度はそれほど速くない。（bowring-inverse は加えて NaN 発生回避が必要⁷。）
karney の誤差は、今回の検討では絶対誤差が常に 10 nm 以下。対蹠点まで安定している。反復型であり計算速度は速くない。処理手順は複雑だが、ライブラリとして使いやすく整備されている。

総合すると、Straight-tunnel-line、Gauss-mid-latitude あたりが、簡素な計算法を求めるが多少近似精度を気にする方の常用のための候補に見えますが、短距離向きですので 100 km 以上での誤差は 1 ppm を超え増大していきます。これを気にする方は、Andoyer あたりが次の候補になるかと思います。

誤差評価結果

相対誤差最大値（縦軸）の距離依存性を評価し、距離（横軸）範囲の段階を変えて３つのチャートを描きました。

計算速度評価

formula	mean calculation time [ms]
flatsurface	23.68
hubeny	38.82
spherical	31.44
straightline	24.79
gaussmidlatitude	41.57
bowring_short	37.55
bowring_short_	45.41
bowring_short_beta_	61.89
andoyer	30.54
andoyerlambert	47.42
andoyerlambertthomas	52.63
vincenty	191.15
bowring-inverse	205.84
bowring-inverse_	204.67
karney	519.06

source code

いずれも haversine 型の距離計算を用いて浮動小数点計算誤差を最小化しています（球面余弦定理を用いるのではなく）⁸。

spherical.js

function spherical(lat1, lon1, lat2, lon2){  // in degrees
  const a_adjusted = 6367000;  // Earth radius adjusted (≈ 0.998 * a_GRS80)
  const degree = Math.PI / 180;
  const sin = Math.sin;
  const cos = Math.cos;
  const hypot = Math.hypot;
  const asin2 = x => x < 1 ? 2 * Math.asin(x) : Math.PI;

  const latdiffhalf = (lat1 - lat2) / 2 * degree;
  const londiffhalf = (lon1 - lon2) / 2 * degree;
  const latmid = (lat1 + lat2) / 2 * degree;

  return a_adjusted * asin2(hypot(
    sin(londiffhalf) * cos(latmid),
    cos(londiffhalf) * sin(latdiffhalf)
  ));
}

flatsurface.js

function flatsurface(lat1, lon1, lat2, lon2){  // in degrees
  const a = 6378137;  // GRS80
  const f = 1 / 298.257223563;  // WGS84
  const e2 = f * (2 - f);
  const degree = Math.PI / 180;
  const abs = Math.abs;
  const sin = Math.sin;
  const cos = Math.cos;
  const sqrt = Math.sqrt;
  const hypot = Math.hypot;
  
  const londiff = abs(lat1) == 90 || abs(lat2) == 90 ? 0 :  // at the north and south poles
                 (((lon1 - lon2 + 180) % 360) - 180) * degree;  // across the 180th meridian
  const latdiff = (lat1 - lat2) * degree;
  const latmid = (lat1 + lat2) / 2 * degree;
  const sinlat = sin(latmid);
  const coslat = cos(latmid);
  const n2 = 1 / (1 - e2 * sinlat ** 2);
  const n = sqrt(n2);  // prime vertical radius of curvature
  const m_by_n = (1 - e2) * n2;  // meridian ratio

  return a * n * hypot(
    londiff * coslat, 
    latdiff * m_by_n
  );
}

straightline.js

function straightline(lat1, lon1, lat2, lon2){  // in degrees
  const a = 6378137;  // GRS80
  const f = 1 / 298.257223563;  // WGS84
  const e2 = f * (2 - f);
  const degree = Math.PI / 180;
  const sin = Math.sin;
  const cos = Math.cos;
  const sqrt = Math.sqrt;
  const hypot = Math.hypot;

  const latdiffhalf = (lat1 - lat2) / 2 * degree;
  const londiffhalf = (lon1 - lon2) / 2 * degree;
  const latmid = (lat1 + lat2) / 2 * degree;
  const sinlat = sin(latmid);
  const coslat = cos(latmid);
  const n2 = 1 / (1 - e2 * sinlat ** 2);
  const n = sqrt(n2);  // prime vertical radius of curvature
  const m_by_n = (1 - e2) * n2;  // meridian ratio	

  return 2 * n * a * hypot(
    sin(londiffhalf) * coslat,
    cos(londiffhalf) * latdiffhalf * m_by_n
  );
}

gaussmidlatitude.js

function gaussmidlatitude(lat1, lon1, lat2, lon2){  // in degrees
  const a = 6378137;  // GRS80
  const f = 1 / 298.257223563;  // WGS84
  const e2 = f * (2 - f);
  const degree = Math.PI / 180;
  const sin = Math.sin;
  const cos = Math.cos;
  const sqrt = Math.sqrt;
  const hypot = Math.hypot;
  const asin2 = x => x < 1 ? 2 * Math.asin(x) : Math.PI;

  const latdiffhalf = (lat1 - lat2) / 2 * degree;
  const londiffhalf = (lon1 - lon2) / 2 * degree;
  const latmid = (lat1 + lat2) / 2 * degree;
  const sinlat = sin(latmid);
  const coslat = cos(latmid);
  const n2 = 1 / (1 - e2 * sinlat ** 2);
  const n = sqrt(n2);  // prime vertical radius of curvature
  const m_by_n = (1 - e2) * n2;  // meridian ratio	

  return n * a * asin2(hypot(
    sin(londiffhalf) * coslat,
    cos(londiffhalf) * sin(latdiffhalf * m_by_n)
  ));
}

andoyer.js

function andoyer(lat1, lon1, lat2, lon2){  // in degrees
  const a = 6378137;  // GRS80
  const f = 1 / 298.257223563;  // WGS84
  const degree = Math.PI / 180;
  const sin = Math.sin;
  const cos = Math.cos;
  const hypot = Math.hypot;
  const asin2 = x => x < 1 ? 2 * Math.asin(x) : Math.PI;

  const latmid = (lat1 + lat2) / 2 * degree;
  const latdiffhalf = (lat1 - lat2) / 2 * degree;
  const londiffhalf = (lon1 - lon2) / 2 * degree;

  const sin_lat = sin(latmid);
  const cos_lat = cos(latmid);
  const sin_latdiffhalf = sin(latdiffhalf);
  const cos_latdiffhalf = cos(latdiffhalf);
  const d = asin2(hypot(  // 0-th order distance
    sin(londiffhalf) * cos_lat,
    cos(londiffhalf) * sin_latdiffhalf
  ));
  const c = (3 * sin(d) - d) * (sin_lat * cos_latdiffhalf / cos(d / 2)) ** 2;
  const s = (3 * sin(d) + d) * (cos_lat * sin_latdiffhalf / sin(d / 2)) ** 2;
  const delta_1st = f * (c - s) / 2;  // 1st order
  return a * (d + delta_1st);
}

精度評価用入力データ生成

もしも半球内（10000km以内）に限定して精度評価するのであれば、それほど稠密にデータを生成する必要はなく、下記のようにデータ生成すれば十分と思います。

$ ./generate_geodtest_dat.sh | wc -l
    3056

generate_geodtest_dat.sh

#!/bin/sh
# constants
angle_step=10  # in degree
dist_max=19970326  # in meter # a * pi * (1 - f)
commands_required="GeodSolve" # GeographicLib utility

# functions
geodsolve_f () {
  args=$@  # lat1 lon1 azi1 s12
  GeodSolve -p 10 --input-string "$args"  # calculates lat2 lon2 azi2
}

geodsolve_pole_f () {
  dist=$1
  dist_half=$(awk "BEGIN { print $dist/2 }")
  lat=$(geodsolve_f 90 0 180 $dist_half | cut -d" " -f1)
  echo $lat 0 0 $lat 180 180 $dist
}

geodsolve_loop_f () {
  dist=$1
  range_lat=$2
  range_azi=$3
  for lat in $range_lat; do
    for azi in $range_azi; do
      dest=$(geodsolve_f $lat 0 $azi $dist)
      echo $lat 0 $azi $dest $dist
    done
  done
}

# main
range_azi="$(seq 0 $angle_step 90)"
range_lat="$(seq -90 $angle_step 90)"
for dist_exponent in $(seq 0 7); do  # 1 10 100 1000 10000 100000 1000000 10000000
for dist_significand in 1 3; do
  dist=$((10 ** dist_exponent * dist_significand))
  [ $dist -gt $dist_max ] && dist=$dist_max
  geodsolve_pole_f $dist
  geodsolve_loop_f $dist "$range_lat" "$range_azi"
done
done

例えば sphericalを除くどの式も、赤道上の二点間に対しては式の上では無誤差です。 ↩
なお文献一覧は、"A geodesic bibliography" (GeographicLib) もよくまとまっています。 ↩
この Thomas (1965) は、Forsyth (1895), Andoyer (1932); (Marie Henri Andoyer - Wikipedia), Lambert: Lambert (1942); (Geographical distance - Wikipedia) の先行研究に基づき、一次近似を再整理したもの。Thomas (1970) は、二次近似を正しく求めたもの。 ↩ ↩² ↩³
これらは経緯度変数の微小量（による局所展開形）に由来し、両極に特異点を持つことが大きな誤差の原因となっている。仮に、評価用入力データからそれらを除外すれば短距離精度は改善される（例えば日本国内の二点間距離を計算したい場合には向いていると言えそう）。 ↩
Boost.Geometry がデフォルト計算法として採用。 ↩
parametric latitude 使用による精度向上の反面、計算速度が若干低下します。 ↩
改善しました： https://qiita.com/kkdd/items/ec4b88659d9baa8bc162 ↩
参考：球面大円距離計算式の比較（浮動小数点計算精度） ↩