Numerical Recipes in C [日本語版] 初版

Presentation on theme: "Numerical Recipes in C [日本語版] 初版"— Presentation transcript:

1 Numerical Recipes in C [日本語版] 初版
宮崎大輔 2003年10月30日（木） PBVセミナー

2 Solution of linear algebraic equations
第２章 連立１次方程式の解法 Solution of linear algebraic equations

3 連立方程式 N=Mならば方程式の数と未知数の数が等しい→解が一意に定まる可能性がある
M<N(またはM=Nでも方程式が特異(singular))なら未知数の数より方程式の数の方が少ない→解は存在しないか無数の解が存在 M>N、未知数の数より方程式の数の方が多い→解は存在しない→最小２乗問題ととらえ解を求める

4 Gauss-Jordan(ガウス・ジョルダン)法 (Gauss-Jordan elimination)
複数の右辺ベクトルbについて方程式の解を求め、同時に逆行列A-1も求める Aの2行と対応するb、1の行を入れ替えても、xとYの解は変わらない Aの行を、その行と別の行の線形結合で置き換え、対応するb、1の行も線形結合で置き換えれば、解は変わらない Aの2列と対応するx、Yの行を入れ替えれば、解の行の順序が入れ替わる 上記3操作を組み合わせAを単位行列にすると右辺に解が得られる ピボット選択(pivoting): 0で割るのを防ぐ

5 Gaussの消去法(Gaussian elimination)と 後退代入(backsubstitution)
簡単に解ける（後退代入） 一般的には

6 LU分解 LU分解はCrout(クラウト)のアルゴリズム(Crout’s algorithm)を使う
Ax=bを解くのが簡単になる

7 連立方程式の解の反復改良(iterated improvement)
このxを求めたい 未知の誤差δxを含む解x+δxが得られたとする δxを求めれば得られた解x+δxから引いて真の解xが求まる ① ② ③ ①を②に代入 ②を③に代入 ④

8 ３重対角(tridiagonal)な連立１次方程式
メモリO(N)、計算時間O(N) 通常は計算時間O(N3) アルゴリズムの詳細は省略

9 Vandermonde(バンデルモンド)行列(Vandermonde matrices)
N個の点(xj,yj)に多項式を当てはめ、その多項式の係数ciを求める方程式 N個の点xiと値qjが与えられ、未知の重みwiを求め、N個のモーメントをqjに一致させる問題 計算時間はO(N2) 通常はO(N3) アルゴリズムの詳細は省略

10 Toeplitz(テープリッツ)行列(Toeplitz matrices)
計算時間はO(N2) 通常はO(N3) アルゴリズムの詳細は省略

11 特異値分解 (singular value decomposition, SVD)
M×N行列A M×Nの列直交行列U N×Nの対角行列W（対角成分は非負） N×Nの直交行列V U,Vは正規直交行列 Uの列、Wの要素、Vの列（つまりVTの行）を入れ替えただけ、は同じ分解とみなす もしWのいくつかの要素がたまたま等しければ、 それに対応するU,Vの列をその線形結合で置き換えたものは同じ分解とみなす

12 SVDによる逆行列と連立方程式の解 Aが正方行列(N×N)のとき また、 の解は ただし、wj=0のときは1/wj=0とする

13 零空間(nullspace)と値域(range)
(a) Aが非特異 Aが特異のとき [階数(rank)がN未満] 解は無数にある SVD解は0に最も近い解を吐き出す SVD解に零空間を足したものが解空間になる 零空間とは Ax=0の解空間 wj=0に対応するVの列の任意の線形結合 (b) Aが特異

14 SVDの良さ Aが特異でなくてもwjが小さいときはwjを0にしてSVD解を求めたほうが精度がいい
未知数より方程式が少ない場合は、足りない行に0を入れてN×N行列に対してSVD解を求める 未知数より方程式が多いときは、SVD解は最小２乗解となる

15 Sherman-Morrison(シャーマン・モリソン)公式 (Sherman-Morrison formula)
正方行列Aの逆行列A-1が求まったとする Aに小さな変更をしたとき、A-1を簡単に計算したい その変更が ならOK (u,vは何らかのベクトル) 例 uが基本単位ベクトルei→vをi行に加える vが基本単位ベクトルej→uをj列に加える u,vが基本単位ベクトルei,ejに比例→aijが加わる 計算方法：

16 Woodbury(ウッドベリ)の公式 (Woodbury formula)
修正すべき箇所が多い場合はSherman-Morrisonの公式を何回も計算する必要がある その場合はWoodburyの公式を使う ここで、AはN×N行列 U,VはN×P行列（P<N,通常P≪N） 逆行列の計算はP×P行列だけで済む

17 Interpolation and extrapolation
第３章 補間と補外 Interpolation and extrapolation

18 補間：低次の多項式か高次の多項式か

19 Lagrange(ラグランジュ)の公式 N個の点 を通る次数N-1の補間多項式は
この公式を効率的に計算するNevilleのアルゴリズムというものがあるが、詳細は省略

20 有理関数による補間 多項式ではうまく近似できない関数でも有理関数でうまく近似できることがある m+1個の点 を通る以下の関数を求める
Bulirsch-Stoer(ブリルシュ-ストア)アルゴリズムを使えばよいが、詳細は省略

21 3次スプライン補間 入力 を上の３次多項式で関数を近似する 以下の式と２つの条件からyi’’が求まる
入力 を上の３次多項式で関数を近似する 以下の式と２つの条件からyi’’が求まる 端点（y1’’,yN’’）の片方または両方を0と置く 端点の片方または両方の１階導関数を与えてやる この連立方程式は３重対角行列になるのでO(N)で解ける

22 双１次補間(bilinear interpolation)
2次元以上の補間(ここでは2次元の場合を説明するが、3次元以上でも同様)

23 双３次補間(bicubic interpolation)
滑らかさを上げる高次補間 各格子点で関数 以外に 勾配 クロス導関数 を人間が与え、それを元に補間する

24 双３次スプライン(bicubic spline)
滑らかさを上げる高次補間 双３次補間に似てるが、格子点での導関数の値を計算してしまう まず行ごとにスプライン補間し その新しい列についてスプライン補間する

25 Integration of functions
第４章 関数の積分 Integration of functions

26 閉じた公式と開いた公式

27 閉じたNewton-Cotes(ニュートン・コーツ)公式 (Newton-Cotes formulas)
台形則 Simpson則 Simpsonの3/8則 Bode則

28 一つの区間についての補外型公式

29 (閉じた)拡張公式 拡張台形則 1/N3次の拡張公式 拡張Simpson則 (extended Simpson’s rule)
もう一つの拡張Simpson則(alternative extended Simpson’s rule)

30 (開いた)拡張公式

31 別のタイプの拡張公式 拡張中点則(extended midpoint rule) 半分開いた公式(semi-open formula)

32 拡張台形則によるアルゴリズム さきほどの拡張台形則 ループの数Nが増えるごとに点が２倍になる
それに応じて合計する点の数が増え、精度が高くなる 誤差が十分小さくなったらループ終了、という使い方ができる

33 B2kはBernoulli(ベルヌーイ)数(Bernoulli number)
Euler-Maclaurin(オイラー・マクローリン)総和公式 (Euler-Maclaurin Summation Formula)によるアルゴリズム B2kはBernoulli(ベルヌーイ)数(Bernoulli number) 今、N個の刻みでSNを得、続いて2N個の刻みでS2Nを得たとする で、最初の誤差項が相殺される 残った誤差は1/N4の大きさ=Simpson則 であることに注意すると は1/N2の大きさの誤差項であることが分かる 誤差の級数にはhの偶数ベキの項しかあらわれていないので、その次の項は1/N4の大きさの誤差項である 刻み幅Nを2倍にすると、hは1/2になり、 は1/4になる

34 Romberg積分 さきほどは、刻み幅をNと2Nの値のみを使った
Rombergの方法では、刻み幅を細かくしながらk回拡張台形則を実行した結果をもとに、誤差級数を根こそぎ取り除く Romberg積分の実装にはNevilleのアルゴリズム(Neville’s algorithm)が使えるが、詳細は省略 このように、パラメータhをいろいろな値にしてアルゴリズムを実行し、その結果からh=0での極限の値を求めるという一般的なアイディアをRichardsonの極限遅延接近(Richardson’s deferred approach to the limit)という

35 変格積分 以下のいずれかが成り立つ場合、その積分は「変格」であるという
x=0でのsin(x)/xのように、正しく計算できない x=0でのx-1/2のように特異点がある 上限が∞か下限が－∞ 第2Euler-Maclaurin総和公式(Second Euler-Maclaurin summation formula) をもとにRomberg積分をするとうまくいく

36 Gauss(ガウス)の求積法 (Gaussian quadratures)
Nを与えたとき、多項式f(x)の積分を求めたい が厳密になるように重みwiと分点xiを探すことが目的：等間隔でなくてもいいし、13/12みたいな重みではない W(x)は既知関数：これをかけることによってf(x)の特異点をなくすことができる W(x) Gauss-(ガウス) 1 Legendre(ルジャンドル) (1-x2)-1/2 Chebyshev(チェビシェフ) xce-x Laguerre(ラゲール) e-x2 Hermite(エルミート)

37 Gauss-Legendre(ガウス・ルジャンドル)積分 (Gauss-Legendre integration)
W(x)=1 分点xiと重みwiが求まれば、f(x)の積分が求まる でPNを求める PN=0となるような分点xiを求める 重みwiを により求める

38 多次元積分 多次元積分の場合は下図のように計算すればよい

39 Evaluation of functions
第５章 関数の計算 Evaluation of functions

40 級数と収束性 任意の関数はある点x0の近傍でベキ級数展開できる
Aitken(エイトケン)のデルタ2乗過程(Aitken’s δ2-process)は、級数の収束を加速する Sn-1,Sn,Sn+1という三つの連続する部分和があったら次のように改良できる 交代級数(alternating series)(各項の符号が順々に入れ替わっているもの)はEuler(オイラー)変換(Euler’s transformation)が有効 ここでΔは前進差分演算子(forward difference operator)

41 連分数の計算 連分数 は以下の漸化式を計算すると求まる

42 多項式 多項式の値を求めるとき、次のようにする人はいまい c0+c1×x+c2×x×x+c3×x×x×x+c4×x×x×x×x
次のように計算すると効率的 c0+x×(c1+x×(c2+x×(c3+x×c4)))

43 Clenshaw(クレンショー)の漸化公式 (Clenshaw’s recurrence formula)
のような和を求めたいとする(ckは係数) このFkがある関数α(n,x)とβ(n,x)に対して以下の漸化式に従うものとする ここでyk(k=N,N-1,…,1)を次の漸化式によって定義する この式を最初の式に代入してckを消すと 残るのは

44 2次方程式 ２次方程式 の解は 上の２通りのうち、片方だけを使うと解が不正確になることがある 次のように計算すべき

45 ３次方程式 ３次方程式 の解の求め方 Q3-R2≧0のとき R2-Q3>0のとき

46 Chebyshev近似 Chebyshev(チェビシェフ)多項式(Chebyshev polynomial)

47 Daisuke Miyazaki 2003 Creative Commons Attribbution 4
Daisuke Miyazaki 2003 Creative Commons Attribbution 4.0 International License.