ニュートンの補間公式

 いくつかの点での関数値を与えて,そのような値を持つ関数の別な点での値を推測しようとすることがあります。このようなことを補間と言います。例えば,

x=0で,y=4であり,x=2で,y=8

であったとします。このとき,

x=1のとき,yはいくつか

と言う問題です。勿論実際には,この前提だけでyの値は決定できません。しかし,座標平面上の点(0,4)と(2,8)を直線で結んでその中点から(1,6)と推測するのは自然です。

 これは,この関数を1次関数と見て,y=2x+4と推測したことになります。

 この方法は多項式補間と言われるものの最も簡単な場合で,求める関数を最も次数の低い多項式で与えられる関数と考えたものです。

多項式補間

 多項式補間は次の事実に基づいています。

<x<x<・・・<xd−1<x を相異なるd+1個の数とし,
,y,・・・,yd−1,y,をd+1個の数とする。このとき,

=f(x), i=0,1,・・・,d

となる高々d次の多項式f(x)が唯一つ存在する。

 この定理は高校の数Tの範囲での定理ですが,確認の意味で簡単に証明をしましょう。


一意性(存在したとすると唯1つであること)

 f(x),g(x)がそれぞれ上の条件を満たす高々d次の多項式とします。このとき,F(x)=f(x) - g(x) と置くと,F(x)は高々d次の多項式で,方程式F(x)=0はd+1個の解x(i=0,・・・,d)を持ちます。

 一方,0でない高々d次の多項式は高々d個の解しか持ちません。従って,F(x)は零多項式でなくてはなりません。故に,f(x)=g(x)となります。


 存在

実際にこのような多項式が存在することは,具体的にそのような多項式を作ることで示します。この作り方は色々ありますが,ここではニュートンの補間公式と言われるものを紹介します。

ニュートンの補間公式

(x)を,f(x)=y,i=0,・・・,j を満たす高々j次の多項式とすると,fj+1(x)=y,i=0,・・・,j,j+1 を満たす高々j+1次の多項式fj+1(x)は適当な定数cj+1に対して次の式で与えられる。

j+1(x)=f(x)+cj+1(x−x)(x−x)・・・(x−x

 

 上のようにして決めた,fj+1(x)がfj+1(x)=y,i=0,・・・,jを満たすことは,上式の右辺第2項がx=x,i=0,・・・,jで0になることから直ちに分かります。また fj+1(xj+1)=yj+1 となることは,これを満たすようにcj+1を決めることで得られます。

 それには,上の式の右辺で,(x=xj+1を代入した式=yj+1)の方程式を解いてcj+1を決めれば求めるものが得られます。。

 上の操作を,定数多項式f(x)=yから始めて,j=dまで行えば求めるものが得られます。


 x=−2,x=0,x=2,x=5,y=−1,y=5,y=3,y=20

とします。

(x)=y=−1,

(x)=f(x)+c(x−x
    =−1+c(x+2)=−1+3(x+2),

  (ここで,-1+c(x+2)=y1より,c=3を求める。)

(x)=f(x)+c(x−x)(x−x
    =−1+3(x+2)+c(x+2)x=−1+ 3(x+2) - (x+2)x,

  (ここで,-1+3(2+2)+c(2+2)2=3より,c=−1を求める。)

(x)=f(x)+c(x−x)(x−x)(x−x
    =−1+ 3(x+2) - (x+2)x+c(x+2)x(x−2)
    =−1+ 3(x+2) - (x+2)x+(x+2)x(x−2)/3
  (ここで,-1+3(5+2)-(5+2)5+c(5+2)5(5−2)=3より,c=1/3を求める。)

と求める式が得られます。

ciの計算(差分法)

ニュートンの補間公式は求める多項式を

f(x)=c+c(x−x)+c(x−x)(x−x)+・・・
     +c(x−x)(x−x)・・・(x−xd−1

として,順次c,c,・・・cを求めていく方法です。筆算でこれらを求める場合は上のような方法で良い のですが,プログラムとして表す場合はもう少し体系的な取り扱いが必要です。上の方法を精密化してプログラムとして書けるようにしましょう。

まず,c=yは直ちに得られます。

1.cの計算:上式にx=xを代入することで,

=(y-y)/(x-x

が得られます。

2.cの計算:上式にx=xを代入することで,

={[(y-y)/(x-x)]-[(y-y)/(x-x)]}/(x-x

が得られます。この式はcに比べて複雑ですが,

y’=[(y-y)/(x-x)],y’=[(y-y)/(x-x)]

と置くと,

=(y’-y’)/(x-x

となり,cと良く似た形になります。


 このような形の式を差分と言います。cは一階の差分,cは差分の差分ということで2階の差分と言います。

 cを求めるために差分の定義をしましょう。

差分の定義

f[x]=f(x) と0階の差分を定め,以下帰納的に

f[x,xi+1,・・・,xi+j]={f[xi+1,xi+2,・・・,xi+j]-f[x,xi+1,・・・,xi+j−1]}/(xi+j-x

によって,j階の差分を定めます。

 この定義に従えば,上の式から,

=f[x,x

=f[x,x,x

となることが分かります。そして一般に,

i=0,1,2,・・・,dに対して

=f[x,x,・・・,x

が成立します。この事実を使ってプログラムを書きましょう。

プログラムの概要

 プログラムの概要を考えましょう。ニュートンの補間公式で得られる関数のグラフを表示することにします。

 多項式の次数dとd+1個の数とそれに対応する関数値をプログラムの中に書くことにします。このとき,

  1. 与えた点をグラフに表示し
  2. ニュートン補間した関数のグラフをその上に描く

 さて,以上の処理を実現するプログラムを書いて見ましょう。

ニュートンの補間公式による関数を求める以外は全く標準的です。

基本データの入力部分

 入力に必要なデータが比較的あるので,Input 文による入力ではなく,プログラムに書き込むことにします。データファイルを別に用意し,それから読み込むようにすることも可能ですが,ここでは簡単のため,プログラムに書き込んで使うことにします。

必要なデータは求める多項式の次数と,点と値との対応データです。


Public PDEG
PDEG = 7 :'多項式の次数
'----------------------
'データの個数=PDEG+1
'     x   y  
Data -8,  2
Data -5,  3
Data -3,  1
Data  0,  2
Data  2,  1
Data  5,  3
Data  8, -4
Data  9, 1
'----------------------

のようにすることにします。PDEG はPolynomial Degree のつもりです。上の場合,高々7次の多項式で補間すると指定してします。対応データは Data 文で書きます。対応データの個数は次数+1ですから,今の場合,8個のデータが必要です。 y のデータは同じでも構いませんが,xの値は全て異なる必要があります。

このデータは配列に読み込み使います。配列変数は

Public Nx(50),Ny(50),NCD(50)

と宣言します。Nx はx座標,Nyはy座標 NCD は上の式での Ciです。これを実際のデータに読み込むルーチンは次の通りです。


Sub ReadData
   For i=0 to PDEG
     Read Nx(i), Ny(i)
     NCD(i) = Ny(i)
   Next i
End Sub

 NCD(i)= Ny(i) としているのは,後の計算の初期設定のためです。

グラフの初期描画

 グラフィックの設定のため,グラフ画面の左下と右上の座標をそれぞれ,(Xmin,Ymin),(Xmax,Ymax) として指定することにします。


Xmin = -10 : Xmax = 10
Ymin = -10 : Ymax = 10

 開く画面の大きさはここでは400×400の正方形で開きます。 背景は紺,座標軸は水色にしました。グラフ画面の初期設定は次の通りです。


Sub GraphInit(Xmin,Xmax,Ymin,YMax)
   BackColor = "Navy"
   GScreen(400,400)
   MathGraph On
   Window (Xmin,Ymin)-(Xmax,Ymax)
   CLS 2
   ForeColor = "Cyan"
   LINE (Xmin,   0)-(Xmax,   0)    : 'x 軸
   LINE (   0,Ymin)-(   0,Ymax)    : 'y 軸
End Sub

 次に入力したデータをグラフ画面に表示します。


Sub DrawInitData
   ForeColor = "Yellow"
   For i = 0 to PDEG
      Circle (Nx(i),Ny(i)),0.3,,,,,F
   Next i
   ForeColor = "White"
End Sub

入力データを表す点を半径0.3 で黄色ので描きます。

f(x)の計算

 次にニュートン補間して得られる関数f(x) のグラフを描きましょう。

まず,Ci の計算をします。配列NCD(i)を Ciのこととします。上の差分法の結果を使って,順次階差を計算し,それをNCD(i)に入れることにします。


Sub CalcNCD
   For Ni = 1 to PDEG
      For Nj = PDEG to Ni step -1
         NCD(Nj) = (NCD(Nj) - NCD(Nj-1))/(Nx(Nj) - Nx(Nj-Ni))
      Next Nj
   Next Ni
End Sub

 NiがNi階の差分を計算する部分です。最終的にNCD(Nj)にはNj階の差分の結果が入るようになります。即ち,上のルーチンを実行すると,Ci=NCD(i) となります。

Ciが計算できたら,f(x) の計算はニュートン補間法の定義式をそのまま使います。


Function Newton(x)
   NewX = NCD(0)
   Prod = 1
   For Ni = 1 to PDEG
      Prod = Prod * (x - Nx(Ni-1))
      NewX = NewX + Prod * NCD(Ni)
   Next Ni
   Newton = NewX
End Function

 これで後は,この関数を使ってグラフを描きます。


Sub DrawGraph(Xmin,Xmax)
   X  = Xmin    
   Y  = Newton(X)
   h  = Xmax -Xmin
   FOR I = 1 TO 80
      X1 = X
      Y1 = Y
      X  = Xmin + h*I/80 
      Y  = Newton(X)
      LINE (X1,Y1)-(X,Y)
   NEXT I
End Sub
実行

 以上の部分に,実行部分を書けばプログラムは完成です。ここでは次の部分を加えました。

完成したプログラムNewtonIP.tbtここにあります。


   Cls : Print "Newton の補間公式"
   Call  ReadData
   Call  GraphInit(Xmin,Xmax,Ymin,YMax)
   Call  DrawInitData
   Call  CalcNCD
   Input "スタートします。リターンキーを押して下さい。"; Ch
   Call  DrawGraph(Xmin,Xmax)
END

実際に実行すると,まず,入力データが画面に表示されます。

リターンキーを押すと,これらの8個の点を通る7次式で補間され,次の画面になります。

まとめ

 今回作成したプログラムは, ニュートンの補間公式についてのある程度の知識が必要です。プログラムはあくまでも道具ですから,書こうとする対象についての知識は必須です。良いプログラマーはプログラム自身を知るだけでなく,当然ですが,書こうとする対象についても知識が必要です。

プログラム以外に得意な分野を持ちましょう。