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

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

9th Physics Based Vision Seminar NUMERICAL RECIPES in C [日本語版]ニューメリカルレシピ・イン・シー C言語による数値計算のレシピ • William H. Press, Saul A. Teukolsky William T. Vetterling, Brian P. Flannery • 丹慶勝市，奥村晴彦，佐藤俊郎，小林誠 • 1993 • 技術評論社 • ISBN4-87408-560-1

第16章２点境界値問題 Two point boundary value problems

9th Physics Based Vision Seminar ２点境界値問題(two point boundary value problem) • 微分方程式： • 境界条件：を解が点x1,x2で満たすねらい撃ち法，射撃法 (shooting method) 緩和法 (relaxation method)

9th Physics Based Vision Seminar 点x1の境界条件を満たす任意のベクトルをVとする Vが決まると点x1での関数の値が分かり，それを初期値としてRunge-Kutta法などを使って計算していくと，点x2での関数の値が分かるこのとき，点x2での境界条件が満たされているか調べるもし満たされていたら計算終了もし満たされていなかったらVを変化させ，もう一度Runge-Kuttaを計算しなおす点x2での「ずれのベクトル」をFとする F=0なら境界条件を満たしている F≠0なら境界条件を満たしていない F=0となるようにNewton-Raphson法で解を求める計算は以下の通りねらい撃ち法

9th Physics Based Vision Seminar 適合点へのねらい撃ち(shooting to a fitting point) • x1からx2へと解く代わりに，まずx1からx1とx2の間にある点xfへと解き，次にx2からxfへと（逆向きに）解いていく

9th Physics Based Vision Seminar 常微分方程式を差分方程式(finite difference equations, FDE)に換える例： k=1,2,…,Mの点の格子（メッシュ），各点xk，最初の境界点x1，最後の境界点xM 点xkでの従属変数をykとする点kで成り立つ：点1で成り立つ：点Mで成り立つ： yの初期値を与え，各反復で新しい値をy+Δyとする点kでΔyを計算：点1でΔyを計算：点MでΔyを計算：この連立１次方程式をGaussの消去法と後退代入で解くことができる緩和法(relaxation methods)

9th Physics Based Vision Seminar メッシュ点の自動割当て • 解が急激に変化する所は格子点を多く取り，その他のところは少なくしてよい→動的に格子（メッシュ）を割り当てる(allocate)

9th Physics Based Vision Seminar 内点での境界条件，特異点の処理 • でD=0の場合→特異点 • 物理の問題では，D→0のときN→0となることが多い • 正則条件(regularity condition) • 特異点があるとき以下で解ける • 緩和法 • 適合点へのねらい撃ち法 • 特異点の場所が未知の場合はどうするか？ • 詳しい方法は省略するが，未知の内点xsにおいて制約条件N=0,D=0を加えれば良い

第17章偏微分方程式 Partial differential equations

9th Physics Based Vision Seminar 初期値問題(initial value problem), Cauchy(コーシー)問題(Cauchy problem) 双曲型(hyperbolic) 波動方程式(wave equation) 放物型(parabolic) 拡散方程式(diffusion equation) 境界値問題(boundary value problem) 楕円型(elliptic) Poisson(ポアソン)方程式(Poisson equation) Laplace(ラプラス)方程式(Laplace’s equation) 境界条件 Dirichlet(ディリクレ)条件(Dirichlet condition)：境界点での値 Neumann(ノイマン)条件(Neumann condition)：境界点での勾配偏微分方程式

9th Physics Based Vision Seminar FTCS(Forward Time Centered Space, 前進時間中心空間)表現 • 流束保存方程式(flux-conservative equation)の一例を解く方法を示す→一般の場合にも適用可能 • xをΔx間隔，tをΔt間隔にとった格子を考える • 点(tn,xj)でのuの値u(tn,xj)をujnと表す • 差分化 • 時間： • 空間： • 以上から • これは不安定でうまくいかない

9th Physics Based Vision Seminar Lax法 • を使うとFTCS法は • von Neumann(フォンノイマン)安定解析(von Neumann stability analysis)を使うと，安定であるための条件は • Courant-Friedrichs-Lewy(クーラント・フリードリクス・レヴイ)の安定性基準，Courant条件(Courant condition) • 波動方程式の速度v，時間格子サイズΔt，空間格子サイズΔx，で決まる

9th Physics Based Vision Seminar 風上差分(upwind differencing) • 安定性の条件はCourant条件と同じ

9th Physics Based Vision Seminar 互い違い蛙跳び(staggered leapfrog)法 • 安定性の条件はCourant条件と同じ

9th Physics Based Vision Seminar ２ステップLax-Wendroff(ラックス・ウェンドロフ)(Two-Step Lax-Wendroff)スキーム • Lax法＋互い違い蛙跳び法 • まとめると • 安定性の基準はCourant条件

9th Physics Based Vision Seminar 拡散方程式差分方程式に換えると：これはFTCSスキーム，全陽的(fully explicit) 安定性基準は完全陰的(fully implicit), 後退時間(backward time) ３重対角の連立方程式を解く任意の刻み幅Δtに対し，無条件に安定陽的と陰的のFTCSスキームの平均Crank-Nicholson(クランク・ニコルソン)スキーム(Crank-Nicholson scheme) 任意のΔtに対し安定拡散初期値問題

9th Physics Based Vision Seminar Schrödinger(シュレーディンガー)方程式 • 陰的スキーム • 無条件に安定だがユニタリ(unitary)ではない • 物理学の問題ではと正規化されていないといけない • と書き直す • ケーリーの形式(Cayley’s form)というのを用いると，次式を計算すればよいことが分かる • Hはxについての差分近似で置き換える • 複素数の３重対角の連立方程式を解けば良い • 方法は安定でユニタリ • Crank-Nicholson法

9th Physics Based Vision Seminar 多次元の初期値問題 • 1次元100個の格子点，2次元100×100の格子点，多次元計算時間がかなりかかる • 小さな格子(8×8とか)でまず試してみて，うまくいったら格子サイズを大きくしたらいいだろう • 詳細な定式化は本を見てください

9th Physics Based Vision Seminar 緩和法 • 以下では，境界値問題としての式を解くことを考える • より一般の場合については本文を見てほしい

9th Physics Based Vision Seminar Jacobi(ヤコビ)法(Jacobi’s method) • J×Jの格子の場合 • スペクトル半径(spectral radius)ρsは • 収束率を表し，１回の反復で残差が減る割合である • 1未満でないと収束しない • 0と1の間の値 • 例：ρs=0.9なら反復ごとに残差が0.9倍になる • 誤差が10-p倍に減少するのに必要な反復回数rは • 例：J=400,p=15のときr=1200000

9th Physics Based Vision Seminar Gauss-Seidel(ガウス・ザイデル)法(Gauss-Seidel method) • J×Jの格子の場合 • スペクトル半径(spectral radius)ρsは • 誤差が10-p倍に減少するのに必要な反復回数rは • 例：J=400,p=15のときr=600000

9th Physics Based Vision Seminar SOR法(successive overrelaxation, 逐次過緩和法) • ω:過緩和パラメータ(overrelaxation parameter) • 0<ω<2の場合だけ収束 • 0<ω<1:劣緩和(underrelaxation),Gauss-Seidel法より収束が遅い • 1<ω<2:過緩和(overrelaxation),Gauss-Seidel法より収束が速い • J×Jの格子の場合 • 最適なωは • 例:J=400のときω=1.984 • スペクトル半径(spectral radius)ρSORは • 誤差が10-p倍に減少するのに必要な反復回数rは • 例：J=400,p=15のときr=2000 • Chebyshev(チェビシェフ)加速(Chebyshev acceleration)によりωを反復ごとに変化させることによりさらに収束を速くする

9th Physics Based Vision Seminar 交互方向陰的方法(ADI)(alternating-direction implicit method) • SORの(8log10J)/J倍の演算回数で済む • 例:J=400のとき反復回数は104 • Poisson方程式を解くためにNumerical Recipes in C[日本語版]に載っているプログラムを実行するのに最低限必要な知識 • tridag関数のプログラムにバグがあるので注意 • uが求めたい値(初期値は0にでもしておく) • 物体領域内では，a=-1,b=2,c=-1,d=-1,e=2,f=-1,g=ρ • 物体領域外では，a=0,b=1,c=0,d=0,e=1,f=0,g=0 • jmaxは格子サイズjmax×jmax • epsは収束判定に使う値(例:eps=10-5) • 例:J=400のときα=6.16847e-5,β=3.999938315 • 僕の場合は，smallval=1e-5,α=smallval,β=4-smallvalにした • N=2k,Nはln(4J/π)に近い2の累乗 • 例：J=400のときln(4J/π)は6.23なのでN=2k=23=8にした

9th Physics Based Vision Seminar 以下にADI法のアルゴリズムを示すが，分かりやすく説明するために多少間違ったことが書かれているので詳しくは本文を読むこと 3に戻るただし，rはk=1の場合，以下を使う kが1じゃないときは本文を見てください rはN=2k個の数列を繰り返す例：k=3のとき(N=2k=23=8) 1回目の反復計算ではr1を使う 8回目ではr8を使い，9回目ではr1を使う ADI法のアルゴリズムの概要

9th Physics Based Vision Seminar Multigrid Methods • ADI法よりも高速な手法が提案されている • multigrid method • FMG(full multigrid) method

9th Physics Based Vision Seminar 不完全Choleski(コレスキー)共役勾配法(incomplete Choleski conjugate gradient method, ICCG) • 格子が小さいならAx=bの形に書き換えてICCG法などで解けばよい • 格子が大きいなら記憶容量の関係で不可能 • 例えばJ=400の場合は，格子は400×400=160000となるが，行列Aのサイズが，160000×160000となる．double型が8バイトのとき，行列Aのメモリサイズは，190GBとなり，32ビット計算機のメモリ空間には入りきらない • Aはほとんどの問題で疎行列なので工夫すればなんとかなることもある

9th Physics Based Vision Seminar 高速な方法 • Fourier(フーリエ)変換法(Fourier transform method) • 巡回還元法(cyclic reduction method) • FACR(Fourier Analysis and Cyclic Reduction):上の２つを合わせた方法，詳細は省略 • いずれも，境界がx軸とy軸に平行でなければならない→長方形 • multigrid法のほうが高速の場合もある

9th Physics Based Vision Seminar 時間空間ではなく周波数空間で解く時間空間の場合周波数空間の場合 Jはx方向のサイズ，Lはy方向のサイズ，mはxの位置，nはyの位置，Δはサンプリング間隔，はをFourier変換したもの，はをFourier変換したものをFourier変換してを計算するを左下の式で計算を逆Fourier変換してを求めるこの方法は周期的境界条件に対して有効 Dirichlet境界条件u=0のときはサイン変換を使う Neumannの境界条件∇u=0のときはコサイン変換を使う Fourier変換法

9th Physics Based Vision Seminar 巡回還元(cyclic reduction, CR) • Poisson方程式を差分方程式で書くと • これを行列の形で書くとこれを以下のように記す

9th Physics Based Vision Seminar 巡回還元(cyclic reduction, CR) • を3つ並べるとこれらを変形して足し合わせると以下のような形になる(T(1)はTから計算でき，gj(1)はgj-1とgjとgj+1とTから計算できる)

9th Physics Based Vision Seminar 巡回還元(cyclic reduction, CR) • さきほど求めた式をまた3つ並べるとこれらを変形して足し合わせると以下のような形になる

9th Physics Based Vision Seminar 巡回還元(cyclic reduction, CR) • 格子のサイズが2の累乗だとすると，最終的に中央の行の方程式１個だけに変形される • なお，u0,uJは既知の境界値 • ３重対角アルゴリズムでuJ/2を解くことができる

9th Physics Based Vision Seminar 巡回還元(cyclic reduction, CR) • その一つ前のレベルf-1では，以下の３つの式が出るが，これも３重対角アルゴリズムで解けるこれを最初のレベルまで解くことにより全てのuを得る ←解ける ←解く必要なし ←解ける

終了

9th Physics Based Vision Seminar 次回 • ９月上旬を予定（ずいぶん間が空いてしまうが） • 発表者２名＋宮崎 • 宮崎はウェーブレットをやろうかと思ってる • ７月１日ニューメリカルレシピ • ９月上旬ウェーブレット前半 • １０月中旬ウェーブレット後半 • ２月変分法と有限要素法 • ３月外積展開？

9th Physics Based Vision Seminar ウェーブレットによる信号処理と画像処理 • 中野宏毅，山本鎭男，吉田靖夫 • 1999 • 共立出版 • ISBN4-320-08549-3

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