計算情報学Ｉ
名古屋大学 情報文化学部
自然情報学科３年
第８回

鈴木泰博
情報文化学部・大学院情報科学研究科
複雑系科学専攻

Science ofinteractions

Agenda

計算情報学 Ｉ 第８回
1.
2.
3.

前回の復習
エイトケンの加速法
ステファンセンの方法

4.
5.
6.
7.

線形補間
2次補間
ラグランジュ補間
スプライン補間

Science ofinteractions

逐次代入法
最も原始的な近似解法
f(x)をx = g( x )と変形した方程式を、x k+1 = g( xk) と逐次に代入することにより解く方法

x = cos x
x1 = 1 (rad)x2 = cos x1 = cos 1….

なにしろ簡単なので昔から広く使われてきた。でも、収束が遅い。

x k+1 = cos xk

演習：逐次代入法をプログラミングしこの問題をといてみよ。

Science ofinteractions

区間の縮小の程度を決めるだから解近傍でのg´(τ)が小さいほど収束が早くなる。

逐次代入法の原理|e i+1|=|g´(τi)|・|ei|≦K|ei|, (i =0,1,…)
ここがポイントより、|e i+1|≦K|ei|≦K2|e i-1|≦…≦K i+1|e0|
K<1より、評価のために定数Ｋで変数g´(τi)を上から押さえる

lim k i+1 = 0 ( i→∞),lim |e i+1|=0 ( i→∞）⇒lim xi = α(i→∞）。

y=g(x)y=x

y=x

y=g(x)
x0x2 α x3x1

区間I収束する場合 Ｋ＝| g´( x ) | < 1x2x1 x0 α
区間I発散する場合 Ｋ＝| g´( x ) | > 1

Science ofinteractions

逐次代入法
数値実験の結果
0.73950.50.450.4
0.350.30.250.20.15
0.1
0.05
00.73930.73910.73890.73870.7385値
収束

1 5 9 13 17 21 25 29 33 37 41 45 49

解の周辺に達するまでに約２５回かかった。
そこから、収束するまでも約２５回の反復計算が必要であった(トータル ５１回の反復計算で正解に到達）。
Science ofinteractions

収束の速さの評価
反復式 x k+1 = φ( xk) ①の解（真の値）をx真とする。
また、第k近似解xkの誤差をek = xk – x 真
つまり、xk = x真 + ekとする。
①より、xk+1 = φ (x真 + ek)。これをテイラー展開すると、

xk 1( x真 )( x真 )( x真 )e k2ek

2

φ(x真)=x真、かつ、xk+1 = x真+ ek+1より、
誤差の評価

ek 1
Ce k

( x真 )ek2k

O (e )|f(x 真)|が>1 … 発散
<1 … 収束
この値に収束の速度が依存「1次収束」、もしくは、「収束次数が１」である
Science ofinteractions

ニュートン法•きわめてシンプルだが、収束が早く実用的•他のより複雑な手法と比べてもさほど遜色がない•唯一の難点はf’’(x)の計算が必要なこと

初期値を適当に選ぶ
k = 0,1,2,...,の順に
f ( xk )xk 1 xkf ( xk )接線の勾配 f’(xk)

収束判定
f(xk)
xk+2
xk+1

xk

f ( xk )f ( xk )xk xk 1 より。
Science ofinteractions

ニュートン法
数値実験の結果
0.752
0.75
0.748
0.746
0.744
0.742
0.74
0.738
0.736
0.7340.7320.3
0.25
0.20.15値
収束

0.10.0501234

逐次法では５１回の反復計算が必要であったが、４回の反復で正解に到達した。
１回の反復で既に解の近傍に到達し、すぐに収束している。
Science ofinteractions

収束の速さの評価
反復式 x k+1 = φ( xk) ①の解（真の値）をx真とする。
また、第k近似解xkの誤差をek = xk – x 真
つまり、xk = x真 + ekとする。
①より、xk+1 = φ (x真 + ek)。これをテイラー展開すると、

ek 1( x真 )ek

ek 1( x真 )ek( x)

f ( x)f ( x)

x

ek 1

Cek2

「2次収束」( x真 )
2ek( x真 )
22ek...
2f(x真)=x真、かつ、xk+1 = x真+ ek+1より、...

{ f ( x)}2 f ( x) f ( x)( x) 1{ f ( x)}2より、φ’(x真)=0O (ek3 )Science ofinteractions

逐次法 vs ニュートン法
e-x – x 2 = 0 を逐次法とニュートン法で解く。

xk 1

e xk

xk 1xk

e xk
e xk

xk2
2 xk

Science ofinteractions

逐次法vs ニュートン法
数値実験の結果
ニュートン法(2次収束）収束の比較（log scale)
10.01

13

57

9 11 13 15 17 19 21 23 25

1.2

0.3

1

0.25

0.8

0.2

0.6

0.15

0.4

0.1

0.2

0.05

0値
収束0123450.0001
逐次法
ニュートン法1E-06

逐次法(1次収束）

1E-08
1E-100.45
0.4
0.35
0.3
0.25
0.20.15
0.10.0500.709

1E-120.707
0.705

収束の比較

0.703
0.7010.699
0.697
0.695
1357値
収束

9 11 13 15 17 19 21 23 25
Science ofinteractions収束の加速
Aitkenの加速法(デルタ2乗法)
ニュートン法はよい方法だが” 微分” が必要
微分を使わないで、かつ、収束が速い方法は？
反復{x n }が” 線形収束” のとき、十分に大きいnに対して

xn 1A( xn), xn 2 A( xn 2)としてもOK(αは解)。Aを消去すると、
2xn( xn 1 xn )xn 2 2 xn 1 xn よりynxn( xn 1 xn ) 2xn 2 2 xn 1 xn

Science ofinteractions

Aitkenの加速法
数値実験の結果
Aitkenの方法

x=cos(x) Aikenの方法
0.74

0.006
0.005

0.735

0.004

0.73

x = cos (x)

0.003
0.002

0.725

0.001

0.72

01

収束の比較

10.1 10.01
0.001
0.0001
1E-05
1E-06
1E-07

値
収束

3557

9 11 13 15 17 19 21

9 13 17 21 25 29 33 37 41 45 49

逐次法
ニュートン法Aitkenの方法

1E-08
1E-09
1E-10
Science ofinteractions

ステファンセンの方法
Aitkenの加速法は逐次法より収束が速いがニュートン法より収束が遅い
ステファンセンの方法はAitkenの加速法を援用して微分をせずに2次収束を得られる方法である
逐次法

y n = φ(x n )zn = φ(y)xn 1znステファンセンの方法( znzn 2yn ) 22 yn xn

Aitkenの加速法
Science ofinteractions

ステファンセンの方法
数値実験の結果
ステファンセンの方法

0.9
0.8
0.7
0.6
0.5

値
収束

0.4
0.30.2x = cos (x)

0.1
01

収束の比較

23

10.1 1 5 9 13 17 21 25 29 33 37 41 45 49
0.01
0.001
0.0001
1E-05
1E-06
1E-07
1E-08
1E-09
1E-104

逐次法
ニュートン法Aitkenの方法
ステファンセンの方法
Science ofinteractions

（数値）解析的な背景
これまで紹介してきた反復法には次のような解析的な背景がある。
定理１ 閉区間 Ｉ から Ｉ への関数 g(x)が x, yI ならば、

| g ( x) g ( y) | L | x y |, (0L 1)

要するに区間に 比例して収縮するということ！

を満たすならば、gは I 内にただ一つの解αを持ち、xn+1 = g(xn)はαに収束する。
（このようなｇは縮小写像 contractive mappingと呼ばれる）。
注)一般にこのような条件を満たすとき、ｇはリプシッツ(Lipschitz)条件を満たすといい、Ｌをリプシッツ定数をいう。
証明 解の存在と収束性の証明については、中間値の定理を用いた収束性の証明と同様であるので省略する。
解の唯一性の証明は、もしα以外に解β, β=g(β)があったと仮定すると、|| | g( )g( ) | L ||より、Ｌ＝１となり条件に矛盾することにより得られる。
縮小写像の原理を応用するポイントはリプシッツ定数Ｌ＜１である！
Science ofinteractions解析的な背景（つづき）
系１ 関数ｇが区間 Ｉ 上で微分可能で、|g’(x)|≦L (xはＩの区間にふくまれる)ならば、ｇはリプシッツ条件を満たす。
証明： 区間内の点x, yについて、平均値の定理を用いると、
g(x) – g(y) = g’(ξ）(x – y) … ①となるξが x と y のうちにとることができるので、ξも区間 Ｉ にふくまれる。
①の両辺の絶対値をとると、|g(x) – g(y)| = |g’(ξ)(x – y)| ≦ L| x – y|
よってｇは |g(x) – g(y)| ≦ L |x – y|なのでリプシッツ条件を満たす。

系２ α=g(α)とし、I = [α – d, α + d] とする。I において、ｇがリプシッツ条件をみたし、Lがリプシッツ定数であるときに、任意のxが区間 Ｉ に含まれるならば、g(x)も Ｉ に含まれる。
証明: |g(x) – α|=|g(x) – g(α)|≦L| x – α|≦Ld<d 。
リプシッツ条件を満たすから

ddα

Lはリプシッツ定数(<1)
だからxＩ

以上より、α=g(α)において、αの近くではgは縮小写像になることが、反復法が収束するための条件であることがわかる。

Science ofinteractions

補間•

補間 とは粗い数表しかなかった時代に数表に無い値を近似的に求めるために考えられたものである。非常に精緻に理論が作られている。
•
数値計算のみならず、統計やコンピュータグラフィックスにも広く応用されている技法である。
Science ofinteractions

線形補間
実験データ

y i+1

離散値
T

2345

V

34.5

67.5

補間 が必要ΔyiyyiΔxi? T=2.5の時のVを知りたい
xixx i+1

一番簡単な補間は線形補間
一般に誤差は分点と分点の間で大きく分点に近ければ小さくなる。

Science ofinteractions

線形補間
ｘがxiからx i+1まで
Δｘ= x i+1 - xiだけ変化するとき、従属変数yはΔｙ＝y i+1 – yiだけ変化するので、変化率は
y i+1
Δyi

yyix

yi 1
xi 1

yixi

これと同じ比率で、x,yが変化していくとすると、

Δxi

xi

yi
xi

x i+1

y yi
x xi

yi 1 yixi 1 xi

yixi

yi

yi 1 yi( x xi )xi 1 xi

より、

yyixi

yi 1xi 1

Science ofinteractions2次補間
同様に3点(x0, y0), (x1, y1), (x2, y2)が与えられたとき、これを2次曲線

p ( x)

a0a 1 x a 2 x 2 を用いて補間を行う。この場合 y=P(x)はa2≠0で放物線となる。
y0a0

a 1 x0a 2 x02y1a0

a 1 x1a 2 x12y2a0a 1 x2a 2 x22を用いて未知数a0, a1, a2を求めることができる。

これを解くとP ( x)( x x1 )( x x2 )y0
( x0 x1 )( x0 x2 )( x x0 )( x x2 )y1
( x1 x0 )( x1 x2 )( x x0 )( x x1 )y2
( x2 x0 )( x2 x1 )が得られる。
Science ofinteractions一般の場合
n点の補間

P ( x)

a0a1 x a 2 x2

L an xn

命題: n点を通るn次多項式は唯一である。
証明: P(x)以外に与えられたn+1点を通るn次以下の多項式をQ(x)とする。すると、P(x)とQ(x)の差は、D(x)=P(x)‑Q(x)となり、高々n次の多項式となる。
D(xk)=P(xk)‑Q(xk)=fk – fk = 0, k=0,1,2,…,nとなるので、D(x)はn+1の零点を持つ。しかし、代数学の基本定理より高々n次の多項式はnの零点しか持たない。
よって、恒等的にD(x)=0となるから、P(xk)‑Q(xk)=0より、P(x)=Q(x).
別証明：P(x)の係数行列はn次正方行列となる。この行列はVandermonde
行列式となるため、正則となり係数は一意に決まる。

Science ofinteractions

ラグランジュ補間
一般にＮ＋１の場合関数 lj(x),(j=1,2,3,…N)を以下のように定義する。

l j ( x)( x x0 )( x x1 ) L ( x x j 1 )( x x j 1 ) L ( x xN )(xjx0 )( x jx1 ) L ( x jx j 1 )( x jx j 1 )L( x jxN )ここでlj(x)は、

l j ( x)1 (i0 (ij)j)Nを満たす。そこで、PN(x)を、 PN ( x)y j l j ( x) と定義すると、高々Ｎ次の多項式となり、
j 0

PN ( x j )yjを満たす。この多項式をラグランジュの補間多項式とよぶ。

⇒例を考えるとわかります。
Science ofinteractions

ラグランジュ補間N

PN ( x)y j l j ( x)j 0

m=3の場合p3 ( x)( x x2 )( x x3 )f1
( x1 x2 )( x1 x3 )( x x1 )( x x3 )f2
( x2 x1 )( x2 x3 )( x x1 )( x x2 )f3
( x3 x1 )( x3 x2 )演習 P3(x)が(x 1,f 1), (x 2,f 2),(x 3,f 3)を通ることを確かめよ。

Science ofinteractions

演習
ラグランジュ補間演習: 以下のx kとf kに対する補間多項式を求めよ。

xk 0 1 2
fk 1 0 1
ヒント：m=3の場合のラグランジュ多項式
p3 ( x)( x x2 )( x x3 )f1
( x1 x2 )( x1 x3 )( x x1 )( x x3 )f2
( x2 x1 )( x2 x3 )( x x1 )( x x2 )f3
( x3 x1 )( x3 x2 )Science ofinteractions

解答
ラグランジュ補間
演習: 以下のx kとf kに対する補間多項式を求めよ。
xk012fk

10

1

p3 ( x)( x x2 )( x x3 )f1
( x1 x2 )( x1 x3 )( x x1 )( x x3 )f2
( x2 x1 )( x2 x3 )( x x1 )( x x2 )f3
( x3 x1 )( x3 x2 )より、

( x 1)( x 2)
( x 0)( x 2)
10(0 1)(0 2)(1 0)(1 2)11( x 1)( x 2)x( x 1) ( x 1) 2
22

p3 ( x)( x 0)( x 1)
1
( 2 0)(2 1)

Science ofinteractions

ラグランジュ補間の誤差
1

f ( x) P ( x)(n 1)!

( x) f n 1 ( )nwhere, ( x)k 0( x xk )証明の概略

H ( x)f ( x) P ( x)where, G ( x)( x)G ( x)f ( x) P ( x)
( x)

とおくと、
H ( x) 0
H ( xk ) f ( xk ) P ( xk )( xk )G ( xk )0

H(x)はx, x 0,…,x nのn+2の点で0となる。

Science ofinteractions

ラグランジュ補間の誤差(つづき）
ロルの定理より
H’(z)はn+1の異なる点で0となる
H’’(z)はnの異なる点で0となる
H’’’(z)はn-1の異なる点で0となる。
…
H (n+1) (ξ)=0なるξがx,x 0,x 1,…,x nの区間内に存在する。
H ( x)f ( x) P ( x)( x)G ( x)f ( x) P ( x)where , G ( x)( x)

H (n+1) (ξ) = f (n+1) (ξ) – π (n+1) (ξ)G(x)=0
また、π (n+1) (ξ)=(n+1)!より、G ( x)

1(n 1)!

f( n 1)H ( x)f ( x) P ( x)( )f ( x) P ( x)1
( n 1)!

1(n 1)!

f ( n 1) ( ) ( x)f ( n 1) ( ) ( x)Science ofinteractions0

スプライン補間
ラグランジュ補間の場合はN+1個の点を通る曲線をN次多項式で生成した。
しかし、点の数が多くなると隣り合う点との間で曲線が激しく振動する現象が現れる（自由度が増大することによる）。

区分的にm次のラグランジュ補間を行い、それらを滑らかにつなぐことにより自由度の増大に伴う振動を抑える。
そのため、以下の条件を課したもとでの補間を行う。
1. 曲線y=s(x)は連続であり、点(xj, yj)(j=0,1,2,…N)をすべて通る。
2. 区間の境目、すなわち、x=xj（j=1,2,..N-1)でy=s(x)の一階微分係数および二階微分係数が連続である。

Science ofinteractionsスプライン補間
3次スプライン補間の場合
実用上一番よく使われる３次スプライン補間を考える。その場合に、補間に用いる3次方程式は、S j ( x)a j (x x j )3b j (x xj )2

c j (x xj ) d j ( j0,1,K , N 1)

スプライン補間の条件を具体的に表すと、まず(1)の条件から、
Sj(xj)=yj(j=0,1,…,N-1)
Sj(xj+1)=yj+1(j=0,1,…,N-1)
が導かれる。次に、(2)の条件から

S j (xj 1)

S j 1 (x j 1 ) ( j0,1,2, K, N 2)

S j (xj 1)

S j*1 ( x j 1 ) ( j0,1,2, K, N 2)が得られる。
Science ofinteractions