数式処理ソフト　DERIVE（デライブ） de ドライブ

８２．補間多項式

（1）関数近似

「数式処理ソフト DERIVE（デライブ）では、ここ何回か、常微分方程式の数値解法とその応用を学んできたんじゃが、第81回で予告した「アダムスの陽的公式」に話を進めていくために、予備的な知識が必要になった。
　そこで、今回は、「関数近似」という大きな項目のなかの基礎的な内容である「補間多項式」について、調べて見よう」

「関数近似というのは、特定の関数を他の関数で近似する、ということね」

「そうじゃな。
　手元の「計算機による数値計算法」（一松信・戸川隼人 著：新曜社：1975年2月初版）によれば、
　関数近似とは、『（何らかの形で）与えられた関数 y(x) に対して，それを（なんらかの意味で）よく近似する簡単な関数－具体的には多項式 p(x)を求めること』と説明されている。
　与えられた関数が「式」の形で与えられている場合もあるし、表で、すなわち、とびとびのxに対するyの値のみしか与えられていない場合もある。
　また、常微分方程式や偏微分方程式のような方程式（と初期条件・境界条件）で与えられる未知の関数の場合もある。
　すでにわしらが知っているものとしては、テイラー展開やフーリエ級数、複素関数に対するローラン展開などが関数近似と言える。
　だから、一口に「関数近似」といっても、非常に幅広い分野が含まれることになるのう」

「そのうちで、今回は、補間多項式を扱うのね」

「そうじゃな。
　上で述べたように、常微分方程式の数値解法の多段階法の一つであるアダムスの陽的公式を説明するには、とびとびの（離散的な）xとyの組を近似する補間法の知識が必要になる。
　なお、コンピューターの発達により、「数表」に載っていない点の値を前後の値から計算するという本来の補間の重要性は、著しく、低下した。
　たとえば、「岩波数学辞典」でも第3版までは、三角関数等に対する近似式が掲載されていたんじゃが、第4版（2007年3月）では、削除されている」

「前に取り上げた「数表」（吉田洋一・吉田正夫 編：培風館：1972年4月第19刷）でも、最初は、「常用対数表」から始まっているもの。
　基本的な関数は、その都度、組み込み関数なりを使って計算させる方がカンタン。
　この30年ぐらいの間の変化は、ものすごく大きかったのね」

「ま、そういうことじゃ。
　次節では、「ラグランジュ多項式」を取り上げよう」

（2）ラグランジュの補間多項式

「離散的変数の値に対して関数の値が与えられている場合じゃ。（関数が簡単な式表示を持たない場合）
　今、「計算機による数値計算法」にならって、ラグランジュの補間多項式（以下「ラグランジュ多項式」という）を取り上げよう。
　ここでは、添え字を書くのが不便なので、変数 x、をx(0)、x(1)、・・x(n)、あるいは、x0、x1、x2 などと略記する。
　また、x(n)に対応する関数値 y(x(0))をy(0)、あるいは、y0と略記する。
　さて、自然に思い浮かぶのが各点を通る関数として、xの多項式を当てはめることじゃ。
　x=x(0)～x(n)までの間で、関数の値は、y(0)、・・・、y(n)とする。
　f(x)=Σa(k)x^k、ただし、k=0～n、とすると、
　n+1個の変数と関数の組を与えて、係数 a(k) を決めればよい。
　同書に次のような例題が載っている。

x	0	1	2	3	4	5
y	1.0	0.9	0.8	0.2	0.1	0.0

上の表のデータを元に、f(x)=Σa(k)x^k、のa(k)を求めよ」

「変数 x の値を、ベクトル u ≔ [0, 1, 2, 3, 4, 5]、関数 y の値を、v ≔ [1, 0.9, 0.8, 0.2, 0.1, 0]、で表す。
　また、係数 a(k) を、a ≔ [a0, a1, a2, a3, a4, a5])、で表す。
　近似関数は、f(x)：=∑(a↓(k + 1)·x^k, k, 0, 5)、なので、
　これらの関係は、f(u↓j) = v↓j、j=1～6、である。
　従って、数列の作成により、VECTOR(f(u↓j) = v↓j, j, 1, 6, 1)、として6元連立方程式を表示する。
　次に連立方程式の解法で、a(k)を解くと、
　DERIVEでは、Solve(関係式、a)とすると、（関係式は、VECTOR(f(u↓j) = v↓j, j, 1, 6, 1)の式番号を指定）、
　[a0 = 1 ∧ a1 = － 149/120 ∧ a2 = 35/16 ∧ a3 = － 4/3 ∧ a4 = 5/16 ∧ a5 = － 1/40]、と係数ベクトルが求まる。
　これにより、近似関数は、
　f(x)=－ x^5/40 + 5·x^4/16 － 4·x^3/3 + 35·x^2/16 － 149·x/120 + 1、である」

「ご苦労さん。
　誤解のないように、言うておくと、このデータに次数の高い近似曲線が適切でないことは、ほぼ明らかじゃ。
　「計算機のための数値計算法」では、実験データが与えられたときは、まず、グラフを描いてみることを勧めている。
　
　で、そのグラフが上の図じゃ。
　同書にあるように、自然な近似曲線は、相互の点を直線で結んだものじゃろう。（上図では、破線で表示している）
　ここに、5次曲線を当てはめるのは、不自然だな。
　ちなみに、ともちゃんが求めてくれた近似曲線をオーバーレイすると下図のようになる。
　
　他に、何も根拠がないならば、このデータに5次近似曲線（ピンク）を当てはめるのは無理があるのう。
　一般に近似曲線の次数が大きくなれば、屈曲が多くなり、与えられた点では、指定した値をとるが、それ以外の点では、かけ離れた値をとる関数になり勝ちじゃ。
　ここでは、そういう知識を持った上で、あえて、このような近似曲線を求めてみたわけじゃ。
　ともちゃんが求めたくれた方法では、近似曲線を求めるのに（手計算の頃は）、手間がかかった。
　そこで、登場するのが、「ラグランジュの補間多項式」というわけじゃ。
　最初に言うておくと、求まる多項式は、ともちゃんが求めてくれたものと同じじゃ。
　ただ、連立方程式を立てて、それを解くという手間が省ける。
　Fn(i，n)=Π（j=0～n、ただし、j<>i）（（x－x(j)/(x(i)－x(j)）、として、
　近似関数を、f(x)=Σ（i=0～n）y(i)× Fn(i，n)、と求めるのが「ラグランジュ多項式」の方法じゃ」

「Fn関数は、積の関数だけど、ただし書き、があるのね。
　こういう関数をDERIVEで定義するのは、いつものことながら、迷うわ。
　Basic言語なら、Fn(i，n)は、
　Fn=1
　　for j=0 to n
　　　if j<>i then Fn=Fn*（x－x(j)/(x(i)－x(j)）
　　next j
　のようなループ処理で、まったくOKなんだけどね。
　DERIVEでは、For ループが自由には書きにくい。
　単純なループ処理なら、Iterates関数 または Iterate関数でいいけど、if のような命令をそのままでは、入れられない。
　そこで再帰関数で定義するということになる。
　試しに、Fnn関数を、再帰的に、
　Fnn(i, n):= IF(n = －1, 1, IF(n = i, Fnn(i, n － 1), Fnn(i, n － 1)·((x － u↓(n + 1))/(u↓(i + 1) － u↓(n + 1)))))、と定義すれば、
　一応、目的は達せられる。
　こうして、f(x)=Σ（i=0～n）v↓(i+1)× Fnn(i，n)、から、f(x)=－ x^5/40 + 5·x^4/16 － 4·x^3/3 + 35·x^2/16 － 149·x/120 + 1、となる」

「わしは、中間に、H(i, j):= IF(i = j, 1, (x － u↓(j + 1))/(u↓(i + 1) － u↓(j + 1)))、のような関数を定義した上で、
　Fn(i, n):= ∏(H(i, j), j, 0, n)、と定義してみた。
　これなら、迷いは少ない。（Πは、DERIVEの積の関数）
　再帰処理で迷うのは、逆向きに考えるのを強いられる点じゃないのかな。
　ともちゃんの再帰関数、Fnnの処理を模擬的に書くと、
　もし、n=－1ならば、Fnnを1と置く。
　nが－1でないならば、
　　　もし、n=i ならば、
　　　　　　Fnn(n－1)を呼び出す　
　　　もし、i<>n ならば、
　　　　　　Fnn(n－1)を呼び出して、それに(x － u↓(n + 1))/(u↓(i + 1) － u↓(n + 1))、を掛ける、
　というようになっている。
　今回の場合は、n=5なので、たとえば、i=1の場合、Fnn(1，5)は、
　n がループ処理の場合の j の役割も果たす。
　nは、n=5～0に至るまでは、n=i 以外では、
　Fnn(i，n－1)×（j=nの時の乗数）を計算。
　n=iの場合は、Fnn(i，n－1)×1。
　とはいえ、いずれにしても、Fnn(i，n－1) は不明じゃ。
　そこで、Fnn(i，n－1)の計算が背後で行われる。
　というように、n=0までは、バックグラウンドで一つ手前のFnnを呼び出していく処理が行われる。
　そして、Fnn(i，－1)が呼び出されると、ようやく、Fnn(i，－1)=1、という具体的な値が決められているので、そこから、前へ前と戻っていき、Fnn(i，n)が計算できる、ということになる。
　だから、回りくどい点があるのは、否めないし、なかなかなじめないのも分かるのう」

「デバッグもしにくいしね。
　なので、おじさんのアイデアをお借りして、再帰呼び出しを利用しないで、ラグランジュ多項式を与えるユーザー定義関数、Lg_user（u，v，x）、を作ったよ。
　　Lg_user(u, v, x) ≔
　　　　PROG(n_ ≔ DIMENSION(u) － 1,
　　　　　　　　 Hs ≔ IF(i_ = j_, 1, (x － u↓(j_ + 1))/(u↓(i_ + 1) － u↓(j_ + 1))),
　　　　　　　　 ∑(v↓(k_ + 1)·∏(LIM(Hs, [i_, j_], [k_, l_]), l_, 0, n_), k_, 0, n_))、
　DERIVE画面では、下図のようになっている。
　

　ところで、第69回「方程式の数値解法（4）（ハレー法、DERIVEのPROGとRETURN）」に出てきたんだけど、
　ユーザー定義関数の中で、新しい関数を定義したとき、2つの制限があるのよ。
　一つは、新しい関数の定義で、引数を指定するための括弧を書いてはいけないこと。
　つまり、Lg_userの中では、Hs がそれよ。
　これをHs(i，j):=定義内容、と書いてはいけないの。
　Hs:=定義内容、としなくてはいけない。
　もう一つは、新しく定義した関数、Hs、の呼び出し時にも、Hsに引数を与える括弧を付けてはいけない。
　じゃ、Hsの定義式の変数、i_、j_、には、どうやって値をセットするのか、という疑問が出ると思うんだけど、
　それを可能にするのが極限を求める関数、Lim 関数なの。
　LIM(Hs, [i_, j_], [k_, l_])、とすることで、Hs(k，l) という意味になるのね。
　しかし、Hsの定義には、u や x も書いてある。
　それは、Lim関数で与えなくてもいいのかという新たな疑問がわくと思う。
　これは、ユーザー定義関数を呼び出す際に与える引数が一定の値を持ち変化しない場合は、OK のようね。
　Lg_user([0, 1, 2, 3, 4, 5], [1, 0.9, 0.8, 0.2, 0.1, 0], z)、とすれば、
　(5 － z)·(6·z^4 － 45·z^3 + 95·z^2 － 50·z + 48)/240、を返してくれる。
　また、
　a ≔ [0, 1, 2, 3, 4, 5]、b ≔ [1, 0.9, 0.8, 0.2, 0.1, 0]、と外部で定義していても、
　Lg_user(a, b, z)、と計算させれば、同じ結果を返してくれる」

「ようできましたな。
　ラグランジュ多項式を与えるためのDERIVEの組み込み関数は、ないのかなとも思ったがの。
　ヘルプを見た限りでは、見つけられなかった。
　ともちゃんが書いてくれた、ユーザー定義関数内で新規に関数を定義して使用するときの制限に関する補足は、有り難い。
　ついでに言うておくと、Prog 命令は、カンマで区切られた命令を順次実行するためのものじゃ。
　また、Dimension(ベクトル)は、引数のベクトルの次元を求める関数じゃな。
　次節では、ラグランジュ多項式を差分で表現する「ニュートンの前進補間公式」について調べて見よう」　　

（3）ニュートンの補間公式

「前節のラグランジュ多項式の計算を差分で表現したものじゃ。
　ラグランジュ多項式の場合、分点が増えると、多項式の係数のすべてに影響があるため、最初から、計算し直さないといけない。
　たとえば、Lg_user([0,1,2],[1,0.9,0.8],x)、と3点の場合、
　多項式、1 － x/10、が得られる。
　ここで、分点が一つ増えて、Lg_user([0,1,2,3],[1,0.9,0.8,0.2],x) となると、
　多項式、－ x^3/12 + x^2/4 － 4·x/15 + 1、が得られる。
　このとき、前者の結果を利用することはできずに、新規に計算し直しとなる。
　差分を利用する場合は、計算量が減るのじゃな。
　まあ、現代では、この違いにぴんとこないかもしれん。
　このような多項式の次数は、せいぜい、10以下じゃろうから、その都度計算し直しても、別段のことはない。
　しかし、手計算の時代は、この違いは、甚だしく大きかったに違いない。
　「（前進）差分演算子」をΔで表す。
　aを任意の配列、kを整数としたとき、
　Δa(k)=a(k+1)－a(k)、じゃ。
　ついでに、あと出てくる「後退差分演算子」を▽（ナブラ）、「中心差分演算子」をδであらわす。
　これらは、▽a(k)=a(k)－a(k－1)、及び、δa(k)=a(k＋1/2)－a(k－1/2)、で定義されるものじゃ」

「前進差分演算子と後退差分演算子は、Δa(k－1)=▽a(k)、の関係にあるのね。
　また、以前に出た、「ずらし演算子」、E、を使えば、Δ=E－1、である。
　Δ^m a(k)=（E－1）^m a(k)=Σ（j=0～m）mCj E^j a(k)×(－1)^(m－j)=Σ（j=0～m）mCj ×(－1)^(m－j)a(k+j)、
　これは、Δ^m a(k)=(－1)^m ΣmCj ×(－1)^j ×a(k+j)、を意味する。
　また、
　a(k+m)=E^m a(k)=（1＋Δ）^m a(k)=Σ（j=0～m）mCj Δ^j a(k)、である。
　ここで、mCj は、二項係数。
　さらに、▽=1－E^－1、なので、▽=Δ×E^－1、とも書ける。
　実際、左辺をa(k)に作用させると、▽a(k)=a(k)－a(k－1)、
　また、右辺を作用させると、Δ×E^－1 a(k)=Δa(k－1)=a(k)－a(k－1)、なので、この関係は、成立している。
　更に、上式に右からEを掛けると、▽×E=Δ、でもある。
　ところで、中心差分演算子 δ は、初お目見えという感じなんだけど、どういうものかしら？」

「配列が、分点でのみ定義されている場合は、このままでは意味を持たないが、配列が分点以外でも意味を持つ場合、
　差分間隔をhと書いたとき、δa(x=kh)=a(x＋h/2)－a(x－h/2)、と定義されるものなんじゃな。
　その場合、a(x)が分かっていても、a(x＋h/2)は、そのままでは、求められないので何らかの近似により、求める必要はある。
　なお、δa(2k)=a(2k+1)－a(2k－1)、ではある。
　さて、f(x)=y(0)＋(y(1)－y(0))(x－x(0))/h、と置いてみると、
　f(0)=y(0)、f(1)=y(1)、となる。
　差分を使って書くと、f(x)=y(0)＋Δy(0)×(x－x(0))/h、
　これを拡張して、
　f(x)=Δ^0 y(0)＋Δy(0)×(x－x(0))/h＋Δ^2 y(0)×(x－x(0))(x－x(1))/2h^2＋Δ^3 y(0)×(x－x(0))(x－x(1))(x－x(2))/6h^3＋・・
　　＋Δ^n y(0)×(x－x(0))(x－x(1))(x－x(2))・・(x－x(n－1))/n! h^n、を定義する。
　ただし、Δ^0=1、とする。
　このf(x)は、結果として、前節の「ラグランジュ多項式」に等しい」

「それが、「ニュートンの（前進差分）補間公式」と言われるものなのね。
　k番目の一般項は、k=1～nと考えると、
　Δ^k y(0)×(1/(h^k k !))Π（j=1～n）（x－x(j－1)）、となるので、
　f(x)=y(0)＋Σ（k=1～n）Δ^k y(0)×(1/(h^k k !))Π（j=1～n）（x－x(j－1)）、
　DERIVEで関数を定義してみた。
　NewtonHokan_user(u，v，x)、がそれ。
　NewtonHokan_user(u, v, x):=
　　PROG(n_ ≔ DIMENSION(u) － 1,
　　　　　　LIM(v_↓1 + ∑(∑(COMB(k_, j_)·(－1)^(k_ － j_)·v_↓(j_ + 1), j_, 0, k_)·∏(x_ － u_↓(l_ + 1), l_, 0, k_ － 1)·(1/h^k_·(1/k_!)), k_, 1, n_), [x_, u_, v_], [x, u, v]))、
　この関数では、すでに記載のように、Δ^k y(0)、を、（E－1)^k y(0)、と考えて展開し、
　Δ^k y(0)=Σ（j=0～k）(－1)^(k－j) y(j＋1)、と計算している。
　実際に前節の例題で実験してみると、
　－ x^5/40 + 5·x^4/16 － 4·x^3/3 + 35·x^2/16 － 149·x/120 + 1、となり、一致している。
　f(x)=Δ^0 y(0)＋Δy(0)×(x－x(0))/h＋Δ^2 y(0)×(x－x(0))(x－x(1))/2h^2＋・・を見ると、
　x=x(0)でのy=y(x)のテイラー展開と類似していることには感心するけど、特に便利とは思えない気がする」

「確かに、コンピューターを使う場合（今は専らこうなっているが）、補間多項式を求めるために、上の定義式のように差分を使って計算するのは、面倒だ。
　手計算の時代は、分点が増えても、項を一つ増やすだけで計算を続行できる点が便利じゃったと思われるがの。
　もっとも、非常に多くの計算を行う場合には、差分を使うと乗除算の回数を減らす効果はある。（ただし、差分のための引き算で桁落ちなどが起こりやすいというデメリットはある）
　なお前述の式で、(x－x(0))/h を t と置き、x(j)=x(0)＋j h、を使うと、
　ニュートンの前進差分補間公式を簡略に表記することができ、
　F(x(0)＋h t)=Σ（j=0～n）Δ^j y(0)×（一般化二項係数 t｜j ）、とまとめられる。
　t｜j で表している「一般化二項係数」は、t(t－1)(t－2)・・・(t－j＋1)/j ! 、を意味する。
　特に、t｜0=1、と約束する。（数学では、一般化二項係数は、横線－を挟んで分数のように記載されることが多いが紛らわしいのでここでは縦線｜で表す）
　また、「後退差分演算子」を使った場合は、「ニュートンの後退（差分）補間公式」と呼ばれ、
　F(x(0)＋ht)=Σ（j=0～n）▽＾ｊ y(0)×（t＋j－1｜j) と書くことができる。
　次節では、他の補間公式について、簡単に触れることにしよう」

（4）ガウスの補間公式など

「前にも述べたように、手計算の時代は、「数表」の引数と引数との間の値を補間公式で求める必要があった訳じゃ。
　しかし、現代では、実験から得られるデータは、別として、関数の分点間の値を補間公式で計算する必要性は、低下している。　
　また、実験データは、最小二乗法等により、直線、あるいは、理論から推定できる低次の曲線に当てはめるのが妥当じゃろう。
　従って、以下に記載する、いくつかの補間公式は、「歴史的」、「教育的」には、意味があるが、補間公式としては、出番は少ないじゃろう。
　ただし、数値微分には、役立つと思われるし、数値積分の公式の導出にも役立つ。
　ところで、（3）節で出た、「ニュートンの前進差分補間公式」について、ちょっと、補足しておこう。
　j=0～n、のn＋1個の点を通るラグランジュ多項式には、「基底関数」と呼ばれる、関数が使われている。
　たとえば、n=2では、以下の2次式を見ると、
　y=y0×1＋(1/h)Δy(0)×(x－x0)＋(1/2h^2)Δ^2 y(0)(x－x0)(x－x1)、
　ここで、赤字で示したものが基底関数じゃな。
　注意すべき点は、(x－x2)は、含まれていない点じゃ。（x=x2を代入すると、y2が返るので、［x2，y2］を通ることは間違いない。）
　n＋1個の点を通る曲線は、高々、n次式で、これは、当然じゃが、次に説明するガウスの前進補間公式などでは、基底関数としてどの点を含ませるか多項式が違ってくるので注意じゃ。
　さて、x(－n)～x(n)、の2n＋1個の点を通る補間曲線を考える。
　では、n=1の場合を考えてみて欲しい」

「えーと、ニュートンの前進差分補間公式、F(t)=Σ（j=0～n）Δ^j y(0) （t｜j ）、
　ここで、始点をx(0)ではなくて、x(－1)に移動して、n=2の場合を考えると、t=（x－x(－1)）/h、として、
　F(t)=y(－1)＋Δy(－1) t ＋Δ^2 y(－1) t (t－1)、
　これは、［x(－1)，y(－1)］、［x0，y0］、［x1，y1］を通る、となるんじゃないかしら」

「うん。
　ともちゃんの解答も、間違ってはいないのじゃが、「計算機のための数値計算法」によると、
　「ガウスの前進補間公式」として紹介されているのがこちらだな。
　y0＋Δy0 t ＋Δ^2y(－1) t (t(－1)、が2次式の場合じゃ」

「えー。
　違うじゃない。なんで」

「わしも、最初、戸惑ったんだがな。
　同書に、「分点をx(0)、x(1)、x(－1)、x(2)、x(－2)、の順でとると、・・」と説明されている。
　分点間をつなぐ曲線を求めるのに、順番は、特に関係がない、と思うかもしれん。
　この点をはっきりとさせるために、3つの関数を下記のように定める。
　yA(x)=a0＋a1(x－x0)＋a2(x－x0)(x－x1)、
　yB(x)=a0＋a1(x－x(－1))＋a2(x－x(－1))(x－x0)、
　yC(x)=a0＋a1(x－x(0))＋a2(x－x0)(x－x(－1))、
　これらの3つの関数は、いずれも、［x(－1)，y(－1)］、［x0，y0］、［x1，y1］、を通るとする。
　DERIVEでは、ベクトル u、v を、
　u ≔ [xm2, xm1, x0, xp1, xp2]、v ≔ [ym2, ym1, y0, yp1, yp2]、とする。
　yA(u↓k)=v↓k、
　yB(u↓k)=v↓k、
　yC(u↓k)=v↓k、
　ここで、k=2～4、として、計算させて、xm1=x0－h、xp1=x0＋h、とすれば、
　yA(x)=y0＋(1/h)Δy0 (x－x0)＋(1/2h^2)Δ^2 y(－1)(x－x0)(x－x1)、
　yB(x)=y(－1)＋(1/h)Δy(－1)(x－x(－1))＋(1/2h^2)Δ^2y(－1)(x－x(－1))(x－x0)、
　yC(x)=y0＋(1/h)Δy(－1)(x－x(0))＋(1/2h^2)Δ^2y(－1)(x－x0)(x－x(－1))、
　これを見ると、ともちゃんの解答は、yB(x) じゃな」

「なるほど。
　つまり、同書に、わざわざ、「順番」が書いてあるのは、基底関数の選択の順番のことなのね」

「そういうことじゃ。
　［x(－1)，y(－1)］、［x0，y0］、［x1，y1］を通るというのは、3つとも同じじゃ。
　だから、関数の表現の仕方の違いではある。
　「計算機のための数値計算法」に書かれているのは、［x0，y0］－＞［x1，y1］－＞［x(－1)，y(－1)］、の順に分点をとる。
　そして、基底関数も、この順に、すなわち、
　(x－x0)、(x－x0)(x－x1)、(x－x0)(x－x1)(x－x(－1))、と増やしていくと言うことなんじゃな。
　従って、模式的に書くと、基底関数の取り方は、
　yA(x)は、ゼロ－プラス1－マイナス1－プラス2－マイナス2・・－プラスn、（x－x(－n) は基底関数に含まれない）、
　yB(x)は、マイナスn－マイナスn－1－・・ゼロ－プラス1－プラス2－・・プラスn－1、（x－x(n)は 基底関数に含まれない）、
　yX(x)は、ゼロ－マイナス1－プラス1－マイナス2－プラス2－・・マイナスn、（x－x(n))は、基底関数に含まれない）、となる。
　では、yA(x)の例に準じて、3次式を考えてみて欲しい」

「yA(x)=a0＋a1(x－x0)＋a2(x－x0)(x－x1)＋a3(x－x0)(x－x1)(x－x(－1))、として、
　DERIVEで計算してみると、
　yA(x)=y0＋(1/h)Δy0(x－x0)＋(1/2h^2)Δ^2y(－1)(x－x0)(x－x1)＋(1/6h^3)Δ^3y(－1)(x－x0)(x－x1)(x－x(－1))、
　（x－x0)/h=t、で表すと、
　yA(x)=y0＋Δy0 (t｜1)＋Δ^2y(－1)(t｜2)＋Δ^3y(－1)(t＋1｜3)、となる」

「4次式の場合は、ゼロ－プラス1－マイナス1－プラス2、の順に基底関数を構成する。
　yA(x)=a0＋a1(x－x0)＋a2(x－x0)(x－x1)＋a3(x－x0)(x－x1)(x－x(－1))＋a4(x－x0)(x－x1)(x－x(－1))(x－x2)、
　と置き、［x0，y0］、［x1，y1］、［x(－1)，y(－1)］、［x2，y2］、［x(－2)，y(－2)］を通る条件により、
　[a0 = y0 ∧ a1 = (yp1 － y0)/h ∧ a2 = － (2·y0 － ym1 － yp1)/(2·h^2) ∧ a3 = (3·y0 － ym1 － 3·yp1 + yp2)/(6·h^3) ∧ a4 = (6·y0 － 4·ym1 + ym2 － 4·yp1 + yp2)/(24·h^4)]、が得られる。
　整理すると、
　yA(x)=y0＋Δy(0)（t｜1）＋Δ^2 y(－1)（t｜2）＋Δ^3 y(－1)（t+1｜3）＋Δ^4 y(－2)（t+1｜4）、である。
　なお、ここで、t=（x－x0）/h、である。
　ともちゃん、（2）節の例で確認してみてくれんかの」

「その4次式では、
　yA(x)=y0＋Δy(0)（t｜1）＋Δ^2 y(－1)（t｜2）＋Δ^3 y(－1)（t+1｜3）＋Δ^4 y(－2)（t+1｜4）、となるので、
　下の表を参照して、［x(2)，y(2)］ を中心にして、係数を計算した。

t=(x－x2)/h、h=1、として、
　yA(t)=0.8＋（－0.6）t＋（－0.5）t(t－1)/2＋（1）(t+1)t(t－1)/6＋（1.5）(t+1)t(t－1)(t－2)/24、が得られる。
　下の図では、中心は、［x(2)，y(2)］。横軸は、t。
　

　試しに、x=1.5の時の値を推定（内挿）してみると、t=1.5－2=－0.5として、yA(t=－0.5)≒0.91640625、となる。
　ちなみに、（2）節の4次のラグランジュ多項式からは、x=1.5では、0.91640625、となり、実は、一致するのね。
　ま、こりゃ、当然か。
　違っていれば、それこそ問題だもんね」

「ニュートンの前進補間公式のところで書いたように、差分を利用する場合は、計算量が減るということじゃな。
　でき上がるものは、同じじゃ。
　いずれも、xのn次式を使って、すべての点を通る曲線という条件で探すと、計算の手段は、異なれど、みな同じものになる。
　なお、念のため、言うておくと、すべての点を通らなくてもよいのであれば、違った構成法もある。
　たとえば、目的の関数と近似関数との差の2乗を最小にする「最小二乗法」、差の絶対値を最小にする「チェビシェフの多項式」、
あるいは、対象とする関数の性質に合った直交多項式による近似式など、多数考えられている。
　また、「スプライン関数」を使ってなめらかな近似式を求める方法もある。
　「スプライン関数」については、「数値計算法基礎」（田中敏幸 著：コロナ社：2006年4月初版）などを参照して欲しい。
　しかし、今回の目的が微分方程式の数値解法に利用しようということなので、これらの詳細には、立ち入らないことにしよう。

　さて、「ガウスの前進補間公式」に対して「ガウスの後退補間公式」もある。
　分点を、ゼロ－マイナス1－プラス1－マイナス2－プラス2・・、という順番にとるものじゃ。
　先に挙げた関数の中で、yC(x)=y0＋(1/h)Δy(－1)(x－x(0))＋(1/2h^2)Δ^2y(－1)(x－x0)(x－x(－1))、がそれじゃ。
　これに対して、前進型の場合は、　yA(x)=y0＋(1/h)Δy0 (x－x0)＋(1/2h^2)Δ^2 y(－1)(x－x0)(x－x1)、だ。
　ガウスの後退補間公式の4次式は、
　yC(x)=y0＋Δy(－1)（t｜1）＋Δ^2 y(－1)（t+1｜2）＋Δ^3 y(－2)（t+1｜3）＋Δ^4 y(－3)（t+2｜4）、となる。
　nが偶数の場合、yA(x)とyC(x)の平均をとったものが、「スターリングの補間公式」と言われるものじゃそうだ。
　また、nが奇数の場合は、同様に平均をとったものは、「ベッセルの補間公式」と呼ばれる。
　一般項の表現や補間公式としての利用例は、「計算機のための数値計算法」を参照して欲しい」　　

最終更新日　2016/5/4