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

方法を選定したい方のためには、簡素な計算のうちでは比較的精度が良い Gauss-mid-latitude もしくは Andoyer あたりが候補だろうと思いました（下記のコメント）。
評価のための処理および各計算法は「測地線距離計算式・計算ライブラリの精度評価」（330K INFO）からそっくり拝借して使いました（ありがとうございます）。
なお実装の際の注意点は余弦定理を用いるのではなく haversine 型の距離計算を使う方が精度上良いです¹。

精度評価用入力データ

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

各計算法の紹介²

Spherical: Menelaus of Alexandria (about 100); Great-circle distance - Wikipedia; 完全な球面近似であり回転楕円体面を考慮しない形; 通称 haversine
Flat-surface: Newton, Cassini (1713); 起源探究は断念しました; Geographical distance - Wikipedia; 一次近似（平面近似）; 日本での通称は Hubeny（簡易式）
Hubeny (1954, 1959)
Straight-tunnel-line: Geographical distance - 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: Forsynth (1895), Andoyer (1932); Thomas (1965)³
Andoyer-Lambert: Lambert (1942); Thomas (1965)³; Geographical distance - Wikipedia
Andoyer-Lambert-Thomas: Thomas (1970)³
Bowring: Bowring (1996); used by Geospatial Information Authority of Japan (GSI) based on this method, 距離と方位角の計算　計算式（国土地理院）
Vincenty: Vincentry (1975, 1976); 回転楕円体上の測地線及び航程線の算出について（海洋情報部研究報告）
Karney: Karney (2011); GeographicLib

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

コメント

今回の評価結果に対するコメントおよび注意点は（厳密ではありませんが）、

spherical は最大相対誤差 5 ‰。距離依存性がほぼなく、対蹠点まで安定している。
なおどの計算法も誤差は緯度依存性も持つが今回の評価結果には表現されない（例えば sphericalを除き、赤道上の二点間に対しては式の上では誤差が無くなる）。
flatsurface と hubeny の方法は経緯度変数の微小量（による局所展開形）に由来し、両極に特異点を持つので、二点間測地線が高緯度・極付近を通過する条件で、短距離でも大きい誤差を生じる⁴。
gaussmidlatitude の誤差は 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 の誤差は安定して小さいが、対蹠点直近で誤差が急増する。計算速度はそれほど速くない。
bowring は加えて実装が容易ではなく、また NaN となる計算結果の回避も不明⁶。
karney の誤差は、今回の検討では絶対誤差が常に 10 nm 以下。対蹠点まで安定している。計算速度は速くない。処理手順は複雑だが、ライブラリとして使いやすく整備されている。

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

評価結果

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

各方法のソース

spherical.js

function spherical(lat1, lon1, lat2, lon2){
  "use strict";
  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 lat = (lat1 + lat2) / 2 * degree; // middle

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

flatsurface.js

function flatsurface(lat1, lon1, lat2, lon2){
  "use strict";
  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 lat = (lat1 + lat2) / 2 * degree; // middle
  const sinlat = sin(lat);
  const coslat = cos(lat);
  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
  );
}

gaussmidlatitude.js

function gaussmidlatitude(lat1, lon1, lat2, lon2){
  "use strict";
  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 lat = (lat1 + lat2) / 2 * degree; // middle
  const sinlat = sin(lat);
  const coslat = cos(lat);
  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){
  "use strict";
  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 lat = (lat1 + lat2) / 2 * degree; // middle
  const latdiffhalf = (lat1 - lat2) / 2 * degree;
  const londiffhalf = (lon1 - lon2) / 2 * degree;

  const sin_lat = sin(lat);
  const cos_lat = cos(lat);
  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

#!/bin/sh
# constants
step=10  # in degree
dist_max=19970326  # in meter # a * pi * (1 - f)

# functions; using GeodSolve (GeographicLib Utilities)
geodsolve_f () {
  args=$@  # lat1 lon1 azi1 s12
  GeodSolve -p 9 --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 $step 90)"
range_lat="$(seq -90 $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 law of cosines（余弦の球面法則）を用いると短距離での計算精度が劣化します。世の中の実装をいくつか調べてみますと、少なからず見られました（例えば、Boost.Geometry、Forsynth, Andoyer, Lambert, Thomas 系実装）。 ↩
なお文献一覧は、"A geodesic bibliography" (GeographicLib) もよくまとまっています。 ↩
この Thomas (1965) は、Forsynth (1895), Andoyer (1932); (Marie Henri Andoyer - Wikipedia), Lambert: Lambert (1942); (Geographical distance - Wikipedia) の先行研究に基づき、一次近似を再整理したもの。Thomas (1970) は、二次近似を正しく求めたもの。 ↩ ↩² ↩³
仮に、評価用入力データからそれらを除外すれば短距離精度は改善されるはず。例えば日本国内の二点間距離を計算したい場合には向いていると言えそう。 ↩
parametric latitude 使用による精度向上の反面、計算速度が若干低下します。 ↩ ↩²
これらへ対処可能であれば有望と思われます。 ↩