6. 補間 (1) ラグランジュ補間, エルミート補間, 直交多項式補間

多項式補間

区間 [a, b] 内の n + 1 個の相異なる点 x₀, x₁, …, x_n における関数 f(x) の値 f(x₀), f(x₁), …, f(x_n) がデータとして与えられているとき, これらの点における値がデータに一致するような近似式を f(x) の補間式という. これら n + 1 個の点を標本点あるいは補間点という. 補間式は x₀, …, x_n 以外の点における値を推測するために用いられる. ただし, a = min(x₀, …, x_n), b = max(x₀, …, x_n) である.

x が [a, b] 内にあるときに「補間」(あるいは「内挿」), x が [a, b] の外にあるときに「補外」(あるいは「外挿」) ということがある.

計算が難しい(あるいは, 時間がかかる)関数 f(x) を計算しやすい補間式で代替するといった使い方をする場合は「関数近似」への応用とみることができる.

多項式補間

補間式として多項式を使用するものを多項式補間という.

n 次多項式 p_n(x) を次のように定義する.

  
  pn(x) = a0xn + a1xn-1 +  ... an-1x + an

 
n + 1 個の相異なる点 x₀, x₁, …, x_n における関数値 f(x₀), f(x₁), …, f(x_n) が与えられているとき, 次の連立一次方程式を解くことにより n 次多項式 p_n(x) の係数 a₀, a₁, …, a_n を求めることができる.

  
  Va = f

 
ただし,

  
      ( x0n  x0n-1  ...  1 )       ( a0  )       ( f0  )
  V = ( x1n  x1n-1  ...  1 ),  a = ( a1  ),  f = ( f1  )
      (       ...          )       ( ... )       ( ... )
      ( xnn  xnn-1  ...  1 )       ( an  )       ( fn  )

 
ここで, V は連立一次方程式の係数行列, f は連立一次方程式の右辺ベクトル(各点における関数値 f_i = f(x_i) を表す), a は連立一次方程式の解ベクトル(多項式の係数)となる.

V はファンデルモンド行列と呼ばれる形をしている. ファンデルモンド行列の行列式は次式で求められることが知られている.

  
  det(V) = Π(xi - xj) (0 ≥ i > j ≥ n)

 
x₀, x₁, …, x_n が相異なる点であれば, 上式より det(V) ≠ 0 である. すなわち, 上の連立一次方程式は解を持つ. この n 次多項式 p_n(x) を補間式として採用する.

数値実験 (1)

区間 [-1, 1] において次の指数関数を多項式補間する.

  
  f(x) = 2ex - 1 - 1

 
n = 2 ～ 5 として(すなわち, 2 ～ 5 次多項式を用いて)補間を行う. ただし, 区間 [-1, 1] を n 等分して, 両端を含め等間隔に並んだ n + 1 点における関数値をデータとして使用することにする. (多項式補間としては各標本点が相異なっていればよく等間隔である必要はないが, 数表などを想定すればこのように標本点を選ぶのは自然である).

上の連立一次方程式を解いて多項式の係数を求め, 区間内における多項式の値をプロットする. 横軸は x, 縦軸は関数値である.

多項式の次数が上がると精度も上がっていき, x = 0.65 付近でみると 2 ～ 5 次多項式はそれぞれ 1 ～ 4 桁程度一致している.

次に, 多項式補間ではどんな関数についても精度を上げたいときには次数を上げていけばよいというわけではないことを示す.

数値実験 (2)

区間 [-1, 1] において次の関数(ルンゲの関数という)を多項式補間する.

  
  f(x) = 1/(1 + 25x2)

 
n = 4, 8, 12, 20 として(すなわち, 4, 8, 12, 20 次多項式を用いて)補間を行う. ただし, 区間 [-1, 1] を n 等分して, 両端を含め等間隔に並んだ n + 1 点における関数値をデータとして使用する.

上の連立一次方程式を解いて多項式の係数を求め, 区間内における多項式の値をプロットする. 横軸は x, 縦軸は関数値である.

n が大きくなると両端で発散した(ルンゲの現象という). 実は最大誤差だけを見ると(すべての n の中で) n = 4 のときが最小になりそれ以上だとかえって悪くなる. ただし, 図のように中央付近では次数が上がるとよい近似をしていることがわかるので, この領域だけを使うという手はあるかもしれない.

この問題は標本点を等間隔ではなくチェビシェフ多項式のゼロ点にとると解決することが知られている. 「直交多項式補間」の数値実験 (4) を参照せよ.

ラグランジュ補間

多項式補間を行うには, 前項のように連立一次方程式を解いて多項式の係数を直接求めるのが直観的な方法ではあるが, 実務的には方程式の条件数が大きくなりやすいなどの理由により, 等価な多項式であるラグランジュ補間公式が用いられる.

ラグランジュ補間公式は次式で表される.

  
  f(x) = Σ Lk(x) f(xk)  (k = 0 ～ n)

  ただし
    Lk(x) = Π(x - xi)/(xk - xi)  (i = 0 ～ k - 1, k + 1 ～ n)

 
L_k(x_i) は i = k のとき 1, i ≠ k のとき 0 である. したがって, 各標本点において与えられたデータに一致する. ラグランジュ補間公式は標本点の並んでいる順序に関係なく成り立つ. また, 標本点は等間隔でなくてもよい.

ラグランジュ補間公式は n 次多項式であり, バラバラにして整理するとその係数は連立一次方程式を解いて求めた多項式の係数と一致するはずである. しかし, 補間においては関数値が求まればよくて多項式の係数は計算に直接は必要ないので, この形のまま計算に使われる.

ラグランジュ補間公式の計算は次のようにして行うとよい.

まず, 与えられたデータについて

  
  wk = f(xk)/(Π(xk - xi) (i = 0 ～ k - 1, k + 1 ～ n))

 
を k = 0 ～ n について計算しておく.

次に, 求めたい点 x それぞれについて次のようにして補間値を計算する.

  
  w(x) = Π(x - xk)  (k = 0 ～ n)
  pn(x) = w(x)Σwk/(x - xk) (k = 0 ～ n)

 
ただし, x = x_k (k = 0, … または n) であれば計算せずに p_n(x) = f(x_k) とすればよい.

上の数値実験 (1), (2) についてラグランジュ補間公式により計算すると多項式補間のときと同じ結果が得られる.

差分商

x = x₀, x₁ に対する差分商を次のように定義する.

  
  f[x0, x1] = (f(x1) - f(x0))/(x1 - x0)

 
これは, x₀ と x₁ が近ければ微分 f'(x₀) を表す. あるいは, これは区間 [x₀, x₁] における微分値の平均値とみなすことができる.

同様に, x = x₀, x₁, x₂ に対する 2 階の差分商は次のように定義できる.

  
  f[x0, x1, x2] = (f[x1, x2] - f[x0, x1])/(x2 - x0)

 
すなわち, 2 階の2階の差分商は 1 階の差分商の差分商である.

便宜上 f[x_i] = f(x_i) と表すと, 高階の差分商まで含めた一般形は n + 1 個の標本点 x₀, x₁, …, x_n について次のように表すことができる.

  
  f[xi, ..., xi+k] = (f[xi+1, ..., xi+k] - f[xi, ..., xi+k-1])/(xi+k - xi) (k = 1 ～ n, i = 0 ～ n - k)

 
ここで, 関数 f(x) を差分商を用いて表すことを考える.

x, x₀ に対する 1 階の差分商

  
  f[x, x0] = (f(x0) - f(x)/(x0 - x)

 
を変形すると次のようになる.

  
  f(x) = f(x0) + (x - x0)f[x, x0]

 
さらに, x, x₀, x₁ に対する 2 階の差分商

  
  f[x, x0, x1] = (f[x0, x1] - f[x, x0])/(x1 - x)

 
を f[x, x₀] について解き, 代入すると次のようになる.

  
  f(x) = f(x0) + (x - x0)f[x0, x1] + (x - x0)(x - x1)*f[x, x0, x1]

 
同様の代入を繰り返すと, n + 1 個の標本点 x₀, x₁, …, x_n について次式が得られる.

  
  f(x) = f(x0) + (x - x0)f[x0, x1] + (x - x0)(x - x1)*f[x0, x1, x2]
       + ...
       + (x - x0)(x - x1)...(x - xn-1)*f[x0, x1, ..., xn/sub>]
       + (x - x0)(x - x1)...(x - xn-1)(x - xn)*f[x, x0, x1, ..., xn]

 
ここで, 最後の項は f[x, x₀, x₁, …, x_n] に変数 x を含むため通常は計算することは難しい.

ニュートン補間公式

関数 f(x) を差分商を用いて表した上式において最後の項を 0 とおいた次式をニュートン補間公式という.

  
  f(x) = f[x0] + (x - x0)f[x0, x1] + (x - x0)(x - x1)*f[x0, x1, x2]
       + ...
       + (x - x0)(x - x1)...(x - xn-1)*f[x0, x1, ..., xn]

 
この式は f(x) の近似式になっており, 0 とおいた最後の項は剰余項と呼ばれ近似誤差を表す.

ニュートン補間公式は, ラグランジュ補間公式を 1, (x – x₀), (x – x₀)(x – x₁), ・・・ について展開した式に相当する. すでに求めた補間式に新たに標本点を追加するとき, ラグランジュ補間公式では全部計算し直す必要があるが, ニュートン補間公式では新たな項を追加するだけでよい.

ニュートン補間公式の計算は次のようにして行うとよい.

まず, 与えられたデータについて差分商を下図のように左から右の順に計算しておく. 実際に補間値の計算に必要なのは対角部分だけであるが, 標本点を追加するときに一番下の行が必要になる.

  
  f[x0]
        ＼
  f[x1] ─ f[x0, x1]
        ＼             ＼
  f[x2] ─ f[x1, x2]   ─ f[x0, x1, x2]
        ＼             ＼                   ＼
    ・・・
        ＼             ＼                   ＼   ・・・   ＼
  f[xn] ─ f[xn-1, xn] ─ f[xn-2, xn-1, xn] ─   ・・・   ─ f[x0, ・・・, xn]

 
次に, 求めたい点 x それぞれについて, 多項式を計算する要領でニュートン補間公式により補間値を計算する.

標本点を追加するときは差分商の一番下の行を追加計算し, ニュートン補間公式の最後の項を補間値に加える.

上の数値実験 (1), (2) についてニュートン補間公式により計算すると多項式補間およびラグランジュ補間公式と同じ結果が得られる.

反復補間法

補間により x における関数値を計算するときに, はじめは少ない標本点を使って計算し, 次第に標本点の数を増やして計算を繰り返し目標精度に達したら(標本点を増やす前後の差が許容範囲に入ったら)収束とみなして終了する方法を反復補間法という. ニュートン補間公式は標本点の並び順に関係なく適用できるが, この場合は標本点をソートしておき x に近い順に追加していくのが収束の速さの点で有利である.

応用例としては, 与えられた数表を使って関数値を計算する場合が考えられる. XLPack の Fitlag を参照せよ.

エルミート補間

ラグランジュ補間は標本点において関数値が一致する近似式を与えた. エルミート補間では, 標本点において関数値の他に微分値も一致する近似式を与える.

区間 [a, b] 内の n + 1 個の相異なる点 x₀, x₁, …, x_n における関数 f(x) の値 f(x₀), f(x₁), …, f(x_n) および微分値 f'(x) の値 f'(x₀), f'(x₁), …, f'(x_n) がデータとして与えられているとき, これらの点における関数値および微分値がデータに一致するような近似式を考える.

条件が 2n + 2 個あるので, これを満足するためには 2n + 2 個の係数がある 2n + 1 次多項式 p_2n+1(x) が必要で, 次式を満たす.

  
  p2n+1(xk) = f(xk)  (k = 0 ～ n)
  p2n+1'(xk) = f'(xk)  (k = 0 ～ n)

 
p_2n+1(x) の係数を a₀, a₁, …, a_2n+1 とすると

  
  p2n+1(x) = Σaix2n+1-i (i = 0 ～ 2n + 1)
  p2n+1'(x) = Σ(2n+1-i)aix2n-i (i = 0 ～ 2n + 1)

 
であるから, 次の連立一次方程式を解くことにより近似多項式の係数を求めることができる.

  
  Va = f

 
ただし,

  
      ( x02n+1      x02n        ...  1 )       ( a0    )       ( f0  )
      ( (2n+1)x02n  (2n)x02n-1  ...  0 )       ( a1    )       ( f'0 )
      ( x12n+1      x12n        ...  1 )       ( a2    )       ( f1  )
  V = ( (2n+1)x12n  (2n)x12n-1  ...  0 ),  a = ( a3    ),  f = ( f'1 )
      (              ...               )       (  ...  )       ( ... )
      ( xn2n+1      xn2n        ...  1 )       ( a2n   )       ( fn  )
      ( (2n+1)xn2n  (2n)xn2n-1  ...  0 )       ( a2n+1 )       ( f'n )

 
ここで, V は連立一次方程式の係数行列, f は連立一次方程式の右辺ベクトル(各点における関数値 f_i = f(x_i) および微分値 f’_i = f'(x_i) を表す), a は連立一次方程式の解ベクトル(多項式の係数)となる.

ラグランジュ補間のときと同じように x₀, x₁, …, x_n が相異なる点であればこの連立一次方程式は解を持つ.

多項式補間のときと同様に連立一次方程式を解く計算法は方程式の条件数などの点から好ましくなく, ニュートン補間に似た以下の計算法が用いられる.

n + 1 個の標本点 x₀, x₁, …, x_n があるときに, j = 0 ～ n について

  
  u2j = xj
  u2j+1 = xj

 
とする.

差分商の計算を次のように変形する.

  
  f[ui] = f(ui) (i = 0 ～ 2n+1)
  f[ui, ui+1] = f'(ui) (i = 0, 2, ..., 2n)
  f[ui, ui+1] = (f[ui+1] - f[ui])/(ui+1 - ui) (i = 1, 3, ..., 2n-1)
  f[ui, ..., ui+k] = (f[ui+1, ..., ui+k] - f[ui, ..., ui+k-1])/(ui+k - ui) (k = 2 ～ 2n + 1, i = 0 ～ 2n + 1 - k)

 
関数 f(x) の補間値は次式で求められる.

  
  f(x) = f[u0] + (x - u0)f[u0, u1] + (x - u0)(x - u1)*f[u0, u1, u2]
       + ...
       + (x - u0)(x - u1)...(x - u2n+1)*f[u0, u1, ..., u2n+1]

数値実験 (3)

区間 [-1, 1] においてルンゲの関数をエルミート補間する.

  
  f(x) = 1/(1 + 25x2)

 
n = 4, 8, 12, 20 として n + 1 標本点を用いて補間を行う(多項式の次数はそれぞれ 9, 17, 25, 41 次である). 区間 [-1, 1] を n 等分して, 両端を含め等間隔に並んだ n + 1 点における関数値および微分値をデータとして使用する.

同じ標本点数のときにエルミート補間多項式の次数はラグランジュ補間のときの倍になるため精度はよい.

この時の微分値の補間の状況をプロットしてみると次のようになった.

比較のためにラグランジュ補間の場合の微分値をプロットしてみると次のようになった.

直交多項式補間

直交多項式系

区間 [a, b] および密度(重み)関数 w(x) が与えられているとき, 次式を満たす多項式系 { p_n(x) } を直交多項式系という.

  
  (pj, pk)w = ∫ pj(x)pk(x)w(x)dx [a, b] = 0 (j ≠ k)

 
(p_j, p_k)_w は p_j と p_k の重み付き内積である.

w(x) は [a, b] で連続で有限個の点で 0 になることを除いて w(x) > 0 とする. また, ∫ x^kw(x) dx [a, b] が k = 0, 1, 2, … について有界とする.

j = k のとき内積は正の値をとり, これを λ_k で表すことにする.

展開形の x^k の係数を μ_k で表すことにする. ただし, μ₀ は正の定数とする.

  
  pn(x) = μ0 + μ1x + μ2x2 + ・・・ + μnxn

直交多項式の性質

n 次直交多項式には次のような性質がある.

• 区間 [a, b] 内に n 個の相異なるゼロ点を持ち, すべて単根である.
• 次のような３項漸化式が存在する (p_-1(x) = 0 とする).

    
    pk(x) = (αkx + βk)pk-1(x) - γkpk-2(x)  (k = 1, 2, ... )
    ただし
      αk = μk/μk-1
      βk = -αk(xpk-1 . pk-1)/λk-1
      γk = αk(xpk-1 . pk-2)/λk-2 = μkμk-2λk-1/(μk-12λk-2)
  
  

• 次式(クリストッフェル・ダルブーの恒等式)が成り立つ.

    
    Σ (pk(x)pk(y))/λk (k = 0 ～ n-1) = μn-1/(μnλn-1) (pn(x)pn-1(y) - pn-1(x)pn(y))/(x - y)
  
  

   

代表的な直交多項式

以下のような直交多項式がよく知られている.

これらの直交多項式の関数値は漸化式を使って求めるとよい. それぞれの漸化式およびその導関数の漸化式は次のとおりである.

それぞれの関数形をグラフにすると次のようになる.

直交多項式補間

ラグランジュ補間において, 標本点を等間隔ではなく n + 1 次直交多項式 p_n+1(x) のゼロ点にとることを考える. 直交多項式の性質よりそのゼロ点はすべて単根であるから, このようにとることは可能である.

x₀, x₁, …, x_n が直交多項式 p_n+1(x) のゼロ点であるとする. クリストッフェル・ダルブーの恒等式において x = x_i, y = x_j とおくと, i ≠ j のときは

  
  Σ (pk(xi)pk(xj))/λk (k = 0 ～ n) = 0 (i ≠ j)

 
となり, i = j のときの値を 1/w_i とおくと

  

  1/wi = Σ (pk(xi))2/λk (k = 0 ～ n) > 0

 
である. すなわち, 次式が成り立つ.

  
  wi Σ (pk(xi)pk(xj))/λk (k = 0 ～ n) = δij

 
これは p_n+1(x) のゼロ点を標本点とするラグランジュ補間のラグランジュ補間係数となっている.

  
  f(x) = Σ (wi Σ pk(xi)pk(x)/λk (k = 0 ～ n)) f(xi) (i = 0 ～ n)

 
さらに, これは次のように直交多項式による補間になっている.

  
  f(x) = Σ ckpk(x) (k = 0 ～ n)
  ただし
    ck = 1/λk Σ wipk(xi)f(xi) (i = 0 ～ n)

同じ次数の補間多項式において, 補間の誤差をできるだけ小さくするための標本点のとり方を決める問題を最良近似問題と呼び, どのような関数にも適用できる公式はなく, 関数ごとに求める必要がある.

しかし, 標本点を直交多項式のゼロ点にとる(直交多項式補間する)と最良に近い結果が得られることがわかっていて, 数値積分などに応用されている.

以下にチェビシェフ補間の例を示す.

数値実験 (4)

区間 [-1, 1] において n = 4, 8, 12, 20 として(すなわち, 4, 8, 12, 20 次多項式を用いて)ルンゲの関数の補間を行う. ただし, n + 1 次チェビシェフ多項式のゼロ点における関数値をデータとして使用する. これはチェビシェフ補間と呼ばれる.

  
  xi = cos(π(2i + 1)/2n) (i = 0, 1, ..., n)

チェビシェフ補間もあらゆる関数についてうまくいくわけではないが, 区間 [-1, 1] で連続な 1 階導関数を持つ関数について n → ∞ のとき f(x) に収束することがわかっている.