2010-05-02

■

Jos Stam の "Stable Fluids"という論文を読んだ。CG系CFDの世界では、有名な論文らしい。

主張としては、NS方程式の数値計算で、移流項の計算だけ、Euler的に計算するのでなくLagrange的に計算するとよいよ！というだけのことらしい。Euler移流だとタイムステップを小さくとらないと破綻するので、移流項だけLagrange的に計算するのをsemi-Lagrangian法とか呼ぶらしい。

Introductionに
Our method cannot be found in the computational fluids literature, since it is custom made for computer graphics applications.
って書いてあるんだけど、semi-Lagrange法は確かに工学系CFDではあまり見ない気がするものの、気象・海洋屋さんが使ってるのをたまに見かける。多分、気象だと、数日とか数ヶ月とか数年とかの単位で計算を追う必要があったりするので、タイムステップを大きくできるかは死活問題だったりしそうな気がする。

あとは大体標準的なMAC法なんかでよくやる通りの方法で、非圧縮性NS方程式は圧力に関する時間微分項がないので、圧力の時間発展を直接解くことができないため、代わりにNS方程式
$\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = -\nabla p + \frac{1}{Re} \Delta \mathbf{u} + \mathbf{f}$
の両辺のdivergenceを取ると、連続の式から、拡散項(粘性項)と速度の時間発展項が消えて、圧力に関するPoisson方程式
$\Delta p = \nabla \cdot (-\mathbf{u} \cdot \nabla \mathbf{u}+\mathbf{f})$
が出るので、現在の速度から圧力決定=>時間発展計算とかすればいい。論文はHelmholtz-Hodge分解とか回りくどい説明をしていて、謎(ではないけど)のスカラーポテンシャルを計算することになっている

MAC法だと数値誤差が蓄積して、連続の式が破れるのを防ぐために、圧力Poisson方程式に適切な補正項をいれるけど、論文ではそういう面倒なことは考えていない。そういうアバウトな計算方法について、Our method cannot be found in the computational fluids literatureと言ってるなら、それはそうかもしれないけど、教科書的な手法の毛を抜いたようなものでオリジナリティを主張されても微妙。

2010-05-02

Operadic Grobner bases

operadは、Peter Mayという人が、1970年代に多重ループ空間とやらを調べるために導入した概念であるらしい。しかし、わたしはトポロジーには興味がないので、それがどういうことなのかは知らない。その後20年以上operadが注目されることはなかったようだけども、Ginzburg/KasparovのKoszul duality for Operadsという論文を契機として、代数や数理物理で注目されるようになったらしい。operadには、Lawvere theoryの精密化という側面があって、色々な代数構造を特徴づけることができる。最近のoperadに対する注目は、このへんの話と関係しているらしい。Lavere theoryと違って、operadはsymmetric monoidal categoryで考えるので、その分だけ一般的。

(symmetric monoidal category C上の)operadの簡潔な定義は、Cに値を持つanalytic functorが合成を積として作るmonoidal categoryのmonoid objectというもの(正確には、Cがcolimitを持つ必要がある一方、operadはcolimitの存在を要求しないけど)。analytic functorはpolynomial functorと呼ばれることもあるけど、無限級数なのでanalyticが適切だと思う。数学系ではpolynomialと呼んで、CS系ではanalyticが多い気がする

operadの定義はぱっと見よく分からないけれども、考えていることは単純で、いくつかの演算(和とか積とか微分とかLie括弧とか)と、それらが満たす適当な関係式(結合律とかJacobi則とかLeibnitz則とか)があったとき、各次数で独立な項はどのくらいあるか?というのが基本的な問題。例えば、opを結合律を満たす普通の積とすると、opで作れる3次の項(op/置換/恒等写像の組み合わせで作れる三重"線形"写像$Hom(V^{\otimes 3},V))$の要素)は、
op(op(x,y),z),op(op(y,x),z),op(op(x,z),y),op(op(z,x),y),op(op(y,z),x),op(op(z,y),x),op(x,op(y,z)),op(x,op(z,y)),op(y,op(x,z)),op(y,op(z,x)),op(z,op(x,y)),op(z,op(y,x))
の12個ある。実際には、このうちいくつかは、結合律によって等しい。例えば、op(op(x,y),z)=op(x,op(y,z))とか。で、独立な項は6個であることはちょっと考えれば分かる。2次の項は、op(x,y)とop(y,x)で、4次の項は、xyzwの並べ替えの個数と等しいだけ独立な項がある。n次の項も同様にして、n!個の独立な項があり、各項は、対称群の要素と1:1に対応している。これが、Associativity operadとかいうもの。対応するanalytic functorは、ベクトル空間の圏では"テンソル代数"で、集合の圏で考えれば、リストを与える。

opが、結合律を満たす可換な積の時は、各次数で独立な項は1個だけで、これをCommutativity operadと呼ぶ。対応するanalytic functorは、"対称代数"(多項式環)あるいは"有限多重べき集合"。同様に、opが反対称で、Jacobi則
op(op(x,y),z)+op(op(z,x),y)+op(op(y,z),x)=0
を満たす時、各次数で独立な項を数えるのはしんどいけど、各次数に(n-1)!個の独立な項があって、Lie operadが定義できる。対応するanalytic functorは"自由Lie代数"を定める。足し算が密輸入のごとく紛れ込んでるけど、、operadを考えるときは、多くの場合k-linear symmetric monoidal category上で考える(kは代数的閉体の時が多い)ので、足し算/引き算は、ベクトル空間のものが自然に入っていると考えるっぽい(集合の圏とかは、k-linearでないので、Lie operadは作れない)。演算は、二項だけじゃなくて、0項、1項、3項、etc...で考えてもいいし、複数あってもいい

operadは、数学で重要な代数構造の多くを扱うことができる。現在知られているoperadの殆どの例は、Lodayが、Encyclopedia of types of algebrasとしてまとめている。最近では余代数構造も重要であるので、余代数構造も扱えるようなoperadの一般化として、PROP(あるいは、ほとんど同等の概念としてproperad)というのがある。余代数を扱うというmotivation以外に、"非線形な"関係式を扱いたいというのもある。例として、alternative algebraという代数は、(xx)y=x(xy),x(yy)=(xy)yを満たす代数構造であるが、xやyについて線形でなく、そのままではoperadの枠組みにのらない。

operadで扱えないもっと顕著な例としては、群がある。逆元に関する公理op(x,inv(x))=op(inv(x),x)=xは、xについて線形でない。Lawvere theoryで扱う時は、舞台がcartesian monoidal categoryだったので、自然な対角射diagが存在して、"point-free style"で
$op \circ (inv \otimes id) \circ diag = op \circ (id \otimes inv) \circ diag = e$
と逆元の公理を書けた。基本的に、Lawvere theoryでは、対角射を組み合わせて非線形な代数構造を扱うことができる。

一般のsymmetric monoidal categoryで考えるときは、diagの代わりに余積を取ると、上記の逆元の公理は、Hopf代数の定義の一部となる。一般に、余積が対角射となるという条件は、直接は代数的な関係式では書けないけど、集合の圏での実現を考えれば、群を得る(ベクトル空間の圏で考えれば、普通のHopf代数を得る)。というわけで、operadでは十分でないケースというのも、それなりにありそうではあるけれども、PROP/properadについては、operadほどよく調べられてはいないらしい。通常、余積は対角射でないので、Lawvere theoryに於ける対角射を余積に置き換えて、同じようなことをやっても、得られるモデルはLawvere theoryのそれとは違うはずだけど、どのくらい違うのかとかはよく分からない。

まあ、そんな一般論はよいのだけど、そもそもoperadとか調べようとした理由に、群とか環とかLie環みたいな、まともな代数構造じゃなくて、なんかもっと変態的な代数構造を調べてみたいというのがあった。変態的な代数構造として、とりあえず念頭にあったのは、Nambu bracketとか。あと、Haskellのmonadは、continuation monadは違うらしいけど、大体finitaryで、対応するLawvere theoryがあるとか、リスト=自由モノイドとか、binary treeは可換性のみを満たす二項演算から作れるとか。

#Lie括弧が交換子[X,Y]=XY-YXの抽象化だとすると、南部括弧は南部交換子[A,B,C] = ABC + BCA + CAB - BAC - ACB - CBAの抽象化である3項演算。Lie環と同じように、抽象的に3項演算を考えるのだけど、どうせ"普遍包絡環"にいけば環になるので結局は環論なんじゃね?と思う。この伝でいけば4項演算とか5項演算とか作って、Lie-4 algebraとかLie-5 algebraとかいくらでも作れるけど、だからなんだ。南部代数の正確な定義は
http://arxiv.org/abs/hep-th/9906248
とかに書いてある。

変態的な代数構造とか調べるのに個別に労力なんて割きたくないので、なんかもうちょっと体系的にアレコレ調べられるとよいとか思う。そんで、ある種のoperadについては、"Grobner basis"が計算できて、Haskellによる実装がある。
http://hackage.haskell.org/package/Operads
基本的に、ここで考えてるoperadは、適当な体上のベクトル空間で考えてるので、足し算は勝手に入ってる。monomialという概念は自然に分かるし、monomial同士の"割り算"も直感的には分かる。順序とかは自明でもないかもしれないけど、言われてみれば、どうってことない入れ方がある。あとは、大体普通のBuchbergerアルゴリズムと全く同じ手順を踏むだけ。

Commutativity operadのGrobner basisの計算は、例えば

Prelude Math.Operad> let m = corolla 1 [1,2]
Prelude Math.Operad> let r1 = oet $ nsCompose 1 m m :: FreeOperad Integer
Prelude Math.Operad> let r2 = oet $ nsCompose 2 m m :: FreeOperad Integer
Prelude Math.Operad> let r3 = oet $ shuffleCompose 1 [1,3,2] m m :: FreeOperad Integer
Prelude Math.Operad> [pp x |x<-operadicBuchberger[r1-r2 , r3-r2]]
["\n+1 % 1*m1(m1(1,3),2)\n+(-1) % 1*m1(1,m1(2,3))","\n+1 % 1*m1(m1(1,2),3)\n+(-1) % 1*m1(1,m1(2,3))"]

みたいな感じで、r1-r2とr3-r2で尽きて新しいものはでてこない。ppはpretty printerで、pp r1とかやると、r1がm(m(a1,a2),a3)に対応していることが分かる。項順序は、FreeOperad型の定義が
type FreeOperad a = OperadElement a Rational PathPerm
となっていて、PathPermが順序を決めているらしいけど、こいつは論文Grobner bases for operadsでは、Path-lexicographic orderingと書かれている順序に相当するものであるらしい。

論文Grobner bases for operadsでAntiCom operadと呼ばれているOperadについて、Grobner basisを計算すると、

Prelude Math.Operad> [pp x |x<-operadicBuchberger[r1+r2 , r3-r2]]
["\n+1 % 1*m1(m1(1,3),2)\n+(-1) % 1*m1(1,m1(2,3))","\n+1 % 1*m1(m1(1,2),3)\n+1 % 1*m1(1,m1(2,3))",
"\n+2 % 1*m1(1,m1(2,m1(3,4)))"]

とかなって、論文にも書いてある通り、新しいm(a1,m(a2,m(a3,a4)))を得る。

Lie operadとかを考えたい場合は若干面倒で、shuffleCompose 1 [3,1,2] m mとかやっても、こいつがShuffle Composeでないのと、もう一つはm(a2,a1)みたいなのが作れないので(a1,a2)->m(a2,a1)を新しい演算として定義してやって、うまく関係式を定義してやる必要がある。演算を増やすと、関係式の作り方は色々あるけど、項順序を固定しておけばどれでも同じグレブナー基底が得られるはず。例えば

Prelude Math.Operad> let u = corolla 1 [1,2]
Prelude Math.Operad> let v = corolla 2 [1,2]
Prelude Math.Operad> let r1 = oet $ shuffleCompose 1 [1,2,3] u u :: FreeOperad Integer
Prelude Math.Operad> let r2 = oet $ nsCompose 2 v u :: FreeOperad Integer
Prelude Math.Operad> let r3 = oet $ shuffleCompose 1 [1,3,2] v u :: FreeOperad Integer
Prelude Math.Operad> [pp x |x<-operadicBuchberger[r1+r2+r3, oet(u)+oet(v)]]
["\n+1 % 1*m2(1,2)\n+1 % 1*m1(1,2)","\n+1 % 1*m1(m1(1,2),3)\n+(-1) % 1*m1(m1(1,3),2)\n+(-1) % 1*m1(1,m1(2,3))"]

一応、こいつを使って、与えられた次数の項を全部数え上げるとかいうこともできる。

参考文献:
[1]Are operads algebraic theories?
http://arxiv.org/abs/math.CT/0404016

[2]Koszul duality for Operads
http://arxiv.org/abs/0709.1228

[3]A Poincare-Birkhoff-Witt criterion for Koszul operads
http://arxiv.org/abs/0709.2286

[4]Grobner bases for operads
http://arxiv.org/abs/0812.4069

[5]Implementing Grobner bases for operads
http://arxiv.org/abs/0909.4950

[6](Non-)Koszulity of operads for n-ary algebras, cohomology and deformations
http://arxiv.org/abs/0907.1505

[7]On the NonKoszulity of (2p+1)-ary partially associative Operads
http://arxiv.org/abs/0812.2687

[8]On the invertibility of quantization functors
http://arxiv.org/abs/math/0306212

[9]A Koszul duality for props
http://arxiv.org/abs/math/0411542

2010-05-02

拘束系の数値積分

拘束のあるHamilton系の数値積分について。拘束のある系の運動方程式を出すには、Lagrangeの未定乗数法が知られている。これは、Lagrangianによるものだけど、Legendre変換すれば当然Hamiltonianが出て、あとはHamilton系で考えることもできる。

未定乗数法のいやなところは、3点あって、第一にメカニズムが分からない。変数を機械的に増やして、なんでうまくいくのか。直感的には、拘束条件で縛れば、相空間の次元は、むしろ減るべきなんじゃないのかと思う。第二に、いちいちLagrangianを経由しないと拘束条件がいれられないのも、うっとうしい話。まあ、別にうまくいくからいいんじゃね？とか思わなくもないけど。第三に二点目と関係するけど、そもそもLagrangianが得られないケースはどうすればいいのか謎。正準方程式は"symplecticでないpoisson空間"上でも定義できるけど、symplecticでない場合に、Lagrangianを得る方法はないような気がする。symplecticでない正準方程式の例としては、Lotka-Volterra系とかEulerのコマとか。

例:(半径Rの)球面上に拘束された自由粒子。
未定乗数法の場合。
$L = \frac{m}{2}(\dot{x}^2+\dot{y}^2+\dot{z}^2) - u(x^2+y^2+z^2-R^2)$
というLagrangianを考えて、Euler-Lagrange方程式を出せばいい。uが未定乗数。今の場合は、初期条件とかから、うまくuを決定できる。

Dirac括弧を使う場合。配位空間(x,y,z)に対応する相空間の座標を(q1,q2,q3,p1,p2,p3)とする。運動を球面上に拘束するには、運動量の方にも何かしら拘束がかかっていないといけない。要するに、実は拘束条件が今のままだと不足しているので、全部出す。
$H = \frac{m}{2}(p_1^2+p_2^2+p_3^2)$
がハミルトニアンで、一個目の拘束条件は
$u_1=q_1^2+q_2^2+q_3^2 - R^2=0$
で、
$u_2=\{u_1,H\}=2m(p_1q_1+p_2q_2+p_3q_3)=0$
として、2個目の拘束条件が見つかる。

$\{u_1,u_2\}=4m(q_1^2+q_2^2+q_3^2)=4mR^2$
を使って、Dirac括弧の関係式として、
$[q_i,q_j] = 0 + \frac{1}{(4mR^2)}(\{q_i,u_1\}\{u_2,q_j\} - \{q_i,u_2\}\{u_1,q_j\})=0$

$[p_i,p_j] = 0 + \frac{1}{(4mR^2)}(\{p_i,u_1\}\{u_2,p_j\} - \{p_i,u_2\}\{u_1,p_j\})=\\0-\frac{1}{(4mR^2)}(-2*q_i*2m*p_j + 2*p_i*2m*q_j)=\frac{1}{R^2}(p_iq_j - p_jq_i)$

$[q_i,p_j] = \{q_i,p_j\} - \frac{1}{(4mR^2)}(\{q_i,u_1\}\{u_2,p_j\} - \{q_i,u_2\}\{u_1,p_j\})=\\\delta_{ij} + \frac{1}{(4mR^2)}(0 - 2m*q_i*2q_j)=\delta{ij} - \frac{1}{R^2}q_iq_j$
を得る。あとは、通常のように正準方程式を考えればいい。

手順は分かるけども、一体どうしてこれでよいのかやっぱり謎。Poisson代数Aを与えて、拘束条件から生成されるイデアルIを考えた時、一般にA/IはPoisson代数になるとは限らない。これはIが第二種拘束条件から生成される場合でもそうで、Aから素直にPoisson構造は入らないけど、うまく変形してやれば、Poisson構造が入るということで、このPoisson構造がどうして自然なのか?他にやり方はないのか?とか、そういう点については、よく分からない。