補間法 - ビジネス

５．２グレゴリー・ニュートン（Gregory-Newton）の補間式
（１）導入
次のような点列に対する補間式を考える。
x1
f ( x1 )
x2
f ( x2 )
f ( x)  a1  a2 ( x  x1 )  a2 ( x  x1 )(x  x2 )  
 an1 ( x  x1 )(x  x2 )()(x  xn )

xn
f ( xn )
点列を代入すると，代入の繰返しで係数を求めることができる。
f ( x1 )  a1
f ( x2 )  a1  a2 ( x 2  x1 )
f ( x3 )  a1  a2 ( x3  x1 )  a2 ( x3  x1 )(x3  x2 )

f ( xn 1 )  a1  a2 ( xn 1  x1 )  a2 ( xn 1  x1 )(xn 1  x2 )  
 an 1 ( xn 1  x1 )(xn 1  x2 )()(xn 1  xn )
［参考］
概念的には，多項式補間を次のような三角行列に変換したことと同じ。
 n
 ( xn 1  x j )
 j 1



0

0


0

また
n 1
 ( xn1  x j )
j 1

0
0
0

( xn 1  x1 ) 1 an 1   f ( xn 1 )
     
 



 
 a3    f ( x3 ) 

( x3  x1 )(x3  x2 ) ( x3  x1 ) 1 
 a 2   f ( x2 ) 
 

0
( x2  x1 ) 1 
 a1   f ( x1 ) 

0
0
1

xi  xi  x j
であるから多項式補間より係数の絶対値の比率は一般に小さくなる。
（２）等間隔法①
表データは，独立変数を等間隔に並べていることが多いので，点間の幅を h
とすると，
h  x2  x1  x3  x2    xn1  xn
f ( x1 )  a1
f ( x1  h)  a1  a2 (h)
f ( x1  2h)  a1  a2 (2h)  a3 (2h)(h)
f ( x1  3h)  a1  a2 (3h)  a3 (3h)(2h)  a4 (3h)(2h)(h)

f ( x1  nh)  a1  na2 (nh)  a3 (nh)(n  1)h  a4 (nh)(n  1)h(n  2)h
 an 1 (nh)(n  1)h(n  2)h (h)
（２）等間隔法②
係数を求めると
a1  f ( x1 )
f ( x1  h)  f ( x1 )
a2 
h
f ( x1  2h)  2 f ( x1  h)  f ( x1 )
a3 
2
2h
f ( x1  3h)  3 f ( x1  2h)  3 f ( x1  h)  f ( x1 )
a4 
3
3  2h

（２）等間隔法③
次のような記法を導入
f ( x1 )  f ( x1  h)  f ( x1 )
 f ( x1 )  f ( x1  2h)  2 f ( x1  h)  f ( x1 )
2
 f ( x1 )  f ( x1  3h)  3 f ( x1  2h)  3 f ( x1  h)  f ( x1 )
3

コーヒーブレーク（Pascalの三角形を思い出そう）
さらに次のようにおく
１
x  x1  uh
１
１
１
１
２
５
１
３
３
６
４
１
１
１０
１
４
１０
１
５
１
（２）等間隔法④
整理すると
2 f ( x1 )
f ( x)  p(u )  f ( x1 )  f ( x1 )u 
u (u  1)
2!
3 f ( x1 )

u (u  1)(u  2) 
3!
n
 f ( x1 )

u (u  1)(u  2)()(u  (n  1))
n!
この式をグレゴリ・ニュートンの補間式と呼ぶ
計算すると（ただしＥｘｃｅｌで定義しよう）
f (0)  0.17365 0.0000 0.17365
（３）例
X(度） sin(X)
0
0.0000
2 f (0)  0.34202 2  0.17365 0.0000
10
0.17365
 0.00528
20
0.34202
30
0.50000
3 f (0)  0.50000 3  0.34202
0.64279
 3  0.17365 0.0000 0.00512 40
4 f (0)  0.64279 4  0.50000 6  0.34202
 4  0.17365 0.0000 0.000316
u  ( x  x1 ) h  (22  0) 10  2.2
0.00528
p (u )  0.0000 0.17365u 
u (u  1)
2
0.00512
0.000316

u (u  1)(u  2) 
u (u  1)(u  2)
3 2
4  3 2
 0.374606
（４）Ｅｘｃｅｌでの定義
ＶＢＡでのプログラム
①データ設定とボタンのＣｌｉｃｋイベントハンドラ
Private X(5) As Double
Private FY(5) As Double
Private UX As Double
Sub データ設定()
XX = 0: DX = 10
For i = 1 To 5
X(i) = XX: FY(i) = Sin(XX * 3.1415926 / 180)
XX = XX + DX
Next
UX = (22 - X(1)) / DX
End Sub
Sub ボタン1_Click()
データ設定
R = Gregory(FY, UX, 5)
MsgBox " 結果=" & R
End Sub
ＶＢＡでのプログラム
②組合せ数と階乗の計算
Function Comb(N, R) As Long '組合せ回数の計算
If R = 0 Or R = N Then
Comb = 1
ElseIf R = 1 Then
Comb = N
Else
Comb = Comb(N - 1, R - 1) + Comb(N - 1, R)
End If
End Function
Function Fact(N) As Double '階乗の計算
Dim D As Double
D = 1
For i = 2 To N
D = D * i
Next
Fact = D
End Function
ＶＢＡでのプログラム
③差分計算と補間計算
Function 差分計算(i, FY) As Double
A = 0: Sign = 1
For k = i To 2 Step -1
C = Comb(i - 1, k - 1)
A = A + C * Sign * FY(k)
Sign = -Sign
Next
差分計算 = A + Sign * FY(1)
End Function
Function 補間計算(i, UX) As Double
U = UX: T = 1
For k = 2 To i
T = T * U: U = U - 1
Next
補間計算 = T
End Function
ＶＢＡでのプログラム
④クレゴリ・ニュートンの補間式の計算
Function Gregory(FY, UX, N)
T = 0
For i = 1 To N
A = FY(1): ii=i
If i >= 2 Then A = A = 差分計算(i, FY) * 補間計算(i, UX)/ Fact(ii)
T = T + A
Next
Gregory = T
End Function
［補足］例題のプログラムの問題点
Function Comb2(N, R) As Long
Dim A() As Double
ReDim A(N)
NN = N: RR = R
If NN - RR < RR Then RR = N - R
If RR = 0 Then
Comb2 = 1
ElseIf RR = 1 Then
Comb2 = N
Else
For i = 2 To RR
A(i) = i + 1
Next
For i = 1 To NN - R - 1
A(1) = i + 2
For j = 2 To RR
A(j) = A(j - 1) + A(j)
Next
Next
Comb2 = A(RR)
End If
End Function
組合せ個数を計算するのに，
nC0=nCn=1, nC1=n,
nCr=n-1Cr-1+n-1Cr
の定義どおり処理しているが，
たとえば，10C5を求めるのに9C4と9C5を使い，
さらに9C5を求めるにも8C4を使う。
このようにまったく同じ計算が何回も出てくる。
左のプログラムでは， nC1から計算して，
結果を配列に保存することで，
無駄な計算を省いている。
課題１
以下のデータは，ある斜面の高さを100m間隔で測定した結果である。グレゴリーニ
ュートンのプサンプルログラムを用いて，250m位置の高さを求めよ。
位置(ｍ)
高さ(m)
0
185
100
186
200
188
300
190
400
191
500
192
課題２
高さｈと角度θを示す衛星の軌道データは以下のとおりである。
このデータを用いて，ｔ＝1,200秒の衛星位置（高さとθ）を推定せよ
。
なお，時間tは衛星が軌道に入った後の時間とする。
衛星
t(秒)
高さ h(mile)
θ(度)
0
100
0
300
107
22
600
160
42
900
230
62
高さ
θ
地球
課題３
バリスタ（非線形抵抗器）の電圧・電流を測定したところ以下のようになった。
このバリスタの特性を示す3次補間式を作れ。
E(ボルト)
I(アンペア)
0
0
10
0.36
20
0.49
30
0.62

Download Report