東京大学 The University of Tokyo 池内研究室 Computer Vision Laboratory Numerical Recipes in C [日本語版] 初版 宮崎大輔 2004年7月1日(木) PBVセミナー 9th Physics Based Vision Seminar 2004/July/1 NUMERICAL RECIPES in C [日本語版] ニューメリカルレシピ・イン・シー C言語による数値計算のレシピ • William H. Press, Saul A. Teukolsky William T. Vetterling, Brian P. Flannery • 丹慶勝市,奥村晴彦,佐藤 俊郎,小林誠 • 1993 • 技術評論社 • ISBN4-87408-560-1 2 東京大学 The University of Tokyo 池内研究室 Computer Vision Laboratory 第16章 2点境界値問題 Two point boundary value problems 9th Physics Based Vision Seminar 2004/July/1 2点境界値問題 (two point boundary value problem) • 微分方程式: • 境界条件: を解が点x1,x2で満たす ねらい撃ち法,射撃法 (shooting method) 緩和法 (relaxation method) 4 9th Physics Based Vision Seminar 2004/July/1 ねらい撃ち法 • 点x1の境界条件を満たす任 意のベクトルをVとする • Vが決まると点x1での関数 の値が分かり,それを初期 値としてRunge-Kutta法な どを使って計算していくと, 点x2での関数の値が分かる • このとき,点x2での境界条 件が満たされているか調べ る • 点x2での「ずれのベクトル」をFと する – F=0なら境界条件を満たしている – F≠0なら境界条件を満たしていな い • F=0となるようにNewtonRaphson法で解を求める • 計算は以下の通り – もし満たされていたら計算終 了 – もし満たされていなかったらV を変化させ,もう一度RungeKuttaを計算しなおす 5 9th Physics Based Vision Seminar 2004/July/1 適合点へのねらい撃ち (shooting to a fitting point) • x1からx2へと解く代わりに,まずx1からx1とx2の間に ある点xfへと解き,次にx2からxfへと(逆向きに)解い ていく 6 9th Physics Based Vision Seminar 2004/July/1 緩和法(relaxation methods) • 常微分方程式を差分方程式 (finite difference equations, FDE)に換える 例: • 点kでΔyを計算: 点1でΔyを計算: • k=1,2,…,Mの点の格子(メッシュ), 各点xk,最初の境界点x1,最後の 境界点xM • 点xkでの従属変数をykとする • 点kで成り立つ: 点1で成り立つ: 点Mで成り立つ: • yの初期値を与え,各反復で新し い値をy+Δyとする 点MでΔyを計算: • この連立1次方程式をGaussの 消去法と後退代入で解くことがで きる 7 9th Physics Based Vision Seminar 2004/July/1 メッシュ点の自動割当て • 解が急激に変化する所は格子点を多く取り,その他 のところは少なくしてよい→動的に格子(メッシュ)を 割り当てる(allocate) 8 9th Physics Based Vision Seminar 2004/July/1 内点での境界条件,特異点の処理 • でD=0の場合→特異点 • 物理の問題では,D→0のときN→0となることが多い – 正則条件(regularity condition) • 特異点があるとき以下で解ける – 緩和法 – 適合点へのねらい撃ち法 • 特異点の場所が未知の場合はどうするか? • 詳しい方法は省略するが,未知の内点xsにおいて制 約条件N=0,D=0を加えれば良い 9 東京大学 The University of Tokyo 池内研究室 Computer Vision Laboratory 第17章 偏微分方程式 Partial differential equations 9th Physics Based Vision Seminar 2004/July/1 偏微分方程式 • • 初期値問題(initial value problem), Cauchy(コーシー)問題(Cauchy problem) – 双曲型(hyperbolic) • 波動方程式(wave equation) – 放物型(parabolic) • 拡散方程式(diffusion equation) • 境界条件 – – Dirichlet(ディリクレ)条件(Dirichlet condition):境界点での値 Neumann(ノイマン)条件(Neumann condition):境界点での勾配 境界値問題(boundary value problem) – 楕円型(elliptic) • Poisson(ポアソン)方程式(Poisson equation) • Laplace(ラプラス)方程式(Laplace’s equation) 2u 2u 0 x 2 y 2 11 9th Physics Based Vision Seminar 2004/July/1 FTCS(Forward Time Centered Space, 前進時 間中心空間)表現 • 流束保存方程式(flux-conservative equation)の一 例 を解く方法を示す→一般の場合にも適用 可能 • xをΔx間隔,tをΔt間隔にとった格子を考える • 点(tn,xj)でのuの値u(tn,xj)をujnと表す • 差分化 – 時間: – 空間: • 以上から • これは不安定でうまくいかない 12 9th Physics Based Vision Seminar 2004/July/1 Lax法 • を使うとFTCS法は • von Neumann(フォンノイマン)安定解析(von Neumann stability analysis)を使うと,安定であるた めの条件は – Courant-Friedrichs-Lewy(クーラント・フリードリクス・レヴ イ)の安定性基準,Courant条件(Courant condition) – 波動方程式の速度v,時間格子サイズΔt,空間格子サイズ Δx,で決まる 13 9th Physics Based Vision Seminar 2004/July/1 風上差分(upwind differencing) • 安定性の条件はCourant条件と同じ 14 9th Physics Based Vision Seminar 2004/July/1 互い違い蛙跳び(staggered leapfrog)法 • 安定性の条件はCourant条件と同じ 15 9th Physics Based Vision Seminar 2004/July/1 2ステップLax-Wendroff(ラックス・ウェンドロフ) (Two-Step Lax-Wendroff)スキーム • Lax法+互い違い蛙跳び法 1 vt u u u u u • 2 2x vt u u u u n 1 2 j 1 2 n 1 j n j 1 n j x n j n 1 2 j 1 2 n j 1 n j n 1 2 j 1 2 • まとめると u u vxt 12 u u 2vxt u • 安定性の基準はCourant条件 n 1 j n j n j 1 n j n j 1 12 u u nj n j u nj1 vt n u j u nj1 2x 16 9th Physics Based Vision Seminar 2004/July/1 拡散初期値問題 • 拡散方程式 • 差分方程式に換えると: • 陽的と陰的のFTCSスキームの 平均 これはFTCSスキーム,全陽的 (fully explicit) • 安定性基準は • Crank-Nicholson(クランク・ニコ ルソン)スキーム(CrankNicholson scheme) • 任意のΔtに対し安定 完全陰的(fully implicit), 後退時 間(backward time) • 3重対角の連立方程式を解く • 任意の刻み幅Δtに対し,無条件 に安定 17 9th Physics Based Vision Seminar 2004/July/1 Schrödinger(シュレーディンガー)方程式 • • 陰的スキーム • 無条件に安定だがユニタリ(unitary)ではない • 物理学の問題では と正規化されていないといけな い • と書き直す • ケーリーの形式(Cayley’s form)というのを用いると,次式を 計算すればよいことが分かる – – – – Hはxについての差分近似で置き換える 複素数の3重対角の連立方程式を解けば良い 方法は安定でユニタリ Crank-Nicholson法 18 9th Physics Based Vision Seminar 2004/July/1 多次元の初期値問題 • 1次元100個の格子点,2次元100×100の格子点, 多次元計算時間がかなりかかる • 小さな格子(8×8とか)でまず試してみて,うまくいっ たら格子サイズを大きくしたらいいだろう • 詳細な定式化は本を見てください 19 9th Physics Based Vision Seminar 2004/July/1 緩和法 • 以下では,境界値問題として の式を解くことを考える • より一般の場合については本文を見てほしい 20 9th Physics Based Vision Seminar 2004/July/1 Jacobi(ヤコビ)法(Jacobi’s method) • • J×Jの格子の場合 – スペクトル半径(spectral radius)ρsは • 収束率を表し,1回の反復で残差が減る割合である • 1未満でないと収束しない • 0と1の間の値 • 例:ρs=0.9なら反復ごとに残差が0.9倍になる – 誤差が10-p倍に減少するのに必要な反復回数rは • 例:J=400,p=15のときr=1200000 21 9th Physics Based Vision Seminar 2004/July/1 Gauss-Seidel(ガウス・ザイデル)法 (Gauss-Seidel method) • • J×Jの格子の場合 – スペクトル半径(spectral radius)ρsは – 誤差が10-p倍に減少するのに必要な反復回数rは • 例:J=400,p=15のときr=600000 22 9th Physics Based Vision Seminar 2004/July/1 SOR法 (successive overrelaxation, 逐次過緩和法) • • ω:過緩和パラメータ(overrelaxation parameter) • 0<ω<2の場合だけ収束 1 1 H xn,y1 H xn1, y H xn11, y H xn, y 1 H xn,y11 xy (1 ) H xn, y 4 4 – 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)によりωを反復 ごとに変化させることによりさらに収束を速くする 23 9th Physics Based Vision Seminar 2004/July/1 交互方向陰的方法(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 1 21 cos 21 cos J J • 例: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にした 24 9th Physics Based Vision Seminar 2004/July/1 ADI法のアルゴリズムの概要 • 以下にADI法のアルゴリ ズムを示すが,分かりや すく説明するために多少 間違ったことが書かれて いるので詳しくは本文を読 むこと 1. ux, y 0 2. x, y 0 1 u 3. x, y 2 r u x1, y u x1, y x, y x, y 4. x, y x, y 2rux, y 1 u x, y 1 u x, y 1 x, y 5. u x, y 2r 6. x, y x, y 2rux, y 7. 3に戻る • ただし,rはk=1の場合,以 下を使う 1. 2. r1 r2 2 2 2 2 2 2 2 2 2 2 – kが1じゃないときは本文を見 てください – rはN=2k個の数列を繰り返す – 例:k=3のとき(N=2k=23=8) • 1回目の反復計算ではr1 を使う • 8回目ではr8を使い,9回 目ではr1を使う 25 9th Physics Based Vision Seminar 2004/July/1 Multigrid Methods • ADI法よりも高速な手法が提案されている – multigrid method – FMG(full multigrid) method 26 9th Physics Based Vision Seminar 2004/July/1 不完全Choleski(コレスキー)共役勾配法(incomplete Choleski conjugate gradient method, ICCG) • 格子が小さいならAx=bの形に書き換えてICCG法な どで解けばよい • 格子が大きいなら記憶容量の関係で不可能 – 例えばJ=400の場合は,格子は400×400=160000となる が,行列Aのサイズが,160000×160000となる.double 型が8バイトのとき,行列Aのメモリサイズは,190GBとな り,32ビット計算機のメモリ空間には入りきらない – Aはほとんどの問題で疎行列なので工夫すればなんとか なることもある 27 9th Physics Based Vision Seminar 2004/July/1 高速な方法 • Fourier(フーリエ)変換法(Fourier transform method) • 巡回還元法(cyclic reduction method) • FACR(Fourier Analysis and Cyclic Reduction):上 の2つを合わせた方法,詳細は省略 • いずれも,境界がx軸とy軸に平行でなければならな い→長方形 • multigrid法のほうが高速の場合もある 28 9th Physics Based Vision Seminar 2004/July/1 Fourier変換法 • 時間空間ではなく周波数空 間で解く • 時間空間の場合 • 周波数空間の場合 – Jはx方向のサイズ,Lはy方 向のサイズ,mはxの位置,n はyの位置,Δはサンプリング 間隔,ˆ は をFourier変換し たもの, uˆ は u をFourier変換 したもの 1. をFourier変換して ˆ を計 算する 2. uˆ を左下の式で計算 3. uˆ を逆Fourier変換して u を 求める • この方法は周期的境界条 件に対して有効 • • Dirichlet境界条件u=0のと きはサイン変換を使う Neumannの境界条件 ∇u=0のときはコサイン変 換を使う 29 9th Physics Based Vision Seminar 2004/July/1 巡回還元(cyclic reduction, CR) • 2u 2u Poisson方程式 x 2 y 2 ( x, y) を差分方程式で書くと u j1,l u j1,l u j,l 1 u j,l 1 4u j,l 2 j,l • これを行列の形で書くと 2 u j 1,0 u j ,0 u j 1,0 j ,0 u 2 u u 1 4 1 j 1,l 1 j ,l 1 j 1,l 1 2 j ,l 1 u j 1,l 0 0 1 4 1 0 0 u j ,l u j 1,l j ,l u 2 u u 1 4 1 j 1 , l 1 j , l 1 j 1 , l 1 j , l 1 2 u u u j , L j 1, L j ,L j 1, L これを以下のように記す u j 1 T u j u j 1 g j 2 30 9th Physics Based Vision Seminar 2004/July/1 巡回還元(cyclic reduction, CR) • u j 1 T u j u j 1 g j 2 を3つ並べると u j 2 T u j 1 u j g j 12 u j 1 T u j u j 1 g j 2 u j T u j 1 u j 2 g j 12 これらを変形して足し合わせると以下のような形にな る(T(1)はTから計算でき,gj(1)はgj-1とgjとgj+1とTから 計算できる) u j 2 T(1) u j u j 2 g(j1) 2 31 9th Physics Based Vision Seminar 2004/July/1 巡回還元(cyclic reduction, CR) • さきほど求めた式をまた3つ並べると u j 4 T(1) u j 2 u j g(j1)22 u j 2 T(1) u j u j 2 g(j1) 2 u j T(1) u j 2 u j 4 g(j1) 22 これらを変形して足し合わせると以下のような形にな る u j 4 T(2) u j u j 4 g(j2) 2 32 9th Physics Based Vision Seminar 2004/July/1 巡回還元(cyclic reduction, CR) • 格子のサイズが2の累乗だとすると,最終的に中央 の行の方程式1個だけに変形される • なお,u0,uJは既知の境界値 • 3重対角アルゴリズムでuJ/2を解くことができる u0 T( f ) u J / 2 u J g(J f/ )22 33 9th Physics Based Vision Seminar 2004/July/1 巡回還元(cyclic reduction, CR) • その一つ前のレベルf-1では,以下の3つの式が出る が,これも3重対角アルゴリズムで解ける u0 T( f 1) u J / 4 u J / 2 g(J f/ 21) 2 u J / 4 T( f 1) u J / 2 u3J / 4 g(J f/ 21) 2 ←解ける u J / 2 T( f 1) u3J / 4 u J g3( Jf / 14) 2 ←解ける ←解く必要なし これを最初のレベルまで解くことにより全てのuを得 る 34 東京大学 The University of Tokyo 池内研究室 Computer Vision Laboratory 終了 9th Physics Based Vision Seminar 2004/July/1 次回 • 9月上旬を予定(ずいぶん間が空いてしまうが) • 発表者2名+宮崎 • 宮崎はウェーブレットをやろうかと思ってる – – – – – 7月1日ニューメリカルレシピ 9月上旬ウェーブレット前半 10月中旬ウェーブレット後半 2月変分法と有限要素法 3月外積展開? 36 9th Physics Based Vision Seminar 2004/July/1 ウェーブレットによる信号処理と画像処理 • 中野宏毅,山本鎭男,吉田 靖夫 • 1999 • 共立出版 • ISBN4-320-08549-3 37 東京大学 The University of Tokyo 池内研究室 Computer Vision Laboratory © Daisuke Miyazaki 2004 All rights reserved. http://www.cvl.iis.u-tokyo.ac.jp/
© Copyright 2025 ExpyDoc