|
与えられたN個の点
（ＳａｍｐｌｅＸ［Ｉ］，SampleY［Ｉ］），Ｉ＝１，２，・・・，Ｎ
を通るＮ－１次の多項式をラグランジュ補間（Lagrange interpolation)という。
あらかじめ大域変数として SampleX［１．．N］，SampleY［１．．N］を与えておくと，関数Interpolation(X)は，Xにおけるラグランジュ補間の値を求める。
いろいろなXの値について点（X，Interpolation(X)）をプロットすれば，最初に与えたN個の点を通る滑らかな曲線ができる。
場合によっては，ラグランジュ補間は奇妙に振動することがあるので，心配な場合はスプライン補間を使用する。
const N=5 (* 例えば *) (* 標本点の数 *)
var SampleX,，SampleY : array [1..N] of real; （* 標本点の座標 *)
function Interpolation(X : real) : real;
var I, J： integer；
T, Y： real；
begin
Y := 0;
for I := 1 to N do begin
T = SampleY[I];
for J <> 1 to N do
if J<>I then T := T * ( X - SampleX[J]) / (SampleX[I] -SampleX[J;
Y := Y+T
end：
Interolation := Y
end]
N個の点
（ＳａｍｐｌｅＸ［ｉ］，ＳａｍｐｌｅＹ［ｉ］）， i = 1，２，・・・，Ｎ
が与えられているとする。 但し，SampleX［１］<SampleX[2]<・・・<SampleX[N］
とする。
これらのN点を通り、区分的にXの3次式で表され，N-2個のつなぎ目で２次導関数まで連続であり，両端点I=1，Nで２次導関数が０になるような曲線の方程式を求めるのが，ここでいう３次スプライン補間である。 スプライン（弾性のある薄板）を曲げたような滑らかな曲線が得られるので，こう呼ばれる。
先ず，上述のSampleX［１．．N］，SampleY［１．．N］を与えて，手続きMakeSpline(SampleX，SampleY，Table)を呼び出す。 すると Table[１．．N］に各点での２次導関数を６で割った値が求められる。 SampleX［１．．N］，SampleY［１．．N］の値は変わらない。なお，この手続きの呼び出しは，最初に１回だけ行えば良い。
そうしてから，関数Interpolation（X，SampleX，SampleY，Table）を呼び出すと，与えられたｘ座標の値Xに対応する関数値（ｙ座標の値）が得られる。 このＸの値は，ＳｍａｐｌｅＸ［１］=<Ｘ<=ＳａｍｐｌｅＸ［Ｎ］の範囲になくても良い。
Ｘの値を少しずつ変えてＹ－Ｉｎｔｅｒｐｏｌａｔｉｏｎ（Ｘ，ＳａｍｐｌｅＸ，ＳａｍｐｌｅＹ，Ｔａｂｌｅ）の値を求め，（Ｘ，Ｙ）をプロットすれば最初に与えたＮ個の点を通る滑らかな曲線ができる。
なお，補間によって得られる曲線が周期関数で表されることが分かっている場合には，周期関数編を使うべきである。
const N = 10; (* 標本点の数 *)
type Vector = array [1..N] of real;
(* 補間のための表を作る *)
procedure MakeSpline( var SampleX, SampleY, Table: Vector);
const Yxx1 = 0; YxxN = 0; (* 両端点での第２次導関数 *)
var I: Integer;
Ｔ，Ｕ： real;
Ｈ，Ｄ： Vector;
begin
for I := １ to N - 1 do H[ I ] := SampleX[ I + 1 ] - SampleX[ I ];
for I := 2 to N - 1 do D[ I ] := 2 * (SampleX[ I + 1 ] - SampleX[ I -
1 ]);
U := （SampleY[ I + 1 ] - SampleY[ I ] ) / H[ I ];
for I := 2 to N - 1 do begin
T := (SampleY[ I + 1 ] - SampleY[ I ] ) / H[ I ] ;
Table[ I ] := T - U； Ｕ := Ｔ
end；
Ｔａｂｌｅ［１］ := Ｙｘｘ１ / 6: Table[ N ] := YxxN / 6;
Ｔａｂｌｅ［２］ := Ｔａｂｌｅ［２］ - Ｈ［１］＊Ｔａｂｌｅ［１］；
Ｔａｂｌｅ［Ｎ - １］ := Ｔａｂｌｅ［Ｎ - 1］ - Ｈ［Ｎ - 1] *Table[N];
for I := ２ to N - 2 do begin
T := H[ I ］ / D[ I ];
Ｔａｂｌｅ［ I + 1 ］ := Table［ I + 1 ］ - Ｔａｂｌｅ［ Ｉ ］*Ｔ；
D[ I + 1］ := D［ I +1 ] －H[ I ] * T
end；
for I := N - １ downto 2 do
Table[ I ] := (Table[ I ] - H[ I ]＊Ｔａｂｌｅ［ Ｉ ＋ １ ］）／Ｄ［ Ｉ ］
end;
(* 実際の補間を行う *)
function Interpolation(X: real; var SampleX, SampleY, Table:Vector) : real;
var I，J, K:integer;
T, H:real；
function F(X: real): real;
begin
F := （sqr(X) - 1) * X
end；
begin { Interpolation }
I := １； J := N;
while I < J do begin
K : = ( I + J ） div 2；
if SampleX[ K ] < X then I := k + 1 else J := K
end;
if I > 1 then I := I - 1;
H := SampleX[ i + 1 ] - SampleX[ I ];
T := (X - SampleX[ I ] ) / H;
Interpolation := T * SampleY[ I + 1 ] + ( 1 - T) * SampleY[ I ]
+ sqr( H ) * (F( T ) * Table[ I + 1 ] + F( 1 - T) * Table[ I ])
end;
スプライン補間
（ａ）補間多項式の表現
標本点がｘ１＜ｘ２＜・・・＜ｘｎと番号づけられているとする。Ｌａｇｒａｎｇｅ補間やＨｅｒｍｉｔｅ補間は，ひとつの多項式を用いた区間［ｘ１，ｘｎ］全域での近似であるが，ここでは，各小区間［ｘｋ－１，ｘｋ］で別々の多項式を用いた近似を考える。 一般に，有限個の分点で区切られた各小区分間上でｍ次多項式であるような関数を区分的ｍ次多項式，（ｍ－１）階までの導関数が全域で連続な区分的ｍ次多項式をｍ次スプラインと呼ぶ。 与えられたｎ個のデータ（ｘｋ，ｙｋ）（ｋ＝１，・・・，ｎ）に対して，補間条件
Ｓ（ｘｋ）＝ｙｋ （ｋ＝１，・・・，ｎ）
を満たすｍ次スプラインＳ（ｘ）による補間をｍ次スプライン補間という。
スプライン補間は，１９４６年 Ｉ，Ｊ．Ｓｃｈｏｅｎｂｅｒｇによって提案され，１９６０年代以後いろいろな分野に盛んに応用されるようになった。スプラインの名は，Ｓｃｈｏｅｎｂｅｒｇの論文にも言及されているように，滑らかな曲線を引く簡単な製図具に由来する。 この製図具によって描かれる滑らかな曲線が連続な２階導関数をもつ区分的３次多項式になること（工学の分野では連続梁の理論で古くから知られている事実の一般化として，Ｓｃｈｏｅｎｂｅｒｇは上記の数学的なスプラインの概念を導入したのである。
応用上最も頻繁に使用される（ｍ＝）３次スプライン補間Ｓ（ｘ）について述べる。 先ず，Ｓ（ｘ）の表式を導く。 ｘ＝ｘｋにおけるＳ（ｘ）の導関数値をｄｋとしておくと，Ｈｅｒｍｉｔｅ 補間の結果から，小区間［ｘｋ－１，ｘｋ］において
但し、ｈｋ＝ｘｋ－ｘｋ－１である。 この作り方から，Ｓ（ｘ）の１階導関数は連続である。 ２階導関数の連続性を考える。
|
