コンピュータ以前の数値計算(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,yの偶奇性によって、分解できることに起因する。

一応、行列要素の計算は、以下のコードでできるはず(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の結果のトレースとは趣旨が外れるので、それも省略