Miyazaki-PBV200407
Download
Report
Transcript Miyazaki-PBV200407
東京大学
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
x2 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 n
vt n
n1 2
n
n
u
u
u
u
u
j
1
2
j
1
j
j
1
j
•
2
2x
u nj1 u nj
vt n1 2 n1 2
u
u j 1 2
x j 1 2
• まとめると u u vxt 12 u u 2vxt u
• 安定性の基準はCourant条件
n1
j
n
j
n
j 1
n
j
n
j 1
12 u
u nj
n
j
u nj1
vt n n
u u j 1
2x j
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, y1 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)
21 cos
J
21 cos
J 1
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法のアルゴリズムの概要
•
1.
2.
3.
4.
5.
6.
7.
以下にADI法のアルゴリ
ズムを示すが,分かりや
すく説明するために多少
間違ったことが書かれて
いるので詳しくは本文を読
むこと
ux, y 0
x, y 0
1
u u x, y x, y
2 r x1, y x1, y
x, y x, y 2rux, y
1
u u x, y
ux , y
2 r x, y1 x, y1
x, y x, y 2rux, y
ux , y
• ただし,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を使う
3に戻る
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 をFourier変換
したもの
1. をFourier変換して ̂ を計
算する
2. û を左下の式で計算
3. û を逆Fourier変換して u を
求める
• この方法は周期的境界条
件に対して有効
•
•
Dirichlet境界条件u=0のと
きはサイン変換を使う
Neumannの境界条件
∇u=0のときはコサイン変
換を使う
29
9th Physics Based Vision Seminar
2004/July/1
巡回還元(cyclic reduction, CR)
•
2u 2u
Poisson方程式 x2 y2 (x, y)
を差分方程式で書くと u j1,l u j1,l u j,l 1 u j,l 1 4u j,l 2 j,l
• これを行列の形で書くと
u j 1,0
u j ,0 u j 1,0 2 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 ,l 1
j 1,l 1
j ,l 1 j 1,l 1
2
u
u
u
j ,L j 1,L
j ,L
j 1,L
これを以下のように記す u j1 T u j u j1 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 ) uJ / 2 uJ g(Jf/ )22
33
9th Physics Based Vision Seminar
2004/July/1
巡回還元(cyclic reduction, CR)
• その一つ前のレベルf-1では,以下の3つの式が出る
が,これも3重対角アルゴリズムで解ける
u0 T( f 1) uJ / 4 uJ / 2 g(Jf/ 21)2
uJ / 4 T( f 1) uJ / 2 u3J / 4 g(Jf/ 21)2
←解ける
uJ / 2 T( f 1) u3J / 4 uJ 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/