2022年度プログラミング演習A・B

第13回:補間関数1

1. 補間とは

複数の点(データ)が与えられたときに、それら全てを通るような関数を求めることを「点の間を補う」という意味で補間(interpolation)と言い、そのような関数を補間関数(interpolating function)と言います。

記号で書けば、$xy$ 平面上の $x$ 座標が相異なる有限個の点 $(x_0,y_0),\ldots,(x_n,y_n)$ に対して、

$$ g(x_i)=y_i\quad (\forall i=0,\ldots,n) $$

を満足する関数 $g$ が補間関数です。

特に、ある関数 $f$ 上の点 $(x_0,f(x_0)),\ldots,(x_n,f(x_n))$ を考える場合には、

$$ g(x_i)=f(x_i)\quad (\forall i=0,\ldots,n) $$

を満足する関数 $g$ を $f$ の節点 $x_0,\ldots,x_n$ における補間関数と言います。

関数の補間は、「節点における値が一致する」という制約の下での関数の近似(approximation)であると解釈することもできます。

2. 多項式による補間

上では補間に用いる関数 $g$ は一般の関数という想定で説明をしましたが、多項式による補間を考えることが最も基本的であり、よく使用されています。

関数 $f$ 上の相異なる $n+1$ 個の点 $(x_0,f(x_0)),\ldots,(x_n,f(x_n))$ が与えられたとき、これら全てを通る、すなわち $f$ の節点 $x_0,\ldots,x_n$ における補間関数となるような $n$ 次以下の多項式 $p$ が一意に存在します。このような多項式を補間多項式(interpolating polynomial)と言います。

$2$ 個の点によって $1$ 次関数(直線)が定まり、$3$ 個の点によって $2$ 次関数(放物線)が定まることを思い出すと理解しやすいです。ただし、$3$ 個の点が直線上に並ぶ場合のように、必ずしも $n+1$ 個の点からちょうど $n$ 次の多項式($n$ 次関数)が定まるわけではないことに注意してください。

補間多項式 $p$ を求めるために

$$ p(x)=\sum_{j=0}^n c_j x^j=c_0+c_1 x+c_2 x^2+\cdots+c_n x^n $$

とすれば、$n+1$ 個の点 $(x_0,f(x_0)),\ldots,(x_n,f(x_n))$ を通るという条件より

$$ c_0+c_1 x_i+c_2 x_i^2\cdots+c_n x_i^n=f(x_i)\quad (\forall i=0,\ldots,n) $$

が成り立つので、

$$ A= \begin{pmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^n\\ 1 & x_1 & x_1^2 & \cdots & x_1^n\\ \vdots & \vdots & \vdots & \vdots & \vdots\\ 1 & x_n & x_n^2 & \cdots & x_n^n \end{pmatrix} ,\quad c= \begin{pmatrix} c_0\\ c_1\\ c_2\\ \vdots\\ c_n \end{pmatrix} ,\quad b= \begin{pmatrix} f(x_0)\\ f(x_1)\\ f(x_2)\\ \vdots\\ f(x_n) \end{pmatrix} $$

として連立一次方程式

$$ Ac=b $$

が得られます。

このとき、行列 $A$ をヴァンデルモンド行列(Vandermonde matrix)と呼び、$x_0,\ldots,x_n$ が相異なることから $\det A\neq 0$ です。

したがって、連立一次方程式 $Ac=b$ は唯一の解を持ち、実際に解くことによって補間多項式 $p$ の係数を求めることができると分かります。

演習1(オプション)

$n$ を自然数とし、関数

$$ f(x)=\sin\left(2\pi x^2\right) $$

の節点

$$ x_i=\frac{i}{n}\quad (i=0,\ldots,n) $$

における補間多項式を考えます。

$n=3,5,7$ の各場合に連立一次方程式を解くことによって補間多項式の係数を求め、$0\leq x\leq 1$ の範囲で元の関数 $f(x)$ と $3$ 個の補間多項式のグラフを描画してください。

連立一次方程式を解く際には、演算子\を使用してよいです。

3. ラグランジュ補間

上で説明した方法では、補間多項式を求める際に連立一次方程式を解く必要があります。

別の方法として、ラグランジュ補間(Lagrange interpolation)を用いれば、より直接的に補間多項式を求めることができます。

3.1. 方法

相異なる $n+1$ 個の節点 $x_0,\ldots,x_n$ に対して、$n+1$ 個の $n$ 次多項式

$$ L_j(x)=\prod_{k=0\\ k\neq j}^n\frac{x-x_k}{x_j-x_k}=\frac{(x-x_0)\cdots(x-x_{j-1})(x-x_{j+1})\cdots(x-x_n)}{(x_j-x_0)\cdots(x_j-x_{j-1})(x_j-x_{j+1})\cdots(x_j-x_n)}\quad (j=0,\ldots,n) $$

ラグランジュ基底多項式(Lagrange basis polynomial)と言います。

ラグランジュ基底多項式の節点での値について

$$ L_j(x_i)=\delta_{ij}= \begin{cases} 1\quad (i=j)\\ 0\quad (i\neq j) \end{cases} $$

が成り立つので($\delta_{ij}$ はクロネッカーのデルタ)、

$$ p(x)=\sum_{j=0}^n f(x_j)L_j(x) $$

とおけば、$p$ は関数 $f$ の節点 $x_0,\ldots,x_n$ における補間多項式となることが分かります。

まとめると、ラグランジュ補間では、ラグランジュ基底多項式に対して関数値を係数とするような線形結合を取ることによって補間多項式を求めることができます。

3.2. 例

例として、$n=2$(節点が $3$ 個)の場合を考えます。

このとき、ラグランジュ基底多項式は

$$ L_0(x)=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)},\quad L_1(x)=\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)},\quad L_2(x)=\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} $$

の $3$ 個です。節点が

$$ x_0=0,\quad x_1=\frac{1}{2},\quad x_2=1 $$

の場合に $0\leq x\leq 1$ の範囲でラグランジュ基底多項式のグラフを描画してみると、次のようになります。

関数 $f$ の節点 $x_0,x_1,x_2$ における補間多項式は

$$ p(x)=f(x_0)L_0(x)+f(x_1)L_1(x)+f(x_2)L_2(x) $$

と与えられます。仮に上記の節点に対して

$$ f(x_0)=2,\quad f(x_1)=-1,\quad f(x_2)=5 $$

を満たすとすると、補間多項式のグラフは次のようになります。

演習2

節点の配列nodeに対してj番目のラグランジュ基底多項式のxでの値を返す関数Lagrange(node,j,x)を定義した上で、節点

$$ x_0=0,\quad x_1=\frac{1}{4},\quad x_2=\frac{1}{2},\quad x_3=\frac{3}{4},\quad x_4=1 $$

に対する全てのラグランジュ基底多項式の $0\leq x\leq 1$ の範囲のグラフを描画してください。

なお、関数の引数であるjの範囲は、0length(node)-11length(node)のどちらで考えてもよいです。

(どうしても難しければ、関数を定義せずにグラフの描画を行ってください。)

演習3

$n$ を自然数とし、関数

$$ f(x)=\sin\left(2\pi x^2\right) $$

の節点

$$ x_i=\frac{i}{n}\quad (i=0,\ldots,n) $$

における補間多項式を考えます。

$n=4,6$ の各場合にラグランジュ補間を用いて補間多項式を求め、$0\leq x\leq 1$ の範囲で元の関数 $f(x)$ と $2$ 個の補間多項式のグラフを描画してください。

第13回レポート課題

演習1~演習3に取り組んでください。