チェビシェフ多項式(Chebyshev polynomials)でフィットさせた係数から値を復元

区間  \displaystyle  [t_s, t_e ] の  \displaystyle  t に対して、  \displaystyle  f(t) をチェビシェフ多項式でフィットさせた係数  \displaystyle  c_0, c_1, c_2, \cdots, c_N が与えられているとします。

係数から、以下の手順で  \displaystyle  f(t) を計算できます。

手順

 \displaystyle  t 区間  \displaystyle  [-1, 1] に正規化したものを  \displaystyle  x とします。

 \displaystyle
x = 2\frac{t - t_s}{t_e - t_s} - 1

チェビシェフ多項式(Chebyshev polynomials)  \displaystyle  T_0, T_1, \cdots, T_N を次の漸化式に従って計算します。

 \displaystyle
\begin{align}
T_0(x) &= 1 \\
T_1(x) &= x \\
T_{k+1}(x) &= 2xT_k(x) - T_{k-1}(x)
\end{align}

最後、  \displaystyle  f(t) は次の式で計算できます。

 \displaystyle
\begin{align}
f(t) &\approx \sum_{k=0}^{N} c_k T_k \\
\end{align}

サンプルコード

以下は係数  \displaystyle  c_i (cs) から  \displaystyle  f(t) を計算するScalaでのサンプルコードです。xは正規化済みとします。高速化のためScalaらしくないコードです。

  def calcChebyshev(cs: IndexedSeq[Double], x: Double): Double = {
    // x は -1.0 <= x <= +1.0
    val n = cs.size;
    if (n == 1) {
      cs(0);
    } else {
      val x2 = 2.0 * x;
      var result: Double = cs(0) + cs(1) * x;
      var i: Int = 2;
      var t0: Double = 1.0;
      var t1: Double = x;
      while (i < n) {
        val t2 = x2 * t1 - t0;
        result += cs(i) * t2;
        t0 = t1;
        t1 = t2;
        i = i + 1;
      }
      result;
    }
  }

JPLのデータを使うときに必要な計算式です。

JPL Planetary and Lunar Ephemerides