Web教材一覧数値解析
(BOK大区分:1 基礎理論、中区分:2 アルゴリズムとプログラミング、中区分:2 アルゴリズム)

補間(ラグランジュの補間公式、ニュートンの補間公式)

キーワード

数値解析、補間、ラグランジュの補間公式、ニュートンの補間公式、スプライン補間公式


補間とは

右図のように、n+1個の点(黒点)(x0,y0),(x1,y1),(x2,y2),・・・,(xn,yn)が与えられているとき、その点を通るなめらかな曲線をひき、x0~xnの間の値xにおける点(赤点)のyの値を求めることを補間といいます。
 n+1個の点を通る曲線は多様ですが、最も単純なのはn次の多項式、
   y=c0+c1x+c2x2+・・・+anxn
です。
 すなわち、補間とは、この係数 cr を算出してn次多項式 f(x) を決定し、与えられた x=x* における f(x*) を求めることです。

連立方程式による係数の決定法

n+1個の点を通るn次多項式の係数は、連立方程式により求めることができます。点 (xk,yk) を通ることは、
  yk=c0 + c1x + c2x2 + … + cnxn  (k=0~n)
で表せます。

4つの点(0,1), (2,5), (3,8), (5,10)がを通る3次式ならば、
(0, 1): 1 = c0 + 0*c1 + 0*c2 +  0*c3
(2, 5): 5 = c0 + 2*c1 + 4*c2 +  8*c3
(3, 8): 8 = c0 + 3*c1 + 9*c2 + 27*c3
(5,10): 10 = c0 + 5*c1 + 25*c2 + 125*c3

を解くことにより、
   c0 = 1, c1 = 0.133, c2 = 1.333, c0 = -0.2
となるので、3次式は、   y = 1 + 0.133x + 1.333x2 - 0.2x3
となります。

しかし、この方法では、点の数(次元)が大きくなると、解くための計算量が大きくなるだけでなく、誤差が累積する欠点があります。
 それを回避する代表的な方法に、ラグランジュの補間公式ニュートンの補間公式があります。

補間での留意点

n+1個の点をn次式で補間するとき、nが大きくなると、右図(赤線、点A)のように、実際とは異なるのではないかと思われる補間値になることがあります。特に、x* が x0~xn の外側の場合(外挿という)には、そうなることが多いのです。

このような欠点を回避するには、あえて全体を一つの式で表すのではなく、補間すべき周辺の数点だけを用いて、低次元の式(青線)で補間を行う(B点)ことがよく行われます。その代表的な方法にスプライン補間公式があります。しかし、この方法は数学的にも高度ですし、計算も複雑になるので、ここでは省略します。→参照:「スプライン補間」

なお、補間と似たような手法に最小二乗法があります。最小二乗法は、「すべての点を通る」という補間の原則を放棄して、「すべての点に近い」低次元の多項式で近似することにより、極端な補間値を避ける方法だということもできます。


1次式での補間

ラグランジュの補間多項式の前に、2点(x0,y0),(x1,y1)が与えられ、1次式により、x0~x1の間にあるxに対するyを求める場合を考えます。

2点(x0,y0),(x1,y1)を通る直線は、
   y= (y1-y0)/(x1-x0) * (x-x0) + y0
です。これを変形すると、次の式が得られます。
   y= (x-x1)/(x0-x1)*y0 + (x-x0)/(x1-x0)*y1

例えば、2点を(1,3)と(4,7)とし、x=3とすると、
  (x-x1)/(x0-x1)*y0=(3-4)/(1-4)*1= 1/3
  (x-x0)/(x1-x0)*y1=(3-1)/(4-1)*7=14/3
ですから、
  y= 1/3 + 14/3 = 5
になります。


ラグランジュの補間多項式

n次式に一般化して、
   y=c0y0+c1y1+c2y2+・・・+cnyn
としたときの、ckは、次の式になります。
  ck = (x-x0)/(xk-x0) * (x-x1)/(xk-x1) * … * (x-xk-1)/(xk-xk-1)
          ((x-xk)/(xk-xk)の項がないことに注意)
    * (x-xk+1)/(xk-xk+1) * … * (x-xn)/(xk-xn)

例えば2点、1次式のときはn=1ですから、上の式でx2以上のものがある項は不要です。それで、
  ck = (x-x0)/(xk-x0) * (x-x1)/(xk-x1)
となり、「(x-xk)/(xk-xk)の項がない」ことから
  c0 = (x-x1)/(x0-x1)
  c1 = (x-x0)/(x1-x0)
となるので、1次式、
  y = (x-x1)/(x0-x1)*y0 + (x-x0)/(x1-x0)*y1
が得られます。これは、前述の「1次式での補間」で示した結果と一致します。

数値例

4つの点(0,1), (2,5), (3,8), (5,10)が与えられ、それを通る3次式を求めます。
  x0=0,  x1=2,  x2=3,  x3=5
  y0=1,  y1=5,  y2=8,  y3=10

c0 =         (x-x1)/(x0-x1) * (x-x2)/(x0-x2) * (x-x3)/(x0-x3)
c1 = (x-x0)/(x1-x0)         * (x-x2)/(x1-x2) * (x-x3)/(x1-x3)
c2 = (x-x0)/(x2-x0) * (x-x1)/(x2-x1)         * (x-x3)/(x2-x3)
c3 = (x-x0)/(x3-x0) * (x-x1)/(x3-x1) * (x-x2)/(x3-x2)

y = c0y0 + c1y1 + c2y2 + c3y3

x=4のときは、
c0y0 = [   2/(-2) * 1/(-3) * (-1)/(-5)]* 1 = 0.067
c1y1 = [4/2     * 1/(-1) * (-1)/(-3)]* 5 = -3.333
c2y2 = [4/3 * 2/1       * (-1)/(-2)]* 8 = 10.667
c3y3 = [4/5 * 2/3  * 1/2        ]*10 = 2.667

y = f(4) = 0.067 - 3.333 + 10.667 + 2.667 = 10.07


ニュートンの補間公式

n+1個の点(x0,y0), (x1,y1), ・・・,(xn,yn)を通るn次式 f(x)を、
   f(x)=c0 + c1(x-x0) + c2(x-x0)(x-x1)
       + ・・・ + cn(x-x0)(x-x1)・・・(x-xn-1
の形で表します。
 そして、
   f(x0)=y0、f(x1)=y1、・・・、f(xn+1)=yn+1
の関係から、cの値を求めて、n次式 f(x) を決定します。

例示

n=3とし、4つの点を(0,1), (2,5), (3,8), (5,10)としましょう。
 上の公式は次のようになります。
   f(x)=c0 + c1(x-0) + c2(x-0)(x-2) + c3(x-0)(x-2)(x-3))
 そして、f(0)=y0=1, f(2)=y1=5, f(3)=y2=8, f(5)=y3=10 です。

これより、
   f(x)=1 + 2*(x-0) + 1/3*(x-0)(x-2) + (-1/5)*(x-0)(x-2)(x-3)
が得られます。

これから、x=4のときを求めれば
   f(4)=1 + 2*(4-0) + 1/3*(4-0)(4-2) + (-1/5)*(4-0)(4-2)(4-3)= 10.07
となります。
 これは、ラグランジュの補間公式での結果と一致します。それは当然のことで、両者の数値例での見かけの式は異なりますが、 y=c0 + c1x + c2x2 + c3x3 の形式に整理すれば、「連立方程式」で求めた
  y = 1 + 0.133x + 1.333x2 - 0.2x3
になっているのです。

差分による説明

rは、次の差分表を作ると、簡単にに計算できます。

n次式の係数 cr は、cr=dr,0になります。
 そして、dr,k は、次の式で計算できます。
  dr,k=(dr-1,k+1-dr-1,k)/(xr+k-xk)   ただし、r=1~n  k=0~n-r  d0,k=yk

ここで、d1,0は、 (x0,y0) と (x1,y1) の傾きですから、f(x) を微分したものになります(階差分商といいます)。d1,1、d1,2 も同様です。
 また、d2,0 は、(x0,d1,0) と (x2,d1,1) の傾きだといえるので、d2,0は、f(x) を2回微分したものだといえます。d3,0は、3回微分したものになります。
 すなわち、dr,k は、点 (xk,yk)を起点とするr階の階差分商であり、ニュートンの補間法の公式は、階差分商を係数にした多項式だといえます。


ラグランジュの補間とニュートンの補間の比較


計算プログラム

次数n=(データの組数は0~nのn+1組)
i 0123456789
x[i]
y[i]
補間値:xmin xmax キザミ幅 dx=
   

数値解析へ