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

第14回:補間関数2

0. 前回の補足

前回の授業では、与えられた関数の $n+1$ 個の節点における $n$ 次以下の補間多項式が一意に定まることを説明した上で、補間多項式を求めるために「連立一次方程式を解く方法」と「ラグランジュ補間を用いる方法」を扱いました。授業では扱っていませんが、他にも「ニュートン補間を用いる方法」が知られています。

前回は説明しませんでしたが、上記の補間多項式を、求める方法に依らずに(ラグランジュ基底多項式の線形結合として表現しない場合にも)ラグランジュ補間多項式(Lagrange interpolating polynomial)と呼ぶことが多いです。

特に、今回の授業では別の条件の下で別の補間多項式を考えることになるため、区別する意味でも前回のものをラグランジュ補間多項式と呼ぶことにします。

1. エルミート補間

1.1. 定義

ラグランジュ補間多項式は「元の関数と節点における関数値が一致する」という条件の下で定まりますが、条件を追加して「元の関数と節点における関数値および微分値が一致する」という状況を考えるのがエルミート補間(Hermite interpolation)です。

微分可能な関数 $f$ と相異なる $n+1$ 個の節点 $x_0,\ldots,x_n$ が与えられたとき、

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

および

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

を満足する $2n+1$ 次以下の多項式 $p$ が一意に存在します。これをエルミート補間多項式(Hermite interpolating polynomial)と言います。

ここで、$'$ は関数の微分を表しており、$2n+1$ 次以下の多項式(係数が $2n+2$ 個以下)になるのは条件が $2n+2$ 個あるからです。

なお、エルミート補間において、より一般には高次までの微分の値が一致するという条件を課す場合もあります。

1.2. 連立一次方程式を解いて求める方法

エルミート補間多項式 $p$ を求めるために

$$ p(x)=\sum_{j=0}^{2n+1} c_j x^j $$

とすれば、

$$ p'(x)=\sum_{j=1}^{2n+1} c_j \cdot j x^{j-1} $$

と条件より

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

および

$$ c_1+c_2\cdot 2x_i+\cdots+c_{2n}\cdot 2n x_i^{2n-1}+c_{2n+1}\cdot (2n+1)x_i^{2n}=f'(x_i)\quad (\forall i=0,\ldots,n) $$

が成り立つので、

$$ A= \begin{pmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^{2n} & x_0^{2n+1}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ 1 & x_n & x_n^2 & \cdots & x_n^{2n} & x_n^{2n+1}\\ 0 & 1 & 2x_0 & \cdots & 2nx_0^{2n-1} & (2n+1)x_0^{2n}\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ 0 & 1 & 2x_n & \cdots & 2nx_n^{2n-1} & (2n+1)x_n^{2n} \end{pmatrix} ,\quad c= \begin{pmatrix} c_0\\ c_1\\ c_2\\ \vdots\\ c_{2n}\\ c_{2n+1} \end{pmatrix} ,\quad b= \begin{pmatrix} f(x_0)\\ \vdots\\ f(x_n)\\ f'(x_0)\\ \vdots\\ f'(x_n) \end{pmatrix} $$

として連立一次方程式

$$ Ac=b $$

が得られます。

この連立一次方程式を解けば係数ベクトル $c$ を求めることができ、それによってエルミート補間多項式 $p$ が求まります。

上の式は一般の $n$ について書いているため分かりにくいですが、例えば $n=1$ の場合を考えると

$$ p(x)=c_0+c_1 x+c_2 x^2+c_3 x^3, $$$$ A= \begin{pmatrix} 1 & x_0 & x_0^2 & x_0^3\\ 1 & x_1 & x_1^2 & x_1^3\\ 0 & 1 & 2x_0 & 3x_0^2\\ 0 & 1 & 2x_1 & 3x_1^2 \end{pmatrix} ,\quad c= \begin{pmatrix} c_0\\ c_1\\ c_2\\ c_3 \end{pmatrix} ,\quad b= \begin{pmatrix} f(x_0)\\ f(x_1)\\ f'(x_0)\\ f'(x_1) \end{pmatrix}, $$$$ Ac=b $$

となります。

演習1

連立一次方程式を解くことによって、関数

$$ f(x)=e^x $$

の節点

$$ x_0=0,\quad x_1=1 $$

におけるエルミート補間多項式を求め、$-2\leq x\leq 3$ の範囲で元の関数 $f(x)$ とエルミート補間多項式のグラフを描画してください。

なお、連立一次方程式を解く際には演算子\を使用し、

$$ f'(x)=e^x $$

は既知のこととしてよいです。

1.3. 基底多項式を用いて求める方法

エルミート補間多項式 $p$ を求めるために

$$ p(x)=\sum_{j=0}^n f(x_j)H_j(x)+\sum_{j=0}^n f'(x_j)K_j(x) $$

として、$2n+1$ 次の多項式 $H_0(x),\ldots,H_n(x),K_0(x),\ldots,K_n(x)$ の節点での関数値および微分値について

$$ H_j(x_i)=\delta_{ij}= \begin{cases} 1\quad (i=j)\\ 0\quad (i\neq j) \end{cases},\quad K_j(x_i)=0, $$$$ H_j'(x_i)=0,\quad K_j'(x_i)=\delta_{ij}= \begin{cases} 1\quad (i=j)\\ 0\quad (i\neq j) \end{cases} $$

が成り立つと仮定します。

このとき、$p$ はエルミート補間多項式の条件

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

および

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

を満足することが容易に確認できます。$H_0(x),\ldots,H_n(x),K_0(x),\ldots,K_n(x)$ をエルミート基底多項式(Hermite basis polynomial)と呼びます。

あとは、エルミート基底多項式の具体的な式が分かればエルミート補間多項式を求めることができます。証明はしませんが、エルミート基底多項式はラグランジュ基底多項式

$$ 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) $$

を用いて

$$ H_j(x)=\left(1-2L_j'(x_j)(x-x_j)\right)\left(L_j(x)\right)^2\quad (j=0,\ldots,n), $$$$ K_j(x)=(x-x_j)\left(L_j(x)\right)^2\quad (j=0,\ldots,n) $$

と表されることが知られています。

上の式は一般の $n$ について書いているため分かりにくいですが、例えば $n=1$ の場合を考えると

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

より

$$ p(x)=f(x_0)H_0(x)+f(x_1)H_1(x)+f'(x_0)K_0(x)+f'(x_1)K_1(x), $$$$ H_0(x)=\left(1-2\frac{x-x_0}{x_0-x_1}\right)\left(\frac{x-x_1}{x_0-x_1}\right)^2,\quad H_1(x)=\left(1-2\frac{x-x_1}{x_1-x_0}\right)\left(\frac{x-x_0}{x_1-x_0}\right)^2, $$$$ K_0(x)=(x-x_0)\left(\frac{x-x_1}{x_0-x_1}\right)^2,\quad K_1(x)=(x-x_1)\left(\frac{x-x_0}{x_1-x_0}\right)^2 $$

となります。

演習2

エルミート基底多項式を用いることによって、関数

$$ f(x)=e^x $$

の節点

$$ x_0=0,\quad x_1=1 $$

におけるエルミート補間多項式を求め、$-2\leq x\leq 3$ の範囲で元の関数 $f(x)$ とエルミート補間多項式のグラフを描画してください。

なお、

$$ f'(x)=e^x $$

は既知のこととしてよいです。

2. 一般の関数による補間

ここまでは多項式による補間について考えてきましたが、一般の関数による補間を行うことも可能です。

2.1. 方法

関数 $f$ の $n+1$ 個の節点 $x_0,\ldots,x_n$ における補間関数となる、すなわち

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

を満足する関数 $g$ を、$n+1$ 個の関数 $g_0,\ldots,g_n$ の線形結合として与えることを考えます。

もし $n+1$ 個の $n+1$ 次元ベクトル

$$ \begin{pmatrix} g_0(x_0)\\ g_0(x_1)\\ \vdots\\ g_0(x_n) \end{pmatrix},\quad \begin{pmatrix} g_1(x_0)\\ g_1(x_1)\\ \vdots\\ g_1(x_n) \end{pmatrix},\quad \ldots,\quad \begin{pmatrix} g_n(x_0)\\ g_n(x_1)\\ \vdots\\ g_n(x_n) \end{pmatrix} $$

が線形独立ならば、そのような補間関数は一意に定まります。このとき、$g_0,\ldots,g_n$ を基底関数(basis function)と呼びます。

実際、

$$ g(x)=\sum_{j=0}^n c_j g_j(x) $$

と表せば、補間関数の条件より

$$ c_0 g_0(x_i)+c_1 g_1(x_i)+\cdots+c_n g_n(x_i)=f(x_i)\quad (\forall i=0,\ldots,n) $$

が成り立つので、

$$ A= \begin{pmatrix} g_0(x_0) & g_1(x_0) & \cdots & g_n(x_0)\\ g_0(x_1) & g_1(x_1) & \cdots & g_n(x_1)\\ \vdots & \vdots & \vdots & \vdots\\ g_0(x_n) & g_1(x_n) & \cdots & g_n(x_n) \end{pmatrix} ,\quad c= \begin{pmatrix} c_0\\ c_1\\ \vdots\\ c_n \end{pmatrix} ,\quad b= \begin{pmatrix} f(x_0)\\ f(x_1)\\ \vdots\\ f(x_n) \end{pmatrix} $$

として連立一次方程式

$$ Ac=b $$

が得られます。

仮定より $A$ の列ベクトルは線形独立であるため、この連立一次方程式を解くことによって係数ベクトル $c$ が一意に定まる、すなわち補間関数 $g$ が一意に定まることが分かります。

2.2. 三角関数による補間

関数 $f$ と区間 $[0,1]$ 内の $2n+1$ 個の節点 $x_0,\ldots,x_{2n}$ が与えられたとき、基底関数を

$$ g_0(x)=1, $$$$ g_j(x)=\sin(j\pi x)\quad (j=1,\ldots,n), $$$$ g_{n+j}(x)=\cos(j\pi x)\quad (j=1,\ldots,n) $$

のように取れば補間関数が一意に定まります。

これは三角関数による補間であり、フーリエ級数の有限個の項を考えていると解釈することもできます。

演習3(オプション)

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

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

と節点

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

に対して三角関数による補間を考えます。

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

第14回レポート課題

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