スプライン補間 関数近似 (2)

関数近似 (2)

スプライン補間

•

補間とは, 関数

が未知であるが, その

個の点

における値

が知

られているとき,

以外の点における

の値を推測しようとすることをいう. 1

次元の補間問題において, min{x

となっているときに補

間, そうでないときに補外と呼ぶ

•

もっとも基礎的な補間は線形補間.

•

これ以外に, 多項式による補間もよく用いら れる.

•

多項式による関数の近似や補間を考えること

に根拠を与えているのが, Weierstrass の多項

式近似定理である.

Weierstrassの多項式近似定理

関数

が区

間

において連続であれば,

に一様収束す

る多項式の列が存在する.

W. Rudin, Principles of Mathematical Analysis, 3/e, McGraw-Hill, 1976

「一様収束」とは,一様ノルムの意味で収束という意味.また,一様ノル ムの定義は,kfk∞= supt∈[a,b]|f(t)|.前回紹介したessential supremum と同じ記号は,意味がほぼ同じだから. このあたりの議論には深入りし

この定理は次のような形で一般化される.

合とする.

1 f, g∈ A,α∈C_ならf+g∈ A,f g∈ A,αf ∈ A 2 x₁, x₂ ∈K,x₁ 6=x₂なら∃f ∈ A, f(x₁)6=f(x₂) 3 ∀x∈K,∃g∈ A,g(x)6= 0.

このとき,Aの閉包はK上の連続関数全体を含む.

W. Rudin, Principles of Mathematical Analysis, 3/e, McGraw-Hill, 1976

• Fourier

解析の基礎を与えるのが, Stone - Weier-

の定理である.

•

三角関数による近似

級数展開) は, 単 位円上で定義された関数

の多項式 による展開と解釈することができる. これを 導いたのも

である.

竹之内,フーリエ展開,秀潤社, 1978.

•

無限次元空間における線形写像のことを線形 作用素と呼ぶことが多い.

• (V,k · kV),(W,k · kW)

をノルム空間,

を

から

への線形作用素とする. このとき,

の作用素ノルムは次のように定義される:

kAk= sup

x∈V,kxkV≤1

kAxkW (岩波数学辞典).

• V, W

におけるノルムの取り方には, 様々な ものがあり得る. これらの取り方に応じて作 用素ノルムは変わる.

•

作用素ノルムのことを誘導ノルム

norm)

と呼ぶこともある

S. Skogestad and I. Postlethwaite, Multivarible feedback control, Wiley, 1996.

• m

行

列の行列は,

次元のベクトルと 解釈することもできるし,

次元数空間から

次元数空間への線形写像

と解釈 することもできる.

•

したがって, 行列

を

行

列の行列とした

とき,

のベクトルとしてのノルムと作用素

としてのノルムが定義できる.

• A

のベクトルとしてのノルムの定義は,

k=1

Pn

l=1|a_kl|^p)^1/p (1≤ p <∞) kAk∞= maxk∈{1,...,m},l∈{1....,n}|akl|

• A

の作用素としてのノルムは

• (V,k · k_p), (W,k · k_p)に対し,V からW への線形 作用素の作用素ノルムをlp誘導ノルムと呼ぶ

S. Skogestad and I. Postlethwaite, Multivarible feedback control, Wiley, 1996.

• (V,k · k_p), (W,k · k_q)としたとき, V からW への 線形作用素の作用素ノルムは(p, q)−誘導ノルム とでも呼ぶべきであろうが,このような用語は定 着していない.

• Scilab

の組み込み関数

は, 行列

に対し て以下のように振舞う.

norm(A) l₂

誘導ノルムを返す

norm(A,1) l₁

誘導ノルムを返す

誘導ノルムを返す

ノルムを返す

(ベクトルとしてのl₂ノルム)

• ややこしいことに, Scilab(およびMATLAB)の オンラインマニュアルは,l_p誘導ノルムのことを 単にl_pノルムあるいはpノルムと呼んでいる. な お,norm(A,’fro’)はAをベクトルと見たときの l₂ノルムを返すが, これはFrobeniusノルムと呼 ばれている.

• 行列のノルムについてはD. S. Bernstein, Matrix Mathematics, 2/e, Princeton, 2009を参照.

• (x_i, y_i)i∈{1,=...,n}

を

個の実数の対とし, これ らのすべての点を通る

次の多項式を求 める, という問題を考える. ただし, これら の実数の対には重複がなく, かつ

から

までの中には値が同じものはないと仮定する.

(xi, y_i)i∈{1,=...,n}

がある関数に対応すると考え

れば, このような仮定は妥当である.

• n−1

次の多項式には係数が

次の項から 第

次の項まで)

個あるから, 多項式の 係数を未知数と解釈すると,

個の未知数に 関する

個の方程式を解くことになり, 解は

(大抵の場合には)

存在する筈である.

•

ここで, 以下の多項式を考える.

p_i(x) = Q

j<i(x−x_j)Q

j>i(x−x_j) Q

j<i(xi−x_j)Q

j>i(xi−x_j)

• i = 1のときにはj < iとなるjは存在しないが, この場合に はQ

j<i(x−xj) = 1,Q

j<i(xi−xj) = 1と定義する. 同様に,

i=nのときにはj > iとなるjは存在しないが, この場合には Q

j>i(x−xj) = 1,Q

j>i(xi−xj) = 1と定義する. これらは,記法

を簡潔にするための約束ごとである.

• p_i(x)

の分子を見ると,

であれば,

となることがわかる.

•

一方,

とすると,

の分母と分子が等し くなるから,

であることがわかる.

•

したがって,

は,

で,

なら

となる,

次の多項式である.

• p_i(x)

の性質を使い,

l(x) =y₁p₁(x) +y₂p₂(x) +· · ·+y_np_n(x)

とすると,

は

を通 る

次の多項式になっていることがわか

る. これが

補間多項式である.

• Lagrange

補間多項式による補間のことを

補間という.

• Lagrange

補間多項式の次数は与えられた標

本点の数から決まる.

個の標本点が与えら

れたときには, 次数は

次になる. そし

て, Lagrange 補間多項式 は与えられた標本

点に対して一意的である.

•

以下の議論では, 区間

で定義された連続

関数

を

補間によって近似すると

いう問題を考える.

• L_k

を, 区間

において

個の標本点を定 め, 関数

から

個の標本点における関数値

を用いて

補間多項式を定める作用

素とする.

• L_k

は線形作用素であることが示せる.

• f

から

によって作られる

補間多 項式を

と書く.

•

問題となるのは, 作用素の系列

が与

えられたとき,

が

を大きく取れば

に

近付くかどうか, ということである.

• 関数の定義域を区間[a, b]とし, (L_k)_k∈Nを,区間

[a, b]におけるあるk個の標本点の取り方に対応

したLagrange補間の系列とすると, (L_k)_k∈Nに よって一様近似できない連続関数が存在する.

• 逆に, 連続関数f を固定した場合には, 標本点の 取り方をうまく定めることで,それを一様近似す るような列(L_k)_k∈Nを作ることができる.

• fをn回以上連続微分可能な関数とする. Lnをn 個の標本点(x₁, . . . , x_n)に対応したLagrange補 間多項式を作る作用素とし, W =Q_n

k=1(x−x_k) とする. このとき, 以下が成り立つ[杉原,室田]:

kf −L_nfk_∞≤ kf⁽ⁿ⁾k_∞kWk∞

n!

• f

が解析関数であっても, 標本点の取り方を 等間隔としてしまうと,

をいくら大きく取っ ても

が

に近付かないことがある. これ

を

の現象という. 標本点の取り方を工

夫すれば, この問題を改善することができる

室田].

• Lagrange

補間には, 数値計算の誤差に弱いと いう問題がある.

•

展開に

Q

j<i(x−xj)Q

j>i(x−xj) Q

j<i(xi−xj)Q

j>i(xi−xj)

のかわ

りに

を使うこ

とにより, この問題は緩和される.

• f_n(x) =c₀+c₁(x−x₁) +c₂(x−x₁)(x−x₂) +

· · ·+c_n−1(x−x₁)· · ·(x−x_n−1)

とおく.

• f(x₁) = c₀,f(x₂) =c₀+c₁(x₂−x₁),f(x₃) = c₀+c₁(x2−x₁) +c₃(x3−x₁)(x3−x+2), . . . , f(xn) =c₀+Pn−1

i=1 c_iQi

j=1(xn−x_j)

である.

•

したがって,

を満た す

が一意的に定まる. この係数 を用いた補間を

補間という.

• Newton

補間と

補間の計算のしかた

は異なるが, 標本点を固定したとき, Newton

補間と

補間によって得られる多項

式は同一である.

• Lagrange

補間

補間) で得ら

れた多項式は, 標本点の取り方によっては, 標

本点間における関数値の脈動が激しい多項式

となる. しかし,

個の標本点を通る

次

の多項式は一意的に定まるため, 多項式の次

数をこのままにした場合には, これを改善す

る余地はない.

• Lagrange

多項式より高い次数の多項式を用 いることで, この「脈動」を抑えることを考 える.

•

具体的には, 標本点において, 関数

と関数

値およびその

階の導関数の値を一致させる

ことを考える.

•

区間

において少なくとも

階微分可能

な関数

を多項式によって近似したい.

個の相異なる点

において,

の値

と, その導関数

の

値

が与えられているものと

する.

• 2n

個の方程式を解くということは, 未知数が

個必要である. すなわち, 2

次の多項

式

個) を使えば,

目的は達成される筈である.

• 以下,斎藤,数値解析入門,東京大学出版会, 2012に基づいてHermite多項式の求め 方を説明する. [杉原,室田]にはVandermonde行列を用いた導出法が記載さている.

• Lagrange

補間で用いた

Q

j<i(x−xj)Q

j>i(x−xj) Q

j<i(xi−xj)Q

j>i(xi−xj) (i = 1, . . . , n)

を再び用いる.

• h_i(x) = (p_i(x))²(1−2p^′_i(x_i)(x−x_i)),

k_i(x) = (p_i(x))²(x−x_i)

とする

•

見にくいが, (1

において,

の係数

は

の点

における 値であり, 定数である. よって,

は

次の多項式である.

•

定義から,

も

次の多項式である.

•

代入して計算すると次の事実が確認できる.

h_i(x_j) =

(1, j=i,

0, j6=i,, h^′_i(x_j) = 0,

k_i(x_j) = 0, k_i^′(x_j) =

(1, j=i, 0, j6=i,

•

したがって, 以下のようにすると

において

およびその導関数と値が一致する

次多項式が得られる. これが

補間多項式である.

q_2n−1(x) =

n

X

i=1

f(xi)hi(x) +

n

X

i=1

f^′(xi)ki(x)

• Hermitte

補間多項式を求めることを

補間という

室田]).

• Hermite

補間多項式の近似誤差は, Lagrange

補間多項式と同様の方法により評価できるが,

詳細は略す

室田]).

• e^−x2cos 2x

および

を区間

に おいて

補間および

補間し た例を次ページに示す.

•

標本点は

点で, 等間隔とした.

• e^−x2cos 2x

については, 区間

の外部で

関数と補間多項式を比較した.

-0.2 0 0.2 0.4 0.6 0.8

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

y

x

Hermite f(x_k)

-600 -400 -200 0 200 400 600 800 1000

-3 -2 -1 0 1 2 3

y

x

Hermite

-4 -2 0 2 4 6 8 10

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

y

x

Hermite f(x_k)

• 標本点では関数とその補間多項式の値は一致する.

• Lagrange補間の, 区間の境界の近くで関数値か

らの乖離が発生するという欠点はHermite補間 では緩和されている.

• Lagrange補間, Hermite補間とも,補外に関する 性能は悪い.

• 関数の脈動と比較して標本点が少ない場合には, 関数とその補間多項式の値との標本点以外の点に おける値の乖離が大きくなる.

• とはいえ,数値計算の誤差の問題があるため,標本 点を無闇に多くするわけにはいかない(特にHer- mite補間多項式).

• 標本点以外の点における近似精度については,線 形補間の方が良いこともある.

• Lagrange補間とHermite補間には, 多項式の次 数が高いという問題と,標本点以外で関数値と補 間多項式の値に大きな乖離が発生することがある という問題があった.

• 線形補間にはこのような問題はないが,一方で微 分可能でないという問題がある.

• 線形補間の良さを踏襲しつつ補間関数を微分可能 にしようというのがスプライン補間の考え方.

•

線形補間では, 区間

における補間直 線と区間

は別物であった.

•

スプライン補間では, 区間

と区間

とで, 相異なる多項式を用いて補間 をおこなう.

•

以下の議論の典拠は

室田], [斎藤].

• n

個の標本点

があれ

ば

補間と異なり, 大小順に並べら

れている必要がある), 各区間で

個, 全体で

個の補間多項式が作られる.

•

有限個の点

が与えられ, 区間

において

が

次の多項式となるとき,

は区分的

次多 項式という.

• (m−1)

階までの導関数が全領域で連続とな

る区分的

次多項式を,

次スプラインと

いう.

• n

個の標本

が与えられ ているとき,

における値が

となるよう な

次スプラインによる補間を

次スプラ イン補間という.

•

応用上は, 3 次のスプラインがよく用いられ

る.

• S(x)

を, これから構成しようとしている

次 の区分的

次多項式とする.

•

スプライン補間は補間であるから,

でなければならない.

•

開区間

では

は多項式だから,

は連続である.

•

よって, 区間の端点における

の右微分と

左部分が一致すれば,

によってスプライ

ン補間が達成される. このようにするために

は,

とし,

を未知数とする方程

式を解けばよい

•

区間

においてスプライン補間の

「部品」を 構成するには, 2 個の 標本点

に関する

補間

多項式 を用いる.

のとき

補

間多項式の次数は

だから, 次数は

合っている.

• どの区間について計算しても同じことなので, {(x1, y1),(x2, y2)}が与えられているものとし,区 間[x₁, x₂]におけるHermite補間を考える.

• 2点に関するLagrange補間はn −2 = 1だか ら線形補間であり, p₁(x) = (x−x₂)/(x₁−x₂), p1(x) = (x −x1)/(x2 −x1)とすると, l(x) = y₁p₁(x) +y₂p₂(x)である.

• δ = x₂ −x₁ とおきh₁, h₂, k₁, k₂ を求めると, h₁(x) = _δ¹3(x −x₂)²(δ + 2(x −x₁)), h₂(x) =

1

δ³(x−x₁)²(δ−2(x−x₂)),k₁(x) = _δ¹2(x−x₁)(x−

x2)²,k2(x) = _δ¹2(x−x1)²(x−x2)となる.

• y₁, y₂, d₁,d₂を用いてHermite補間多項式を求 めると,y1h1(x) +y2h2(x) +d1k1(x) +d2k2(x)と なる. そのx₁における右微分は d₁,x₂における 左微分はd₂である.

• 区間[x₂, x₃]以降でも同様に操作をおこなう. 区 間[x₂, x₃]の左端x₂ における補間多項式の右微 分はd₂で, [x₁, x₂]の右端x₂における補間多項式 の左微分と一致するから, 1階の導関数は連続で ある. d₁, d₂, . . . , d_nはまだ決まっていない.

• m= 3だから, 2階の導関数も全領域で連続でな ければならない. ここからd_kを定める.

• この補間多項式のx₁における右2階微分とx₂に おける左2階微分はそれぞれ6^(y2_δ^−y₂ ¹⁾−4^d_δ¹−2^d_δ²,

−6^y2_δ^−y2 ¹ + 2^d_δ¹ + 4^d_δ².

• δ₁ =x₂−x₁,δ₂ =x₃−x₂とおき, [x₂, x₃]にお いて同様の計算をした後に,x₂における左右から の2階導関数の値を等号で結んで整理すると,

d₁ δ₁+2

1 δ₁ + 1

δ₂

d₂+d₃ δ₂ = 3

y₂−y₁

δ₁² +y₃−y₂ δ₂²

• δ_k=x_k+1−x_k (1≤k≤n−1)とおいて全区間 について同様の計算をおこない, γ_k = 1/δ_k−1 + 1/δ_k+1 (2≤k≤n−1),µ_k = 1/δ_k (1≤k≤n), ν_k= (y_k−y_k−1)/δ_k−1² + (y_k+1−y_k)/δ_k² (2≤k≤ n−1) とすると,これらの式を次ページように連 立1次方程式として纏めることができる.







µ₁ γ₂ µ₂ . .. ... . ..

µ_n−1 γ_n−1 µ_n−1











 d₁

... d_n





=





 ν₂

... ν_n−1







• この連立1次方程式は,n個の変数に対して方程 式がn−2個なので, 解は一意ではない.

• 解を一意にするために, 上記に式を2個追加して 解を求めることが普通である.

• 補間したい関数fが既知で微分可能のとき,f^′(x₁) = S^′(x1+0),f^′(xn) =S^′(xn−0)とする方法がある.

• これ以外に,S^′′(x₁+ 0) =S^′′(x_n−0) = 0とする 方法もある(教科書にはこちらだけが書かれてい る).

• 上記において, +0は右極限, −0は左極限をあら わす.

• Scilabでスプライン補間をおこなう関数は色々.

• 関数interpは,x= (x₁, . . . , x_n)における関数値 y= (y₁, . . . , y_n)とその導関数値dy = (y₁^′, . . . , y^′_n) および補間関数値を計算したい点のデータz = (z₁, . . . , z_m)を受け取り,zにおけるスプライン補 間関数値を返す.

• 次ページに例を示す.

スプライン補間をおこない, 0から0から3まで0.1刻 みで関数値を評価する:

x=[0 %pi/4 %pi/2 3*%pi/4 %pi];

y=sin(x);

dy=cos(x);

z=0:0.1:3;

val=interp(z,x,y,dy);

補間点における関数値を返すが, xとyのペアを2行n 列の行列にして与える点と(導関数は不要), 「0.1刻み で関数値を計算」などのように刻みを与える点, 返却値 は2行m列の行列で, 第1列にx軸上で取られた点の 値,第2列に補間関数値が入る点が異なる.

たとえば, 先ほどと同じx,yについて, x軸上で0.9刻 みで補間関数値を計算したいときには次のようにする.

x=[0 %pi/4 %pi/2 3*%pi/4 %pi];

y=sin(x);

val=smooth([x;y],0.9);

valの内容は以下のようになる(端点も含まれる):

val =

0. 0.9 1.8 2.7 3.1415927

0. 0.7817805 0.9724687 0.4333687 1.225D-16

• 関数splinは,与えられた標本点からスプライン 関数を生成し,その標本点におけるスプライン関 数の導関数を返す.

• ScilabのオンラインマニュアルのInterpolation の項を見ると,これ以外にもスプライン補間に対 応する関数があるのがわかるが,この講義ではこ れ以上述べない.