2020-12-13

Fermi推定の推定

Fermiが書いた文書My Observations During the Explosion at Trinity on July 16, 1945では、1945年7月に行われたトリニティ実験に於ける原爆のエネルギー推定を行っているが、どういう計算をしたのか詳細は書かれていない。以下に、数値情報が含まれている部分を抜粋する。

On the morning of the 16th of July, I was stationed at the Base Camp at Trinity in a position about ten miles from the site of the explosion.

After a few seconds the rising flames lost their brightness and ...

About 40 seconds after the explosion the air blast reached me. I tried to estimate its strength by dropping from about six feet small pieces of paper before, during and after the passage of the blast wave. Since at the time, there was no wind I could observe very distinctly and actually measure the displacement of the pieces of paper that were in the process of falling while the blast was passing. The shift was about 2 1/2 meters, which, at the time, I estimated to correspond to the blast that would be produced by ten thousand tons of T.N.T.

Fermiがどういう計算をしたか疑問に思う人はいるようで、How did Enrico Fermi calculate the classical Fermi Problem?で、回答している人もいるけど、(私が判断する限り)正しいかどうかは疑わしい。

そもそも、Fermiが計算しようとした量は、原理的には以下のようなものだと思われる。

$E = \displaystyle \int_{0}^{R(t)} \left( \dfrac{1}{2}\rho u^2 + \dfrac{p-p_{0}}{\gamma-1} \right) 2\pi r^2 dr$

$r$ は、爆発地点からの距離。 $R=R(t)$ は、時刻 $t$ に於いて爆風が到達した爆発地点からの距離。 $\rho(r,t) , u(r,t) ,p(r,t) , p_{0} , \gamma$ は、それぞれ、大気の密度、流速、圧力、爆発前の大気圧、比熱比。

爆発は、半球状に広がり、地表からの反射があるので、正確には球対称性はないと思われるが、ここでは方向依存性は無視した。 $E$ は時刻 $t$ ごとに計算されるが、エネルギー保存則によって一定と考えて差し支えないだろう(実際は、粘性や地表との摩擦によって減衰するだろうが)。爆発がなければ、被積分関数は恒等的に0(流速は0で、圧力は至る所大気圧と等しい)で、どこまで積分しても0になる。

上の式には、輻射に関する項が含まれてない。今日、物理のいくつかの分野では、輻射を伴う流体力学を使うことがあって、輻射流体力学(radiation hydrodynamics)と呼ばれることもある。輻射流体シミュレーションの基礎という解説記事には『現在，世界的に見ると，輻射流体コードと呼ばれるものは宇宙物理と慣性核融合のコミュニティを中心として多数存在する．』と書いてる。

輻射流体力学の研究が本格的に進展したのは戦後のことだと思うが、Betheらによる報告書BLAST WAVE(PDF)でも、輻射に関する考察は、かなり行われている。しかし、私の知る限り、輻射の影響を簡単に見積もる方法は確立されてない。報告書の中で、Betheは、

In the Trinity test it was measured that about 15 to 20 percen t of the total energy was emitted in the form of radiation which could penetrate the air to a distance of several miles.

と書いている。このことの妥当性はともかく、今の問題はFermiの推定方法で、輻射に関する見積もりをできたとは考えにくいので、輻射を無視することにする。

G.I. Taylorは、論文The formation of a blast wave by a very intense explosion. - II. The atomic explosion of 1945で、輻射を無視した爆発エネルギーの大きさを見積もり、それは上式で与えられるものだった(正確には、上式は半球上のみの伝播を考えているのに対して、G.I.Taylorは全方位に伝播する球対称解を考えているので、上式の2倍)。

Taylorは、今日、Sedov-Taylor自己相似解と呼ばれる理論解で仮定される相似則が(少なくとも最初の数十ミリ秒の間は)実際に成立していることを、写真の解析から示し、爆発エネルギーを、 $\rho(r,t) , u(r,t) ,p(r,t)$ の理論解に基づいて計算した結果は、上エネルギーを16800TNTトン相当だとした。

Sedov-Taylor自己相似解では、爆発エネルギー $E$ は、

$E \propto \dfrac{R^5 \rho_0}{t^2}$

となり、比例係数は、比熱比に依存するが、大体1くらいである(比例係数は無次元の量になる)。 $\rho_0$ は標準状態での大気の密度。

G.I. Taylorの論文TABLE Iのデータを使うと、、例えば、t=62(msec)の時、R=185.0(m)とあり、空気の密度を1.293(kg/m^3)とすると、

E ≒ (185**5)*1.293/(62/1000.0)**2/4.184e9 = 17421(tons of TNT)

などと計算でき、G.I.Taylorの見積もりに近い値が出る。

長崎や広島の原子爆弾は、高度600メートル付近で爆発させたようなので、トリニティ実験も、同程度の高さで爆発させたとすれば、圧力波が600メートルほど広がるまでは、球対称解は十分良い近似だったと考えられ、R=185.0(m)で、この理論を当てはめるのは問題ないように思う。

G.I.Taylorとvon Neumannは、1941年には、Sedov-Taylor自己相似解を独立に得ていたようだが、公開されたのは、いずれも終戦後らしい。Fermiは、マンハッタン計画に参加していたので、このような結果は知ってたと思われる。が、Fermiのいた地点で、この式を適用しても法外な値が出てくる。実際、t=40+few secondsで、R=16(km)を代入すると、E=1.3億TNTトンというような値になってしまう。それに、これだと、わざわざ紙を落とした意味がない。

Fermiが紙を落として計測しようとした量は、多分、衝撃波が通り過ぎた直後の風速じゃないかと思う。高さ6フィート(1.8m)の地点から、紙片を落として、落下地点が2.5mずれたとあるので、落下時の空気抵抗を無視すれば、落下に要する時間が0.6秒で、平均風速4(m/s)程度と計算される。正しい風速の計算は、紙片の重量や形状に依存して厄介だが、軽い物体なので、殆ど瞬時に風速まで加速され、一定速度に達したと見なしたのだろう。

爆風が通過した直後の風速は、Rankine–Hugoniotの理論で、よくparticle velocityと呼ばれてる量である。衝撃波面後方の圧力、密度、内部エネルギーを $p_1,\rho_1,e_1$ として、衝撃波の速度を $u_s$ 、particle velocityを $u_p$ とすると

$\rho_0 u_s = \rho_1 (u_s-u_p)$

$p_1-p_0 = \rho_0 u_s u_p$

$e_1-e_0 = \dfrac{1}{2}(p_1+p_0)(\dfrac{1}{\rho_0} - \dfrac{1}{\rho_1})$

が成り立つ。 $\rho_0,p_0,e_0$ は、標準状態での大気の密度、圧力、内部エネルギー。ちょっと頑張ると、以下のような式が得られる。

$\dfrac{\rho_1}{\rho_0} = \dfrac{(\gamma-1)p_0 + (\gamma+1)p_1}{(\gamma+1)p_0+(\gamma-1)p_1} = \dfrac{2\gamma p_0 + (\gamma+1)(p_1-p_0)}{2\gamma p_0+(\gamma-1)(p_1-p_0)}$

$u_p = \dfrac{c_0 p_1}{\gamma p_0} \dfrac{1}{\sqrt{1+\dfrac{\gamma+1}{2\gamma}\dfrac{p_1-p_0}{p_0}}}$

$c_0$ は音速。物理の教科書では、あんまり見ない気がする式だけど、 $(\Delta p_{S}) = p_1-p_0$ は最大静止過圧に相当する量になっている。で、最大動圧 $q = \dfrac{1}{2} \rho_1 u_p^2$ が

$q = \dfrac{ (\Delta p_{S})^2 }{2\gamma p_0 + (\gamma-1)(\Delta p_{S})}$

と書ける。Fermiが生身で観測して平気だったということは、 $(\Delta p_{S}) << p_0$ だったと思っていい。そうすると、 $\rho_1 \approx \rho_0$ で、 $u_p$ が4(m/s)だった時、 $q$ は約10(N/m^2)で、 $(\Delta p_{S})$ は1670(N/m^2)くらい。0.02気圧以下なので、この程度の爆風は、人体に損傷を与えない

Taylorのエネルギー式に戻ると、 $(\Delta p_{S})$ は $r$ が $R$ より少し小さい地点での $p-p_0$ に相当する量になっている。爆風が通り過ぎると、比較的素早く大気圧に戻ると予想されるので $r=R$ に近い領域で積分すれば十分で、また、Taylorのエネルギーは第二項の方が支配的になる(これについては、Taylorの論文The formation of a blast wave by a very intense explosion I. Theoretical discussionの式(31)と(32)を参照)。両者を勘案して、積分領域は狭いので、積分領域の $p(r,t)-p_0$ を一定値 $(\Delta p_{S})$ に置き換えてみると、以下のようになる。

$E \approx \dfrac{(\Delta p_{S})}{\gamma-1} \displaystyle \int_{R-\epsilon}^{R} 2\pi r^2 dr \approx \dfrac{2 \pi (\Delta p_{S}) R^2}{\gamma-1} \epsilon$

$\epsilon$ は適当なカットオフで、 $r \lt R-\epsilon$ では、圧力 $p$ は、十分、大気圧 $p_0$ に等しいと見なせるくらいの数値に取る。 $R$ は16(km)という大きな距離なので、 $\epsilon \lt\lt R$ としている。

$\epsilon$ は、どれくらいに取るべきだろうか。オーダー計算なので、1000(m)なのか100(m)なのか、10(m)なのかという話。そもそも、紙片が落下するのに0.6秒くらいかかってて、その間、風が吹いてたなら、その間に爆風が伝播した距離340(m/s)×0.6(s)= 204(m)くらいに取るべきという気もする(衝撃波の伝播速度 $u_s$ はほぼ音速に等しい)。

一方、実際の $p(r,t)-p_{0}$ は、 $r$ が $r=R$ から少し離れるだけで、急激に小さくなる(通常、過圧は、時間の指数関数オーダーで減衰するとされる)が、 $(\Delta p_{S})$ をそのまま掛けてるので、 $\epsilon$ は小さめの10(m)とかにするのが正しそうにも思える。 $\epsilon$ を10(m)とした場合、 $E$ は、およそ16000(tons of TNT)になる。

あるいは、以下のように考えてみることもできる。

$p(r) \sim p_0 + (\Delta p_{S}) \exp(-B(r-R)/\epsilon)$

という形を仮定してみる(衝撃波が通り過ぎた時、圧力は一旦上昇した後、大気圧以下まで下がってから戻るので、これは、あんまり正確ではなく積分を過大評価しそうであるが)と、主要な項は

$\dfrac{ 2 \pi (\Delta p_{S})}{\gamma-1} \displaystyle \int_{R-\epsilon}^{R} r^2 e^{-B(r-R)/\epsilon}dr \sim \dfrac{2 \pi (\Delta p_{S}) R^2}{\gamma-1} \dfrac{\epsilon}{B}$

と書ける。 $B$ は減衰の早さを表すが、 $r\lt R-\epsilon$ では、 $p(r)$ は $p_0$ に十分近くなっているとしていた。 $(p(R-\epsilon)-p_0)/(p(R)-p_0) = e^{-B}$ であるから、 $B=1$ では、不十分だろう。 $B=5$ なら、比率は1%以下になる。

所詮、オーダー計算なので、あまり深く悩んでも仕方ない。 $\epsilon$ が100(m)で、 $B$ を10にすれば、 $E$ は、およそ16000(tons of TNT)になる。

G.I.Taylorの見積もりを知ってるから、それらしい数字を選べるけれど、正直なところ、一点に於ける最大静止過圧が分かっても、積分を正しく見積もるのは、結構難しいように思える。爆発は、半球状に広がり、方向依存性がないとしたのも、妥当な近似かどうかは疑わしい。

マンハッタン計画では、当然、原子爆弾の威力についても、爆風の伝播についても、事前に、かなりの検討が行われたはずだから、Fermiは、それらの知識を動員した上で、見積もりを行ったのかもしれない。何らかの実験データに基づいた数値をFermiが使ったということがあれば、見積もりの詳細を書かなかったことも理解できる。

Betheらによる報告書BLAST WAVEのChapter 5は、爆発地点から遠く離れた場所での爆風の伝播を考察していて、式(5.144)及び(5.144a)に、数値計算(IBM calculationなどと書かれてるが、マンハッタン計画では、電子計算機以前のIBMの"計算機"が使用されてたらしいので、それのことだろう)の結果に基づいて

$\dfrac{(\Delta p_{S})}{p_{0}} = \dfrac{260.7}{R\sqrt{\ln(R/260.1)}}$

$\dfrac{(\Delta p_{S})}{p_{0}} = \dfrac{278.6}{R\sqrt{\ln(R/129.1)}}$

などの(経験)式が書かれている。Rの単位はメートルで、10000~13500TNTトン相当のエネルギーが爆発によって解放されたと仮定した計算(エネルギーは保存せず増加していったらしい)だったようだ。この数値計算も、圧力波伝播の方向依存性はないものと仮定した計算になっている。当時の計算機の性能では、空間次元を減らすことは、とても重要だったのだろう。

Fermiのいた地点が、R=16000(m)であるから、どちらの式を使っても、最大静止過圧/大気圧比 $\frac{(\Delta p_{S})}{p_{0}}$ は、0.008程度。一方、既に書いたように、最大静止過圧は、風速から見積もることができ、 $u_p$ が4(m/s)とした時、 $(\Delta p_{S})$ は1670(N/m^2)くらいで、 $(\Delta p_{S})/p_{0}$ は、0.0167程度。数値計算とは2倍くらいずれているが、爆発のエネルギーが数値計算と同程度のオーダーだったと考えても差し支えないだろう。

Fermiは、単に、この結果を知っていて、10000TNTトン相当と結論付けたのかもしれない。

2020-06-14

Tales of Tetranitrogenia

高校化学(≒量子力学以前の化学)の知識だと、窒素4つからなる分子は、以下のように、2つの構造異性体が考えられる。

周期表で、窒素の下にあるリン元素の場合は、リン2分子からなる二リンと、正四面体構造を取る白リン(黃リン)があるらしい。リンの場合は、二リンより正四面体構造を取る方が安定で、二リンは高い反応性を持つとWikipediaには書いてある。理屈上は、立方体構造を取るN8とP8も可能と思われるけど、どっちも合成されてないっぽい(N8の方は、オクタアザキュバンと呼ばれる)。正多面体構造だと、正十二面体も可能そうだけど、合成されてるのか調べてない

窒素では、N4が観測されたとか合成されたという話を聞かないし、きっと不安定で、すぐに2つの窒素分子に解離するに違いないのだろうと思ったけど、実際どうなのか。ということが、ふと気になって、検索したら

Thermochemistry of Tetrazete and Tetraazatetrahedrane: A High-Level Computational Study
https://doi.org/10.1021/jp952026w

という論文が1996年に出ていた。論文によると、上の2つの構造以外に、N4分子としては、open-chain structureもあるよと書いてあって、正四面体構造を取る方が、TetraazaTetrahedrane、もう一個がtetrazeteと呼ばれているので、以下では、そう呼称することにする

TABLE Iの計算結果によると、窒素原子4つからなる構造の中では、窒素分子2つが、エネルギー的に一番安定らしい。Abstractにも、At the G2 level, the enthalpies of formation,∆Hf298, of tetraazatetrahedrane, tetrazete, and the open-chain tetranitrogen are 732.5±8.0, 746.5±7.6, and 686.6±7.6 (kJ /mol), respectivelyなどと書いてある。このエネルギー差は、窒素一重結合、窒素二重結合が、三重結合に比べると、全然弱いからだと、彼らの説明には書かれている。

AbstractにあるG2 levelというのは、
Quantum chemistry composite methods
https://en.wikipedia.org/wiki/Quantum_chemistry_composite_methods
にあるGaussian-2 methodのこと。まぁ、なんか胡散臭い方法ではあるけど、G2の論文が1991年で、G3の論文が1998年なので、当時は新しい手法だったっぽい。

Gamess(2019/09/30のソースからビルドしたもの)のexam43.inpに、G3MP2があるので、参考にして、G3MP2で計算しておく。Tetrazeteの方は、以下のinputで計算した

 $contrl scftyp=rhf coord=zmt runtyp=g3mp2 $end  
 $system timlim=2 mwords=4 $end
 $scf    dirscf=.true. $end
 $data
Tetrazete...G3(MP2,CCSD(T))
C1
N
N 1 r1
N 2 r2 1 a1
N 3 r1 2 a1 1 0.0

r1=1.287
r2=1.542
a1=90.0
 $end

一応、出力の結合を確認しておくと

          -------------------------------
          BOND ORDER AND VALENCE ANALYSIS     BOND ORDER THRESHOLD=0.050
          -------------------------------

                   BOND                       BOND                       BOND
  ATOM PAIR DIST  ORDER      ATOM PAIR DIST  ORDER      ATOM PAIR DIST  ORDER
    1   2  1.287  2.368        1   4  1.542  1.128        2   3  1.542  1.128
    3   4  1.287  2.368

とか出ているので、二重結合が、1.287Å、一重結合が1.542Åで、まぁ、良さそう

    ----------------------------------------------------------------
                   SUMMARY OF G3(MP2) CALCULATIONS
    ----------------------------------------------------------------
    MP2/6-31G(D)    =  -218.202935   CCSD(T)/6-31G(D) =  -218.248194
    MP2/G3MP2LARGE  =  -218.400052   BASIS CONTRIBUT  =    -0.197116
    ZPE(HF/6-31G(D))=     0.015581   ZPE SCALE FACTOR =     0.892900
    HLC             =    -0.091700   FREE ENERGY      =    -0.007307
    THERMAL ENERGY  =     0.020515   THERMAL ENTHALPY =     0.021459
    ----------------------------------------------------------------
    E(G3(MP2)) @ 0K =  -218.521429   E(G3(MP2)) @298K =  -218.518365
    H(G3(MP2))      =  -218.517420   G(G3(MP2))       =  -218.546187
    ----------------------------------------------------------------

    ----------------------------------------------------------------
          HEAT OF FORMATION   (0K):     185.12 KCAL/MOL
          HEAT OF FORMATION (298K):     183.47 KCAL/MOL
    ----------------------------------------------------------------
          HEATS OF FORMATIONS BASED ON NIST DATABASE FROM 
          COMPUTATIONAL CHEMISTRY COMPARISON AND BENCHMARK DATABASE
    ----------------------------------------------------------------

298Kに於ける生成エンタルピーが、768.15kJ/mol相当

tetraazatetrahedraneは

 $contrl scftyp=rhf runtyp=g3mp2 $end  
 $system timlim=2 mwords=2 $end
 $scf    dirscf=.true. $end
 $data
Tetraazatetrahedrane...G3(MP2,CCSD(T))
Td

X
N  7.0  0.85332     0.85332     0.85332
 $end

でやると、生成エンタルピーを計算してくれなかった。ダミー原子のせい？

    ----------------------------------------------------------------
                   SUMMARY OF G3(MP2) CALCULATIONS
    ----------------------------------------------------------------
    MP2/6-31G(D)    =  -218.222646   CCSD(T)/6-31G(D) =  -218.245052
    MP2/G3MP2LARGE  =  -218.421910   BASIS CONTRIBUT  =    -0.199265
    ZPE(HF/6-31G(D))=     0.014131   ZPE SCALE FACTOR =     0.892900
    HLC             =    -0.091700   FREE ENERGY      =    -0.006329
    THERMAL ENERGY  =     0.018831   THERMAL ENTHALPY =     0.019776
    ----------------------------------------------------------------
    E(G3(MP2)) @ 0K =  -218.521885   E(G3(MP2)) @298K =  -218.518880
    H(G3(MP2))      =  -218.517936   G(G3(MP2))       =  -218.544040
    ----------------------------------------------------------------

 NO INTERNALLY STORED ATOMIC THERMOCHEMICAL DATA FOR ATOM Z=  1

面倒なので、座標書いて計算してやると、298Kに於ける生成エンタルピーは、766.26(kJ/mol)らしい。

論文の値Tetraazatetrahedrane732.5±8.0(kJ/mol), tetrazete746.5±7.6(kJ/mol)と、オーダーは合ってる

生成エンタルピー的には、窒素分子2つというのが安定なので、窒素であふれる地球大気中に、N4分子が全然見当たらなくても不思議ではないのだけど、人工的に合成できたとして、どれくらいの時間で解離するのかという問題は、活性化エネルギーが分からないといけない。

一般的に信じられているところでは、遷移状態=ポテンシャルエネルギー曲面上の鞍点ということになっているけど、窒素4原子でも、ポテンシャルエネルギー曲面は、6次元なので、全探索するのは、非効率的だし、かといって、遷移状態の構造最適化をやるにも、ある程度、どのへんに鞍点があるかアタリを付けてやらないと、収束してくれなかったり、関係ないところに収束するかもしれない。

遷移状態探索の常套手段を知らないのだけど、とりあえず、Tetraazatetrahedraneの解離を調べるとして、4つの窒素原子の座標を(-A/2,B/2,0),(-A/2,-B/2,0),(A/2,0,B/2),(A/2,0,-B/2)として、2つのパラメータA,Bを動して、一点計算をやれば、おおまかなポテンシャルエネルギー曲面の形状が掴めると期待する。このパラメータは、AとBをうまくとれば、正四面体にできるし、一方、Aを非常に大きく取れば、N2が遠く離れた解離状態も再現できると思われる。

x座標が等しい窒素原子同士の距離はBであるけど、これを大きくすることは、あまり興味がない。これらの原子間の結合は、Tetraazatetrahedraneでは、一重結合で、冒頭の論文によれば、原子間距離は1.478Åとされていて、Aが大きくなるにつれて、三重結合に近付くと思われる。窒素分子N2では、結合距離が1.098Åくらいとされてるので、Bの範囲としては、1.098Å〜1.478Åを見れば十分だろう。少し広めに、1.0〜1.5Åとする

そんなことを考えて、一点計算を繰り返して得られたポテンシャルエネルギー曲面が以下(計算は、RHF/6311G**で行った)。

f:id:m-a-o:20200531182915p:plain

なんか断崖絶壁みたいになっているけど、Tetraazatetrahedraneは、A=1.0451,B=1.478のあたりで、高台の窪地みたいな位置にある。これを元に、遷移状態を計算しようとしたけど、かなり初期値をうまく選ばないと、明らかに関係ないとこに収束してしまう。

 $CONTRL SCFTYP=RHF RUNTYP=SADPOINT UNITS=ANGS $END
 $SYSTEM TIMLIM=10 $END
 $BASIS GBASIS=N31 NGAUSS=6 $END
 $STATPT OPTTOL=1.0E-6 NSTEP=500 HESS=CALC HSSEND=.T. $END
 $DATA
N4 ... TS search
C1
N  7 -0.75170436 0.63485445  0.000
N  7 -0.75170436  -0.63485445  0.000
N  7 0.75170436  0.0000  0.63485445
N  7 0.75170436  0.0000  -0.63485445
 $END

で、とりあえず、鞍点には収束した。

TOTAL ENERGY = -217.2507936836

だそうだ。

 COORDINATES OF ALL ATOMS ARE (ANGS)
   ATOM   CHARGE       X              Y              Z
 ------------------------------------------------------------
 N           7.0  -0.5134742719   0.7564702157   0.1903549687
 N           7.0  -0.4320829219  -0.7401136741  -0.1328901129
 N           7.0   0.3912730128   0.1115247255   0.8406026425
 N           7.0   0.5542841810  -0.1278812671  -0.8980674975

     FREQUENCIES IN CM**-1, IR INTENSITIES IN DEBYE**2/AMU-ANGSTROM**2,
     REDUCED MASSES IN AMU.

                            1           2           3           4           5
       FREQUENCY:       667.84 I      1.51        1.17        0.49        0.14

で、667.84(/cm)に虚の振動数が一個出てる。こうして得た座標を初期値として、6311G**基底とかで再度、鞍点探索すると、鞍点に到達できる

6311G**で計算した場合、TOTAL ENERGY = -217.5073394100となった。

そんで、一般的に、ポテンシャルエネルギー曲面に鞍点は複数存在するので、得られた鞍点が求める反応の遷移状態か確認するためにIRC計算というものをやるのが標準的な手続きらしい。出力してくれるのは、途中及び最終的な位置座標と、エネルギーだけど、FORWARD=.T.の方は、窒素原子間距離を計算すると、全てのペアで距離が約1.3918Åになっているので、Tetraazatetrahedraneにたどり着いたっぽい雰囲気。論文と結合距離が結構違うけど、基底関数が違うせいだろう。

FORWARD=.F.の方は、どこに向かってるのか分からないところで終了してしまった。4つの原子が平面上にあって、2つ目の窒素原子が、残りの3原子の重心近くにいるような配置。何となく、2つのN2に解離していきそうな雰囲気もあるけど謎。。。

よく分からんけど、これが求める遷移状態だったとして、始状態で一応構造最適化してエネルギーは、-217.5869661854(Hartree)だと言ってるので、Tetraazatetrahedraneと鞍点でのエネルギーとの差を単純にkJ/molに変換すると、210kJ/mol。予想に反してでかい。室温で、すぐ解離する程度のエネルギーだろうと思ってたけど、意外と、室温で安定するのかもなぁという感じになった。

tetrazeteもやろうと思ってたけど、鞍点に収束させるのに、手間取ったので終了

という計算をした後、改めて探してみると、2000年の以下の論文では、窒素気体中で電気放電を行い超低温(6.2~35K)に急冷すると、tetraazatetrahedraneが合成できた(かも)と主張してる。論文では、tetrahedral tetrazetesと呼ばれている。

Tetrazete (N4). Can it be prepared and observed?
https://www.sciencedirect.com/science/article/abs/pii/S000926140000926X

根拠がIRスペクトルしかないっぽいので若干微妙だけど、合成の方法は、フラーレンとかを最初に合成したのと同じような方法っぽい(?)

2020-06-14

コンピュータ以前の数値計算(2)20世紀初頭の計算需要

コンピュータ以前の数値計算(1) 三角関数表小史 -
の続き。

18世紀や19世紀のヨーロッパでは、三角関数表のように、事前に数表を作っておいて、必要なら補間して参照するとかいうのが数値計算の典型的な運用だったようである。1938年になっても、アメリカでは、(雇用創出の一環として)Mathematical Tables Projectという大規模な数表作成プロジェクトがあったらしい。計算尺みたいな計算補助用のアナログツールも、表という形ではないけど、事前に行った特殊関数の計算に帰着させるという点では、同種の方法と言える

一方で、大規模な連立一次方程式を解くとか、微分方程式の解法なんかは、一般には、(少数変数の)特殊関数の計算に帰着できない。微分方程式の場合、有理関数や代数関数の積分、あるいは有名な特殊関数を使って解が書ける状況では、数表を作って置くという方法は有効だっただろうし、それで、初期には解が具体的に書ける問題が追求されたという面もあったのかもしれない。

ヨーロッパで、物理現象を微分方程式で記述するというスタイルが普及してきたのは18世紀頃だけど、古い所では、Eulerが常微分方程式の数値積分を行っている(今でもEuler法は最初に習う)し、解を適当な特殊関数で書いたり近似できるとは限らないという疑念も初期からあったのだろう。

電子計算機ができた直後には、数表を作るために電子計算機が使われたりもしたようだけど、電子計算機開発直前の様子を見ると、新しい計算機の開発者たちは、単純な特殊関数の計算に帰着できない問題を解くために、試行錯誤を行っていたように見える。

20世紀初頭の応用数学書

20世紀の初め頃、実用的な問題への利用を視野に入れつつ応用数学を解説した本が、数学者によって書かれるようになった。最初に、それらの本を眺めて、当時の雰囲気を概観する。

元々、19世紀頃は、まだ分野の専門分化がそれほど進んでなくて、現在は数学者に分類される人でも、物理や工学を研究するのは普通だったと言っていい。当時は、理論物理の中心は、力学だった(光学はあったし、電磁気学や熱力学は研究され始めてたけど、力学で理解できるという期待もあったと思われる)ろうし、力学の理論は数学の一部と見なす人も、いたっぽい。また、物理学者とされる人も、工業的な問題に根差した題材を扱うことは、しばしばあって、数学、物理、工学の垣根は低かったようである。

そんな時代ではあっても、応用数学とか実用数学とか工業数学みたいな呼称は存在して、使われてたっぽい。

イギリスのFrancis Bashforth(1819-1912)は、弾道学に関する1868年の論文

On the resistance of the air to the motion of elongated projectiles having variously formed heads
https://royalsocietypublishing.org/doi/abs/10.1098/rstl.1868.0015

で、Professor of Applied Mathematicsという肩書や用語を使っている。

In the spring of 1864, when I was appointed Professor of Applied Mathematics to the Advanced Class of Artillery Officers at Woolwich,...

尤も、どういう分野を指して"応用数学"と呼んでるのかは定かでない。また、所属が王立陸軍士官学校なので、当時、どれくらい研究機関としての機能があったのかは分からない。更に、現代の基準からすれば、Bashforthの仕事は、数学者らしいものとは言えないかもしれない。

常微分方程式の数値解法であるAdams-Bashforth法という名前は、この人にちなむ。Adamsの方は、John Couch Adams(1819-1892)という人で、海王星の存在と位置を理論的に予測したことで知られる(計算結果は、Airyらに渡されたが、すぐに探索を開始しなかったので、海王星の発見はドイツでなされたようである)。Bashforthは弾道学の人だけど、Adams-Bashforth法は、以下の文献で使われた手法に由来するっぽい

An attempt to test the theories of capillary action by comparing the theoretical and measured forms of drops of fluid.
With an explanation of the method of integration employed in constucting the tables which give the theoretical forms of such drops
https://archive.org/details/attempttest00bashrich/page/n7/mode/2up

Adamsファミリーには、Adams-Moulton法というのもあり、Forest Ray Moultonは、天文学者として知られている人だけど、1926年に、New methods in exterior ballisticsという本を出版している。

それはともかく、少なくとも、AdamsやBashforthは数値計算の専門家ではないし、19世紀には、数値計算の専門家みたいな人は多分いなかったと思われる(計算手という仕事はあったっぽいけど)。自分自身で計算するか他人に計算させるかはともかく、数値計算自体は、数学者や科学者だろうと技術者だろうと、文字通り誰でも行うものだったので、分野として独立させる意味がなかったのだと思われる

1903年のScienceには、
ON THE FOUNDATIONS OF MATHEMATICS
http://doi.org/10.1126/science.17.428.401
という記事が載っていて、pure and applied mathematicsについて、度々言及されている。"This separation between pure mathematics and applied mathematics is grievous even in the domain of elementary mathematics."とか"The fundamental problem is that of the urtification of pure and applied mathematics. "などと書いてるので、文脈をすっとばしても、純粋数学と応用数学に大きな乖離があると、当時の人々が認識していたらしきことは伝わる。イギリスでは、工学者のジョン・ペリーが、数学教育の改革を訴え、一方、数学研究としては、Peanoによる算術の公理化や、Hilbertの"幾何学基礎論"が出たような時期で、そういったことが影響していたようである。

ドイツでは、Felix Klein(1849-1925)の尽力もあって、1904年に、ゲッティンゲン大学で、いくつかの研究所が新設された時、応用数学の責任者として、Carl Rungeが招聘された。この時、Prandtlも、応用力学部門の責任者として招かれた。Prandtlの境界層の理論は、1904年にハイデルベルグで開かれたICM/国際数学者会議で発表されている。

応用数学が何を指すかは曖昧であるものの、現代だと、テーマを近似解法に限定しないのが一般的だと思う(例えば、グラフ理論とかは、応用数学に分類されることもあると思う)けど、この時の"応用数学"は、"数学的問題の数値的及び図式的な手法による近似解法の研究を行う分野"と位置づけられていたよう。上のScienceの記事には

As a basis of union of the pure mathematicians and the applied mathematicians, Klein has throughout emphasized the importance of a clear understanding of the relations between those two parts of mathematics which are conveniently called 'mathematics of precision' and 'mathematics of approximation,' and ...

という記述があって、応用数学=近似の数学という認識が、Kleinにあったらしきことが窺える。

数値計算を専門的に扱う部門が、ドイツで登場したのは、多分これが初めて。他の欧米諸国でも、数値計算を専門的に扱う部門というのは存在してなかったんじゃないかと思うけど、確証はない

Rungeは、Weierstrassの弟子で、Runge-Kutta法のプロトタイプ版を、1895年に論文で発表していたし、近似解法を専門的に研究した最初期の人だと思う。

【補足】Runge-Kutta法(時々、Runge-Kutta-Heunとする方が適切と書かれている)の歴史の日本語で読める簡単な説明は、例えば、以下の文献の前半にある。
Runge-Kutta 法-その過去, 現在, 未来-
https://doi.org/10.11429/emath1996.1998.Spring-Meeting_93

かつては、数値計算だけでなく、図式解法も、数値計算より精度は低いものの高速で有用とされていたらしい。Rungeは、1909〜10年に、コロンビア大学で講義を行い、講義録がGraphical methodsという本になっている。
Graphical methods
https://archive.org/details/graphicalmethods00rung/page/n5/mode/2up

そのRungeの下で学位を取得した人に、Horst von Sandenという人がいて、1913年にドイツ語で、数値計算と図式計算に関する本を書いていて、英訳が以下で公開されてる
Practical mathematical analysis. With examples by the translator, H. Levy
https://archive.org/details/practicalmathema00sanduoft

日本では、1928年に翻訳らしきものが出てる(オープンアクセスにはなってない)。国会図書館で適当に検索する限り、これ以前に、類書は見当たらない
実用解析学 : 数値計算、図計算、機械計算ノ概念
https://iss.ndl.go.jp/books/R100000002-I000000777141-00

Sandenの本には、(総ページ数が200ページにも満たないのに)約30ページに渡って、計算尺や"計算機"の使い方が説明されている。現在では実用性のない記述だけど、当時の状況を知るには役立つ。

計算機の例としては、Burkhardt's calculating machineを選ぶと書いてる。現在、Burkhardt arithmometersと呼ばれてるものと同じだろうと思うけど、当時は、他にも、数字を入力してハンドルを回すと、足し算や掛け算ができる(つまり、人力を動力とする)機械が沢山あったそうだ

最も基本的な計算尺は、対数目盛りのついた定規を組み合わせたような道具で、掛け算や割り算の近似計算に広く使われたアナログツールだったらしい。計算尺で値を読む時は、目盛りを使うので、精度には限界がある。現在の定規は最小目盛りが通常1mmで、目盛り自体の幅が、0.2~0.3mmくらいらしい。視力的には、5m先から、1.45mmの隙間が(片目で)見えることが視力1.0の定義なので、例えば、30cm地点からなら、0.1mmを区別できるかもしれない。

仮に、0.1mmごとに目盛りが打ってあっても、長さ25cmの計算尺で、高々2500個以下の目盛りしか打てない。ビット数にして、11bitか、あるいは、対数目盛りであることを考慮すれば、それ未満。手に持つ都合上、計算尺の長さは数十cmのものが普通で、どうしても精度が必要な場合には、長さ1mを超える長大な計算尺が使うこともあったらしい。現代で言うと、半精度浮動小数点数の仮数部の幅が11bitなので、よくて、その程度の精度だったということになる。基本的に、単純な乗算や除算の補助というのが主目的っぽいので、その程度でも問題なかったのだろう

計算速度に関しては、1850年代に存在した初期のArithmometerで、8桁の数同士の乗算に18秒という記録があるようなので、有効桁数2〜3桁で十分という状況なら、計算尺の方が速かっただろう。

イギリスでは、E.T.Whittakerが、1913年から、数値計算と応用数学に関する講義を開始した。1913〜23年までの間に行った講義がまとめられて、1924年に以下の本が出版された(と、前書きに書いてある)。

The Calculus of Observations. A Treatise on Numerical Mathematics
https://archive.org/details/calculusofobserv031400mbp/page/n6

アメリカでは、James B. Scaboroughという人が、1930年に、数値計算に関する書籍の初版を出版している

Numerical Mathematical Analysis
https://archive.org/details/in.ernet.dli.2015.149743/page/n1

表紙によれば、Scaboroughは、アメリカ海軍兵学校(U.S. Naval academy)の数学の名誉教授(?)(professor emeritus)だと書いてある(他の人の記述では、associate professorとなっている)。

WhittakerやScaboroughの本は、Sandenの本と違い、記述が数学的内容に限定されていて、計算尺や計算機の使い方の説明もなければ、図式計算に関する記述も少ない。けど、Whittakerの本のChapter VIの冒頭を見ると、方程式の解法として、
(α)literal methods
(β)Numerical or computer's methods
(γ)Graphical methods
(δ)Nomographic methods
(ε)Mechanical methods
が挙げられている。

(γ)のノモグラフというのは、特定の関数(典型的には、二変数関数)の計算に特化したアナログ計算ツールで、まぁ、そういうのがあったらしい。画像検索とかすると、目盛りが3つあって、2つの入力値を結ぶ直線を引いて、第三の目盛りとの交点の数値を読み取るというタイプのものが多いっぽい。計算尺同様、精度は高くないだろうが、高速で、また高等数学教育を受けてない人でも扱えるということで重宝されていたっぽい。現代だと、近似解法は、数値計算とほぼ同義だけど、当時は、まだそうなってなかったらしいことが窺える。

他に、(Krylov部分空間法で名前が残っている)ロシア海軍の造船技術者Alexei Krylov(1863-1945)は、1906年に"О приближённых вычислениях"(On approximate caluculations)という講義録を出版している。英語訳や日本語訳はないようなので、内容は確認できてないけど、数学や天文学で知られていた近似計算法を扱っているらしい。

Sanden、Whittaker、Scaboroughの本を全部読むのは長すぎるから、数値線形代数に関する記述に限定して見ると、Whittakerの本では、Chapter Vに"Determinants and Linear Equations"という章があるけど、5ページくらいしかなく、Cramerの公式による方法(行列式は、普通に消去法で計算する想定らしい)が説明されている。また、Sandenの本のChapter IX §1が、連立線形方程式の解法に関する記述となっている。ページ数は同じく5ページくらい。

連立一次方程式の消去法による解法は、紀元前後に成立した九章算術に記述があるらしいけど、現在、Gaussの名前を冠して呼ばれることがある(単に"消去法"と呼ぶと、複数の語義がありうるので、以下では、Gaussの消去法と呼ぶことにする)。Gauss以前のヨーロッパでは、1750年にCramerがCramerの公式を発表し、1772年にLaplaceが行列式の余因子展開を与えているけど、行列式をどういう方法で計算するにしても、Cramerの公式による連立一次方程式の解法は、アルゴリズムとしては実用的でない

Scaboroughの本は、最終章が"Numerical Solution Of Simultaneous Linear Equations"となっていて、数十ページほどある。内容は、solution by determinants/solution by successive elimination of the unknowns/solution by inversion of matrices/solution by iterationとなっている。Gaussの消去法を含み、適切なピボット選択が必須である旨も注意されていて、当時の記述としては、実用性のあるものだった。

【補足】Gaussの消去法をLDU分解として再定式化したのは、Alan Turingの1948年の論文のようである。
ROUNDING-OFF ERRORS IN MATRIX PROCESSES
https://academic.oup.com/qjmam/article/1/1/287/1883483
この前年には、von NeumannとGoldstineが、"Numerical inverting of matrices of high order"という論文を書いているので、比較すると面白い

一番新しいScaboroughの本は、連立一次方程式の数値解法に割かれるページ数が増大しているけど、どの本も、連立一次方程式の解法のみがあって、固有値問題の記述は存在しない。eigenvalue/eigenvectorという用語は、1904年、Hilbertが積分作用素の研究で導入したeigenwert/eigenvektorに由来するけど、固有値の計算自体は、もっと古くから見られる。

二次形式や定数係数線形常方程式系の問題は、固有値問題に帰着でき、18〜19世紀には現れている。物理では、二次形式は、ポテンシャルエネルギーを極小点近傍で二次の項まで展開した時に現れる(つまり、二次形式を与える正定値実対称行列は、極小点に於けるHessianの値)。物質科学の人は、調和近似と呼んでることもあるけど、簡単な割に、古典的な弾性論から、量子力学や量子化学でも、有用なので頻繁に使われる。このようなポテンシャルエネルギーで記述される古典力学系は、連成振動系なので、定数係数線形常方程式系となる。

また、18世紀に、慣性主軸の決定は、Eulerによって取り上げられた。実対称行列の固有値が実数であることを証明したのは、Cauchyらしく、彼の動機の一つは、二次曲面の理論にあったらしい。定数係数線形常微分方程式系に対する関心は、1743年のD'Alambertの本"Traité de dynamique"で、連成振動の考察を行ったことに始まる。また、天体力学に於ける永年方程式からも、定数係数線形常微分方程式系の解析が必要になった。
cf)Lagrange and the secular equation
https://link.springer.com/article/10.1007/s40329-014-0051-3

19世紀でも、特性方程式の根の積が行列式に等しいというようなことは理解されていたけど、現在、英語圏でeigenvalueという名前が普及してるのを見ると、それ以前は広く使われる名前もなく、取り立てて、重要視もされてなかったのかもしれない。英語圏で、eigenvalueという単語の初出は、1927年の

Eigenvalues and Whittaker's Function
https://www.nature.com/articles/120117a0

らしい。Abstractを読むと、水素原子のシュレディンガー方程式の解法自体は、Whittaker-Watsonの教科書に載ってるみたいなことが書いてるので、量子力学が動機といえる。

数値計算の本ではないけど、Courant-Hilbertの有名な本"Methoden der mathematischen Physik"の第1巻は、1924年に出版されていて、この本では、固有値・固有ベクトルが、至る所に出てくる。

Methoden der mathematischen Physik (online reproduction of 1924 German edition)
https://gdz.sub.uni-goettingen.de/id/PPN38067226X

Courantは、ゲッティンゲン時代(1933年に、ドイツを脱出している)に、線形代数の体系的な講義を初めて行ったらしい。Courant-Hilbertの本は、量子力学と関係ないけど、固有値問題が市民権を得る、もうひとつの契機になったかもしれない。

線形代数の起原は、他にも、1844年のGrassmannの本"Die lineale Ausdehnungslehre ein neuer Zweig der Mathematik"にあると書いてある場合もある。これは、Peanoが、ベクトル空間の公理化などを行った時の動機がGrassmannの本にあったためらしい。Grassmann自身の動機は、(現代の言葉では)外積代数とEuclid空間上の微分形式の理論の構築だったっぽい。この方向での萌芽は、Jacobi行列式とかにも見られる。

超大雑把な話として、外積代数を支配するのは一般線形群だし、量子力学を支配するのはユニタリー群で、二次形式を支配するのは直交群なので、このへんの違いから、線形代数の雑多な起原が生じているのだと思われる。

1930年代以降、多くの人が、線形代数という枠組みを認識し始め、線形代数が、理工学系学生の必修になるのは、世界的に1960年代か、それ以降のようだけど、あまり、ちゃんとした資料はない。

線形代数のことはともかく、数値計算への関心と教育の変化が、20世紀初頭に出現していたということが見て取れる。

おまけ:統計学の成立

上に書いたような数値計算は、物理や工業分野の問題を想定している印象がある。現在では、応用数学の一分野と見做されることもある統計学は、かなり独立した起原を持っている。

statisticsという用語は、ドイツ語のStatistikに由来するとされていて、Statistikという語は、Gottfried Achenwallという人の1749年の著書の出版以降、広く使われるようになったと言われている。語源は、stateと同じで、政治学に近い分野だったらしいけど、定性的な記述を行うもので、現在の定量的な"統計学"とは全く異なるものだったようである。

英語のstatisticsという用語は、イギリスのJohn Sinclairという人が、ドイツ語のStatistikの訳として作ったらしく、1791-99年に、Statistical Accounts of Scotlandという統計データの収集を行った。Sinclairは、この用語について、 "the idea I annex to the term is an inquiry into the state of a country, for the purpose of ascertaining the quantum of happiness enjoyed by its inhabitants"などと書いていて、なんか功利主義っぽいけど、関係あるのかは不明。

cf)統計表の概念史
http://doi.org/10.14992/00011558

19世紀ヨーロッパでは、"統計"が流行ったそうなのだけど、印象的な出来事として、
1)"medical statistics"の出現
2)欧米諸国で、国家機関として、統計局が設立される
3)欧米諸国に於ける統計学会の設立
を挙げることができる。

Google scholarで検索すると、英語圏では、遅くとも1820年代には、"medical statistics"という用語が使われてたっぽい。フランスでは、1840年に、Casimir Broussaisという人による論文
De la Statistique appliquée à la pathologie et à la thérapeutique
https://books.google.co.jp/books?id=h6ZdAAAAcAAJ
が出版されている。また、1833年頃、スイスの医者っぽいF. Zschokkeなる人が、"Beiträge zur medizinischen Statistik und Epidemiologie des Bezirkes Aarau"(Contributions to medical statistics and epidemiology of the Aarau district)という論文(?)を書いているようである(?)。このへんの流れは、よく分からないけど、公衆衛生のような、政治学と医学に跨る問題が取り上げられた結果、元々は純粋に社会科学であったstatisticsに、一部の医療関係者が接近したということかもしれない。

統計局としては、ザクセン王立統計局と、プロイセン統計局は、エルンスト・エンゲルが長官を務めたことで有名。当時は、まだ統計的手法は殆ど何もなかったと言っていよく、統計局の仕事は、データを集めることだったようである。このような国家機関が出現するということは、"統計"というものが、当初の叙述的なものから定量的な性格へ変化していたということなのだろう。

日本では、1871年に"統計司"という部署が設置され、この時に、"統計"という訳語が使われている。明治初期の日本では、statisticsの訳語として、「国勢学」「知国学」「国誌学」「政表学」などが提案されたこともあり、これらは、元々の社会科学的な内容を引きずった名前になっている。統計学会は、1834年に、イギリスで、ロンドン統計学会が設立された後、欧米各国でも同様の学会ができたようである。

とにかく、そんな感じで、19世紀ヨーロッパでは、"統計学"という分野が広く認知されるに至ったようである。

Karl Pearsonが1892年に出版したThe grammar of scienceは、自然科学全般を扱っているが、statisticsという用語が沢山出てくる。Pearsonに続いて、20世紀初頭のイギリスでは、進化論、遺伝学、優生学や心理学など、従来、定量的なアプローチが採られてなかった分野での利用を動機として(※)、数理統計学が確立した。応用数学としての統計学でイメージされる内容の起原は、このへんにあるけども、数値計算法という観点からは、それほど興味を引くものがない。

※)例えば、Ronald Fisherは、統計的手法を開発しつつ、遺伝学や優生学の研究を行い、因子分析の提案者であるCharles Spearmanは、心理学者だったそうだ

【余談】いわゆる"統計学"とはずれるけど、statistical mechanicsという用語を最初に使った人は、調べた限りでは、Gibbsかもしれない。MaxwellやBoltzmannの理論は、今でも、気体分子運動論と呼ばれてたりするけど、彼らは、統計力学とは呼んでないっぽい。Gibbsは、1884年に、"On the fundamental formula of statistical mechanics, with applications to astronomy and thermodynamics"という短い報告を書き、1902年に"Elementary Principles of Statistical Mechanics"という本を出版している。

コンピュータ以前の計算道具

※)電子計算機以前のいくつかの"計算機"のハード面について、以下のブログで、解説されている。

http://parametron.blogspot.com/

タグとしては、微分解析機/古い計算機/Maxwellの面積計/面積計を使う調和解析器/手回し機械式計算機/Zuseの計算機/九元連立方程式求解機

既に見たように、20世紀初頭の応用数学書を読むと、当時、様々な計算補助のツールが存在していたことが窺える。Horst von Sandenの本には、計算尺や"計算機"の使い方が書かれていたし、Whittakerの本には、ノモグラフというツールへの言及がある。計算補助の道具は、昔から世界各地に存在していたけれど、19世紀のヨーロッパでは、計算のための道具が、特に色々開発されたようだ。

計算尺は、対数とほぼ同じ歴史を持つ。1908年の"The American Mathematical Monthly"という謎雑誌の記事

Notes on the History of the Slide Rule
https://doi.org/10.1080/00029890.1908.11997404

を読むと、20世紀初頭のアメリカで、直線の計算尺は"Gunter slide rule"と呼ばれていて、Gunterは、イギリスのEdmund Gunter(1581-1626)という人らしいけど、彼は1620年に本を書いただけで実際に制作はしなかったそうだ。その後、ニュートンの時代には、代数方程式を解くのに使われていたと書いてある。

Arithmometerは、1851年のロンドン万博で展示され、商業生産を開始したデジタル"計算機"。開発者は、フランスのCharles Xavier Thomas de Colmarという人で、保険会社を経営していたらしいので、当初の利用目的も、自然科学や工学とは異なるけど、後には、理工学系の人も利用していたらしい。1914年には、アメリカで、キーボード入力式のモンロー計算機というものが発表、販売開始された。

また、ノモグラフは、1884年に発明されたらしい。汎用性は低いけど、使い方が簡単で、誰でも扱える道具になっている。発明者は、フランスの技術者Philbert Maurice d'Ocagneという人らしい。比較的最近、以下のような本が出ている(読んではいないけど)
The History and Development of Nomography
https://books.google.co.jp/books/about/The_History_and_Development_of_Nomograph.html?id=V0A_elrmNBYC

Hilbertの第13問題(一般7次方程式を2変数の関数だけで解くことの不可能性)の意図は、私は最初見た時全く理解できなかったけど、ノモグラフが背景にあったようである

アメリカや西ヨーロッパ諸国では、そろばん/アバカスは、殆ど使われてなかったようだけど、日本や、中国、ロシアでは、そろばん/アバカスが使われていた(それ以外の地域で、どうだったのかは、よく分からない)。日本の場合は、Wikipediaの"そろばん"の項に、

日本では昭和中期くらいまでは、銀行の事務職や経理の職に就くにはそろばんによる計算（珠算）を標準以上にこなせることが採用されるための必須条件だった

とか

なお、この時代、手動式アナログ計算器としては計算尺があり、理系の人間はそちらも使いこなした

とか書いてある。理工学系の人が、どの程度、そろばんを使ってたのかは不明だけど、高精度でなくていいなら、計算尺の方が高速ではあるだろう。

中国のそろばんは、suanpan/算盤と、現地では呼ばれてるっぽいのだけど、20世紀途中まで、学校教育で教えられていたようなので、事情は、日本と大差なかったかもしれないが、詳細不明。ロシア〜ソ連でも、1990年代まで、学校教育でアバカスの使い方を教えていたらしい。

元々、西ヨーロッパ方面でも、16世紀くらいまでは使われてたっぽいのだけど、衰退した理由は謎。1886年のNatureに載ってる記事では、アラビア数字による記数法や小数点の誕生に原因を求めてるけど、根拠のある議論にも見えない

The Abacus in Europe and the East
https://www.nature.com/articles/034093a0

他に、プラニメータという面積測定機械も19世紀前半に登場し、これは積分計算機の一種と思える。国境や領地境界の形状は単純な図形ではないので、地図上の特定領域の面積計測や、pV線図から蒸気機関の行った仕事を計算するのに使われていたらしい。

19世紀後半、いくつかのアナログコンピューターの提案がされている。1876年に、物理学者William Thomsonと、兄James Thomsonは、線形常微分方程式のソルバーについて、複数の論文で記述した。
On an integrating machine having a new kinematic principle
https://royalsocietypublishing.org/doi/abs/10.1098/rspl.1875.0033

Mechanical integration of the general linear differential equation of any order with variable coefficients
https://royalsocietypublishing.org/doi/abs/10.1098/rspl.1875.0036

機械的な部分を考案したのは、前者のJames Thomsonの論文で、William Thomsonは、それを組み合わせて、任意階数の線形常微分方程式を解く方法を示唆したようである。前者の論文を読むと、冒頭で、MorinのDynamometerとかSangのプラニメータという道具への言及がされている。論文タイトルの"new kinematic principle"とは、これらの機械とは異なる機構によるということらしい。

※)James Thomsonの論文には、1851年のロンドン万博で、Maxwellが、Sangのプラニメータを目にしたことが書いてある

また、William Thomsonは、1878年に、連立一次方程式を解くアナログ計算機を提案した。
On a machine for the solution of simultaneous linear equations
https://royalsocietypublishing.org/doi/10.1098/rspl.1878.0099

これらのアナログコンピュータが、すぐに作られたわけではなく、実際にアナログコンピュータの開発が増えるのは20世紀以降っぽい。網羅的に調べてはいないので、確信はないけど、アナログコンピュータの開発は、1930年以降のものが多いように見える。

1920年代後半に、アメリカのVannevar Bushは、微分方程式を解くための機械を開発し始め、1927年の論文では、continuous integraphと呼んでいるけれど、
A continuous integraph
https://doi.org/10.1016/S0016-0032(27)90097-0

1931年の論文で、微分解析機(differential analyzer)と命名している。その後、Bushの方式に倣って、微分解析機が何台か作られた。

The differential analyzer. A new machine for solving differential equations
https://doi.org/10.1016/S0016-0032(31)90616-9

イギリスのDouglas Hartreeも、Bushから微分解析機について学び、作成している。Hartreeは、化学/物性物理/原子核物理では、Hatree-Fock近似で有名だけど、第一次大戦時には、弾道計算に従事していて、数値計算の専門家という感じの人である。

1934年に、Bushは、微分解析機を弾道計算に使うことを提案し、1935年に、アバディーン性能試験所に研究部門が設立され(数年後に、Ballistic Research Labとなった)、微分解析機は、そこで弾道計算に使われたそうだ。

連立一次方程式求解機は、MITのJohn Wilburが1936年に、以下の論文を書いている。
THE MECHANICAL SOLUTION OF SIMULTANEOUS EQUATIONS
https://doi.org/10.1016/S0016-0032(36)91387-X

日本でも模倣品が作られ、以下に解説がある

情報処理技術遺産 : 九元連立方程式求解機
https://ipsj.ixsq.nii.ac.jp/ej/?action=pages_view_main&active_action=repository_view_main_item_detail&item_id=67195&item_no=1&page_id=13&block_id=8

1930年代に作られたアナログコンピュータは、他にも例えば、イギリスのMallock machineとか、ソ連のWater integratorなどは、どちらも微分方程式の求解用だったらしい。変わったものとして、Fermiがモンテカルロ法の計算用に開発したアナログコンピュータがあったらしい(FERMIACという呼称があるっぽいけど、多分、ENIACにちなんで、ずっと後から付けられたものだろう)

19世紀から20世紀初頭に、計算器具の増加や、アナログコンピュータの開発という流れがあったのを見ると、計算に対する需要が社会全体で増大したのだろうと推測される。

今の所、エネルギーなしで動く計算機はなく、先進国では、一人あたりエネルギー消費は、減少し始めている。電力化率は伸びているけど、計算機のエネルギー効率が際限なく向上しない限り、"社会全体の計算量"も、いずれは、頭打ちになるはず

1930年代のデジタル計算機とアプリケーション

機械式のデジタル計算機であるArithmometerは、20世紀初頭には沢山存在してたけど、1930年代後半になると、電気機械式、電子式のデジタル計算機が作られるようになった。何があったのかよく分からないけど、同時多発的に、多くの人が、似たようなことを考えたっぽい。

多分、1930年頃でも、アナログツールは速いけど低精度、デジタル計算機は遅いけど高精度みたいな住み分けがあったんじゃないかと思うけど、常微分方程式の時間積分や(手計算では厳しいサイズの)連立一次方程式の数値解法など、数値誤差の蓄積が問題になりうるケースが増えてきた結果、(中間値を読み取らなくていい)アナログ計算機で頑張るか、高速なデジタル計算機を作るかという選択に、多くの人が直面したのかもしれない。

アメリカでは、George Stibitzという人が、Model Kという、リレーを使ったデジタル計算機を1936年に開発してるらしい。少なくとも、当初は、二進法の足し算を行うだけの簡単なものだったっぽい。Shannonの修論"A Symbolic Analysis of Relay and Switching Circuits"(謝辞に、Vannevar Bushが含まれる)が1937年で、それに影響を受けたものというわけではなさそう。真剣に使おうというより、おもちゃ的な位置付けだったのかもしれない。

ドイツのKonrad Zuseは、1935年か36年頃、(数値の内部表現として、二進法を採用した)機械式のデジタル計算機Z1の開発を開始した。Zuseは、1936年に、ドイツで特許"Verfahren zur selbsttätigen Durchführung von Rechnungen mit Hilfe von Rechenmaschinen"を申請しているらしいけど、現在、オンラインでは内容を確認できなかった。

cf)現物は、既に存在しないらしいし、特許文献も確認できていない。日本語での解説が、以下にある。
Zuse Z1の演算機構
https://www.iij-ii.co.jp/members/ew/zusez1.2.pdf

続いて、Zuseは、1938〜40年頃にZ2の開発を行い、1941年にZ3を完成させた。Z2やZ3では、一部で、リレーを用いていて、また、Z3はプログラム可能だったとされている。Zuseの協力者だったHelmut Schreyerは、1930年代には、電子式の計算機の可能性をZuseに提案し、自身でも、小さなものを試作したこともあったようではある。

Z3の計算速度は、加算一回一秒、乗算一回三秒とかだったようで、終戦までに完成したZuse Z4も、Z3と比較して、せいぜい数倍速いとかいう程度だったよう。初期のArithmometer(8桁の数同士の乗算で18秒)と比較しても、劇的に速いわけでもない。後のENIACは、毎秒5000回の加減算と、400回程度の乗算(※)ができたようなので、ピーク性能では、ENIACの1000倍以上、遅かったことになる。

※)1947年のNeumannとGoldstineの論文"Numerical inverting of matrices of high order"のfootnote10に、"Fully automatic electronic computing machines"(念頭にあったのは、ENIACや開発中のEDVACだと思われる)の乗算の性能について記述がある。

Zuse自身のキャリア（土木工学を専攻し、1935年に卒業後、ヘンシェル航空機製造会社でエンジニアとして勤務）のために、想定アプリケーションとしては、構造力学の比重が大きいように見える。

コンラート・ツーゼにおける計算機観
https://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/210112?locale=ja

によれば、Zuseの1940年の報告書について

5.2節で例として挙げられているのは航空機建造における外郭構造計算（Berechnungvon Schalenkonstruktionen im Flugzeugbau）と気象計算（Wetterrechnung）である．ツーゼは航空機の構造設計や気象予報のような分野は理論的には大きく進歩したにもかかわらず，その計算があまりにも膨大であるという理由で発展が妨げられてきたと述べたうえで，自らの計算機を使えばこうした問題が（解析的にではなく）計算的に解を見出だせるようになる可能性があるとしている

とあり、5.3節でも

例えばモーターの構造や鋳物の正確な剛性計算，振動計算，複雑な航空力学的問題，原子論などといった数多くの領域が，これまで不十分な仕方でしか計算的に開拓されてこなかった．ここにおいてもこの機械は費用のかかる実験を省略し，これまでに予期し得なかったほどの学術的進展を促す可能性がある

のように書いているらしい。機械/航空工学的な問題への適用は、Zuse自身の実体験に基づくものかもしれない。気象予測は、Richardsonの1922年の本に影響されたものだろう。

アメリカのJohn Attasoffは、1936年にアナログコンピュータの開発(詳細はよく分からない)を行った後、1937年頃から、連立一次方程式を解くことを目的とした、(プログラム可能でない)電子式デジタル計算機の開発を考えるようになった。院生のClifford Berryと共同作業だったので、Attanasoff-Berry computerとかいう名前で呼ばれている。電子式だと、数値の内部表現に二進法を採用する方が簡単だと思うので、AttanasoffとBerryのコンピュータも、そうなってたようである

Atanasoff-Berry Computerの最初の計算は1942年に披露され、最大で29変数の線形連立方程式を解くことができたらしい。Wikipediaには、1秒間に30回の加減算ができたと書いてるので、演算速度自体は、電子式でないZuseのコンピュータよりは速い。

1940年に書かれたAttanasoffの報告

Computing Machine for the Solution of large Systems of Linear Algebraic Equations
https://link.springer.com/chapter/10.1007/978-3-642-61812-3_24
[PDF] http://jva.cs.iastate.edu/img/Computing%20machine.pdf

には、想定アプリケーションとして、以下のように書かれている。

The following list indicates the range of problems in which the solution of systems of linear algebraic equations constitutes an essential part of the mathematical difficulty.
1. Multiple correlation
2. Curve fitting
3. Method of least squares
4. Vibration problems including the vibrational Raman effect
5. Electrical circuit analysis
6. Analysis of elastic structures
7. Approximate solution of many problems of elasticity
8. Approximate solution of problems of quantum mechanics
9. Perturbation theories of mechanics, astronomy and the quantum theory

Atanasoffは、電気工学で学士号、数学で修士号、理論物理で博士号を取ったそうで、その経歴を反映してか、技術的なものより、基礎科学への応用の方が詳細な印象を受ける。

理工学的な問題でなくても、初等的な問題は、しばしば、連立一次方程式に帰着する。一例として、Leontiefは1941年に"The structure of American economy"という本を出版していて、44セクターの産業連関表を作ったけど、44x44の逆行列は計算できないので、10セクターにまで整理したらしい(本は、オープンアクセスにはなってないので、確認はしてない)。

上のAtanasoffの報告には、"Since an expert computer requires about eight hours to solve a full set of eight equations in eight unknowns"という記述がある。computerが計算機を指すとすると遅すぎる(※)ので、人間/計算手を指していたと思われる

(※)8変数8個の連立一次方程式を解くのに必要な算術演算は、除算が36回と、乗算及び減算が196回ずつなので、除算を無視すると、乗算や減算一回に平均一分かかる場合、演算そのものに必要な時間は6.5時間。最初のArithmometerの計算速度(8桁の数同士の乗算を18秒)から考えても、計算機だったら、8時間は遅すぎるだろう

あまり関係ないけど、デジタル信号処理は、1940年代に現れたと言っていいと思う。二次大戦中、アメリカとイギリスの間で、暗号化された音声通信を行うために、(1940年8月に)SIGSALYプロジェクトが開始され、SAISALYで、デジタル音声通信が実現した。標本化定理が、当時、知られていたのか不明(※)だし、量子化の方も、明確な提案者というのは、いないよう。とりあえず、ハードとして、A/D変換器とD/A変換器が作られたのは、この時らしい

cf)以下の文献によれば、標本化定理を述べたのは、1915年のWhittakerの論文が最古だそうで、定理としては存在していたけど、一般に(特に、技術者に)知られていたかは分からない。また、Nyquistは、1920年代に標本化定理を発見していたと書かれてることもあるけど、Nyquistの論文に、明示的に標本化定理は書かれていないようである。NyquistとShannonは、共にSIGSALYプロジェクトに参加していて、Shannonは戦後に標本化定理の論文を発表しているので、SIGSALY参加者たちは、既知の事実とは見なしてなかったんじゃないかと推測される
Interpolation and Sampling: E.T. Whittaker, K. Ogura and Their Followers
https://link.springer.com/article/10.1007/s00041-010-9131-8

アメリカのENIACは、スポンサーがBallistic Research Labで、砲外弾道学が応用と考えられていた(実際に最初に行った計算は、核爆弾開発のためのモンテカルロ・シミュレーションだったらしい）

1930年代前半に、微分解析機の調査にアメリカに行ったHartreeは、ENIACも視察に行ったようで、1946年に、Natureで紹介記事を書いている。

The Eniac, an Electronic Computing Machine
https://www.nature.com/articles/158500a0

常微分方程式の数値解法と砲外弾道学

常微分方程式の数値解法の重要な応用の一つは天文学だろうけど、19世紀の天文学について書くのは辛いので、もう一つの応用であり、かつ、文献で詳しく言及されることが少ない弾道学について、概要を見ることにする。

単に、弾道学(ballistics)というと、弾道計算をイメージするけど、現在では、これは砲外弾道学(exterior ballistics)と呼ばれる。他に、砲内弾道学とか終末弾道学のような区分が存在する。

ヨーロッパに於ける砲外弾道学は、16世紀イタリアのタルタリアに始まるとされている(他の地域、特にオスマン帝国で、同種の試みがなかったのかは不明)。彼の本を見ても、放物軌道は見当たらない
Nova scientia
https://archive.org/details/novascientia00tart/page/n10

放物運動は、ガリレイ(タルタリアの孫弟子になるらしい)の『新科学対話』で提示されたものらしいけど、実際には、空気抵抗のせいで、放物軌道を描かないことも、多くの人に気付かれていたようである。そもそも、タルタリアの本にある挿絵とかを見ても、水平方向の速度が失速しているように描かれている。もし、ガリレオが一般的に言われてるように、先入観を持たず観測事実をありのままに記述するような人物だったなら、放物軌道は却下されていたかもしれない。

Ballisticsという用語は、フランスのMersenneによる本"Ballistica et Acontismologia"(1644年)に由来するそうだ(語源としては、"I throw"を意味するギリシャ語βάλλωに遡るらしい)。この頃は、砲外弾道学/砲内弾道学/終末弾道学という区分はなく、弾道学=砲外弾道学だったのだろう。

少し時代は下って、イギリスのBenjamin Robins(1707-1751)は、1742年にnew principles of gunneryという本を書いている。
New principles of gunnery
https://archive.org/details/newprinciplesgu00wilsgoog/page/n7

この本は、Eulerが、ドイツ語訳を行ったことでも知られている。Eulerは、砲外弾道学に関する論文を複数書いていて、以下のものがよく引用されている。
Recherches sur la véritable Courbe que décrivent les Corps jettés dans l'air, ou dans mn autre fluide quelconque
http://eulerarchive.maa.org/pages/E217.html

19世紀前半は、砲外弾道学が、あんまり盛り上がってなかったように見えるけど、19世紀後半になると、砲外弾道学の研究は再び盛んになった。直接の契機は、弾丸形状の変化にあるっぽい。クリミア戦争や南北戦争(日本だと、戊辰戦争)の頃には、ミニエ銃やエンフィールド銃というものが使われ、小銃の弾丸はミニエ弾と呼ばれるものに変わっていった。同時期に、大砲の砲弾形状にも変化があり、1858年にイギリス軍で採用されたアームストロング砲の砲弾は、ミニエ弾と類似の形状(当然、大きさは違う)だったらしい。

Francis Bashforth(1819-1912)の1868年の論文

On the resistance of the air to the motion of elongated projectiles having variously formed heads
https://royalsocietypublishing.org/doi/abs/10.1098/rstl.1868.0015

は、タイトルから、球状でない弾丸の空気抵抗に関する議論であることが分かる。この論文の冒頭で過去の研究に触れていて、18世紀は、多くの人が名前を挙げられているけど、19世紀前半は、散発的な実験があっただけという印象。

19世紀終盤には、多分、火薬の変化があったせいで(※)弾速も速くなった。これは、砲外弾道学研究のもうひとつの動機になったと思われる。火薬兵器の登場から500年以上に渡って、弾丸初速は、音速を少し上回る程度で、上Bashforthの論文で扱われてる速度も、最大で、1500ft/s≒460m/s程度に留まっている。これが19世紀末になると、初速がマッハ2を超える大砲も作られるようになっているので、Bashforthの測定は不十分なものとなったはず。

※)1845年に、スイスの化学者がニトロセルロースの合成法を発見し、ついでに、1847年頃、イタリアの化学者が、ニトログリセリンの合成法を公表している。1880年代になると、ニトロセルロースをベースとして、最初の無煙火薬が実用化された。また、ピクリン酸も、ほぼ同時期に火薬として利用されるようになった。米西戦争や日露戦争の頃には、各国が無煙火薬を使用し、日露戦争終盤の日本海海戦では、日本海軍は、ピクリン酸を主成分とする下瀬火薬も使用したとされている

1870〜1900年までに、各国で、砲外弾道学の本が出版されている。例えば、1872年に、Nicolas Mayevskiという人が

Traité de balistique extérieure
https://books.google.co.jp/books?id=jE_zAAAAMAAJ

という本を書いている(ロシアの人らしいけど、何故かフランス語?)。この人については、情報が少なく、詳細な経歴は分からない。

日本には、明治時代に、イタリアからScipione Braccialiniという人がやってきて、「砲外弾道学」(1894年出版)という本を残している。国立国会図書館で、公開されているので、見ることができる
砲外弾道学. 1: https://dl.ndl.go.jp/info:ndljp/pid/844786
砲外弾道学. 2: https://dl.ndl.go.jp/info:ndljp/pid/844787
砲外弾道学. 3: https://dl.ndl.go.jp/info:ndljp/pid/844788
砲外弾道学. 4: https://dl.ndl.go.jp/info:ndljp/pid/844789
砲外弾道学. 5: https://dl.ndl.go.jp/info:ndljp/pid/844790
砲外弾道学. 6: https://dl.ndl.go.jp/info:ndljp/pid/844791
砲外弾道学. 附表: https://dl.ndl.go.jp/info:ndljp/pid/844792

縦書きで、全部合わせると1000ページ近い。

cf)1000ページも読むのは辛いので、概要を掴むのに、以下の文献とか読むと良さそうなのだけど、入手出来てない。
砲外弾道学テキストに見る明治の微分方程式解析法 : 運動方程式の解析法に関する技術数学
https://ci.nii.ac.jp/naid/10024452473

アメリカでは、James M. Ingalls(1837-1927)という人が、弾道学関係の本を何冊も書いていて、有名。弾道学研究者でも、大抵は、どこそこのprofessorとかいう肩書だけど、Wikipediaで、IngallsはAmerican soldierと紹介されている。Ingallsは、彼が計算して作った表で知られ、Google booksで公開されている。

Ingalls' Ballistic Tables
https://books.google.co.jp/books?id=nOOgAAAAMAAJ

この本のTABLE Iは、元々は、1893年に計算したものだそうだけど、今日、Siacci's methodと呼ばれる(近似)手法で使われる補助関数の値になっている。

Siacci's methodの考案者は、イタリアのFrancesco Siacci(1839-1907)という人で、1875〜1892年まで、トリノ大学で、力学の教授をやっていた。Siacciの後任はVito Volterraで、Volterra積分方程式やLotka-Volterra方程式などで知られる。Siacciは、砲外弾道学の本も書いているけど、弾道学以外の業績もあるらしい

Siacci's method自体は、多分、以下の2つの理由で、歴史的価値しかないと思う

(1)Siacci's methodは、平面上を運動する質点系に使う方法だけど、精度を追求する場合、弾丸を、3次元空間中を飛行する剛体として扱いたい
(2)そこまで精度を追求しない場合でも、常微分方程式を積分する標準的なアルゴリズムと、普通のコンピュータを使えば、数値積分できる

しかし、19世紀末〜20世紀初頭の人が、どういう計算をしていたか知るには、丁度いい題材かもしれない。

当時の弾道計算では、弾丸に働く力は、重力と空気抵抗のみで、平面内を動くとして、空気抵抗の大きさは、弾丸の速度のみの関数とするモデル化が一般的に使われていた。よく知られている通り、物体が音速に比べて十分低速であれば、空気抵抗の大きさは、速度の二乗に比例すると置く近似が有効だけど、銃の登場初期から、弾速は、しばしば音速を超えていたので、空気抵抗の速度依存性は、そんなに単純な関数では書けない。これは、Benjamin Robinsによって、すでに知られていた。

当時は、弾丸の速度と空気抵抗力の関係を測定した後、何らかの関数でフィッティングすることが行われた。しばしば、Mayevskiの公式と呼ばれるフィッティングは、速度区間を、いくつかに分割して、各区間では、 $A v^n$ の形になっているもので、Ingallsの本では、例えば、速度が0~790(ft/sec)の区間では
$R = Cr = 10^{5.668914-10} v^2$
で、790~970(ft/s)の区間では
$R = Cr = 10^{2.7734430-10} v^3$
などという式が書かれている。速度 $v$ の単位は、全部(feet/sec)。Cはballistic coefficientと呼ばれる量で、弾丸の質量を、弾丸直径の2乗で割った量(しばしば、形状などに応じて、補正係数が掛けられる)。 $r$ が弾丸が受ける抗力で、 $R=Cr$ が、弾丸速度のみに依存する関数になるとしている。

[注意]Ingallsの本に書かれてる式は、OCRのミスか、元々ある間違いなのか知らないけど、係数Aがおかしいので、ここでは、本来意図されたはずの形に修正してある

本では、関数は、3600ft/s≒1100m/s(〜マッハ3)くらいまで用意されている。もうひとつ、Siacciのempirical formulaとして、形は複雑だけど、単一の式で、全速度域をカバーする代数的な関数も書かれてる(式の具体的な形は、本の最初の方を参照)。

とにかく、R=R(V)として何か関数が与えられると、考える微分方程式は
$\dfrac{dV_x}{dt} = -\dfrac{R(V)}{C V}V_x$
$\dfrac{dV_y}{dt} = -\dfrac{R(V)}{C V}V_y - g$
となる。xが進行方向で、yが鉛直方向。 $C$ は、ballistic coefficientで、 $g$ は重力加速度。 $V=\sqrt{V_x^2+V_y^2}$ で、物体が受ける空気抵抗の大きさは $R(V)/C$ になっている。

Ingallsの本のTABLE Iの補助関数の定義は、
$S(u) = \displaystyle \int_{u}^{u_{max}} \dfrac{V}{R(V)} dV$
$I(u) = \displaystyle \int_{u}^{u_{max}} \dfrac{2g}{V R(V)} dV$
$A(u) = \displaystyle \int_{u}^{u_{max}} \dfrac{V \cdot I(V)}{R(V)} dV$
$T(u) = \displaystyle \int_{u}^{u_{max}} \dfrac{dV}{R(V)}$
となる。 $u_{max}$ は、任意に決めていいのだけど、Ingallsの本では、3600(feet/sec)とされている。Mayevski型の式は、この手の積分を計算するには、とても都合が良い。

この補助関数の出処というか使い方だけど
$(V_x,V_y) = (V\cos \theta , V \sin \theta)$ として、 $V,\theta$ に対する方程式に書き換えると
$\dfrac{dV}{dt} = -\dfrac{R(V)}{C} - g \sin \theta$
$V \dfrac{d\theta}{dt} = -g \sin \theta$
となって、特に、[\theta \approx 0]の時は、
$\dot{V} = -\dfrac{R(V)}{C}$
となる。この微分方程式を積分して、時間tを速度Vの関数として書く場合、
$t = C \left( T(V) - T(V_0) \right)$
となる。

Ingallsの本を見ると、ちょっと違う式を書いてあって、これはpseudo-velocityと呼ばれる量を使ってる(通常、pseudo-velocityの使用はSiacci's methodの一部として説明されるので、Siacciが始めたのかもしれない)からだけど、pseudo-velocityは、普通の速度と定数倍違うだけなので、特に必要ない。当時の人が、どういう気持ちで、わざわざpseudo-velocityを導入したのか謎。

pseudo-velocityを使わず、普通の速度を使う場合、位置 $x$ は
$x = x_0 + C \cos \phi (S(V) - S(V_0)$
によって計算できる。同様に、 $y$ や $\tan \theta$ も、補助関数を使って、速度の関数として書ける。Siacci's methodは、大体、こういうもの

第一次大戦では、参戦国は、砲外弾道学の研究に、人材を投入したようだけど、フランスでは、著名な数学者も、かなり砲外弾道学研究に参加したらしい。何人かの人は報告を書いている。例えば、Hadamardの報告は

Rapport sur les travaux examinés et retenus par la Commission de Balistique pendant la durée de la guerre
https://zbmath.org/?q=an:47.0727.01

で読める。

Siacciは、 $V$ を従属変数にするような変換を考えたけど、 $\theta$ を従属変数とする微分方程式を考えると
$\dfrac{dx}{d\theta} = -\dfrac{V^2}{g}$
$\dfrac{dy}{d\theta} = -\dfrac{V^2}{g} \tan \theta$
$\\dfrac{d(V \cos \theta)}{d\theta} = \dfrac{V R(V)}{gC}$
が得られ、最後の方程式から、 $V=V(\theta)$ が得られれば、 $x=x(\theta)$ と $y=y(\theta)$ の解は、定積分の計算に帰着する。

その性質から、この最後の方程式は、principal equationと呼ばれることもある。フランスの数学者Jules Drachは、更に、変数 $\tau = \sin(\theta)$ を導入して、 $V$ と $\tau$ を変数にして
$\dfrac{d(\log V)}{d\tau} = \dfrac{R(V) + \tau}{1-\tau^2}$
の形の方程式を調べている。Drachの1920年の論文は、95ページもある上に、フランス語なので、詳細を追うには辛すぎる。数値計算の研究というわけではないし、これが何かの役に立ったかは疑わしい

L'équation différentielle de la balistique extérieure et son intégration par quadratures
http://www.numdam.org/item/ASENS_1920_3_37__1_0/

イギリスでは、(物理方面の人には、Hartee-Fock近似で知られる)Douglas Hartreeの1920年の論文が、Natureに載っている(昔のNatureは、時々、弾道学の話題も載ったようで、Bashforthも1884年に寄稿している)。
Ballistic Calculations
https://www.nature.com/articles/106152a0

数学的内容はないけれど、高い場所からの砲撃(※)を行うようになって、高さに依存する空気の密度の変化を考慮しないといけなくなったとか実際的なことが書いてあって面白い。また、"To perform this integration a step-by-step method is employed"とか書いてる。空気の密度の高度依存性を考慮するとかしだすと、Siacci's methodが有効でなくなり、一次大戦頃から、汎用的な常微分方程式の数値解法に移行したっぽい(?)。このような手法は、天文学では以前から行われていたが、砲外弾道学では、大変と思われて、使われてなかったと書いてある

(※)"when guns began to be used at higher elevations, and the muzzle velocities of howitzers were increased, these approximate solutions became unsatisfactory"とあり、高所からの砲撃が行われるようになった理由は、説明がない(大砲は重いので、航空機に載せたわけではないだろう)。多分、射程が伸びて、目視での着弾確認が難しく(海岸線から水平線までの距離が4〜5kmなので、射程が10kmとかになると、着弾の確認も難しいと思われる)、高所から撃つことが増えたとか、角度を付けると、弾丸が相当高くまで打ち上がるようになったとか推測できるけど、詳細は分からない

弾丸を質点として扱うモデルが長く使われていたが、3次元空間を移動する剛体として、弾丸を取り扱おうというのは誰でも考える。このような方向での研究は、1921年のイギリスで論文が出ている

The aerodynamics of a spinning shell
https://royalsocietypublishing.org/doi/10.1098/rsta.1921.0010

20世紀後半の砲外弾道学の論文を見ると、上記論文の次に、しばしば引用されてるのは1946年の論文

On the motion of a spinning shell
https://doi.org/10.1090/qam/17019

で、時間の断絶があるので、電子計算機が使えるようになるまでは、あまり盛んに使われるモデルでもなかったっぽい。

このように、砲外弾道学に於いては、第一次大戦以降、モデルの精密化が進み、常微分方程式の数値解法が必要になっていったことが分かる。

おまけ。Siacciの空気抵抗則とMayevskiの空気抵抗則で、弾道が、どの程度変わるか、直接、微分方程式を数値積分して確認してみる。Ingallsの本の"examples of the use of the tables"の最初の例では、"Take the 6-inch gun, M.V.=2,600foot-seconds, firing a 108-pound long-pointed projectile. We will take C=4.894 and ..."という記述があり、これを例題として使う。スペック的には、丁度1900年頃開発された

BL 6-inch Mk VII naval gun
https://en.wikipedia.org/wiki/BL_6-inch_Mk_VII_naval_gun

とかが近そう？弾丸重量/弾丸直径の二乗は、108/(6*6)=3.0(lb/in^2)だけど、補正されたCの値4.894を、使うことにする

また、初速を800(m/s)、角度を15度とする。計算によると、着弾時の速度は、約310(m/s)≒1100(km/hr)で、予測される着弾地点は、空気抵抗則によって、180メートルほど違う。割と大きい誤差という気がするけど、放物軌道との違いを見ると、弾道計算の補助なしで勘で撃った場合、何キロも離れたところに落ちても不思議ではなさそうなので、この程度の精度でも、十分役立ったのかもしれない

f:id:m-a-o:20200531181540p:plain

画像は、以下のコードで描画したもの

import math
import matplotlib.pyplot as plt
import numpy as np

# the units are [ft・lb/(sec^2・in^2)]
def siacci(v):
  a = 0.284746
  b = -224.221
  c = 0.234396
  d = -223.754
  e = 209.043
  f = 0.019161
  g = -984.261
  h = 371.
  i = 656.174
  return ( a*v+b+math.sqrt((c*v+d)**2+e) + f*v*(v+g)/(h+(v/i)**10) )*0.896


def mayevski(v):
  A = [math.pow(10, 7.6090480-10) , math.pow(10,7.0961978-10) , math.pow(10 , 6.1192596-10) , 
       math.pow(10, 2.9809023-10), math.pow(10,6.8018712-20), math.pow(10,2.7734430-10),math.pow(10,5.6698914-10)]
  N = [1.55 , 1.7 , 2 , 3 , 5, 3 , 2]
  if v>2600:idx = 0
  elif v>1800:idx=1
  elif v>1370:idx=2
  elif v>1230:idx=3
  elif v>970:idx=4
  elif v>790:idx=5
  else:idx=6
  return A[idx]*math.pow(v,N[idx])


# units of C are [kg/m^2]  g=9.8(m/s)
#-- vars = (x,y,Vx,Vy)
def f(vars ,t , Cinv , g , type):
  def ft2meter(L):return L*0.3048
  def inch2meter(L):return L*0.0254
  def lb2kg(M):return M*0.453592
  def meter2ft(L):return L/0.3048
  Vx , Vy = vars[2],vars[3]
  V = math.sqrt(Vx*Vx+Vy*Vy)
  if type==0:  #-- siacci air resistance law
     R = Cinv*siacci(meter2ft(V))*((lb2kg(1.0)*ft2meter(1.0))/(inch2meter(1.0)**2))
  else:        #--mayevski air resistance law
     R = Cinv*mayevski(meter2ft(V))*((lb2kg(1.0)*ft2meter(1.0))/(inch2meter(1.0)**2))
  return [Vx , Vy , -R*Vx/V , -R*Vy/V - g]


if __name__=="__main__":
  from scipy.integrate import odeint
  V0 = 800.0    #--(m/s)
  phi0 = 15.0*math.pi/180  #--(rad)
  C = 4.894*0.453592/(0.0254**2) 
  times = np.linspace(0, 33.0, 10001)
  #--siacci
  sols = odeint(f , [0 , 0 ,V0*math.cos(phi0) , V0*math.sin(phi0)] , times , args=(1/C , 9.80665 , 0))
  plt.plot(sols[:,0] , sols[: ,1] , color='blue' , label='siacci ARL')
  #-- mayevski
  sols2 = odeint(f , [0 , 0 ,V0*math.cos(phi0) , V0*math.sin(phi0)] , times , args=(1/C , 9.80665 , 1))
  plt.plot(sols2[:,0] , sols2[: ,1] , color='red',label='mayevski ARL')
  #-- no ARL
  sols3 = odeint(f , [0 , 0 ,V0*math.cos(phi0) , V0*math.sin(phi0)] , times , args=(0.0 , 9.80665 , 0))
  plt.plot(sols3[:,0] , sols3[: ,1] , color='green',label='Parabola')
  plt.xlabel('x(meter)')
  plt.grid()
  plt.legend()
  plt.show()

偏微分方程式の数値解法と構造物の計算

現代の工学では、盛んに偏微分方程式の数値計算が行われるけど、20世紀初頭には、そんなことは無理だったわけで、最初に、当時の技術者はどうしていたのか、概略を見る。

建物、道路、橋、ダム、堤防、鉄道、船、車、航空機、機械などの構造物設計に関わる計算は、工学では対象ごとに細かく分野が分かれているけど、個別に調べるのはしんどいので、解析法の大きな括りとして(勝手に)以下の3つに分ける。

(1)静力学的解析
(2)連成振動系による振動解析
(3)弾性体の偏微分方程式に基づく振動解析

Attanasoffが列挙しているapplicationで言うと、

4. Vibration problems including the vibrational Raman effect
6. Analysis of elastic structures
7. Approximate solution of many problems of elasticity

あたりに相当するもの。ラマン効果は、量子力学的な問題だけど、Attanasoffの意図では、多分、古典的な振動論(振動工学とかの対象になるもの)と、両方、念頭にあったんじゃないかと思う。

【補足】分子振動の解析の手法として、初期には、GF methodというものが使われ、アメリカの化学者Edgar Bright Wilson(1908-1992)が1939年に書いた論文

A Method of Obtaining the Expanded Secular Equation for the Vibration Frequencies of a Molecule
https://doi.org/10.1063/1.1750363

あたりが起原と思われる。

静力学的解析(と勝手に呼んでるもの)は、力やモーメントの釣り合いの計算で、素朴に考えると、連立一次方程式に帰着する。建築物なんかの構造物は、激しく変形しないのが普通で、こういう計算は特に有益と思われる。対照的なのは、ロボットや、生体力学の問題で、静力学では済まないことが多いと思われるけど、19世紀には、殆ど扱われてないっぽい

現代の建築関係の教科書を見ると、トラス(truss)構造というのが、説明されてたりする。トラス構造自体は、古代ローマからあったそうだけど、19世紀に金属トラス橋が現れ、その静力学的な解析手法は、19世紀に研究された。研究した人の中には、MaxwellやMöbiusのような人も含まれる。

トラス構造は、木や鉄の棒を端っこ同士で接合して、三角形を単位として作られる骨組の一種だとか書いてある。鉄道橋以外に、エッフェル塔、東京タワー、スカイツリーなどをよく見ると、三角形を単位として組まれているのがわかるけど、あれのことらしい。東京タワーの設計は、計算尺を使って行われたらしいので、電子計算機以前の技術によって設計されたと言える。

19世紀に行われた、最も単純な解析では、棒にかかる力を軸力(引張or圧縮。曲げや捩れは無視する)に限定して、節点ごとに力の釣り合いを考えるというタイプの方法だったらしい。何も考えずに、方程式を書くと、手で解くには変数の数が多くなると思われるし、変数の数と方程式の数が一致するという保証もないので、19世紀には、こうした課題に対応するための手法が色々作られたようだし、20世紀初頭には、また別の発展もあったようだけど、建築学固有の話すぎるので省略。

弾性体の振動は、物理では古くから研究がある。工学分野では、19世紀終わり頃〜20世紀になると、蒸気タービンの出現によって、機械類の振動問題が関心を集めるようになった。(ラバールノズルで有名な)スウェーデンのde Lavalや、イギリスのCharles Algernon Parsonsが、1880年代に蒸気タービンを考案し、1904年時点で、de Lavalの蒸気タービンが10000〜30000rpmで動作していたようである

cf)10000〜30000rpmという記述は、The Electrical World and Engineer, 第 43 巻(1904年刊?)の1086ページにあった
https://books.google.co.jp/books?id=geVQAAAAYAAJ&pg=PA1086&lpg=PA1086

Parsonsは、1897年に蒸気タービンを搭載した実験艇のデモを行ったそうで、この場合も数万rpmで動作したようである。それ以前の"rotor dynamics"の研究では、危険速度を超えると、回転軸の振れ回り運動が不安定になると結論されていたらしく、数万rpmという回転速度で安定動作できるというのは、予想されてなかったことのようだ。

静止弾性体の微小振動について、連続体モデルで方程式を書いて、有限要素法などの近似を使うと、連成振動系と同じ形の常微分方程式系に帰着する。回転する剛体の運動方程式は、Euler topやLagrange topとして知られ、無条件には非線形になってしまうが、回転速度を一定として軸対称性を仮定すると、線形な運動方程式を得ることが出来る。同様に、軸対称な回転する弾性体で、軸周りの回転速度が一定の場合の運動方程式は、gyroscopic matrixを追加するだけで、定数係数の2階線形常微分方程式系に帰着できる。

定数係数の2階線形常微分方程式系は、数学としては面白みのない問題ではあるけれど、自由度が大きいと、(例外的な場合を除いて)コンピュータなしでの計算は困難になるので、コンピュータのない時代は、ごく少数のモードのみを取り扱うようにモデル化されていたっぽい。現代でも工学の教科書を見ると、1自由度や2自由度での議論が最初に載ってるのを見ることが出来る

土木工学に於ける振動の研究は、地震への応答に動機があって始まったようである。1906年のサンフランシスコ地震や、1923年の関東大震災などがあって、日本や欧米でも、この手の研究がされるようになった。

以下の"二木造洋風建築物の振動験測結果報告"の引用文献は、国内のものばかりだけど、1920年代のものが多い
https://www.jstage.jst.go.jp/article/zisin1929/4/11/4_11_657/_article/-char/ja/

引用文献の中に、末廣恭二という人のものがあるけど、Wikipediaによると、最初は造船工学者として出発し、三菱造船所の研究所の所長になり、その後、地震研究所初代所長になったそうで、こういう研究が当時は新しいものだったのだろうことが推察できる。

寺田寅彦の追悼記事に"大正十二年関東大震災以前から既に地震学に興味をもっていたが、大震災の惨害を体験した動機から、地震に対する特殊の研究機関の必要を痛感"、"地震学に興味をもつようになったのは船舶の振動に関する研究から自然に各種構造物の振動に関する問題の方に心を引かれるようになったためであるらしい。"などと書かれている。

cf)工学博士末広恭二君
https://www.aozora.gr.jp/cards/000042/files/43076_18496.html

また、違った種類の現象として、1920年代には、航空機のフラッタが研究されるようになった。翼に働く力には流体力学的な力を考慮する必要があるので、1920年代初頭、薄翼理論が作られて、ある程度、揚力、モーメントの計算ができるようになったため、こういう計算が始まったのだと思われる。このへんの話題も、現代の工学の教科書を見れば載ってるし、数値計算としては、別に面白くないので略。

弾性体の振動は、通常、偏微分方程式で記述される。しばしば境界値問題の形で定式化され、Dirichlet境界条件とかNeumann境界条件とかに名前が残ってるDirichlet(1805-1859)やCarl Neumann(1832-1925)は、19世紀の人であることからも、問題の定式化が確立したのは19世紀だろうと思う

偏微分方程式の数値解法を試みた論文が、20世紀初頭に、いくつか出ている。10x10の行列の計算にも苦労する時代に、実用的な計算ができたのかは謎だけど、調べた限り、電子計算機の登場以前には、偏微分方程式の数値計算が広く行われたとは言えないようである。

当時の論文を読むと、biharmonic equationあるいはGermain-Lagrange方程式と呼ばれる4階の偏微分方程式が扱われている。太さや厚さを無視した弦や膜の(振幅の小さい)振動については、2階の波動方程式で記述されるけど、構造材に使われる梁や板のように、太さや厚さを無視できない場合には、Euler-Bernoulli beam theoryとか、Kirchhoff–Love plate theoryというものに基づいた4階の微分方程式を使い、Kirchhoff–Love plate theoryの特殊ケースに於ける支配方程式が、Germain-Lagrange方程式。細かろうが薄かろうが3次元物体であることには変わりないので、原理的には、3次元領域での(2階の)偏微分方程式を使って記述できるけど、昔の人は、空間次元を減らして扱いやすくする道を選んだようである

1910年頃までに、偏微分方程式の数値解法として、差分近似による方法と、解を適当な関数の線形結合で書く方法が試みられた。微分を差分で近似するのは、常微分方程式でも、散々やってるから、それ自体は目新しくはない。Carl Rungeの1908年の論文

Über eine Methode, die partielle Differentialgleichung Δu=constant numerisch zu integrieren
Z. Math. Phys.56, 225–232

は、差分法による数値計算を試みたっぽいのだけど、ネット上では、Abstractも確認できなかった。

Richardsonは、1910年頃、偏微分方程式の差分法による数値計算について論文を書いている。
On the approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam
https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1910.0020

The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam
https://royalsocietypublishing.org/doi/abs/10.1098/rsta.1911.0009

後者の論文は、一次元の拡散方程式や二次元のラプラス方程式などの例を扱った後、biharmonic equationでの計算を行っている。冒頭の記述によると、複雑な形状の物体では解析解を期待できないので、数値解が必要と考え、例としてダムの計算を扱っている。このような差分近似による計算は、電子計算機以前には、実用化には至らなかったように思う。

Richardsonは、その後、Vilhelm Bjerknesの1904年の論文に書かれた着想に従って、現在、primitive equationと呼ばれる方程式群を書き下して、数値計算を試みた。これは数値天気予報の最初の試みとして知られている。Richardsonの試みは、1922年に出版された

Weather prediction by numerical process
https://archive.org/details/weatherpredictio00richrich/

Richardsonの計算をトレースしてみようかと思ったけど、ページ数多いので、ちょっと辛かった。そのうち、気が向いたらやってみたい

cf)Bjerknesの論文の英訳が、以下にあるけど、数式は一つもない。性質を調べるのが難しい方程式を書いてみることに、意味を感じなかったのかもしれない
The problem of weather prediction, considered from the viewpoints of mechanics and physics
https://www.schweizerbart.de/papers/metz/detail/18/74383/The_problem_of_weather_prediction_considered_from_the_viewpoints_of_mechanics_and_physics

偏微分方程式を差分近似すれば数値計算できることは、殆ど明らかとはいえ、数学的には自明とは言えない。1928年には、現在、CFL条件で有名なCourant-Friedrichs-Lewyの論文が出版された。原文は、ドイツ語だけど、英訳が存在している
On the partial difference equations of mathematical physics
[PDF]https://web.stanford.edu/class/cme324/classics/courant-friedrichs-lewy.pdf

この論文の目的は、古典的な線形偏微分方程式で、差分法による解が、ちゃんと微分方程式の解に収束することを示すという点にある。ここでは、biharmonic equationと、波動方程式が扱われている。Courantたちの研究が、数値計算での差分法の利用を念頭に置いたものだったかは不明だけど、これによって、偏微分方程式の解の存在証明ができた。

Nowhere do we assume the existence of the solution to the differential equation problem - on the contrary, we obtain a simple existence proof by using the limiting process.

Our method of proof may be extended without difficulty to cover boundary value and eigenvalue problems for arbitrary linear elliptic differential equations and initial value problems for arbitrary linear hyperbolic differential equation

常微分方程式では、Cauchyの折れ線近似(polygonal approximation)による解の構成が知られていて、折れ線近似は、アルゴリズムとしては、Euler法と同じものなので、Euler-Cauchy polygonal approximationと呼ばれてたりする(多分、解の存在証明に利用しようと試みた最初の人が、Cauchyだったのだろう)けど、差分近似の極限として解を構成する発想は同根と言える。

また、Courant-Friedrichs-Lewyの論文のSection3は、"Connections with the problem of the random walk"となっており、偏微分方程式の境界値問題に対するモンテカルロ法を示唆していると見ることができる。

In addition to the main part of the paper, we append an elementary algebraic discussion of the connection of the boundary value problem of elliptic equations with the random walk problem arising in statistics.

と書いてるので、オマケ的位置付けだったっぽいし、彼らが、数値計算に使うことを想定していたかは分からない。少なくとも、数値計算例はない

解を適当な関数の線形結合で書くという発想は、Laplaceが球面調和関数展開をした(※)り、Fourierが三角級数を使って熱方程式の解を書いたあたりに、萌芽が見られる。

※)Théorie des attractions des sphéroïdes et de la figure des planètes
http://sites.mathdoc.fr/cgi-bin/oetoc?id=OE_LAPLACE__10

ドイツのWalther Ritzは、1909年頃、微分方程式が変分構造を持つ時、解を有限個の既知関数の線形結合で近似し、変分汎関数に放り込んで停留条件を課して、線形結合の係数を決めるという方法を提案した。線形偏微分方程式に対しては、Ritz法は、線形代数の計算に帰着する

[Ritz1908] Über eine neue Methode zur Lösung gewisser Variationsprobleme der mathematischen Physik
https://doi.org/10.1515/crll.1909.135.1

[Ritz1909] Theorie der Transversalschwingungen einer quadratischen Platte mit freien Rändern
https://doi.org/10.1002/andp.19093330403

[Ritz1908]を読むと、Hilbertへの言及があって、Hilbertが当時研究していた変分法の直接解法を知っていて、この方法を考えたのだろう。Ritzは、一時期、ゲッティンゲン大学に在籍して、Hilbertに会ったことがあるっぽい。

1915年にGalerkinは、変分汎関数を経由することなく、(現在で言うところの)重み付き残差法によって係数を決定する方法を提案して、Ritzの計算法を、変分構造を持たない微分方程式でも使えるように拡張した。Galerkinのアルゴリズムは、何らかの汎関数を最小化する関数を探しているわけではないけど、変分問題に帰着できる場合には、Ritz法と同じ方程式に帰着する。[Ritz1908]の式(41)に、ひっそりと、Galerkin法から得られるのと同じ式が出ている。

(収束先が、どのような関数空間に含まれるのか気にしないならば)形式的には、任意の線形偏微分方程式の境界値問題を、Galerkin法で解くことができ、差分法や有限要素法と違って、(基底関数をどう選ぶかは考える必要があるけど)非有界領域でも使える。

Galerkinの論文は、原文はロシア語で書かれていて、英訳も存在するらしいのだけど、オンラインでは、英語版もロシア語版も公開されてないようなので、内容の確認はできなかった。

Rods and plates : series in some questions of elastic equilibrium of rods and plates
https://books.google.co.jp/books/about/Rods_and_Plates.html?id=gXLCHAAACAAJ

代わりとして、以下を参照した。
The history of Galerkin's method and its role in M. V. Keldysh's work
http://www.mathnet.ru/php/archive.phtml?wshow=paper&jrnid=ipmp&paperid=1929&option_lang=eng

Ritz法は、Rayleigh-Ritz法と呼ばれることがあり、Rayleigh自身が、Ritzと同じ方法を著書The theory of soundの中で書いていると主張する論文

On the calculation of Chladni's figures for a square plate
https://doi.org/10.1080/14786440808637121

を書いたために、こういう名前が付いてるっぽい。タイトルにあるChladni's figures(クラドニ図形で検索すれば、絵が沢山出てくる。数学的には、固有モードの節線)の計算は、[Ritz1909]で見られるもので、Germain-Lagrange方程式は、元々は、この図形を説明するために導入された方程式だった。Rayleighは、Ritzの数値解を、practically complete solution of the problem of Chladni's figures on square platesと評している一方、"As has been said, the general method of approximation is very skillfully applied, but I am surprised that Ritz should have regarded the method itself as new. "とも書いている

Rayleighが論文の中で触れている箇所(§168,175,228や§88, 89, 90, 91, 182, 209, 210, 265及び、vol ii appendix Aを見ろと書いている)を見ても、Ritzと同じ方法であると解するのは、難しいと思われる。例えば、

The theory of sound §182
https://archive.org/details/theorysound06raylgoog/page/n249/mode/2up

では、振動の形を仮定した上で、ポテンシャルエネルギーの最大値と運動エネルギーの最大値が等しいと置いて、振動数を決定している。

TimoshenkoのVariational Problems in Engineeringの初版は1928年に書かれたものであるけど、この方法の解説があり、Ritz's methodは、further development of Rayleigh's methodだと書いている。

cf)Vibration Problems In Engineering(第二版)370ページ §62はRitz methodというタイトルになっている
https://archive.org/details/vibrationproblem031611mbp/page/n385/mode/2up

ただ、この本来の"Rayleigh法"は、計算量が少ないので、昔は、よく使われたようである。一方、Ritz法は、計算量的な問題で、電子計算機の発展までは、それほど出番はなかったらしい。

The historical bases of the Rayleigh and Ritz methods
https://www.sciencedirect.com/science/article/abs/pii/S0022460X05000362

には、以下のように書かれている。

Although the Rayleigh method is used frequently, the Ritz method has found tremendous usage during the past three decades in obtaining accurate frequencies and mode shapes for the vibrations of continuous systems, especially for problems not amenable to exact solution of the differential equations. This is especially because of the increasing capability of digital computers to set up and solve the frequency determinants arising with the method. Even before that, the writer found 15 publications that used the Ritz method to solve classical rectangular plate vibration problems prior to 1966.

Galerkin近似も、Ritz近似と呼ぶ方が適切だと思うけど、Ritzさんは割と不遇である(1909年に死んでしまい、積極的に宣伝してくれる弟子とかがいなかったせいかもしれない)

同様の発想を、非線形偏微分方程式へ適用したのは、多分、数学者が最初っぽい。非圧縮性Navier-Stokes方程式の解が、適当な有限個の関数の線形結合で近似できるとして、非圧縮性Navier-Stokes方程式の弱形式に放り込むと、線形結合の係数は、時間に依存する関数で、この係数は非線形常微分方程式を満たす(この構成は、NS方程式の近似解を構成するアルゴリズムを与えていると見ることもできる)。この近似解の極限は、Leray-Hopfの弱解と呼ばれるものになる。

差分法であれ何であれ、第二次世界大戦前に、Navier-Stokes方程式を、直接数値計算で解いてみようと試みた人は、おそらくいなかったと思われる。そもそも、19世紀〜20世紀初頭の流体力学に於いては、ポテンシャル流の理論が中心にあって、Navier-Stokes方程式は、それほど基本的なポジションにあったわけではないと思われる。考えてみれば、古典力学での摩擦や、音波や電磁波の減衰は、しばしば無視されるか、後付で現象論的に扱われる方が、今でも普通なので、粘性がオマケ扱いでも一貫性はある。

このことは、Lambの教科書Hydrodynamicsの構成にも反映されている。この本は、1895年に初版が出版され、それから半世紀の間に、6版まで出たようであるけど、Viscosityは、かなり後ろの方のChapterで扱われている(初版だと、Chapter11)

Hydrodynamics
https://books.google.co.jp/books/about/Hydrodynamics.html?id=h-FJAAAAMAAJ

一応、ViscosityのChapterで、Navier-Stokes方程式は導かれているけど、"Navier-Stokes方程式"とは呼ばれておらず、普及した名称がなかったらしい(NavierとPoissonが導き、教科書で使用した導出法は、Saint VenantとStokesによるものだとは書かれている)ところを見ても、方程式自体が広く使われてたわけでもなかったのだろう。

Leray-Hopfの弱解を最初に導いたとされるLerayの1930年代の論文では、"equations de Navier"と書かれている。Google scholarで検索すると、1940年代には、Navier-Stokes方程式という名称が、ちらほら出現しているっぽいので、名称が固定したのも、この時期なのだろう。Landau-Lifschitzのfluid mechanicsは、初版が1959年のようだけど、Navier-Stokes方程式という名称が使われ、方程式の導出は、かなり最初のほうで行われている。

流体力学のことはともかく、以上を総括すると、結局のところ、電子計算機以前の偏微分方程式の数値計算は、proof of concept的な位置付けに留まっていたというのが妥当な評価だと思われる。

Ritzによる正方形板のクラドニ図形の計算

最後に、[Ritz1909]の結果を、論文に従って追ってみる。後から気付いたけど、2012年に書かれたレビュー論文

Chladni Figures and the Tacoma Bridge: Motivating PDE Eigenvalue Problems via Vibrating Plates
https://doi.org/10.1137/10081931X

には、Ritzの結果のトレース、Ritzの方法で、次元を大きくした時の結果、有限差分法との比較などが書かれていて、オープンアクセスになってるので、こっちを読んだ方がいい。

正方形板を考えるので、領域の座標を、 $(x,y) \in [-1,-1] \times [-1,1]$ とする。方程式や変分汎関数は、texを打つのが面倒なので省略。境界条件は、自由端とする。具体的な式は省略

試行関数について、Ritzは、以下のように取った。

自然数 $m$ が偶数の時は $\tan k_{m} + \tanh k_{m}=0$ で、奇数の時は $\tan k_m - \tanh k_{m}=0$ であるような正実数 $k_{m}$ を小さい方から順にとる。tex:k_{0}=k_{1}=0]で、Ritzの論文では、 $k_{2} = 2.3650$ としている。 $k_{m}$ が大きくなると、 $\tanh k_{m}$ は、ほぼ1になり、 $\tan k_{m} \approx \pm 1$ となるように周期的に取ればよくなるから、Ritzの論文では、 $k_{4}=\dfrac{7 \pi}{4},k_{3} = \dfrac{11 \pi}{4}$ と近似してるっぽい

同様に、 $k_{3}=3.92660$ で、 $k_{5} = \dfrac{9\pi}{4},k_{7}=\dfrac{13\pi}{4}$ など。

$u_{0}(x)=\dfrac{1}{\sqrt{2}} [tex:u_{1}(x)=\sqrt{\dfrac{3}{2}}x$
で、 $m$ が2以上の偶数の時
$u_{m}(x) = \dfrac{\cosh k_m \cos k_m x + \cos k_m \cosh k_m x}{\sqrt{\cosh^2 k_{m} + \cos^2 k_{m}}}$
で、自然数 $m$ が3以上の奇数の時は、
$u_{m}(x) = \dfrac{\sinh k_m \sin k_m x + \sin k_m \sinh k_m x}{\sqrt{\sinh^2 k_{m} - \sin^2 k_{m}}}$
として、基底関数を $u_{m}(x)u_{m}(y)$ という形に選ぶ。

論文の式(20)(21)などにある通り、関数 $u_{m}$ は
$\dfrac{d^4 u_{m}}{dx^4} = k_{m}^4 u_{m}$
$\dfrac{d^2 u_m}{dx^2} = \dfrac{d^3 u_{m}}{dx^3} = 0 (x=\pm 1)$
を満たす。係数は、正規直交基底となるように選ばれている。

論文では、 $m,n$ を、0から6まで取っているので、単純に考えると、49個の基底がある。現代なら49x49の行列の計算は、どうってことないけど、当時の計算力的には厳しい。実際に行列要素を計算してみると、多くの成分が0になり、単純な要素の置換によって、ブロック対角形になっていることが分かる。このように次元が落ちるのは、x,yの入れ替えに関する対称性や、x,ｙの偶奇性によって、分解できることに起因する。

一応、行列要素の計算は、以下のコードでできるはず(sympyの計算が、物凄く遅くて、数時間かかった)。論文に従って、Poisson比μ=0.225に固定しておく

import sympy
import numpy as np

x,y=sympy.symbols('x y')
k2,k3,k4,k5,k6 = 2.3650,3.92660,7*math.pi/4,9*math.pi/4,11*math.pi/4   #--values used by Ritz

u0 = (x**0)/math.sqrt(2)
u1 = math.sqrt(1.5)*x
u2 = (math.cosh(k2) * sympy.cos( k2*x ) + math.cos(k2)*sympy.cosh(k2*x)) / math.sqrt(math.cosh(k2)**2 + math.cos(k2)**2)
u3 = (math.sinh(k3) * sympy.sin( k3*x ) + math.sin(k3)*sympy.sinh(k3*x)) / math.sqrt(math.sinh(k3)**2 - math.sin(k3)**2)
u4 = (math.cosh(k4) * sympy.cos( k4*x ) + math.cos(k4)*sympy.cosh(k4*x)) / math.sqrt(math.cosh(k4)**2 + math.cos(k4)**2)
u5 = (math.sinh(k5) * sympy.sin( k5*x ) + math.sin(k5)*sympy.sinh(k5*x)) / math.sqrt(math.sinh(k5)**2 - math.sin(k5)**2)
u6 = (math.cosh(k6) * sympy.cos( k6*x ) + math.cos(k6)*sympy.cosh(k6*x)) / math.sqrt(math.cosh(k6)**2 + math.cos(k6)**2)


Basis = []
for u_m in [u0,u1,u2,u3,u4,u5,u6]:
  for u_n in [u0,u1,u2,u3,u4,u5,u6]:
     Basis.append( u_m*u_n.subs(x,y))
 

D=np.zeros( (len(Basis) , len(Basis)) )
mu = 0.225   #--Poisson ratio

for n,w1 in enumerate(Basis):
    for m,w2 in enumerate(Basis):
       V1 = sympy.diff(sympy.diff(w1,x),x)*sympy.diff(sympy.diff(w2,x),x)
       V2 = sympy.diff(sympy.diff(w1,y),y)*sympy.diff(sympy.diff(w2,y),y)
       V3 = 2*mu*sympy.diff(sympy.diff(w1,x),x)*sympy.diff(sympy.diff(w2,y),y)
       V4 = 2*(1-mu)*sympy.diff(sympy.diff(w1,x),y)*sympy.diff(sympy.diff(w2,x),y)
       D[n,m] = sympy.integrate(sympy.integrate(V1+V2+V3+V4 , (x,-1,1)), (y,-1,1))

最も次元が大きいブロックは、(m,n)=(0,2),(0,4),(0,6),(2,0),(2,2),(2,4),(2,6),(4,0),(4,2),(4,4),(4,6),(6,0),(6,2),(6,4),(6,6)の組み合わせで、15次元。これでも、手計算でやるには、少し次元が大きい。Ritzの論文では、 $u_{m}(x)u_{n}(y)\pm u_{n}(x)u_{m}(y)$ に基底を取り替えたりすることで、更に次元を落としている。

例として、論文の式(54)を見ると、6次元まで落として計算している。対角成分だけ比較すると、

最初は、18.95であるが、これは、D[8,8]=13.9499...で一致している
次は、411.8であるが、これと比べるのは、(D[22,22]+D[10,10]+D[22,10]+D[10,22])/2=406.1...で、ややずれている
3つ目は、1686となってるけど、D[24,24]=1684.4....なので、これも少しずれている
4つ目は、2945で、(D[36,36]+D[12,12]+D[36,12]+D[12,36])/2=2945.4...なので、まぁいい
5つ目は、6303で、(D[38,38]+D[26,26]+D[38,26]+D[26,38])/2=6325.4...なので、オーダーは合ってる
最後は、13674で、D[40,40]=13672.2...で、これも概ね合っている。

Ritzの論文の式(54)の行列は、以下の通り。当時は、補助的な計算ツールがあったとはいえ、基本は手計算なので、固有値を決めるのは割と大変だったと思われる。少なくとも、現代の線形代数のテストで、この行列の固有値を計算しろとか言われたら絶望する

matrix([[   13.95,   -32.08,    18.6 ,    32.08,   -37.2 ,    18.6 ],
        [  -16.04,   411.8 ,  -120.  ,  -133.6 ,   166.8 ,   140.  ],
        [   18.6 ,  -240.  ,  1686.  ,  -218.  , -1134.  ,   330.  ],
        [   16.04,  -133.6 ,   109.  ,  2945.  ,  -424.  ,   179.  ],
        [  -18.6 ,   166.8 ,  -567.  ,  -424.  ,  6303.  , -1437.  ],
        [   18.6 ,   280.  ,   330.  ,   358.  , -2874.  , 13674.  ]])

固有値を、コンピュータで数値計算すると、[12.4712661, 378.039661, 1582.08757, 2886.81821, 5943.14563, 14231.1877]となった。論文だと、最小固有値は、12.43だと書かれている。

この場合のRitzの計算結果は、論文末尾のA. Lösungen, die in x und y ungerade und symmetrisch sind.にまとめられていて、そこに計算されたクラドニ図形も載っている

一方、上のコードで得られた計算値から、同じ行列を計算すると、以下のようになった

matrix([[   13.95      ,   -32.21607048,    18.59991392,    32.21616838,  -37.19994089,    18.60002697],
        [  -16.10803524,   406.11892027,  -119.9245119 ,  -133.5507308 ,  172.08041016,   -83.55409784],
        [   18.59991392,  -239.84902379,  1684.45812633,   217.98225565,  -1136.62502909,   329.53550565],
        [   16.10808419,  -133.5507308 ,   108.99112782,  2945.46763315,  -427.12669121,   179.24264235],
        [  -18.59997045,   172.08041016,  -568.31251454,  -427.12669121,  6325.43274432, -1441.53670456],
        [   18.60002697,  -167.10819567,   329.53550565,   358.4852847 ,  -2883.07340912, 13672.21548558]])

固有値は、[12.4880464, 379.339436, 1559.26831, 2899.93081, 5961.30651, 14235.3098]となった。Ritzが何を使って計算したかは不明だけど、上位2桁くらいは合ってるっぽい

クラドニ図形の確認には、固有関数を見る必要があるけど、Ritzの結果をトレースできてそうなので、省略。また、それっぽい絵が出るというだけでは、モデルの妥当性の確認には不十分だけど、Ritzの結果のトレースとは趣旨が外れるので、それも省略

2020-06-14

調和代数幾何解析(調和解析と不変式論)

GoodmanとWallachの本

Symmetry, Representations, and Invariants
https://www.amazon.co.jp/dp/1441927298

の12章を眺めていて、ふと、ここに書いてあるのは、調和解析+代数幾何 = 調和代数幾何解析だなぁと思った(勝手に命名)。まぁ、この本の範囲では、"解析学"は、全然出てこないのだけど

【参考文献】同じく、GoodmanとWallachの本

Representations and Invariants of the Classical Groups
https://www.amazon.co.jp/Representations-Invariants-Encyclopedia-Mathematics-Applications/dp/0521663482

の12章あたりにも同等の内容があるっぽい。以下のPDFにも、同様の内容がある

An Algebraic Group Approach to Compact Symmetric Spaces
https://sites.math.rutgers.edu/~goodman/pub/symspace

Harmonic Analysis on Compact Symmetric Spaces : the Legacy of Élie Cartan and Hermann Weyl
https://sites.math.rutgers.edu/~goodman/pub/weyl_goodman.pdf

コンパクト対称空間 $G/K$ 上の二乗可積分関数の空間 $L^2(G/K)$ を $G$ 作用で既約分解するという問題は、Cartanが1929年に解いたらしい。この問題が、GoodmanとWallachの本の12章のそもそもの主題である。Cartanの議論には瑕疵があったらしく、1970年にHelgasonがギャップを埋めたということで、Cartan-Helgasonの定理と呼ばれていることがある。古典的な例として、 $L^2(S^1)$ の $SO(2)$ 作用による分解はFourier級数で、 $L^2(S^2)$ の $SO(3)$ 作用による分解は球面調和関数展開で記述できる。この定理は、対称空間を舞台とする調和解析の基本と言える。

Cartan-Helgasonの定理には、Peter-Weylの定理と同様に、解析的なバージョンと代数的なバージョンが存在する。

Peter-Weylの定理の場合は、コンパクトLie群 $G$ に対して、 $L^2(G)$ の既約分解を記述するのが解析的バージョンで、一方、コンパクトLie群には、実代数群の構造が一意に入る(※)ので、複素化して、群上のアフィン座標環を( $G$ の複素化群の作用で)既約分解するのが代数的バージョン。

※)代数群側から始めるなら、複素簡約(代数)群 $G_{\mathbf{C}} \subset GL(n,\mathbf{C})$ には、 $G=G_{\mathbf{C}} \cap U(n)$ として、コンパクト実型が取れると言ってもいい。

個人の好みは別として、Peter-Weylの定理の解析バージョンは、関数空間がHilbert空間だから積構造は捨てられているのに対して、代数バージョンの方は、可換環なので、よりrichな情報を持っている。球面調和関数の積とか議論することを考えると、代数的な方が、面倒事がない。

簡単な例として、 $SU(2)$ の場合を考える。 $SU(2)$ は、幾何学的には $S^3$ と同一視出来るので、Peter-Weylの定理の適用対象であると同時に、Cartan-Helgasonの定理の適用対象にもなる。Peter-Weylの定理で考える場合は、 $SU(2)$ の作用で分解し、Cartan-Helgasonの定理で考える場合は、例えば、 $SO(4)$ などの作用で既約分解する。 $L^2(SU(2))$ の基底として、Wigner D-matrixを取ることができ、 $L^2(S^3)$ の基底として、hyperspherical harmonicsを取ることができる。

$SU(2)$ を複素化した群は、 $SL(2,\mathbf{C})$ で、座標環は $\mathbf{C}[a,b,c,d]/(ad-bc-1)$ と書ける。それで、
$\left(\begin{array}{cc} a & b \\ c & d \end{array} \right) = \left( \begin{array}{cc} x_1 + i x_2 & -x_3 + i x_4 \\ x_3+i x_4 & x_1 - i x_2 \end{array} \right)$
と変数変換を行えば、
$\mathbf{C} [a,b,c,d]/(ad-bc-1)\simeq \mathbf{C}[x_1,x_2,x_3,x_4]/(x_1^2+x_2^2+x_3^2+x_4^2-1)$
となって、後者は、3次元球面上の座標環の複素化と思える。

$\mathbf{C}[a,b,c,d]/(ad-bc-1)$ は代数群の座標環なので、Peter-Weylの定理を考えることができる。一次の元は、 $a,b,c,d$ と4つあって、これは、二次元既約表現の行列要素を与える。二次の元は、 $a^2,b^2,c^2,d^2,ab,ac,ad,bc,bd,cd$ の10個あるけど、関係式 $ad-bc=1$ によって、独立な元の数は9個。適当な基底を取って、三次元既約表現の行列要素を見れば、これらの元が全部出ることが分かる。

GoodmanとWallachの本には、複素簡約群とその閉部分群(それぞれ $G,K$ とする)に対して、以下のような環が定義されている
$\mathcal{R}(G/K) = \{f \in \mathbf{C}[G] : f(gk) = f(g) \verb|, for all | k \in K,g \in G\}$
Goodman-Wallachは、 $K$ の自明表現の誘導表現として、この環を定義している。誘導表現という視点は、GやKが有限群であるような状況で議論を拡張するには有用であるけど、ここでは、 $\mathcal{R}(G/K)$ が $K$ 作用による不変式環であるという点に着目してみる。

また簡単な例として、 $\mathcal{R}(SL_{2}/GL_{1})$ を考える。 $\mathbf{C}[SL_{2}] = \mathbf{C}[a,b,c,d]/(ad-bc-1)$ と座標を取って、
$\left(\begin{array}{cc} a & b \\ c & d \end{array} \right) \left( \begin{array}{cc} t & 0 \\ 0 & t^{-1} \end{array} \right) = \left( \begin{array}{cc} at & bt^{-1} \\ ct & dt^{-1} \end{array} \right)$
なので、 $ab ,ad, bc, bd$ は、 $K=GL_{1}$ 作用で不変である。 $K=GL_{1}$ 作用で不変な元は、これらの元で生成されることは明らか。4つあるけど、条件 $ad-bc=1$ によって、独立な元は3つ以下。

更に $ad+bc , ab-cd , ab+cd$ に対して
$(ad+bc)^2 + (ab-cd)^2 - (ab+cd)^2 = (ad - bc)^2 = 1$
という関係式があって、独立な元は2つとなる。この二次の関係式と球面の関係は明らかだけど、上でやったのと同様に
$(a,b,c,d) = (x_1 + i x_2 , -x_3 + i x_4 , x_3+i x_4 , x_1 - i x_2)$
と変数変換を行うと
$ad+bc = x_1^2+x_2^2-x_3^2-x_4^2$
$ab-cd = 2(x_1x_3+x_2x_4)$
$ab+cd = 2i(x_1x_4-x_2x_3)$
なので、 $x_1,x_2,x_3,x_4$ を実数に制限すると、 $SU(2) \simeq S^3$ から $S^2$ へのHopf写像を与えていることが分かる

cf)Hopf fibration
https://vertexoperator.github.io/2018/11/20/elementary_Hopf_fibration.html

以上のような計算によって、
$\mathcal{R}(SL_{2}/GL_{1}) \simeq \mathbf{C}[y_0,y_1,y_2]/(y_0^2+y_1^2+y_2^2-1)$
が分かり、これは $SU(2)/U(1) \simeq S^2$ という同型の複素化に相当する。

【超余談】
三角形の合同条件と不変式論
https://m-a-o.hatenablog.com/entry/20130104/p1
の計算は、 $\mathcal{R}(GL_{2}/SO(2))$ を決定する問題になっている。三角形のモジュライ空間自体は、鏡映で移るものを区別するなら、 $GL(2,\mathbf{R})/SO(2)$ で、同一視する場合は、 $GL(2,\mathbf{R})/O(2)$ になり、いずれにせよコンパクトではないんだけど、複素化すると、上のHopf fibrationの計算と殆ど同じことになる。

$\mathbf{C}[GL_{2}] \simeq \mathbf{C}[a,b,c,d,1/(ad-bc)]$ で、 $ad-bc$ は、 $SO(2)$ 作用で不変(三角形のモジュライという観点からは、面積が回転不変ということ)なので、結局、主要な仕事は、 $\mathbf{C}[a,b,c,d]^{SO_2}$ の決定になる。この不変式環の生成元は、 $a^2+b^2,c^2+d^2,ac+bd$ となり、 $(a^2+b^2)(c^2+d^2) = (ac+bd)^2 +(ac-bd)^2$ という関係式が成立する。

ところで、平面上の三角形は平面上の並進自由度も持つことを考えると、三角形のモジュライ空間は、 $GA(2)/SE(2)$ あるいは、 $GA(2)/E(2)$ でもある( $(GA(n)$ はアフィン変換群で、 $SE(n),E(n)$ はEuclid変換群)けど、 $GA(n)/SE(n) \simeq GL(n,\mathbf{R})/SO(n)$ 、 $GA(n)/E(n) \simeq GL(n,\mathbf{R})/O(n)$ なので、簡約群の計算で済む。Riemann計量は、各点で $GL(n,\mathbf{R})/O(n)$ に値を取る場であるが、 $GA(n)/E(n)$ に値を取る場と思うこともできる。一般相対論の場合、直交群はローレンツ群に、Euclid群はPoinare群に変更すれば、同じ同型があって、Einsteinによる定式化は前者に基づくけど、後者に基づく理論は、Einstein-Cartan理論とか呼ばれる。このCartanは、Cartan-Helgasonの定理のCartanと同じ人。

複素簡約群GとKに対して、G/Kがアフィン代数多様体であるというのは、松島の定理の簡単な部分である。考えてみると、コンパクト対称空間を複素化すると、アフィンになるというのは、ちょっと面白いことだと思う。複素射影空間や複素グラスマン多様体を、複素化するというのは、ぱっと見奇妙な感じがしないでもないけど、仕方ない。

【松島の定理】松島の定理は、複素簡約群Gに対して、G/Kがアフィンであることと、Kが複素簡約群であることが同値という定理。この定理は、1960年に発表され、色んな証明があるようである。2000年以降でも、以下の論文を見つけた
Invariant Ideals and Matsushima's Criterion
https://arxiv.org/abs/math/0506430

The Kempf-Ness theorem and Invariant Theory
https://arxiv.org/abs/math/0605756

別に定理を持ちださなくても、その座標環の計算は、上で見たように、不変式の計算で原理的には出来る。 $SL_{n}/SO_{n}$ や $SL_{2n}/Sp_{n}$ のような場合は、極分解や、そのvariationで、行列の分解に帰着して計算できるけど、一般的には、それだけで済むのか分からない。21世紀になって、DerksenやKemperによって、具体的なアルゴリズムも与えられている(※)。群が少し大きくなると、計算が終わらないと思うけど。

※)計算不変式論の使い方:合同変換群の多項式不変量の計算事例
https://m-a-o.hatenablog.com/entry/20131227/p2

【行列の分解による不変式の計算】例えば、non-oriented Grassmann多様体の場合でも、極分解と似たような方法で、不変式を構成できる。実Grassmann多様体 $O(p+q)/(O(p)\times O(q))$ を例に取ると、行列をブロック化して
$\begin{pmatrix} A & B\\ C & D \end{pmatrix} \begin{pmatrix} P & 0 \\ 0 & Q \end{pmatrix} = \begin{pmatrix} AP & BQ \\ CP & DQ \end{pmatrix}$
なので、 $AA^{T}$ や $CA^{T},CC^{T},\cdots$ なんかの行列成分に不変式が出る。問題は、これらで、不変式の生成元が尽きるのかということだけど、ちゃんと数えてないので不明。まぁ、何となく、他にないやろうという気はするけど

そういうわけで、 $\mathcal{R}(G/K)$ は、自明でない環構造を持ってるのだけど、Cartan-Helgasonの定理の証明には、環構造を具体的に知る必要は全然ないし、環構造があることすら重要ではない。これは、Cartan-Helgasonの定理の解析的バージョンが $L^2$ 関数空間の分解として書かれることを思えば、当然ではある。GoodmanとWallachの本にある戦略は、 $\mathbf{C}[G]$ に対する(代数的)Peter-Weylの定理から、 $G$ の有限次元既約表現を $K$ に制限して、 $K$ の一次元表現が何個出るか見ればいいよね的なものである

まぁしかし、せっかく環構造があるのだから、もうちょっと、そこに注目してみてもいいのじゃないかと漠然と思う。別に、それで何か新しい物が出てくるとか強気に主張できないのが、苦しいとこではあるけれど。

逆に、Cartan-Helgasonの定理による既約分解から、何とかして環構造を復元できるかということを考えると、よく分からない。球面調和関数の空間は、対称トレースレステンソルの空間と同一視できるけど、EastwoodのCartan productという(結合的かつ可換な)積構造を入れることができ、この環は、(複素化した)球面上の座標環とは異なる。

The Cartan Product
https://doi.org/10.36045/bbms/1110205624

尤も、この場合、違いは、 $r^2=1$ とするか、 $r^2=0$ とするかなので、この2つの環は一方から他方へ連続変形できる。他の場合は、どうなってるのかというのは謎

ところで、Cartanは、Schoutenと共同で、絶対平行性を持つRiemann多様体の分類を調べる中で、曲率テンソルの共変微分が0になる空間(今では局所対称空間と呼ばれる)に関心を持って、対称空間を調べ始めたようだ。例えば、定曲率空間は、局所対称空間なので、定曲率空間を大幅に一般化したと思うことができる。Schoutenとの論文が1926年出版で、この時点では調和解析の片鱗もない。

元々の動機は、今となっては割とどうでもいいものになってしまったけれど、その後、調和解析をやるだけなら、もうちょっと条件の緩い空間でもいいんじゃねと思ったのか何なのか、20世紀中盤〜後半にかけて、
・弱対称空間(Selberg 1956)
・Gelfand pair(Gelfand 1960s?)
・spherical pair(Vinberg&Kimel'feld 1978, Krämer 1979)
・spherical homogeneous space(Brion 1986)
などの概念が現れた。これらは同じようなもので、別に、互いの仕事を認識してなかったというわけでもないようだけど、なんか色々な名前が付いててめんどくさい。

Gelfand pairの心は、不変微分作用素環が可換というものだと思う。Selbergは、弱対称空間が、この性質を持つことを示したそうだ(論文は見てないけど)。Selbergの論文は、Poisson和公式の一般化/Selberg trace formulaにあったそうだ。多分、Selbergの仕事が先っぽいけど、Gelfand pairは、有限群とかも対象になる点で、より一般的なものではある。このへんは、解析学らしさが感じられる話。

一方、spherical pairやspherical homogeneous spaceとかは、なんか関数空間を、G作用で既約分解すると、multiplicity-freeになるという感じのもので、Cartan-Helgasonの定理から、対称空間が、この性質を持つことが分かる。

spherical homogeneous spaceという名前を、最初に使ったのは、Brionという人だと思われる。
Quelques proprietes des espaces homogenes spheriques
https://link.springer.com/article/10.1007/BF01168684

もうちょっと一般化したspherical varietyという用語も、この人が命名したのが、多分、最初っぽい
Spherical Varieties an Introduction
https://link.springer.com/chapter/10.1007/978-1-4612-3702-0_3

私は全然知らなかったけど、spherical varietyの研究は(論文が沢山出てくるので)結構、流行ってるっぽい。このへんの話は、調和解析も表現論も離れて、また違う動機の元で進んでるように見えるので、よく分からない。

で、不変微分作用素環の可換性と、関数空間をG作用で既約分解したらmultiplicity-freeになるのが同値という事実は、色んなバリエーションがある。今の話に関係しそうなとこだと
Multiplicity-free spaces
https://doi.org/10.4310/jdg/1214438422

Compact weakly symmetric spaces and spherical pairs
https://arxiv.org/abs/math/9808039

Weakly symmetric spaces and spherical varieties
https://link.springer.com/article/10.1007/BF01236659

とりあえず、不変式環を見ると、(原理的には)計算できる！と思って安心する

2020-01-12

Siegelモジュラー形式環と計算機代数

asir

Rungeという人の1993年の論文

On Siegel modular forms. Part I.
https://doi.org/10.1515/crll.1993.436.57

に、Siegelモジュラー形式環の構造が、有限群の不変式の計算と、環の商体内での正規化(整閉包)によって計算できるという結果が書いてある(Theorem3.12)。有限群が作用する環は、テータ関数(テータ零値)の多項式たちで生成され、単純な多項式環というわけではない。偶数ウェイトのSiegelモジュラー形式環は、もう少し、簡単な計算でいける(Corollary3.17)。

Siegelモジュラー形式は1930年代終わり頃、Siegelが導入した。Siegelモジュラー形式に関する概略が、以下の論説の1〜2節などに書かれてる。
コンパクト化の今昔
https://www.jstage.jst.go.jp/article/sugaku1947/51/2/51_2_129/_article/-char/ja/

歴史的に見ると、(19世紀頃の)初期のモジュラー形式論は、個々のモジュラー形式を理解しようという特殊関数論的(というか古典解析的)な動機が強かったように見える。一方で、Siegelモジュラー形式環を調べる動機の一つは、離散群による商空間という扱い難い対象を理解しようというところにあって、解析的側面は薄れてて、代数幾何学的な話になっている。

【モジュライ空間史】On the early history of moduli and Teichmüller spaces
https://arxiv.org/abs/1602.07208
を読んだところ、Siegelは、1935年の論文で、Torelli空間とSiegel上半平面、Siegel基本領域を定義し、Torelli空間から、Siegel上半平面への複素構造を保つ写像を定義したということが書いてあった。Torelli空間は、おそらく、Weilによる命名で、今では、タイヒミュラー空間のTorelli群による商と定義されることが多い。タイヒミュラー空間の論文は、この少し後に書かれたらしいけど、TeichmüllerがSiegelの仕事を知ってたのかは謎。論文は読んでたとしても、当時、Teichmüllerは、ナチの反ユダヤ主義を支持し、Siegelは、ユダヤ人だったのかは知らないけど、1940年にアメリカに亡命してるので、仲良く交流したとは思えない。それより前の時代だと、Schottkyなんかは、純粋に解析的なアプローチを取っていたようだし、1913年にTorelliの定理を示したTorelliも、周期行列の集合に良い幾何構造が入るという認識は持ってなかったっぽい。1920年代に、対称空間の分類とかが行われてて、Siegelに、その影響があったかは不明。モジュライ空間を厳密に定式化するという点で、Siegelは大きな貢献があったようである

私は、コンパクト化の議論とか読んでも、そんなんええから、Siegelモジュラー形式環の生成元と関係式を書いてくれへんかなという気持ちになる。現在、種数3までのSiegelモジュラー形式環の構造が決定できているっぽい。種数1のSiegelモジュラー形式環は、19世紀には、既に知られていた通り、4次と6次のEisenstein級数で生成される。この環と、有限群の不変式環との間に自然な同型が存在するというのは、誰が最初に気付いたのか分からないけど、有限群の不変式環とある種のモジュラー形式環に自然な同型があるという現象は、1970年頃には知られてたっぽい(もっと古い起源があるのかは不明)。

例えば、1972年のフランス語の論文

Polynômes des poids de certains codes et fonctions thêta de certains réseaux
http://www.numdam.org/item/ASENS_1972_4_5_1_157_0/

には、位数192の群の不変式環が、 $\mathbf{C}[E_4,\Delta]$ と自然に対応するということが書かれてるように見える(フランス語は読めない)。

そういう話とは関係なく、Siegelモジュラー形式環を頑張って決めようとした人もおり、種数2のSiegelモジュラー形式環は、井草(1962)
On Siegel modular forms of genus two
https://www.jstor.org/stable/2372812

種数3のSiegelモジュラー形式環は、露峰(1986)の結果がある。
On Siegel Modular Forms of Degree Three
https://www.jstor.org/stable/2374517?seq=1

種数3の場合は、計算が大変なようで、生成元がminimal setかどうか分かってなかったらしいけど、2019年に、minimal setを決めたと主張している論文が出ているよう(内容は未精査)。19個の生成元と55個の関係式があるらしいので、結果をトレースするだけでも、手計算は辛そう

Siegel modular forms of degree three and invariants of ternary quartics
https://arxiv.org/abs/1907.07431

Rungeの結果を使って、モジュラー形式環を決定するには、テータ零値の間に存在する代数的関係式を決める必要があるけど、一般的な方法は知られてないっぽい。 $m=(m',m'')\in \{0,1\}^g \times \{0,1\}^g$ に対して
$\theta_m(z,\tau) = \displaystyle \sum_{x \in \mathbf{Z}^g} \exp 2\pi i\left(\dfrac{1}{2} {}^{t}[x+\dfrac{1}{2}m'] \tau [x+\dfrac{1}{2}m'] + \langle x+ \dfrac{1}{2}m',z+\dfrac{1}{2}m'' \rangle \right)$
とすると、zについて偶関数であるものが、 $2^{g-1}(2^g+1)$ 個ある。z=0の時しか考えないので、oddな方は、どうでもいい。Rungeの論文では、テータ関数の二乗が、
$f_a(\tau) = \theta_{(a,0)}(0,2\tau) = \displaystyle \sum_{x \in \mathbf{z}^g} \exp 2\pi i \left( {}^{t}[x+\dfrac{a}{2}] \tau[x+\dfrac{a}{2}] \right)$
を使って書けることから、この形のテータ零値の作る環を中心に考える。この形のテータ零値は、どれも0でなく、 $2^g$ 個の線形独立な関数となっている(Rungeの論文p59)

テータ関数やテータ零値には、沢山の恒等式があって、関係式の間の関係式が必要なほど、ややこしい。必要なテータ零値の間の関係式も、知られている関係式から得られるものもある。Rungeの論文§6では、種数3の場合に、既知のテータ関係式から、 $f_a$ の関係式を得て、Siegelモジュラー形式環の構造を決定している(Theorem6.2)

cf)A remark on a theorem of Runge
https://link.springer.com/article/10.1007/s000130050220
[PDF]https://www.mathi.uni-heidelberg.de/~freitag/preprints/runge.pdf
の式(2)も参照

もっと高い種数では、既知のテータ恒等式に帰着できない関係式が出るのか、よく分からない。

A theta relation in genus 4
https://projecteuclid.org/euclid.nmj/1114631553
では、種数4の場合に、一つの関係式を見つけて、"It is a natural question, whether our relationis a consequenceof Riemann's theta relations. We could not decide this."と書いてる(2001年付けだから、少し前だけど、凄く古いわけではない)。全部が既知のテータ恒等式に帰着できるなら、原理的には、変数消去の計算で、全ての関係式を見つけることができる。実際上は、evenなテータ関数が $2^{g-1}(2^g+1)$ 個あって、(τではなく)2τの関数である $f_a$ が、 $2^g$ 個あるので、単純に考えると、g=4で152変数、g=5で560変数の多項式の変数消去計算をやらないといけないことになって、現在のコンピュータとアルゴリズムで計算できるかは怪しいけど

【参考文献】arXiv:1401.5368やarXiv:1510.02699には、色々なテータ恒等式の発見に、誰が寄与したかという話が少し書いてある

【テータ恒等式とPlucker関係式】2x4行列の2x2小行列式 $D_{ij} (1 \leq i \lt j \leq 4)$ で生成される環は、多項式環を、一個のPlucker関係式で生成されるイデアルで割ったものになる。
$\mathbf{C}[D_{12},D_{13},D_{14},D_{23},D_{24},D_{34}] \simeq \mathbf{C}[X_{12},X_{13},X_{14},X_{23},X_{24},X_{34}]/(X_{12}X_{34}-X_{13}X_{24}+X_{14}X_{23})$
小行列式 $D_{ij}$ に対して、4x4対角行列 $h \in H_4 \subset GL(4)$ の作用を
$\rho(h)(D_{ij})(x) =D_{ij}(xh)$
で定める。 $\rho(h)(D_{12}D_{34}) = det(h) D_{12}D_{34}$ etc.などが成立するけど、 $h \in H_4$ の作用で、 $det(h)$ の多項式倍だけ変化する元の全体は環をなし、
$\mathbf{C}[D_{12}D_{34} , D_{13}D_{24} , D_{14}D_{23}] \simeq \mathbf{C}[y_0,y_1,y_2]/(y_0-y_1+y_2)$
のようになる。一般の群と指標に対して、この手の環は、半不変式環と呼ばれてることがある。一方、種数1のテータ零値には
$\theta_{00}^4 = \theta_{10}^4 + \theta_{01}^4$
という関係式があるので、環同型 $\mathbf{C}[\theta_{00}^4,\theta_{10}^4,\theta_{01}^4] \simeq \mathbf{C}[D_{12}D_{34} , D_{13}D_{24} , D_{14}D_{23}]$ が成立する。

代数的に書くと、大した意味はなさそうだけど、これは偶然でなく、19世紀から知られる超幾何関数と楕円関数の不思議な繋がりが背景にある。 $\mathbf{C}[\theta_{00}^4,\theta_{01}^4,\theta_{10}^4]$ は、モジュラー群 $\Gamma(2)$ のモジュラー形式環 $M(\Gamma(2))$ になっている。このモジュラー形式環は、楕円曲線の単位元以外の等分点 $e_1,e_2,e_3$ を使うと、
$M(\Gamma(2)) = \mathbf{C}[e_1 , e_2 , e_3]$
のようにも書ける。テータ零値で書くと
$e_1=-\dfrac{1}{3}(\theta_{00}^4 + \theta_{10}^4)$
$e_2=\dfrac{1}{3}(\theta_{00}^4 + \theta_{01}^4)$
$e_3=-\dfrac{1}{3}(\theta_{01}^4 - \theta_{10}^4)$
で、 $e_1+e_2+e_3=0$ が成立する。 $SL(2,\mathbf{Z})/\Gamma(2) \simeq S_3$ は、 $e_1,e_2,e_3$ の置換として作用し、 $M(SL(2,\mathbf{Z})) \simeq M(\Gamma(2))^{S_3}$ で、 $M(SL(2,\mathbf{Z}))$ の生成元として、 $e_1e_2+e_2e_3+e_3e_1$ と $e_1e_2e_3$ が取れて、それぞれ、定数倍を除いてEisenstein級数と一致する( $E_4=-3(e_1e_2+e_2e_3+e_3e_1) , E_6 = \dfrac{27}{2}e_1e_2e_3$ )(逆に、 $M(\Gamma(2))$ は $M(SL(2,\mathbf{Z}))$ 上の自由加群になってて、この構造を具体的に書いてみるのも面白い)。一方、 $\mathbf{C}[D_{12},D_{13},D_{14},D_{23},D_{24},D_{34}]$ のProjは、Grassmann多様体Gr(2,4)で、これを $H_4$ で"簡約"した空間が、Gauss超幾何関数の住む空間と思える。

"私説超幾何関数"(ISBN:978-4-320-01576-0)という本を読むと、種数2でも類似の同型が作れると書いてある。今度は、Grassmann多様体Gr(3,6)を考える。3x6行列の3x3小行列式 $D_{ijk}(1 \leq i < j < k \leq 6)$ は、20個あって、代数的トーラスのdetに関する半不変式環Z(3,6)は、10個の生成元と5個の独立な関係式からなる。種数2のテータ零値は10個あるけど、テータ零値の2次式を10個うまく選ぶと、それらの二乗から生成される環(適当なモジュラー群のモジュラー形式環となるらしい)と半不変式環Z(3,6)とが同型になるようだ。種数3以上に一般化する話は知られてないっぽいし、根源的な理由も分からないので、モジュラー形式というより、テータ零値の恒等式(の一部)に、何か良い構造があるという可能性もあるんじゃないかと勝手に思ってる

それ以外のステップ。有限群の不変式環を計算するアルゴリズムは以前に実装した

有限群の不変式環の生成元を計算するアルゴリズム
 https://m-a-o.hatenablog.com/entry/20130316/p2

asir_misc/finvar.rr
https://github.com/vertexoperator/asir_misc/blob/master/finvar.rr

にも、Risa/Asirで実装したものが置いてある。しかしまぁ、過去の私の記述によれば、"有限群の不変式の計算は、SingularやMagamaにも実装されてるらしいので、計算するだけなら、それ使うのが一番確実かもしれない"らしいので、Singular使えばいい気がする。

正規化の方は、Risa/Asirには、ないので、使いたければ自分で実装するしかない(と思う)けど、Singularには、既に実装されている

D.4.18 normal_lib
https://www.singular.uni-kl.de/Manual/4-0-3/sing_1269.htm#SEC1344

つまり、総合的に見て、Singular使えばいいのでは。一応、Singularで使ってる正規化アルゴリズムは、以下の論文に書かれているようだ。

Normalization of Rings
https://arxiv.org/abs/0904.3561

以前測定した限りでは、グレブナー基底の計算は、Singularの方が、Risa/Asirより、速い(数倍くらい違うこともあった気がする)。ただ、Risa/Asirは、何故か、代数的数のサポートが割と手厚い気がするので捨てがたい。あと、Singularは、(Mac以外)インストールめんどくさい

種数1の場合は、容易に計算できるし、Rungeの論文では、Remark3.6の直後に、結果が書いてある。種数1の場合は、群がShephard-Toddの(既約擬鏡映)群No.8というものになる

cf)List of irreducible complex reflection groups
https://en.wikipedia.org/wiki/Complex_reflection_group#List_of_irreducible_complex_reflection_groups

【用語について】正方行列Aが擬鏡映とは、Aが有限位数で、rank(I-A)が1であることを指す。英語だと、pseudo-reflection。位数が2の時は、ある超平面に関する鏡映になる。同じ概念に対して、complex reflection,unitary reflectionとかいう呼び方をされることもある。擬鏡映は、条件から、1でない固有値を一つだけ持ち、それは、１のべき根になってるから、こういう名前が付いてるのだと思う。+1でない固有値が実であれば、-1しかありえず、Aは鏡映になり、普通の鏡映は、実鏡映ということになる。個人的には、擬鏡映という名前が一番わかり易い。擬鏡映で生成される群が、擬鏡映群

私の実装でも

Z2 = newalg(x^2+1);
T=newmat(2,2,[[1,1],[1,-1]])*(1+Z2)/2;
D=newmat(2,2,[[1,0],[0,Z2]]);

finvar(grp([T,D]),[f0,f1]);

とかやって、一瞬で正しい答えが出た。
$f_0 = \theta_{00}(2\tau) , f_1 = \theta_{10}(2\tau)$
で、
$f_0^2 = \dfrac{\theta_{00}^2 + \theta_{01}^2}{2}$
$f_1^2 = \dfrac{\theta_{00}^2 - \theta_{01}^2}{2}$
が成立し、 2つの独立な不変式は
$g_4 = f_0^8+14f_1^4 f_0^4+f_1^8$
$g_6 = f_0^{12} - 33f_0^4f_1^4(f_0^4+f_1^4)+ f_1^{12}$
となる。Eisenstein級数は
$E_4 = \dfrac{1}{2} ( \theta_{00}^8 + \theta_{10}^8 + \theta_{01}^8)$
$E_6 = \dfrac{1}{2} (-3\theta_{10}^8(\theta_{00}^4+\theta_{01}^4) + \theta_{00}^{12} + \theta_{01}^{12})$
なので、両者が一致することは、単なる計算問題

種数2の場合。群 $H_2$ の生成元は、整係数Symplectic群の生成元から来ていて、Jから来る方はいいけど、整係数対称行列は、無限個ある。しかし
$S = \left( \begin{matrix} a & b \\ b & c \end{matrix} \right)$
と置くと
$D_S = diag(1 , i^c , i^a , i^{a+2b+c})$
になって、a,b,cは、0~3の範囲を舐めれば十分。有限個の生成元が決まったので、群の要素を列挙しようと思ったけど、pythonのsetに相当するものがないRisa/Asirの実装では、まる一日かけても、46080個の元を列挙するのは無理だった。仕方ないので、pythonで列挙して、ファイルから群の要素を読み込むという方法で対応(計算時間は、夜寝ている間に終わる程度)。

Rungeの論文には、これは擬鏡映群だと書いてあるけど、群に含まれる擬鏡映をコード書いて列挙してみると60個ある。位数は全て2で、どうやら、Shephard-Toddの群No.31っぽい(ちゃんとチェックはしてないけど、擬鏡映を列挙できてるので、比較すれば、力ずくで確認できるはず)。Molienの公式でHilbert級数を計算する(これはRisa/Asirで実装してある。昼寝してたら終わってた)と、
$\dfrac{1}{(1-t^8)(1-t^{12})(1-t^{20})(1-t^{24})}$
になっている。なので、この群を、4変数多項式環に作用させた不変式環は、8,12,20,24次の代数的に独立な不変式で生成される(8*12*20*24=46080で、群の位数が得られる)。具体的な式は、Rungeの論文p77冒頭に記載されているけど、Reynolds作用素の計算を、コンピュータにやってもらえばできる。論文にも書いてあるけど、20次と24次の不変式の空間は、一次元ではないので、定数倍を除いても、論文と同一の式にならない可能性はある。というか、私はならなかった。

モジュラー形式としては、ウェイト4,6,10,12に対応し、(Corollary3.17によって)種数2、偶数ウェイトのSiegelモジュラー形式環は、この4つの元で生成される。奇数ウェイトのモジュラー形式として、ウェイト35のものがあり、この5つの元が、種数2のSiegelモジュラー形式環を生成する

種数3の場合は、Rungeの群の位数が約3.7億なので、全部要素を列挙して計算するというのは望みがない。多分、私のPCでは、メモリに載らない。詳細は、
On Siegel modular forms. II
https://projecteuclid.org/euclid.nmj/1118775400
の方にあって、群の構造を調べて何とかしてるようだ。

現在のアルゴリズムは、Molienの公式やReynolds作用素の直接計算に依存してるので、位数の大きな群に対して使うには、どうしても限界がある。けど、種数が大きくなるごとに試行錯誤してたらキリがないので、超でかい有限群の不変式と戦えるアルゴリズムが欲しい。次元が上がると、今度はグレブナー基底の計算が辛くなるという二重苦が待っている

2020-01-12

解析力学から見た古典/量子Kepler系の解法

古典Kepler系の解析解
https://vertexoperator.github.io/2016/12/01/kepler2d_classical.html

2次元量子Kepler系の束縛状態
https://vertexoperator.github.io/2016/12/01/kepler2d_quantum.html

で、2次元のKepler系を、古典力学と量子力学で解いている。多少見かけは違ってるけど、同じ方法で空間変数の取り直しを行って、調和振動子に帰着(全く同じなので、書いてないけど、正エネルギー/散乱状態の場合は、inverted harmonic oscillatorというものに帰着する)、周波数の形も古典/量子で同じ形になる

ところで、古典Kepler系の方は、時間変数も取り直しているのに対して、量子Kepler系では、時間変数の変換は出てきてない。量子系の方は、時間非依存のシュレディンガー方程式を考えているのだから当然だけど、量子Kepler系のHamiltonianに対して、時間依存シュレディンガー方程式を考えて、時間変数の取り換えを行ってから、(新しい時間変数に対する)定常シュレディンガー方程式を考えると、結局、Hamiltonianは、左からrを掛けるという変換を受けることになる。

"量子Kepler系の束縛状態"の解法2の方で、「量子Kepler系のシュレディンガー方程式を
$r(\hat{H} - E)\Psi = 0$
と置くと・・・」と書いてあるのは、自明な式変形を書いてるのではなく、シュレディンガー方程式をこのように見ることが、古典系に於ける時間変数の取り換えと同等の操作を行うことに相当している。また、spectrum generating algebraに含まれるのは、Hamiltonianそのものではなく、 $r(\hat{H} - E)$ であることからも、この変換の必要性が理解できる。

ただ、これだけだと、ややad hocな印象もある

Duru-Kleinertによる、量子Kepler系の経路積分に関する以下の論文には、"pseudo Hamiltonian"というものが書いてある(論文中の式(10))。彼らは、これ以前にも同じような論文を書いてるけど、そっちには、このような説明はないので全然参考にならない

Quantum Mechanics of H-Atom from Path Integrals
https://doi.org/10.1002/prop.19820300802
[PDF]http://users.physik.fu-berlin.de/~kleinert/83/83.pdf

$r(\hat{H} -E)$ は、Duru-Kleinertが"pseudo Hamiltonian"と呼んでいるものに相当する。

で、Duru-Kleinertの論文の式(15)を見ると、このpseudo Hamiltonianが住んでるのは、拡大相空間であることが読み取れ、pseudo Hamiltonianに対して正準方程式を考えると、普通の運動方程式+時間変数の取り換えの式が出ている。Duru-Kleinertは、拡大相空間とは言ってないけど、昔の私が調べた所によると、1949年にLanczosが導入した比較的新しい概念らしい。昔の私が書いてることを信じるなら、H-Eはextended Hamiltonianと呼ばれることもあるそう

cf)extended phase space
https://formalgroup.tumblr.com/post/105356647155/extended-phase-space

古典力学では、H-E=0は定義みたいなものであり、量子力学では、H-Eを波動関数に作用させると0になるというのは、シュレディンガー方程式そのもの。

拡大相空間を考えるなら、(明らかに、拡大相空間も、自然なsymplectic構造を持ってるので)正準変換も"拡大"されるべきと思うのだけど、

Stackel系の全ての保存量を保つ離散化
https://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/42747

とかいう関係なさそうな文献を見ると、そのような変換を考えた人はいるらしい。Duru-Kleinertのpseudo Hamiltonianは、extended Hamiltonianに、一般化した正準変換を施したものであると理解できる。

DuruとKleinertは、経路積分のために、pseudo Hamiltonianを考えたけど、経路積分自体は、どうでもよく、背景に、こういう構造があることが重要なのだと思う

古典力学では、ハミルトニアンは関数なので、右から関数を掛けるか左から関数を掛けるか気にする必要はないけど、量子力学では、事情が違う。Kepler系では、 $r(\hat{H}-E)$ としたのだけど、 $(\hat{H}-E)r$ でない理由があるだろうかという疑問が残っている。一般には
$f_{l}(\hat{H}-E)f_{r}$
という形の変換が考えられる。この形の変換を取り扱った文献で、一番古いのは、調べた限り

On the generalization of the Duru-Kleinert-propagator transformations
https://link.springer.com/article/10.1007/BF01318171

じゃないかと思う。

Kleinertの本が、Kleinert自身によって公開されているけど、

New Path Integral Formula forSingular Potentials
http://users.physik.fu-berlin.de/~kleinert/b5/psfiles/pthic12.pdf

この形の変換が扱われている。式(12.26)とか。本を読めばわかるけど、新しい時間変数と、時間変数の関係は、 $f_l$ と $f_r$ の積のみで決まる。式(12.37)など。

この形の変換で一番良く見るのは、相似変換で、その場合は、時間変数を変える必要はない。相似変換については、当たり前すぎて、注意を払われることもないけど、それより、ずっと一般的な変換を行っても、ツケを時間変数に押し付けてしまえるというのは、言われれば当たり前だけど、何かに使えることもあるかもしれない。Kepler系の場合は、rを掛けることによって、Coulombポテンシャルの特異性が見かけ上消えるという見方もできる

【おまけ】Kleinertは、1967年に、A.O. Barutと共に
Transition Probabilities of the Hydrogen Atom from Noncompact Dynamical Groups
https://doi.org/10.1103/PhysRev.156.1541

という論文を書いていて、3次元量子Kepler系の束縛状態の空間が、so(4,2)の既約表現空間になってることを、(私が知る限り)初めて明確に書いた人でもある。この論文には、まだ"pseudo Hamiltonian"は書かれていない。もう少し後のBarutの論文では、Hamiltonianそのものより、"pseudo Hamiltonian"(とは呼んでないけど)を重視する書き方をしているものがあったと記憶している

Group Theory and Orbital Fluctuations of the Hydrogen Atom
http://users.physik.fu-berlin.de/~kleinert/kleiner_re222/barut92.html

のAbstractで、Kleinertは、以下のように書いている。

The group theoretic development was triggered by A.O. Barut who suggested to me the search for a dynamical group larger than SO(4). In this way he became partly responsible also for the later and most recent path integral development.

2020-01-12

Chapman-Jouguet理論

物理で定義とか真面目に考え出すと、辛いことになるけど、"法線方向の速度が不連続に変化する波面の伝播"を衝撃波と総称する立場からすると、発熱反応を伴う圧縮衝撃波がデトネーション波ということになる。発熱反応を伴う衝撃波では、膨張波も存在することができ、こっちは、デフラグレーション波と呼ぶ。日本語では、デトネーションは爆轟と訳され、デフラグレーションは爆燃と訳されている。

※)接線方向の速度が不連続に変化する不連続面は、渦層と呼ばれ、衝撃波面とは異なる。同じ不連続面とはいえ、使われる手法とかも全然別のよう。

普通の衝撃波は、化学反応とかは伴わないので、そのへんがデトネーション/デフラグレーションとの違い。デフラグレーションは、不連続波面を持つけど、伝播速度は音速より遅いので、衝撃波のイメージに、そぐわない。まぁ、あまり深く考えず、気相爆轟=超音速で伝わるガス爆発くらいの認識でいいっぽい。

気相爆轟は、19世紀の終わり頃から研究されるようになり、理論的な研究は、Mikelson(1890),Chapman(1899),Jouguet(1905)などが独立に行ったらしい。これらは、同じ結果を与えるらしく、現在では、Chapman-Jouguet理論(以下、CJ理論)として知られる。現在、CJ理論は、衝撃波のRankine-Hugoniot理論を拡張した形で展開される。

どうでもいいけど、CJ理論のChapmanとChapman-Enskog展開のChapmanは別人。Chapman-Enskog展開の人は、Sydney Chapman(1888-1970)という名前で、CJ理論の人は、David Chapman(1869-1958)らしいけど、David Chapmanには、Sydney John Chapman(1871-1951)という弟がいて経済学者だったらしい。ややこしいわ

CJ理論では、波面の前後が、それぞれ未燃ガス、既燃ガスになっていて、反応は一瞬で終了するものとされている。未燃ガス、既燃ガスは共に、理想気体として記述するのが一般的。基本的に、代数的な計算しか出てこないので、難しいわけでないけど、既燃ガスの比熱比は常温のものと違うことに気付かず、計算してみると、値が合わないという罠にはまった

後になって気付いたけど、デトネーションのことは、ランダウ・リフシッツの流体力学の教科書に載っていて、CJ理論が説明されてる。この本も、測定値と理論値の比較とかはないので、助けにはならなかったろうけど。

Fluid Mechanics §121
https://archive.org/details/FluidMechanics/page/n491

罠にはまった時に、Chapmanの原論文を読んだ。

On the rate of explosion in gases
https://www.tandfonline.com/doi/abs/10.1080/14786449908621243

Chapmanの議論は、現在の定式化とは違っていて分かりにくいけど、TABLE IIIに、"Velocity of Explosion"の式として

$V^2 = \dfrac{2RJ}{\mu C_{v}^2} \left( {(m-n)C_p + m C_{v} }C_{p} t_0 + (C_p+C_v)\mathrm{h} \right)$

が書かれていて、様々な混合気体の燃焼反応に対して、Vを計算し、測定値と比較している。単位とかが書かれてないけど、私が解読した限りでは、多分、以下のようになっている。

Chapmanの論文には、多くの混合気体に対する測定値と計算値が載ってるけど、全部やるのは大変なので、水素と酸素がモル比で1:1に混合した気体の燃焼反応(2モルの水素と2モルの酸素を燃焼させて、1モルの水蒸気が発生し、1モルの酸素が反応せず残るという状況)
$2\mathrm{H}_2 + 2\mathrm{O}_2 = 2 \mathrm{H}_{2}\mathrm{O} + \mathrm{O}_2$
を例に取ると、

RJ = 8.314462618(J/mol/K) : 気体定数
μ = 68(g) : 単位量あたりの気体の質量(反応前後で保存されるはずなので、水素2モル+酸素2モルの質量であり、また、水2モルと酸素1モルの質量でもある)
m = 3(mol) : 反応前のモル量(水2モルと酸素1モル)
n = 4(mol) : 反応後のモル量(水素2モルと酸素2モル)
h = 117000(cal) : 反応熱=1molのH2O(g)につき、245kJ弱の計算
t0 = 273.15(K) : 基準温度
Cv = 11.63(cal/mol/K) : 生成ガスの定積モル比熱(表に記載されてる値)
Cp = 13.62(cal/mol/K) : 生成ガスの定圧モル比熱(Cvから、マイヤーの関係式を使って計算)
V = 2341.8(m/s) : Chapman-Jouguet速度(論文では2340と計算されている)

だと思う。

生成ガスの比熱について、論文には、以下のようなexcuseが長々と書いてあって、時代を感じる

A few words are necessary regarding explosions in which water is formed. If the specific heat of steam is taken as 3/2 × specific heat of the diatomic gases, the found rates of explosion fall below the calculated rates when the dilution with inert gas is great, and vice versd when the dilution is small.

It is possible to account for this by two theories. The first theory is that at high temperatures the water is dissociated, whereas at low temperatures the combination of hydrogen and oxygen is complete. The second theory is that the specific heat of steam rises more rapidly with the temperature than the specific heat of the diatomic gas.
(中略)
We are therefore encouraged to test the first theory, i.e. that the specific heat of steam rises more rapidly with the temperature than that of the diatomic gases.

水蒸気の比熱を、2原子分子の1.5倍とすると、うまくいかないので、高温では、水蒸気の比熱は、2原子分子より急速に上昇するのだろうとか。Chapmanが書いている1.5倍というのは、多分、エネルギー等分配則の誤った適用による(比熱が一分子あたり原子数に比例すると考えたんだろう)ものだろう。Chapmanは、もう一つの可能性として、現在でいうところの弱電離プラズマのようなものを考えて却下している。

反応直後のガスは、高圧で、温度は数千度に達するらしいので、常温の場合と違って、水や酸素でも、分子振動の比熱への寄与が無視できなくなる。分子内振動を調和振動子で近似して、独立なモードに分離すると基準振動が得られる。基準振動波数は、赤外分光やラマン分光で測定でき、水素分子の基準振動波数は、ν=4160(/cm)程度。分光業界(?)では、波数の記号は、νの上にチルダを付けたものを使うことが多い気がするけど、以下では、普通のνを使う。回転スペクトルを測るマイクロ波分光はGHzを使ってるっぽいので、赤外/ラマンも、THz単位に揃えろやと思うけど。水の場合、適当に検索した数値によると、1595(/cm),3756(/cm),3657(/cm)。

世の中には合成法が確立してなかったり合成難度が高い化合物もあるけど、その場合でも(低分子なら)割と高速かつ、それなりの精度で、基準振動波数を数値計算することもできる。例えば、現時点で、合成法が確立してない低分子の一つに、テトラヘドランC4H4があり、分子量などからすると、常温常圧で気体の可能性がある(構造異性体の1,2−シクロブタジエンや、1,3-シクロブタジエンの沸点が40度前後という情報がある)

cf)Tetrahedrane
https://en.wikipedia.org/wiki/Tetrahedrane

振動モードが一個の場合、分子振動の分配関数は
$q_{\mathit{vib}} = \left( \dfrac{e^{-\beta c \mathit{h} \nu/2}}{1 - e^{-\beta c \mathit{h} \nu}}　\right)^N$
なので、一個の振動自由度による定積比熱への寄与は
$C_{v,\mathit{vib}} = N k_{B} \dfrac{x^2 e^{x}}{(e^{x}-1)^2}$
と計算できる。但し、 $x = \dfrac{c h \nu}{k_{B}T}$

ν=4160(/cm)、T=300(K)だと、x=19.95で、 $y=\dfrac{ C_{v,\mathit{vib}} }{N k_{B}}$ は、8.62e-7とほぼ無視できる。T=3000(K)だと、x=1.995で、y=0.725になって、無視できない寄与があることが分かる。

十分高い温度では、振動の自由度1つにつき、モル比熱は気体定数Rだけ増加する。振動自由度を高温近似で扱うと、高温の二原子分子理想気体では、定積比熱が3.5R,定圧比熱が4.5Rになり、高温の非直線型三原子分子理想気体では、定積比熱は6R,定圧比熱は7Rになる(はず)。

6R=50(J/mol/K)=12(cal/mol/K)程度で、上のChapmanの計算に出ている11.63(cal/mol/K)に近い値になっている。実際に上の計算で扱ってるのは、水蒸気2モルと酸素1モルの混合気体の定積比熱なので、それぞれ理想気体として、高温近似が成立するならば、定積比熱は、10.27(cal/mol/K)程度という計算。完全に一致とは言えないけど、それなりに合ってそうではある

100年以上前の論文の数値と理論値だけを頼りにするのもどうかと思うので、もうちょっと新しいデータを見ておく(数千度とかで比熱を測定できるのか知らないけど、以下のデータは、測定値というわけではない)

Water@NIST Chemistry Webbook
https://webbook.nist.gov/cgi/cbook.cgi?ID=C7732185&Units=SI&Mask=1&Type=JANAFG&Table=on#JANAFG

を見ると、水の6000(K)に於けるCpの値は、60.59(J/mol/K)となっていて、これは、7R=58.20(J/mol/K)と近い。4100(K)では、丁度、58.20(J/mol/K)となっている。二原子分子気体の例として、H2(g)を見ると

Hydrogen@NIST Chemistry Webbook
https://webbook.nist.gov/cgi/cbook.cgi?ID=C1333740&Units=SI&Mask=1&Type=JANAFG&Table=on#JANAFG

に、3200(K)で、Cp=37.53(J/mol/K)とある。これは、4.5R=37.415(J/mol/K)と、ほぼ一致している。上の基準振動の数値に基づく計算を考慮すると、3200(K)で、この数値に達するのちょっと早くね？という気もするけど。6000(K)に於けるCpの値は、41.97(J/mol/K)となっている。

この手の小さな多原子分子では、3000〜4000(K)で、比熱への振動自由度の寄与を高温近似で扱うことが許されそうにも思える。幸い、デトネーションに於ける生成ガスの温度は、この水準にある

これらのデータが、どれくらい信用できるものなのか分からないけど、並進、回転、分子内振動を単純に考慮しただけのモデルとは少しズレがあるけど、実験値と比較しないことには、何も言えそうにない。

【おまけ】大学一年の時、エネルギー等分配則による気体の比熱の説明を習った際、結構ずれてる気体もあって胡散臭いと思った記憶がある(条件とか分子は忘れた)。ファインマンの教科書の

The Feynmann Lectures on Physics　Vol.I 40-6 The failure of classical physics
http://www.feynmanlectures.caltech.edu/I_40.html#Ch40-S5

を見ると、古典統計力学で気体の比熱を説明できないことは、丁度Chapmanの論文が書かれた時代の物理学者も不可解に思ってたと書いてある。そんな説明を21世紀にもなって教えないで欲しかった。ファインマンは、分子振動を調和近似して量子統計力学で扱えばいいと言ってるけど、計算例は一つもないし、他で見たこともないので、やってみたのが、以下の表。

分子	Cp(J/mol/K)	Fr	Cp/R-2.5-Fr/2	基準振動(/cm)	X	R*(2.5+Fr/2+X)
N2(g)	29.124	2	0.0028	2358.57	0.001477	29.113
O2(g)	29.376	2	0.0331	1580.19	0.028395	29.337
HCl(g)	29.136	2	0.0043	2990.95	0.000112	29.102
HI(g)	29.156	2	0.0067	2309.01	0.001798	29.116
Br2(g)	36.048	2	0.8356	325.32	0.8177	35.899
CO2(g)	37.129	2	0.9656	[667,667,1388,2349]	0.9565	37.053
N2O(g)	38.617	2	1.1446	[589,589,1285,2224]	1.1428	38.603
NO2(g)	36.974	3	0.4470	[750,1320,1617]	0.4653	37.127
O3(g)	39.238	3	0.7192	[719.6 ,1062.1 ,1157.8]	0.6745	38.866
H2S(g)	34.192	3	0.1124	[1183,2615,2626]	0.1098	34.171
NH3(g)	35.652	3	0.2880	[950,1627,1627,3337,3444,3444]	0.2671	35.479
ND3(g)	38.225	3	0.5974	[793,1225,1225,2495,2652,2652]	0.5254	37.626
PH3(g)	37.102	3	0.4623	[992.1,1118.3,1118.3,2322.9,2327.7,2327.7]	0.4656	37.130
BF3(g)	50.446	3	2.0673	[480.7,480.7,696.7,888,1463.3,1463.3]	2.0673	50.446
C2H2(g)	44.095	3	1.3034	[612,729,729,1974,3289,3374]	1.2929	44.077
CH4(g)	35.639	3	0.2864	[1306,1306,1306,1534,1534,2917,3019,3019]	0.2864	35.639
CH3Cl(g)	40.731	3	0.8988	[732,1017,1017,1355,1452,1452,2937,3039,3039]	0.8998	40.739
SF6(g)	96.960	3	7.6616	[346,346,346,524,524,524,...](略)	7.6648	96.987
C3H8(g)	73.60	3	4.8520	[216, 268, 369, 748,...](略)	4.524	70.874

定圧モル比熱Cpの出典は、JANAF熱化学データ表かNIST Chemistry WebBookで、298.15(K) , 0.1(MPa)の値。Frは回転の自由度で、Rは気体定数。ND3は、アンモニアの水素を全部重水素に置換したもの。

表の値Xは、全基準振動波数について、以下の和を計算したもの(温度Tは、298.15K)
$X = \displaystyle \sum_{i} f(\dfrac{c h \nu_i}{k_{B}T})$
$f(x) = \dfrac{x^2 e^{x}}{(e^{x}-1)^2}$
並進と回転運動からのみ比熱への寄与がある理想気体を高温近似で考えると、Cp=(5＋Fr)R/2になる。分子内振動からの比熱への寄与を量子統計力学に基づいて計算し、加えると、Cp=(2.5+Fr/2+X)*Rになる。これを測定値と比較すればいい

基準振動の出典は、以下のどれか
(1)JANAF熱化学データ表
(2)Vibrational Modes in HITRAN2012 Database
https://www.cfa.harvard.edu/hitran/vibrational.html
(3)Tables of molecular vibrational frequencies, consolidated volume I (Takehiko Shimanouchi)
https://nvlpubs.nist.gov/nistpubs/Legacy/NSRDS/nbsnsrds39.pdf
https://archive.org/details/tablesofmolecula39shim

(2)のアセチレンのデータは、縮退度にミスがあるので修正してある。また、六フッ化硫黄の基準振動は(2)から、
[346,346,346,524,524,524,615,615,615,643,643,775,948,948,948]
を使い、27個あるプロパンの基準振動は、(3)に従って、以下の数値を使った
[216,268,369,748,869,922,940,1054,1158,1192,1278,1338,1378,1392,1451,1462,1464,1472,1476,2887,2887,2962,2967,2968,2968,2973,2977]

計算値は、これなら、過去の自分も納得しただろうと思う程度には合ってる。どっちみち、調和近似が完全に正確な結果を与えるはずはないので、多少のずれは、あって然るべきもの

Chapmanの論文にあるCJ速度の式は、実際の数値を見ると、大抵の場合、第二項が卓越しており、更にマイヤーの関係式を使うと
$V^2 \approx \dfrac{2(C_p - C_v)}{\mu C_v^2} (C_p + C_v) h = \dfrac{2(C_p^2-C_v^2)}{C_v^2} \dfrac{h}{\mu}$
と書ける。この近似式は、現代の標準的なCJ理論でも導くことができる。冒頭のランダウ・リフシッツの教科書だと式(121.9)

水素2モルと酸素1モルの爆発の場合、生成ガスは、水蒸気2モルのみで、
h = 242(kJ/mol)
μ = 18(g/mol)
を使うことにする。更に、高温の非直線型三原子分子理想気体に対するCp=7R,Cv=6Rを使うと、CJ速度は
V = 3116(m/s)
となる。

デトネーション伝播の基礎
https://www.jstage.jst.go.jp/article/jcombsj/55/174/55_317/_pdf/-char/ja

の表1を見ると、CJ速度は2841(m/s)で、"既燃ガスの実効的な比熱比"は、1.129となっている。観測されるデトネーション速度も、こんなものらしい

上で使用した反応熱は298K,1atmのもので、0K換算の値を使うべきと思う。0Kに於ける生成エンタルピーは、定圧比熱を、0(K)から298(K)まで積分した値を引けば得られる。面倒なので、(振動の影響を無視した常温での)定圧比熱と基準温度の積で見積もることにすると、H2(g)の場合、定圧比熱が3.5Rとすれば、8.672(kJ/mol)くらい引けばいい。O2(g)とH2O(g)に対して、同様に熱補正を考慮すると、0(K)に於ける反応エンタルピーは、上の値から3.1(kJ/mol)ほど引いたものになる。まぁ影響は小さい

もうひとつ、上の近似式は、反応ガス、生成ガスが共に理想気体と仮定されていて、定圧比熱が温度に依らず一定であることは重要。もうちょっとよく考えると、CJ理論のRankine-Hugoniotの式を使う定式化では、波面前後でのエネルギー保存則が、
$h_1(T_1) + \dfrac{u_1^2}{2} = h_2(T_2) + \dfrac{u_2^2}{2}$
と書かれる。添字の1と2は、反応ガスと生成ガスの物理量を表し、 $h_i,u_i,T_i$ は、比エンタルピーと流速、温度。 $q = h_1(0)-h_2(0)$ が、0Kに於ける(単位質量当たりの)反応熱(Chapmanの論文の記号に従うと、h/μに相当する量)。で、生成ガスが理想気体と仮定すれば
$h_2(T_2) = h_2(0) + \dfrac{\gamma_2}{\gamma_2-1} \dfrac{N k_B T}{\mu}$
となる。 $\gamma_2$ は生成ガスの比熱比。 $\mu,N$ は、反応前後の質量と、反応後の分子数。

$T_2$ は高温であるが、単原子分子気体の場合を除けば、分子内振動の影響で、比熱は全温度領域で一定でありえない。単原子分子気体は反応性が低く、デトネーションの文脈では多原子分気体が、いつでも重要だろう。多原子分気体では、上のエンタルピーの式を、そのままで使うのは、どう考えても良くない

考えられる処方として
$C = \lim_{T_2 \to \infty} \left( h_2(T_2) - h_2(0) - \dfrac{\gamma_2}{\gamma_2-1} \dfrac{N k_B T}{\mu} \right)$
を定義して、高温領域で
$h_2(T_2) \approx (h_2(0) + C) + \dfrac{\gamma_2}{\gamma_2-1} \dfrac{N k_B T}{\mu}$
と近似すれば、 $q = h_1(0)-h_2(0)-C$ と取りなおすだけで、Chapman-Jouguet速度に関する解析的な式を使うことも許されるはず

少し手を抜いて、並進と回転成分については、全領域で高温近似が妥当だとみなし、振動成分からくる比熱のみ、量子統計力学と古典統計力学で、差があると近似する。基準振動モードが一つで、波数νの場合は
$C = \dfrac{N k_B T}{\mu} \displaystyle \int_{0}^{\infty}(f(\dfrac{c h \nu}{k_B T}) - 1)dT$
$f(x)=\dfrac{x^2 e^{x}}{(e^{x}-1)^2}$
で、積分を実行すると
$C = -\dfrac{N}{\mu} \dfrac{c h \nu}{2}$
となる。基準振動モードが複数ある場合は足すだけ。H2O(g)の場合、1595(/cm),3756(/cm),3657(/cm)を使うと

Cμ = -6.02214076e23 * 2.99792458e8 * 6.62607015e-34 *(1595e2+3756e2+3657e2)*0.5 / 1000.0 = -53.88(kJ/mol)

となる。符号が負なので、CJ速度は、より大きく評価されるようになり、ますます測定値と合わなくなる(ダメ)けど、少なくとも、Chapmanの原論文の意図からすると、このような数値を使うのが正しい気がする。Cの値は熱力学データから決まらない(比熱の温度依存性からフィッティングを行えば、基準振動を知ることができるかもしれないけど、普通、そんなことはしない)せいか、普通、教科書を読んでも、こうは書いてない

まだ、生成ガスの比熱比の値が、考慮に入れてない効果によって、7/6より小さいかもしれないという希望は残ってるけれど

Lasing action and the relative populations of vibrationally excited carbon monoxide produced in pulse-discharged carbon disulfide-oxygen-helium mixtures
https://pubs.acs.org/doi/abs/10.1021/j100639a019
では、CS2-O2-He混合ガス火炎中に生成するCOの振動準位の分布が、点火後の経過時間で、どう変化するか計測している。論文Figure9によると、点火後、数十μsまでは、ボルツマン分布と異なり、特定の振動準位にピークがあり、200μs以上かけて、ボルツマン分布に近付くらしい。デトネーションでも、平衡に達するのに、同じような時間スケールが必要だとすると、例えば、速度3000(m/s)で、200μsの間に、60cm進む。これは、平均自由行程より遥かに長いから、化学反応が一瞬で終わって、短時間で平衡に達するという近似の方が問題かもしれない。

【補足】Chapman-Jouguet理論で必要な熱力学パラメータは、未知化合物でも、量子化学の力によって原理的には計算できる。0(K)に於ける生成エンタルピーは、電子エネルギー+分子内振動の零点エネルギーで計算すればいい。電子エネルギーは、量子化学計算で、ほぼ必ず計算する量で、電子系のハミルトニアンの最低固有値に相当する。分子内振動の零点エネルギーは、分子内振動を調和振動子で近似して最低エネルギーを見る。例えば、H2(g)の場合、基準振動波数ν=4160(/cm)とすれば、
0.5*6.02214086e23*2.99792458e8*6.62607004e-34*4160e2/1000 = 24.88(kJ/mol)
くらいの影響がある。実際上の問題としては、反応熱には、複数の分子の電子エネルギーの差が寄与することになるけど、通常、その大きさは電子エネルギーと比較すると非常に小さい。従って、電子エネルギー自体が1%程度の誤差で計算できても、反応熱の誤差は、割と大きくなり得る　