Transcript Miyazaki-PBV200310

Numerical Recipes in C
[日本語版]
初版
宮崎大輔
2003年10月30日（木）
PBVセミナー
第２章 連立１次方程式の解法
Solution of linear algebraic equations
連立方程式
a11x1  a12x2  a13x3   a1N xN  b1
a21x1  a22x2  a23x3   a2 N xN  b2
a31x1  a32x2  a33x3   a3N xN  b3

aM 1 x1  aM 2 x2  aM 3 x3   aMN xN  bM
A x  b
 a11 a12
a
a22
A   21
 


aM 1 aM 2
 a1N 
 b1 
b 
 a2 N 
b 2

  

 
 aMN 
bM 
• N=Mならば方程式の数と未知数の数が等しい→解が一意に
定まる可能性がある
• M<N(またはM=Nでも方程式が特異(singular))なら未知数の
数より方程式の数の方が少ない→解は存在しないか無数の
解が存在
• M>N、未知数の数より方程式の数の方が多い→解は存在し
ない→最小２乗問題ととらえ解を求める
Gauss-Jordan(ガウス・ジョルダン)法
(Gauss-Jordan elimination)
• 複数の右辺ベクトルbについて方程式の解を求め、同時に逆行列A-1も求
める
a11
a
 21
a31

a41
a12
a22
a32
a42
a13
a23
a33
a43
a14   x11   x12   x13   y11
      
a24   x21   x22   x23   y21




a34   x31   x32   x33   y31
       
a44   x41   x42   x43   y41
y12 y13 y14   b11   b12   b13   1
       
y22 y23 y24   b21   b22   b23   0




y32 y33 y34   b31   b32   b33   0
       
y42 y43 y44   b41   b42   b43   0
0
1
0
0
0
0
1
0
0 

0 
0 

1 
A x1  x2  x3  Y  b1  b2  b3 1
A  x1  b1
A  x2  b2
A  x3  b3
A Y  1
1. Aの2行と対応するb、1の行を入れ替えても、xとYの解は変わらない
2. Aの行を、その行と別の行の線形結合で置き換え、対応するb、1の行も線
形結合で置き換えれば、解は変わらない
3. Aの2列と対応するx、Yの行を入れ替えれば、解の行の順序が入れ替わ
る
• 上記3操作を組み合わせAを単位行列にすると右辺に解が得られる
• ピボット選択(pivoting): 0で割るのを防ぐ
Gaussの消去法(Gaussian elimination)と
後退代入(backsubstitution)
 a12
 a13

a11
 0 a a
22
23


 0 0 a33

0
0 0
   x1  b1 
a14
   x2  b2 
a24
    
   x3  b3 
a34
    
   x4  b4 
a44
簡単に解ける（後退代入）

x4  b4 a44
一般的には
1
xi 
aii
1


x3 
b3  x4a34

a33
N


bi   aij x j 
j i 1


LU分解
0
0
1

0
 21 1
31 32 1

 41  42  43
0 11 12 13
0  0 22 23

0  0
0 33
 
1  0
0
0
14  a11
24  a21

34  a31
 
44  a41
a12
a22
a32
a42
a13
a23
a33
a43
a14 
a24 
a34 

a44 
L U  A
Ax=bを解くのが簡単になる
A  x  (L  U)  x  L  (U  x)  b
L y  b
Ux  y
• LU分解はCrout(クラウト)のアルゴリズム(Crout’s
algorithm)を使う
• LU分解を使って逆行列を求めることもできる
• LU分解された行列の行列式は対角要素の積
連立方程式の解の反復改良(iterated
improvement)
A x  b
①
• このxを求めたい
• 未知の誤差δxを含む解x+δxが得られたとする
A  (x  x)  b  b
②
• δxを求めれば得られた解x+δxから引いて真の解x
が求まる
①を②に代入
A  x  b
③
②を③に代入
A  x  A  (x  x)  b
④
３重対角(tridiagonal)な連立１次方程式
 b1 c1 0
a b c
 2 2 2





  u1   r1 
 u  r 

  2   2 
      

 
  
 aN 1 bN 1 cN 1  uN 1  rN 1 
 0
aN bN   uN   rN 
• メモリO(N)、計算時間O(N)
– 通常は計算時間O(N3)
• アルゴリズムの詳細は省略
Vandermonde(バンデルモンド)行列
(Vandermonde matrices)
1 x1

1 x2
 

1 xN
x  x   c1   y1 

x  x   c2   y2 



   
  
2
N 1  
c
xN  xN   N   yN 
2
1
2
2
N 1
1
N 1
2
• N個の点(xj,yj)に多項式を
当てはめ、その多項式の係
数ciを求める方程式
1
 1
 x
x2
 1
 x12
x22


 x1N 1 x2N 1
 1   w1   q1 
 xN   w2   q2 
 xN2    w3    q3 
    

    
 xNN 1  wN  qN 
• N個の点xiと値qjが与えられ、
未知の重みwiを求め、N個
のモーメントをqjに一致させ
る問題
計算時間はO(N2)
通常はO(N3)
アルゴリズムの詳細は省略
Toeplitz(テープリッツ)行列(Toeplitz matrices)
R1 R2
 R0
 R
R0
R1
 1
 R2
R1
R0


RN 2 RN 3 RN 4

 RN 1 RN 2 RN 3
 R N 2 R N 1 

 R N 3 R N 2 
 R N 4 R N 3 



 R0
R1 

 R1
R0 
計算時間はO(N2)
通常はO(N3)
アルゴリズムの詳細は省略
特異値分解
(singular value decomposition, SVD)











A
 
 
 
 
 

 
 
 
 
 
U


 
  w1
 
 
w2
 
 



 
 

 
 
 
wN  


VT








M×N行列A
M×Nの列直交行列U
N×Nの対角行列W（対角成分は非負）
N×Nの直交行列V
U,Vは正規直交行列
Uの列、Wの要素、Vの列（つまりVTの行）を入れ替えただけ、は同じ分解とみなす
もしWのいくつかの要素がたまたま等しければ、
それに対応するU,Vの列をその線形結合で置き換えたものは同じ分解とみなす
SVDによる逆行列と連立方程式の解
• Aが正方行列(N×N)のとき
A1  V [diag(1/ wj )] UT
• また、 A  x  b の解は
x  V [diag(1/ wj )] (U  b)
T
ただし、wj=0のときは1/wj=0とする
零空間(nullspace)と値域(range)
• Aが特異のとき
[階数(rank)がN未満]
• 解は無数にある
• SVD解は0に最も近い解を
吐き出す
• SVD解に零空間を足したも
のが解空間になる
• 零空間とは
(a) Aが非特異
– Ax=0の解空間
– wj=0に対応するVの列の任意
の線形結合
(b) Aが特異
SVDの良さ
• Aが特異でなくてもwjが小さいときはwjを0にしてSVD
解を求めたほうが精度がいい
• 未知数より方程式が少ない場合は、足りない行に0
を入れてN×N行列に対してSVD解を求める
• 未知数より方程式が多いときは、SVD解は最小２乗
解となる
Sherman-Morrison(シャーマン・モリソン)公式
(Sherman-Morrison formula)
•
•
•
•
正方行列Aの逆行列A-1が求まったとする
Aに小さな変更をしたとき、A-1を簡単に計算したい
その変更が A  (A  u  v) ならOK
(u,vは何らかのベクトル)
例
1. uが基本単位ベクトルei→vをi行に加える
2. vが基本単位ベクトルej→uをj列に加える
3. u,vが基本単位ベクトルei,ejに比例→aijが加わる
•
計算方法：
(A1  u)  (v  A1 )
(A  u  v)  A 
1 
1
1
  v  A1  u
Woodbury(ウッドベリ)の公式
(Woodbury formula)
• 修正すべき箇所が多い場合はSherman-Morrison
の公式を何回も計算する必要がある
• その場合はWoodburyの公式を使う
(A  U  VT )1  A1  [A1  U  (1  VT  A1  U)1  VT  A1 ]







U

1
 


 
  1  VT  A1  U  
 
 
 

VT
• ここで、AはN×N行列
• U,VはN×P行列（P<N,通常P≪N）
• 逆行列の計算はP×P行列だけで済む




第３章 補間と補外
Interpolation and extrapolation
補間：低次の多項式か高次の多項式か
Lagrange(ラグランジュ)の公式
• N個の点 y1  f ( x1 ), y2  f ( x2 ), , yN  f ( xN )
を通る次数N-1の補間多項式は
( x  x2 )(x  x3 )( x  xN )
( x  x1 )(x  x3 )( x  xN )
y1 
y2
( x1  x2 )(x1  x3 )( x1  xN )
( x2  x1 )(x2  x3 )( x2  xN )
( x  x1 )(x  x2 )( x  xN 1 )
 
yN
( xN  x1 )(xN  x2 )( xN  xN 1 )
P( x) 
• この公式を効率的に計算するNevilleのアルゴリズム
というものがあるが、詳細は省略
有理関数による補間
• 多項式ではうまく近似できない関数でも有理関数でう
まく近似できることがある
• m+1個の点(x , y ),, (x , y ) を通る以下の関数を求める
i
i
i m
i m
P ( x) p0  p1x   p x
Ri (i1)(im) 

Q ( x) q0  q1x   q x
• Bulirsch-Stoer(ブリルシュ-ストア)アルゴリズムを使
えばよいが、詳細は省略
3次スプライン補間
y  Ayj  Byj 1  Cyj  Dyj1
A
x j 1  x
x j 1  x j
B  1 A 
x  xj
1
1
C  ( A3  A)(x j 1  x j )2 D  (B3  B)(x j 1  x j )2
x j 1  x j
6
6
• 入力 yi  y( xi ) を上の３次多項式で関数を近似する
• 以下の式と２つの条件からyi’’が求まる
x j  x j 1
x j 1  x j 1
x j 1  x j
y j 1  y j y j  y j 1
yj1 
yj 
yj1 

6
3
6
x j 1  x j x j  x j 1
– 端点（y1’’,yN’’）の片方または両方を0と置く
– 端点の片方または両方の１階導関数を与えてやる
• この連立方程式は３重対角行列になるのでO(N)で
解ける
双１次補間(bilinear interpolation)
• 2次元以上の補間(ここでは2次元の場合を説明する
が、3次元以上でも同様)
x1  X1L
t
X1U  X1L
u
x2  X2L
X2U  X2L
y(x1, x2 )  (1 t)(1 u) y(点1)  t(1 u) y(点2)  tuy(点3)  (1 t)uy(点4)
双３次補間(bicubic interpolation)
• 滑らかさを上げる高次補間
• 各格子点で関数 y( x1, x2 ) 以外に
y x1
y x2
勾配
2

y x1x2
クロス導関数
を人間が与え、それを元に補間する
双３次スプライン(bicubic spline)
• 滑らかさを上げる高次補間
• 双３次補間に似てるが、格子点での導関数の値を計
算してしまう
– まず行ごとにスプライン補間し
– その新しい列についてスプライン補間する
第４章 関数の積分
Integration of functions
閉じた公式と開いた公式
閉じたNewton-Cotes(ニュートン・コーツ)公式
(Newton-Cotes formulas)
• 台形則

x2
x1
• Simpson則

x3
x1
1 
1
f ( x)dx  h f1  f 2   O(h3 f )
2 
2
4
1
1
f ( x)dx  h f1  f 2 
3
3
3
• Simpsonの3/8則

x4
x1
• Bode則

x5
x1

f3   O(h5 f ( 4) )

9
9
3 
3
f ( x)dx  h f1  f 2  f3  f 4   O(h5 f ( 4) )
8
8
8 
8
64
24
64
14 
14
f ( x)dx  h f1 
f2 
f3 
f4 
f5   O(h7 f (6) )
45
45
45
45 
 45
一つの区間についての補外型公式

x1
x0

x1
x0

x1
x0

x1
x0
f ( x)dx  h f1  O(h2 f )
1 
3
f ( x)dx  h f1  f 2   O(h3 f )
2 
2
16
5
 23
f ( x)dx  h f1  f 2 
12
12
12

f3   O(h4 f (3) )

59
37
9 
 55
f ( x)dx  h f1 
f2 
f3 
f 4   O(h5 f ( 4) )
24
24
24 
 24
(閉じた)拡張公式
1
f
(
x
)
dx

h
x1
 2 f1  f 2  f3 
 (b  a)3 f  
1 

 f N 1  f N   O
2
2 
N


xN
13
5
f
(
x
)
dx

h
f

x1
12 1 12 f 2  f3  f 4 
13
5 
 1 
 f N 2  f N 1  f N   O 3 
12
12 
N 
xN
4
2
4
1
f
(
x
)
dx

h
f

f

f

x1
 3 1 3 2 3 3 3 f 4 
2
4
1 
 1 
 f N 2  f N 1  f N   O 4 
3
3
3 
N 
xN
•
拡張台形則
•
1/N3次の拡張公式
•
拡張Simpson則
(extended Simpson’s rule)
•
もう一つの拡張Simpson則(alternative extended Simpson’s rule)

xN
x1
59
43
49
17
f ( x)dx  h f1 
f2 
f3 
f 4  f5  f6 
48
48
48
48

49
43
59
17 
 1 
  f N 4 
f N 3 
f N 2 
f N 1 
f N   O 4 
48
48
48
48 
N 
(開いた)拡張公式

xN
x1
3
3

 1 
f ( x)dx  h f 2  f3  f 4   f N 2  f N 1   O 2 
2
2

N 

xN
x1

xN
x1

xN
x1
7
 23
f ( x)dx  h f 2  f3  f 4  f5 
12
12
7
23

 1 
 f N 3  f N 2 
f N 1   O 3 
12
12

N 
13
4
 27
f ( x)dx  h f 2  0  f 4  f5 
12
3
 12
4
13
27

 1 
 f N 4  f N 3  0 
f N 1   O 4 
3
12
12

N 
5
63
49
109
f ( x)dx  h
f 2  f3 
f4 
f5  f6  f7 
48
48
48
48

49
63
5
109

 1 
 f N 5 
f N 4 
f N 3  f N 2 
f N 1   O 4 
48
48
48
48

N 
別のタイプの拡張公式
• 拡張中点則(extended midpoint rule)

xN
x1


 1 
f ( x)dx  h f3 2  f5 2  f7 2   f N 3 2  f N 1 2  O 2 
N 
• 半分開いた公式(semi-open formula)

xN
x1
7
 23
f ( x)dx  h f 2  f3  f 4  f5 
12
12
13
5 
 1 
 f N 2  f N 1  f N   O 3 
12
12 
N 
拡張台形則によるアルゴリズム
•
•
•
•
1 
1
h f1  f 2  f3   f N 1  f N 
2 
2
さきほどの拡張台形則
ループの数Nが増えるごとに点が２倍になる
それに応じて合計する点の数が増え、精度が高くなる
誤差が十分小さくなったらループ終了、という使い方ができる
Euler-Maclaurin(オイラー・マクローリン)総和公式
(Euler-Maclaurin Summation Formula)によるアルゴリズム

xN
x1
1 
1
f ( x)dx  h f1  f 2  f3   f N 1  f N 
2 
2
B2k h2k ( 2k 1)
B2h2

( f N  f1)  
( fN
 f1( 2k 1) )  
2!
(2k )!
• B2kはBernoulli(ベルヌーイ)数(Bernoulli number)
• 今、N個の刻みでSNを得、続いて2N個の刻みでS2N
を得たとする
4
1
• S  3 S2N  3 SN で、最初の誤差項が相殺される
• 残った誤差は1/N4の大きさ=Simpson則
•
•
•
xN  x1
B2h2
であることに注意すると 
( f N  f1) は1/N2の大きさの誤差項で
N
2!
あることが分かる
h
誤差の級数にはhの偶数ベキの項しかあらわれていないので、その次の項は
1/N4の大きさの誤差項である
B2h2
( f N  f1) は1/4になる
刻み幅Nを2倍にすると、hは1/2になり、 
2!
Romberg積分
• さきほどは、刻み幅をNと2Nの値のみを使った
• Rombergの方法では、刻み幅を細かくしながらk回
拡張台形則を実行した結果をもとに、誤差級数を根
こそぎ取り除く
• Romberg積分の実装にはNevilleのアルゴリズム
(Neville’s algorithm)が使えるが、詳細は省略
• このように、パラメータhをいろいろな値にしてアルゴ
リズムを実行し、その結果からh=0での極限の値を
求めるという一般的なアイディアをRichardsonの極
限遅延接近(Richardson’s deferred approach to
the limit)という
変格積分
• 以下のいずれかが成り立つ場合、その積分は「変
格」であるという
– x=0でのsin(x)/xのように、正しく計算できない
– x=0でのx-1/2のように特異点がある
– 上限が∞か下限が－∞
• 第2Euler-Maclaurin総和公式(Second EulerMaclaurin summation formula)

xN
x1

f ( x)dx  h f3 2  f5 2  f 7 2   f N 3 2  f N 1 2

B2k h2k
B2h2

( f N  f1)  
(1  22k 1 )( f N( 2k 1)  f1( 2k 1) )
4
(2k )!
をもとにRomberg積分をするとうまくいく
Gauss(ガウス)の求積法
(Gaussian quadratures)
• Nを与えたとき、多項式f(x)の積分を求めたい
N
b
 W ( x) f ( x)dx   w f ( x )
a
i 1
i
i
• が厳密になるように重みwiと分点xiを探すことが目
的：等間隔でなくてもいいし、13/12みたいな重みで
はない
• W(x)は既知関数：これをかけることによってf(x)の特
異点をなくすことができる
W(x)
Gauss-(ガウス)
1
Legendre(ルジャンドル)
(1-x2)-1/2
Chebyshev(チェビシェフ)
xce-x
Laguerre(ラゲール)
e-x2
Hermite(エルミート)
Gauss-Legendre(ガウス・ルジャンドル)積分
(Gauss-Legendre integration)
•
•
•
•
•
N
wi f ( xi )
W(x)=1
x f ( x)dx  
i 1
分点xiと重みwiが求まれば、f(x)の積分が求まる
(i 1)Pi 1  (2i 1) xPi  iPi1 でPNを求める
PN=0となるような分点xiを求める
2
重みwiを wi 
により求める
2
2
x2
1
(1  xi )[PN ( xi )]
多次元積分
• 多次元積分の場合は下図のように計算すればよい
第５章 関数の計算
Evaluation of functions
級数と収束性
• 任意の関数はある点x0の近傍でベキ級数展開できる

f ( x)   ak ( x  x0 )k
k 0
• Aitken(エイトケン)のデルタ2乗過程(Aitken’s δ2-process)は、
級数の収束を加速する
• Sn-1,Sn,Sn+1という三つの連続する部分和があったら次のよう
(Sn1  Sn )2
に改良できる

Sn  Sn1 
Sn1  2Sn  Sn1
• 交代級数(alternating series)(各項の符号が順々に入れ替
わっているもの)はEuler(オイラー)変換(Euler’s

transformation)が有効  s
()s s
() u
s 0
s
 u0  u1  u2  un1  
s 0
s 1
2
[ un ]
• ここでΔは前進差分演算子(forward difference operator)
un  un1  un
2un  un2  2un1  un
3un  un3  3un2  3un1  un
連分数の計算
• 連分数
f ( x)  b0 
b1 
a1
a2
b2 
a3
b3 
a4
b4 
a5
b5  
• は以下の漸化式を計算すると求まる
An
fn 
Bn
A1  1
B1  0
A0  b0
B0  1
Aj  bj Aj 1  a j Aj 2
Bj  bj Bj 1  a j Bj 2
j  1,2,, n
多項式
• 多項式の値を求めるとき、次のようにする人はいま
い
c0+c1×x+c2×x×x+c3×x×x×x+c4×x×x×x×x
• 次のように計算すると効率的
c0+x×(c1+x×(c2+x×(c3+x×c4)))
Clenshaw(クレンショー)の漸化公式
(Clenshaw’s recurrence formula)
N
•
f ( x)   ck Fk ( x)
のような和を
求めたいとする(ckは係数)
• このFkがある関数α(n,x)と
β(n,x)に対して以下の漸化
式に従うものとする
Fn1 ( x)   (n, x)Fn ( x)   (n, x)Fn1 ( x)
k 0
• ここでyk(k=N,N-1,…,1)を次
の漸化式によって定義する
yN 2  yN 1  0
yk   (k , x) yk 1   (k 1, x) yk 2  ck
(k  N, N 1,,1)
• この式を最初の式に代入し
てckを消すと
f ( x)  
 [ y8   (8, x) y9   (9, x) y10]F8 ( x)
 [ y7   (7, x) y8   (8, x) y9 ]F7 ( x)
 [ y6   (6, x) y7   (7, x) y8 ]F6 ( x)
 [ y5   (5, x) y6   (6, x) y7 ]F5 ( x)

 [ y2   (2, x) y3   (3, x) y4 ]F2 ( x)
 [ y1   (1, x) y2   (2, x) y3 ]F1 ( x)
 [c0   (1, x) y2   (1, x) y2 ]F0 ( x)
• 残るのは
f ( x)   (1, x)F0 ( x) y2  F1 ( x) y1  F0 ( x)c0
2次方程式
• ２次方程式
の解は
ax2  bx  c  0
 b  b2  4ac
x
2a
x
2c
 b  b2  4ac
• 上の２通りのうち、片方だけを使うと解が不正確にな
ることがある
• 次のように計算すべき

1
q   b  sgn(b) b2  4ac
2

x1 
q
a
x2 
c
q
３次方程式
• ３次方程式
x3  a1x2  a2 x  a3  0 の解の求め方
2a13  9a1a2  27a3
a12  3a2
Q
R
9
54
• Q3-R2≧0のとき

  arccos R Q3
  a
x1  2 Q cos   1
 3 3

   2  a1
x2  2 Q cos

 3  3
   4  a1
x3  2 Q cos

 3  3
• R2-Q3>0のとき


13
x1   sgn(R) R2  Q3  | R | 




 a
 1
13
R2  Q3  | R |  3
Q

Chebyshev近似
• Chebyshev(チェビシェフ)多項式(Chebyshev polynomial)
T0 ( x)  1
T1 ( x)  x
T2 ( x)  2x 2 1
T3 ( x)  4x3  3x
T4 ( x)  8x 4  8x 2  1

Tn1 ( x)  2xTn ( x)  Tn1 ( x) n  1
• Chebyshev多項式で
関数を近似
N
 1
f ( x)   ckTk ( x)  c1
 k 1
 2
2 N    k 1 2   jk 1 2 
c j   f cos
 cos

N k 1  
N
N
 

Tn ( x)  cos(n arccosx)
(c) Daisuke Miyazaki 2003
All rights reserved.
http://www.cvl.iis.u-tokyo.ac.jp/