中世の数理科学

Wikipediaによれば、mathematicsの語源は、ギリシャ語のmáthēmaで、knowledge, study, learningを意味してたそうだ。polymath(博学)という単語には、原義が残っている。そして、ラテン語や英語に於いて、1700年頃まで、mathematicsという単語は、占星術(時々、天文学)を指していたと書かれてるけど、これは正しくない。中世ヨーロッパに於いて、「数学」という単語が、どういう分野を指してたかは、それほど固定的でもない。

以下、ヨーロッパと書いてる時、西方のラテン語文化圏を念頭に置いている。また、中世は、地域ごとに、時代区分が微妙に異なるけど、以下では、西暦800〜1600年を指す時代区分のつもりで使っている。

アラビア数学とQuadrivium

12〜13世紀のカスティーリャ王国には、トレド翻訳学派という、アラビア語文献をヨーロッパの言語(ラテン語、カスティーリャ語)に翻訳する仕事をしていた学者がいた。トレド翻訳学派には、ユダヤ人が多かったとも言われ、カスティーリャ王国では、1492年にアルハンブラ勅令が出されるまで、ユダヤ系学者が活動してた節がある。

トレド翻訳学派の一人にDominicus Gundissalinus(1115~1190)という人がいて、1140年に書いたDe divisione philosophiaeの中で、"数学は、算術、幾何、音楽、天文学、scientia de aspectibus(光学)、scientia de ponderibus(science of weights)、scientia de ingeniis(science of devices)の7つのartsを含んでいる"と書いている

Mathematica quoque universalis est, quia sub ea continentur septem artes, quae sunt arithmetica, geometria, musica et astrologia, scientia de aspectibus, scientia de ponderibus, scientia de ingeniis.

ラテン語のscientiaは英語のscienceと同根の語だけど、中世には、知識とか学問のようなニュアンスで、結構広い意味だったようである。格言「知識は力なり」の原文は、scientia est potentiaらしい。

astrologiaは、直訳では、英語のastrology/占星術だけど、中世ヨーロッパの大学では、リベラル・アーツという括りで、文法学・修辞学・論理学の3学と、Quadrivium(算術、幾何、天文学、音楽)という区分をよく使っていたことを踏まえると、天文学と解する方が適切に思える。Gundissalinusの分類では、Quadriviumより広い分野を含んでるけど、この分類は、多分、アラビア科学の学問論から来ている。Gundissalinusは、トレド翻訳学派の学者だったから、別に意外ではない。


アラビアの学者al Farabi(872?〜950?)がKitáb al-Ihsa' al-'ulum("Book of enumeration of the sciences")のChapter3で、"数学"('ilm al-Ta'alim)として分類している分野は、算術、幾何、天文学、音楽以外に、後のGundissalinusが、scientia de aspectibus、scientia de ponderibus、scientia de ingeniisと呼んだ分野に相当するものがあったらしい。

アラビア語Ta'alimは、instruction/teaching/educationなどの意味だそうで、ギリシャ語のmáthēmaと似た感じっぽい(現在のアラビア語で、「数学」を意味するのは別の単語)。ラテン語にしろ、アラビア語にしろ、この分野を指すのに、あまり適切な名前とは思えないんだけど、Quadriviumにしたって、内容を反映した命名ではないので、大した理由はなかったのかもしれない。al Farabiは、これらの分野を更にtheoreticalなものとpracticalなものに分けたとよく書いてあるけど、詳細は知らない。例えば、practicalな算術は、商業とかで使われる計算ということらしい。

Quadrivium自体は、少なくとも、古代ローマ帝国末期まで遡れるらしいけど、al Farabiの分類は、どっから来たものか、よく分からない。


ついでに、De divisione philosophiaeでは、"自然科学"(scientia naturalis)として、

  • scientia de medicina(薬学)
  • scientia de iudiciis(占術?)
  • scientia de nigromantia secundam physicam(science of necromancy!?)
  • scientia de imaginibus(science of image magic?)
  • scientia de agricultura(農学)
  • scientia de navigatione(航海術)
  • scientia de speculis(反射光学?)
  • scientia de alquimia(錬金術)

の8つを挙げている。

自然科学に分類されてるscientia de speculisと、数学に分類されてるscientia de aspectibusは、現代では、どちらも"光学"の一部に相当するが、多分、ユークリッドのCatopriticaとOpticaに対応したものと思われる。Catopriticaは、鏡で反射する光の理論で、Opticaは、視覚が目から出ている視線によって生じるという考えをベースにした視覚の幾何学のようなもの。aspectibusは、aspectusの複数形で、aspectusはvisionを意味するらしいから、scientia de aspectibusは、直訳では、science of visionsとなる。

但し、プトレマイオスの本"Optica"では、両方を扱ってるっぽいし、「光学の父」と呼ばれるIbn al-Haitham(965〜1040)は、視覚が物体から目に到達した光によって生じると理解してたようなので、この区別はGundissalinusの不見識という気もする。Gundissalinusの"自然科学"には、完全にオカルトな分野も混じってるけど、いくつかは、当時のアラビア科学でも時代遅れなものだった可能性もある。


この時代には、類似の分類は色々あって、Gundissalinusのも、その一つに過ぎない。別の例として、Hugh of Saint Victor(1096〜1141)がDidascalicon de studio legendiの中で述べてる分類では、数学は、以下のように、算術、音楽、幾何、天文学の4つの知識に分かれるみたいなことが書いてあって、Quadriviumと同一と思ってよさそう。

mathematica dividitur in quattuor scientias. prima est arithmetica, quae tractat de numero, id est, de quantitate discreta per se. secunda est musica, quae tractat de proportione, id est, de quantitate discreta ad aliquid. tertia est geometria, quae tractat de spatio, id est, de quantitate continua immobili. quarta est astronomia, quae tractat de motu, id est, de quantitate continua mobili.

ここでは、4つの分野が、それぞれ異なる種類のquantitate(quantity/量)を扱うとして特徴付けされている。al Farabiの分類も、基本的には、「何らかの量を扱うのが数学だ」という思想上にあると思われる。

12世紀に於いて、mathematicaは、「定量科学」くらいのニュアンスで認識されてたのだろう。現代では、数学と実験科学は遠い印象があるけど、scientia de aspectibusの研究者としてIbn al-Haithamは、光学実験をしてたようなので、アラビア科学では、"数学"だから実験しないということは、別になかったと思われる。そもそも、アラビア天文学が、観測と計算が合うことを重視した分野だし、中世ヨーロッパの天文学者も、しばしば天文観測機器を自分で作っている。

多分、中世ヨーロッパの大学に於いて、数学とはQuadriviumのことでしたというのは、そんなに間違ってないけど、大学の外では、もっと緩いイメージで捉えられてたかもしれない。とりあえず、以下で、Gundissalinusが「数学」に分類した分野を個別に概観していく。

但し、幾何学と算術は、現代の意味でも数学の一部と見なされてるので、特に触れない。また、音楽は、よく分からないので、やっぱりpass。

数理天文学

天文学は、Quadriviumの一つではあるけど、西暦1200年頃のヨーロッパは、数学の下地が何もない状態なので、数理天文学ができるようになるのは長い時間が必要だった。江戸時代の日本も、(商人の会計や財政管理はあったにしても)数学の下地がほぼない状態から、約250年かけて和算が発展し、平方根や代数方程式の数値解法、三角関数、初歩の微分と積分(極値計算法と求積法)、対数と、それらの計算法を、測量や天文学、建築(規矩術)などに応用する術を習得した。

現代の数学教育でも、小学校から(殆どが1700年までに発見された結果しか含んでいない)高校数学を学習するのに約10年かけてることを考えると、知識が整理されてない時代に、知識源にアクセスするには言語の壁がある状態で、数百年単位の時間が必要だったとしても、そういうものなのかもしれないという気がする。


以下では、中世ヨーロッパにおける、数理天文学の受容状態を簡単に追跡する。

トレド翻訳学派に、クレモナのGerard(1114〜1187)という人がいる。多くのアラビア語文献を翻訳したが、1175年に、プトレマイオスのアルマゲストのラテン語翻訳を行い、写本しかない時代にしては、広く回覧されたと言われる。13世紀初頭には、パリ大学にも写本があったようだ。

SacroboscoのJohannes(1195〜1256)という人は、その生涯は曖昧だが、オクスフォード大学で学んだ後、1220年代初頭からパリ大学でQuadriviumを教えるようになったと伝えられる。1230年頃、De sphaera mundi("On the sphere the World"、天球論)という本を書いて、この本は、中世ヨーロッパの大学で、天文学の標準的な"教科書"となり、1472年には、フェラーラで印刷されたらしい。ガリレオの時代でも使われてたようだ。内容の詳細は知らないけど、プトレマイオスや中世アラビア天文学の初歩的な解説を与えていたらしい。三角関数のような内容は含んでなかったぽいので、この本で、十分な計算ができるようにはならなかったと思われる。

13〜14世紀のヨーロッパで、数理天文学をやってた人たちは、情報があんまりない。13世紀のカスティーリャ王国では、アルフォンソ天文表が作られており、1270年代に完成した。計算に関わったのは、トレド翻訳学派の学者Yehuda ben MosheやIsaac ben sidという人らしい。

Levi ben Gershon(1288〜1344)は、おそらくフランス出身のユダヤ系学者らしいが、生涯は殆ど何も分かってないそうだ。1342年にローマ教皇クレメンス6世の依頼で、Sefer Milhamot Ha-Shem("The Wars of the Lord")という"科学/哲学史"の本の一部をラテン語に翻訳したとか言われている。1342年に、De sinibus, chordis et arcubus("On sine, chord, and arcs")という三角法の本を書き、ヤコブの杖という天体測量機器の発明者とされることもある。

Georg von Peuerbach(1423〜1461)は、ウィーン大学(1365年創立)で学び、天文測量機器の作成や三角関数表の計算を行い、アルフォンソ天文表の修正を試み、そして1450年代に世界最初の天文学教授となった。彼は、パドヴァ大学やボローニャ大学、フェラーラ大学などで教えたこともあるようだ。Peuerbachは唐突に出現したわけではなく、彼が師事したJohannes von Gmunden(1380~1442)は、ウィーン天文学派の創始者とされ、本格的に数理天文学を開始しようと努力していたようだ。

1514年に、Viri Mathematici quos inclytum Viennense gymnasium ordine celebres habuitという本が出版されていて、ウィーン天文学派の数学者の伝記が記載されているそうである。16世紀初頭には、天文学者が、mathematici(mathematician)と見られていたことが分かる。

Peuerbachの生徒に、レギオモンタヌス(1436~1476)がいて、1464年に 三角法の著書De Triangulis omnimodusを著し、1474年に天文表Ephemerisも出版した。レギオモンタヌスの著書は、丁度、ヨーロッパで商業出版が開始した時期だったこともあり、かなり広く読まれ、影響が大きかったと思われる。

レギオモンタヌスは、1475年、教皇シクストゥス4世から、ユリウス暦の改暦相談のため、ローマに呼ばれている。グレゴリオ暦の運用を開始したのが1582年だから、改暦に至るまで100年以上かかってる。ユリウス暦への批判自体は、もっと古くからあったようである。江戸時代の日本で、西洋天文学ベースの改暦を試みた際も、何度か失敗して長い時間がかかってるので、そういうものらしい。

カスティーリャ王国出身のAbraham Zacuto(1452〜1515)は、サラマンカ大学で、天文学を勉強し、1492年に、アルハンブラ勅令でスペインを追放された後、ポルトガル王ジョアン2世に仕えた。1496年には、ポルトガルでもユダヤ人追放令が出され、その後の足跡は、詳しくは分かってないっぽい。

Zacutoは、1478年頃、Almanach perpetuumという本をヘブライ語で完成させた。これには、太陽の黄道座標系での予測位置以外に、赤緯表も含まれていたっぽい。Zacutoが、この本を書いた理由は知らない。レギオモンタヌスが、同時期に書いて1490年に出版されたTabulae Directionum profectionumqueにも、赤緯表が含まれてるらしいが、どれが何の表なのか、よく分からない。黄道座標での位置を赤道座標に変換するのは、別に難しくない気がする(プトレマイオスの「アルマゲスト」でも出来てた話だと思う)けど、それ以上の何かがあったのか、それとも、当時は、座標を変換するのも簡単とは思われてなかったのか、不明。

ポルトガル王ジョアン2世が集めた天文学者/数学者たち(しばしば、Junta dos Mathematicos/数学委員会という呼称が使われるが、当時から、そのように呼ばれてたかは不明)は、Zacutoの表をベースにしつつ、1484年にRegimento do estrolabio e do quadranteという、簡単な数表を含む航海マニュアルにまとめ、Regimento das léguas、Regimento do norte、Regimento do Solと呼ばれる遠洋航海の手法が記載されているらしい(1914年版がオープンアクセスになってるけど、殆ど何もわからない)。もう少し後に、Regimento do Cruzeiro do Sulというのも追加され、これらの手法は、江戸時代の「元和航海書」(1618年)にも、天文表と共に解説されているっぽい。

詳細はともかく、15世紀末に、ジョアン2世の雇用した数学者たちが、新しい遠洋航海術を発明して、100年以上に渡って使われたことは確からしい。コロンブスは、アメリカ大陸到達時に、レギオモンタヌスの出版したEphemerisを携行していたという逸話が残っている(これは本当かどうかは分からないが、何らかの天文表は持ってたのだろう)。ヴァスコ・ダ・ガマも、何らかの天文表を使っていたと信じられているが、どの本かについては、多少の議論があるようだ(cf. Nautical Tables for Vasco da Gama, 1497–1500?)。

そんなわけで、15世紀には、ヨーロッパも、数理天文学を習得し、航海術などに応用することができるようになった。

scientia de aspectibus

ヨーロッパで光学の研究が盛んになるのは17世紀で、理論方面では、フェルマー、スネル、デカルト、ホイヘンスなどの光学研究は現在まで知られている。望遠鏡や顕微鏡のような光学機器も、この時期に作られた。ヨーロッパでは、13世紀末に老眼鏡が作られたと言われ、1600年頃のヨーロッパでは、近視/遠視用共に眼鏡は一般的なものになっていたらしい。ガラス職人のレンズ加工技術も高水準になってたんだろう。

ヨーロッパの光学研究先駆者に、13世紀のポーランド人Witelo(あるいはVitello)という人がいる。彼は、1260年頃、パドヴァ大学で学び、1270年代にPerspectivaという本を書いたそうで、この本は、Ibn al-HaithamのKitāb al-Manāẓir(The Book of Optics)に基づいて書かれたと考えられている。1572年に、Friedrich Risner(1533〜1580)が、Opticae thesaurusという本を書いて、そこにal-Haithamの本のラテン語訳も含まれてたようである。この本は、ヨーロッパで広く読まれ、光学研究ブームの切っ掛けとなったっぽい。Keplerは、1604年に、"Ad Vitelonem paralipomena”という本を書いた

記録上では、最古の望遠鏡は、オランダ人Hans Lipperheyが1608年に特許を申請したものとされる。Keplerの本は、それ以前のもので、内容的には、光学というより視覚論のようである。天体観測への利用を想定してたわけではないっぽい。尤も、ティコ・ブラーエは、Camera obscuraを利用してたし、大気の屈折率の高度依存性によって、星の視高度がずれること(大気差というわかりにくい名前が付いてる)にも気付いてて、補正を行ってたらしいので、視覚論への興味も、天体観測の問題から生じてた可能性はある(が、本当かどうかは分からない)。

13世紀に生きてたらしいJordanus de Nemore(1225?〜1260?、Jordanus Nemorarius)は、その生涯が殆ど何も分かってない人で、Liber de Speculisという光学の著書があるそうだが、未出版。

Jordanusは、Liber de ponderibusという"静力学"の著書で最も知られる。また、Jordanusは、De plana speraという本で、ステレオ射影に関する命題を幾何学的に証明したようだが、動機は、天文観測機器アストロラーベにあったらしい(cf.The Mathematics of the Astrolabe and Its History)。Jordanusは、算術なんかの教科書も書いてるようだが、扱ってるテーマは、Quadriviumを超えて、al FarabiやGundissalinusが"数学"とした領域に及んでいて、アラビア科学の"数学"から影響を受けたのでないかと思われるが証拠はない。

いずれにせよ、Jordanusの光学研究が、ヨーロッパに影響を与えたことはなく、中世の間、ヨーロッパに光学研究は殆どなかったように思える。おかげで、ヨーロッパの光学研究が、アラビア科学に起源を持つことは、かなりはっきりしている。

その後、光学の知識は、望遠鏡の改良のために、かなり長い間、天文学者の教養の一部のようなものとなった。

scientia de ponderibus

ponderibusは、pondusの複数形で、pondusはweightなどの意味らしい。ポンドは、重さの単位なのでわかりやすい。従って、scientia de ponderibus('ilm al-athqal)は、直訳では、science of weightsとなる。

The Arabic Transformation of Mechanics: The Birth of the Science of Weightsに、該当するアラビア語文献が32冊挙げられている。6冊は、古代のギリシャ語文献の翻訳で、1600年以降の文献を、少なくとも8冊含み、19世紀の本まである。ギリシャ語文献は、アリストテレスのMechanica、HeronのMechanica、PappusのMechanicaを含むので、scientia de ponderibusは、Μηχανικά/Mechanicaの直接の子孫と言っていいだろう。

Euclidの"on balance"と"on heaviness and lightness"という著作もあるらしい(原題は不明)。これらの本は、梃子の原理だったり滑車だったりを扱ってて、ユークリッドの他の著書と同様、公理的なスタイルで書かれてたようで、後者は、9つの公準と6つの定理が含まれていたとも書いてある。

Mechanicaは機械学と訳されることが多いので、scientia de ponderibusも"機械学"と訳すのが妥当っぽい。ただ、機械と言っても、梃子、天秤、滑車とかなので、現在イメージする機械とは、ちょっとずれてる。mechanicaは、英語では、mechanicsになるけど、classical mechanicsやfluid mechanicsのように、力学の意味で使う今と違って、昔は、機械学の意味で使われてる。現在、機械工学は、英語だとmechanical engineeringなので、mechanicsから派生したんだと推測できる。


mechanicsの使用例として、ガリレオ(1564〜1642)が書いたRaccolta di Quelle cognizioni che à perfetto Cavaliero and Soldato si richieggono, the quali hanno dependenza dalle scienze matematiche ("完璧な騎士と兵士に求められる数理科学の知識")という文書を見る。これは、1608年、パドヴァに作られたAcadémie Deliaの教育カリキュラムに関するガリレオの提案書らしい。この文書では、算術、幾何、体積計測、機械学(scienze mecaniche)、コンパスや計測機器の使い方に加えて、砲術、軍事建築/築城術(Architettura militare)などが挙げられている。

ここのscienze mecanicheは、後に続く文章を見ると、

Cognizione delle scienze mecaniche, non solo intorno alle loro ragioni e fondamenti communi, quanto intorno a molte machine ed instrumenti particolari, insieme con la resoluzione di moltissime questioni e problemi da essa cognizione mecanica dependenti.

のように、"machine ed instrumenti"(machines and instruments)というフレーズが出てきて、様々な機械や機器に精通して、トラブルに対処できなければならないという感じのことを書いているので、"機械学"と訳すべきだろう。

ついでに、ガリレオがArchitettura militareと書いてるのは、16~19世紀に存在した築城術(fortification)という分野を指していると思われる。15~19世紀には、ヨーロッパを中心として、銃火器に対抗するために生まれた星形要塞というのが流行っていて、日本では、函館五稜郭に見られる。最古の星形要塞は、イスタンブールにあるYedikule要塞(1458年建設)とされる。

築城術は星形要塞の外形形状設計術みたいな分野っぽい。少なくとも、中世の時点では、土木工学的要素は殆どないけど、築城術は、数学の応用と見なされていたらしく、1736年の本The Young Trigonometer's Compleat Guideにも、築城術への応用が書かれている。この本は、適当に検索して見つけた古い本で、別に有名とかではないと思う。ガリレオ自身も、fortificationに関する原稿を書いている(例えばTrattato di fortificazione)。

実際のところ、中世の築城術自体には、(当時の水準で言っても)数学や物理として高度なことは何もないように見える。ただ、ウィトルウィウスのDe architectura(BC30頃?)には、幾何学、日時計、天文学、機械に関する記述があって、建築を監督する者は幅広い数学的素養を持つべきという思想が存在したそうである。築城術を数学の一分野のように見なす発想も、そういうところから来てるのかもしれない。


また、mechanicsの17世紀英語の文例として、John Wallis(1616〜1703)が、1640年代のロンドンにあった非公式の集まりについて書いてる文章には

Our business was ( precluding matters of theology and state affairs ) to discourse and consider of Philosophical Enquiries , and such as related thereunto ; as Physick , Anatomy , Geometry , Astronomy , Navigation , Staticks , Magnetics , Chymicks , Mechanicks , and Natural Experiments ;

とある(出典:The Royal Society 1660-1940)ようで、Staticksが静力学に対して、Mechanicksは機械学だろう。

scientia de aspectibusと違って、アラビア科学のscientia de ponderibusのヨーロッパへの影響は、はっきりしてないが、少なくとも、"mechanics"とscientia de ponderibusは、どっちも、古代ギリシャ/オリエントのMechanicaに由来する同一ルーツの分野と思われる。そして、1600年頃のヨーロッパで、mechanicsを数理科学の一つに位置付ける場合があり、Gundissalinusも、scientia de ponderibusを数学の一分野としていた。


mechanicsのニュアンスが変化した要因は、EulerのMechanica sive motus scientia analytice exposita (1736)が大きいかもしれない。その前後の力学の本を見てみると、Jakob Hermannは、 Phoronomia, sive De viribus et motibus corporum solidorum et fluidorum libri duo (1716)を書き、d'Alembertは、Traité de dynamique (1743)を書いている。Lagrangeは、Mécanique analytique (1788)を書いた。phoronomyは、kinematicsと同じ意味だそうで、kinematicsは、運動を直接扱い、運動の原因(古典力学なら力や慣性)は扱わないというニュアンスがあるとされている。

多分、dynamicsとmechanicsの差は、あんまりない。流体力学は、昔はhydrodynamicsだった(Bernoulliの1738年の本が最初で、Lambの1895年の本のタイトルにも見られる)けど、最近は、fluid mechanicsが一般的な気がする。量子力学は、quantum mechanicsなのに、量子電磁力学や量子色力学が、quantum electrodynamicsやquantum chromodynamicsなのは、どういう経緯があったのか謎だけど、後者の2つは、相互作用の原因に迫ってるから(?)。古典力学だと、重力やバネの力、摩擦力なんかを、現象論的に導入するが、その正体には関知しないので、dynamicsよりmechanicsの方が適切というニュアンスが現在はあるかもしれない。だとしても、この使い分けも、18世紀や19世紀に存在したものではないと思う。

scientia de ingeniis

scientia de ingeniis('ilm al-hiyal)は、該当文献が不明で、実態が一番よく分からない。

Jordanusのアストロラーベとステレオ射影に関する本De plana speraを、scientia de ingeniisに分類してる人もいる。Jordanusは、そう考えてたのかもしれないが、アラビア語文献ではないので、al Farabiの考えとは違うかもしれない。

イスラム法学では、hiyalとは、"本来なら違法にしか達成できない目的を法律的な工夫で合法に達成する手段"を指すらしい。現代だと節税などが例になるそうだ。おそらく、原義は"工夫"とか、その程度の意味だったと思われるが、そのhiyalと同じ単語っぽい(?)。また、ラテン語ingeniisは、ingenium(辞書によると、machine,deviceの意味があったらしい。engineと同根ぽい)の複数形で、ingeniumに"人"を意味する接尾辞torと複合すると、ingeniatorになる。ingeniatorは、英語のengineerなので、scientia de ingeniisは、engineering scienceつまり工学としてもいいだろう。


単語としてはそうでも、実際の内容は違う可能性もある。NOTES ON MECHANICS AND MATHEMATICS IN TORRICELLI AS PHYSICS MATHEMATICS RELATIONSHIPS IN THE HISTORY OF SCIENCEには、"The science of devices refers to applications of mathematics to practical use and to machine construction"と書いてある。

タイトルにingeniisが含まれるラテン語写本には、

  • De ingeniis spiritualibus
  • De decem ingeniis seu indicationibus curandorum morborum (1300年頃)
  • De aquarum conductibus et ingeniis erigendis

などがあって、1つ目は、ビザンチウムのフィロンによる著作Mechanike syntaxisのPneumaticaの一部の翻訳らしい。pneumaticsは空気圧工学とか訳されるらしいが、当時は流体で駆動する"機械"(水車やポンプ、水オルガンなど)全般を指してて、現代の水理学に相当する内容も含んでた可能性もある。2つ目は、何か"医療機器"(?)の話っぽいけど、中世のラテン語に確信が持てない。3つ目は、アレクサンドリアのHeronのPneumatica(アイオロスの球という蒸気機関の一種の記載があるらしい本)のラテン語訳っぽい。訳者はトレド翻訳学派のWillem van Moerbekeらしい。

タイトルにhiyalが含まれるアラビア語文献としては、

  • Kitáb al-Ḥiyal(一般に"The Book of Ingenious devices"と訳される)
  • Kitáb fi ma'rifat al-hiyal al-handasiyya(一般に"The Book of Knowledge of Ingenious Mechanical Devices"と訳される)

があり、前者は、9〜10世紀のバグダッドにいたBanū Mūsā兄弟の著書。当時知られていた多くの発明を記載しているそうである。彼らは、数学や天文学の著書も残してるらしい。後者は、12世紀のアラブ人Al-Jazariによる著書で、カムシャフトやクランクシャフトを記載した最初の文献とされる。この中には、水力機械の記述もあるらしいので、pneumaticaの影響もあったのかもしれない。

本のタイトルを見る限り、このへんの分野が、scientia de ingeniisという気もするけど不明。


オスマン帝国の数学者/天文学者Taqi al-din(1526〜1585)が書いた機械式時計の本が2冊あって、一冊は、al-Kawākib al-durriyya fī waḍ ҁ al-bankāmāt al-dawriyya(1556or1559?)で、もう一冊は成立年不詳のKitab al-Ṭuruq al-saniyya fī al-ālāt al-rūḥāniyyaである。当時のオスマン帝国は商業出版を行ってない(活版印刷はイスラムの戒律に反するという主張があったとか言われている)ので写本だと思うけど、後者はオープンアクセスになってるのを見つけた。私には微塵も読めないので、本物かどうかすら分からない。

Wikipediaによれば、後者には、Banū Mūsā兄弟とAl-Jazariへの言及があり、時計のgeometrical-mechanical structureを強調したと書いてある。つまり、Taqi al-dinは、自身の機械式時計の研究開発、Banū Mūsā兄弟の本、Al-Jazariの仕事を同分野と位置付け、数学的なものと見なしていた可能性が高い。Taqi al-dinの時計は、最小単位が5秒だったそうで、天文学者でもあったTaqi al-dinは、機械式時計を天文観測に利用した。

イスラム世界では、13世紀末に、muwaqqitという礼拝の時間を管理する仕事が誕生し、天文学の知識を持つ者が選ばれたそうである。それに関連して、ʿilm al-mīqāt(science of timekeeping)という分野が発生したらしい。どういう内容なのか知らないけど、初期の研究者として、Shams al-Din al-Khalili(1320〜1380)とIbn al-Shatir(1304〜1375)という人がいる。ヨーロッパでも、時計学/horologyという分野名が、かつては存在していたけど、ʿilm al-mīqātとの関連性は不明。Taqi al-dinの時計研究には、そういう方向からの影響もあるのかもしれないが、これも不明。

仮に、水車、オートマトン、機械式時計などの研究が、scientia de ingeniisだとすると、こっちも"機械学"に見える。個人的な感覚では、scientia de ponderibusとscientia de ingeniisの違いは、"機械"の複雑さという感じがする。scientia de ponderibusで扱う"機械"(天秤、滑車、梃子など)は、原始的な機械要素が多く、理論的解析はしやすいが、機械の印象が薄い(レバーとか投石機になると、機械感がなくもないけど)。scientia de ponderibusがmechanicsだとすれば、scientia de ingeniisは、mechanical engineering/機械工学に近い分野かもしれない。artes mechanicae(織布、建築、軍事・狩猟、鍛冶・冶金などの技術で、リベラル・アーツと対比される)との区別が曖昧でもある。


Taqi al-dinの一冊目の本では、ヨーロッパの機械式時計コレクションを調査した旨が述べられているそうなので、16世紀時点でヨーロッパの機械技術は高かったと推測される。イギリス人Richard of Wallingford(1292~1336)という人が書いたTractatus Horologii Astronomici (1327)は、機械式時計に関する(多分、世界で)最初期の文献らしい。彼は、オクスフォード大学で学んだ修道僧だったそうだ。

中世ヨーロッパでは、多くの天文時計が作られたが、1579年に、Jost Bürgi(1558~1632)が秒針を持つ時計を作ったとされる。Bürgiが、どこで教育を受けたのか、全く分かってないらしい。Bürgiは、Napierと独立に対数を発見したとも言われ、後にKeplerの計算助手になった。また、2015年になって、三角関数を計算するBürgiのアルゴリズムの詳細と証明を書いた論文が出ている(Jost Bürgi's Method for Calculating Sines (arXiv:1510.03180))。

オランダのホイヘンス(1629〜1695)は、機械式時計に対して理論的なアプローチを行い成功した。ホイヘンスは、1673年に、Horologium oscillatoriumという本を出している。タイトルは、振り子時計を意味する。ホイヘンス以前は、日差15分程度くらいの誤差があったらしいので、大体、1%くらいのずれがあったようだ(cf.機械式時計の歴史 1)

ついでに、当時の時計の分解能では、重力加速度を、自由落下時間によって高精度に決定するのは難しかったので、振り子の周期と長さから間接的に重力加速度を決定することが試みられた。最初に、この測定をしたのは、フランスのメルセンヌっぽいけど、単振り子では、振れ角に依存して周期が変わることも実験的に発見していたようだ。ホイヘンスは、サイクロイド振り子を使って、周期が二秒になる振り子(seconds pendulumと呼ばれる)の長さを測定し、この本に数値が書いてあるらしい(しかし、どこに書いてあるのか発見できなかった。ラテン語力がほしい)。

現代の単位では、この長さは、約1メートルで、18世紀末フランスのメートル法制定時に、これを基準値とする案は強く推された(cf.Why g〜π2、Why does the meter beat the second? (arXiv:physics/0412078)の特に§3と4)。また、ニュートンは、月の加速度と、(ホイヘンスの測定値から算出した)重力加速度の比率が、地球半径と地球・月間距離の比の二乗に近いことを確認した。

補足)逆二乗則の推測は、等速円運動のkinematicsのみでもできる(現在の知識だと、等速円運動の軌道を座標で書いて二回微分すると、遠心加速度が分かる)。地球と月の距離が、地球半径の約60倍であることは、プトレマイオスの時代には視差の観測によって知られていた。地球半径も、紀元前から多くの推定があったが、フランスのJean Picardの1669年の測量結果では、(現在の単位に換算して約)6372(km)だった。月の公転周期は約27.3日なので、仮に完全な等速円運動をしていたとすると、遠心加速度は、(60 * 6372e3) * (2pi / (27.386400))**2=0.0027129(m/s2)と見積もられ、9.8(m/s2)は、その3612倍になる。これは、60の二乗に近い。この簡単な試算から逆二乗則に到達した人は複数いたかもしれない。また、等速円運動の時は、Keplerの第三法則から逆二乗則を導くことも容易(角速度は、公転周期に反比例し、遠心加速度は、軌道半径と角速度の二乗の積だから、Keplerの第三法則が正しければ遠心加速度は軌道半径の二乗に反比例する)。ニュートンは、この結果を楕円軌道の場合に拡張していた。重力加速度と地球半径の高精度の測定は、惑星運動と自由落下が同一の原因で生じている可能性を検証できるようにした。大雑把に言って、ここまでが1700年頃の到達点だけど、等速円運動近似に基づく試算は、十分とは言えず、もっと精密なものが求められただろう。けど、精密な検証をしようとすると、月の軌道が、地球と太陽の両方から無視できない影響を受けてるせいで、難易度が一気にあがる。この困難は、ニュートンもオイラーも解決できなかった。

ホイヘンスの本と同時期の1674年に、デンマークの「数学者」Ole Roemerは、エピサイクロイド歯車に関する研究を発表した(文献未確認)。これ以降、歯車(の歯形曲線)を研究した「数学者」は多く、その中には、オイラーも含まれる。


16世紀〜17世紀前半の時点で、ヨーロッパに、機械に詳しい「数学者」が結構いたことは確かだけど、彼らが複雑な機械も数学的に解析できると思ってたのかは分からない(統一見解があったかも疑わしい)。当時は、複数分野に手を出す人が沢山いたから、単に機械にも詳しい「数学者」だったかもしれない。例えば、Gerolamo Cardanoは、三次・四次方程式の研究で知られる「数学者」だけど、本職は医者で、Cardan joint(Hooke's jointなど、色んな名前がある)というリンク機構に名前を残していて、錬金術に関する著書もある。

1695年に出版されたフランス人LahireによるTraité de mecaniqueという本の4ページには、"nouvelles machines qui ont été inventées par de tres-sçavans Mathematiciens,& par de tres-habiles ouvriers"(有能な数学者と職人によって発明された新しい機械)というフレーズがあって、17世紀末には、数学者を、機械の発明にも関与する人々と見なす場合があったっぽい。これはフランス語の本で、ヨーロッパ内でも、国によって認識の差異があった可能性は否定できない。


江戸時代の日本では細川頼直(??~1796)という人が、天文学と暦学を学んだ「天文学者」であり、『機巧図彙』(1796)という機械の本を書いている。この中には、時計や、お茶くみ人形のようなオートマトンが記述されているらしい。江戸時代のからくりも、機械式時計で使われていた機構を応用して発展したととか検索すると書いてある。第四章 実学としての和算|江戸の数学には、竹内度道(1780〜1840)という和算家が水車を設計し大工を指導してたとか書いてある。機械に関わる和算家というのが、数学者が機械に強いというヨーロッパの風潮が日本にも伝播しただけなのか、日本で独自に発展したものなのかは分からないけど、どっちだとしても、面白い事象ではある。

scientia de ingeniisは不明点が多いけど、結局は、al Farabiが10世紀に"数学"だと認識してたものを、ヨーロッパは、分裂した形で受容して、バラバラに発展させた後、再統合したというだけの話かもしれない。そして、19世紀には、また分裂した。

おまけ:16世紀の力学研究者

16世紀のmechanicsは、天文学と無縁の分野だったので、大学のカリキュラムと関係ない。そのせいか、16世紀の力学研究者は、その研究によって大学の職を得ていたわけではないし、独学の人も多い(従って「数学者」かどうかの定義が曖昧になる)。中世ヨーロッパの天文学者が、大体、大学で天文学を勉強してるのと比べると、対照的。

16世紀の力学研究者として、William Whewellの1840年の本History of the inductive sciencesに載ってる人として、タルタリア(1500?~1557)、Gerolamo Cardano(1501~1576)、Giambattista Benedetti(1530〜1590)、Michael Varro(1542〜1586)、Guidobaldo del Monte(1545~1607)、Simon Stevin(1548〜1620)、ガリレオ(1564~1642)などがいる。多分、何らかの著作が残ってる人を挙げていると思う。

ガリレオを除けば、大学で"数学"を教えた人は、いないと思う。こういう"野良数学者"が、昔からいて商業出版によって可視化されるようになったのか、この時期に実際に増えたのかは分からない。似た現象として、中世ヨーロッパの有名な魔術師や錬金術師も、この時期の人が多い。WikipediaにEuropean alchemistsというリストがあるけど、誕生年で分類すると、12世紀:2人、13世紀:5~7人,14世紀:2人,15世紀:9人,16世紀:27人,17世紀:21人、18世紀:6人となっている。

ガリレオは、大学を中退して、(トスカーナ大公2代目に仕えた)宮廷数学者のOstilio Ricci(この人はタルタルアの弟子)に数学を習ったとされる。1586年に書いた論文は、La Billancetta("little balance"/小天秤)で、アルキメデスの原理の実践的利用法みたいなお話。

これが、大学で数学を教える仕事に結び付くのか不思議だけど、彼は推薦状を書いてもらうために、偉い人に面会して、イエズス会士Christopher Clavius、Guidobaldo、その伝手でトスカーナ大公(3代目)などのコネを得たらしい。最終的に、トスカーナ大公の推薦でピサ大学で教えることになったが、当時、そんな求職方法が、一般的だったわけではないだろう。

Claviusは、イエズス会の中心的教育機関ローマ学院で、幾何学や天文学を教えてた学者で権威とされていた。詳細は知らないけど、イエズス会の数学教育を改革したとか言われている(イエズス会学事規定1599年版)。ガリレオが異端審問にかけられた時、Claviusは既に死んでたけど、生きてたら、宗教裁判を回避できたかは定かでない。Claviusには、Opera mathematicaという著書があって、読んだわけではないけど、内容は、幾何、算術と代数、天文学の話題で構成されているっぽい。アストロラーベに関する著書Astrolabiumは、漢訳『渾蓋通憲図説』が存在するという事実から、よく知られている。

Guidobaldoは、何者かよく分からないけど、Wikipediaによれば、裕福な家の出身で、1564年のBattle of Saint Gotthardに参加した後、Federico Commandino(1509~1575、ヘレニズム時代のギリシャ語文献を何冊もラテン語翻訳した)に師事して、数学を学んだらしい。1588年からは、トスカーナ大公国の"築城監督官"の役職を得ていた。ガリレオがピサ大学をクビになった時も、Guidobaldoの援助があって、パドヴァ大学の職を得たらしい。

そんなわけで、ガリレオも野良数学者みたいな人に師事して、苦労して、大学で仕事を得たっぽい。

Simon Stevinの力学上の仕事は、1586年に出版したDe beghinselen der weeghconst(The principles of weighing)とDe Beghinselen des Waterwichts(The principles of water balance)の2冊の本。Stevinも、独学だったそうだが、やはり誰に何を学んだかは分かっていない。

A first mathematics curriculum : Stevin’s Instruction for engineers (1600)によると、1600年、Simon Stevinは、ライデン大学内で、"Duytsche Mathematique"というコースを開くよう要請されたそうだ。

名前は、数学だけど、目的は、military engineerの教育。カリキュラムは、fortificationがメインのようで、それに必要な算術や幾何学、測量も含まれる。機械の使い方なんかは、記述がないっぽい(?)し、砲術も含まれてない。建築機械については、実地で学ぶ感じだったのかもしれない。デカルトやOtto von Guericke(1602~1686)も、ここで勉強したことがあるらしい。

おまけ2:中世以降の砲工兵教育

17世紀にはヨーロッパ諸国で、砲兵隊、工兵隊が専門分化し、17〜18世紀には、西ユーラシア諸国で、砲工兵学校が作られた。17世紀にも学校は存在したが、正式名称も定かでないものが多い。18世紀のものを、いくつか挙げておく。

  • 1720年のスペインで、Real Academia Militar de Matemáticas y Fortificación(直訳は、Royal military academy of mathematics and fortification)
  • 1739年にサルデーニャ王国のトリノで、Regie Scuole di Artiglieria e Fortificazione(直訳では、Royal School of artillery and fortification)
  • 1741年のイギリスで、WoolwichのRoyal Military Academy(王立陸軍士官学校)
  • 1748年のフランスで、École royale du génie de Mézières(メジエール工兵学校)
  • 1773年のオスマン帝国で、Mühendishâne-i Bahrî-i Hümâyun(Imperial Naval Engineering School)
  • 1795年のオスマン帝国で、Mühendishâne-i Berrî-i Hümâyun(Imperial School of Military Engineering)

Regie Scuole di Artiglieria e Fortificazioneは、若い頃のLagrangeが教師をやってたことがあり、メジエール工兵学校では、Mongeが教えていたことがある。Mongeは、Monge-Ampere方程式に名前が残ってるけど、エコール・ポリテクニークの教育カリキュラム決定に大きく関与した人でもある。

The Teaching of Mathematics at the Royal Military Academy: Evolution in Continuityに、1792年のWoolwich Royal military academyのカリキュラムが列挙されているが、初等的な内容以外に、力学(mechanics)、流率法(fluxions,微積分)、水理学(hydrostatics and hydraulics)、pneumatics、砲術(Gunnery)などが含まれている。

ここのmechanicsは、力学とする方が適切だと思う。詳細は、以下のように書いてある。

Mechanics: Including Motions equable and variable; Forces constant, variable and percussive; Gravity, Sound and Distances; Inclined Planes; Projectiles; Pratical Gunnery; Pendulums; Centres of Gravity; Percussion, Oscillation, and Gyration; Ballistic Pendulum, &c.;

出典は、1892年のRecords of the Royal Military Academy, 1741-1892という本らしいけど、1792年の文言を、どれくらい正確に収録してるのかは知らない。機械工学の方は、1846年には、イギリスで、スチーブンソンが発起人となってInstitution of Mechanical Engineers(英国機械学会)が作られているので、19世紀には、mechanicsとmechanical engineeringが区別された。

18世紀当時は、水理学や弾道学も、「数学者」の研究対象だった。Johann Bernoulliは、1743年に"Hydraulics"という本を出している。弾道学も、例えば、Eulerが研究を行っている。Eulerは、論文を書く以外に、イギリスの砲工兵Benjamin　Robinsによる弾道学の本を翻訳して、詳細な注釈を付けている。BernoulliやEulerが関与していることから推測される通り、これらの分野は、当時としては、かなり数学的なものではある。

‘Mixed’ mathematics in engineering education in Spain: Pedro Lucuce's course at the Barcelona Royal Military Academy of Mathematics in the eighteenth centuryは、スペインのRoyal military academyのカリキュラムについて書いている。おそらく18世紀半ば頃のもので、Curso Mathematico para la Instruccio ́n de los Militares(軍事教育の数学コース)は、大きく8つの主題から構成され、算術、ユークリッド幾何、実用幾何(平面三角法や対数、測量など)、築城術、砲術、cosmography(celestial sphere/天球と書いてるのは、天文学?他に、地理学や海図作成など)、静力学(staticsと書いてるが、内容は、"the motion of heavy bodies"、機械学/machinery、水理学、光学など)、土木建築になっている。各主題は、概ね同程度の分量だったっぽい。まだ微積分はないけど、オイラーのMechanicaが出て10年足らずの時期のカリキュラムなので、まだ十分に、その有用性が知れ渡ってはいなかったのだろう。

mixed mathematicsは、応用数学の古い呼び方で、1800年頃までは、よく使われたらしい。The Evolution of the Term "Mixed Mathematics" には、mixed mathematicsという用語を最初に使ったのは、フランシス・ベーコンだと書いてある。

18世紀は、理論物理学が形成され、そのための数学的道具として、微分方程式論が開発された時代でもある。 一方、d'Alembertは、幾何、算術、代数、微積分学を純粋数学に分類し、力学、天文学、光学、確率論をmixed mathematicsに分類したそうだ。これは、確率論を除けば、現代の認識に近い。確率論は、現在は、数学の一分野と思われてるけど、ヒルベルトの第6問題「物理学の公理化」で、物理学と確率論の公理化をあげてるので、1900年頃でも、微妙な立ち位置にいたのだろう。

経験から完全に切り離して、機械的な論理操作のみで数学ができるという発想が、いつから存在してたのか分からない。少なくとも、算術の公理は19世紀末になるまで存在してない。純粋数学は、日常感覚に任せてた部分を整理する方向に進み、mixed mathematicsに分類されてた分野は、自然科学と融合して、物理学になった。この仕切り直しが、本当に良かったのかは分からないけど、とりあえず、名前を変えてくれなかったせいで、昔の話をする時は、ややこしくなった。

遺伝子の起源

生物が持つ遺伝子は、翻訳されて、数百〜数千アミノ酸のタンパク質になるのが一般的である。数百〜数千残基の長大なタンパク質or遺伝子が、ランダムにアミノ酸やヌクレオチドの重合を繰り返しただけで出現する確率は、普通に考えると極めて低い。ヒトの遺伝子で最も巨大なものの一つはtitin遺伝子で、363個のエクソンで構成され、最大のisoformは36000残基近くになるようである。

一般的に、生物の持つタンパク質は、複数の共通するドメインからなることが知られていて、ヒトの遺伝子数は2〜3万とされている(選択的スプライシングがあるので翻訳されるタンパク質の種類は10万種類くらいになると考えられている)が、1000個程度のドメインの組み合わせでバリエーションを増やしているらしい。また、真核生物、真正細菌、古細菌間でも、共通するドメインの利用が確認されているとされる。

cf)Evidence for PDZ domains in bacteria, yeast, and plants.
https://dx.doi.org/10.1002%2Fpro.5560060225

SMART: Domain List
http://smart.embl-heidelberg.de/smart/domain_table.cgi
を見ると、1400以上のドメインが列挙されている。

多くのドメインは、数十〜数百残基程度のアミノ酸からなり(例として、PDZドメインは、90残基弱で、ホメオドメインは60残基程度、SAMドメインは70残基程度、SH2ドメインは100残基程度、pairedドメインは100残基以上、MAMドメインは170残基程度...)、固有の立体構造を取り、単体でも何らかの機能を有し、ドメインの組み合わせによって、新規遺伝子が作られることもあったと考えられている。が、これでも、まだちょっと長い。

例えば、１Lの水中で、毎秒1mM相当の新規タンパクが生成する(つまり、毎秒6.02×10^20個の新規タンパクが生成される)ような(法外な)イベントがあったとしても、一億年かけても、最大で10^36個オーダーのタンパク質しか試せない。現在の海水(その量は、10^21リットルのオーダーだそうだ)を全部使っても、20の100乗よりは大分小さく、ランダムな重合によって、100残基のタンパクを全部網羅するには到底及ばない。

というわけで、生物が遺伝子ガチャに成功している理由が説明されるべきで、以下のように、いくつかの可能性が考えられる。
①全ての遺伝子は、単一の遺伝子に由来する
②ランダムに生成したアミノ酸配列でも高確率で何らかの活性を持っている:例えば、ドメインの配列も高度に保存された(連続してるとは限らない)短い配列と、それ以外の部分からなることが多々あり、保存された配列以外の部分は、割と何でもいいとか
③ドメインより小さい構造単位が存在する


現在、多くの遺伝子が、遺伝子重複によって、既存の遺伝子を雛形にして生成し、重複遺伝子は元になった遺伝子と独立な変異を経ることで、しばしば機能を変化させてきたと考えられている。ヒトに多くの相同遺伝子(パラログ)が含まれていることが分かっているし、大腸菌や酵母も、パラログ遺伝子を持っていて、これらの事実は遺伝子重複によって説明できる。

1970年に出版された大野進の本Evolution by gene duplicationには、“In a strict sense, nothing in evolution is created de novo. Each new gene must have arisen from an already existing gene.”と書かれていて、これは①の考え方に相当する。当時は知られてなかった遺伝子の水平伝播も、今ではありふれた現象と考えられてるけど、既存の遺伝子を雛形として新しい遺伝子を獲得するという点では共通している。

進化の歴史上、このような方法によってのみ遺伝子の獲得が起きたとすれば、確率的な困難は、最初の遺伝子が誕生した時に一回だけ存在したことになり、一回だけなら、奇跡が起きることもあるだろう。但し、このような奇跡が必要だったなら、地球以外の星に生命体が存在する可能性は極めて低いことになる。


遺伝子重複は、古くから予想されていたことではあるけど、最近は、個別の遺伝子について、かなり詳細な進化の軌跡を辿ることができるようになっている。

遺伝子重複の顕著な例は、ヒトのLオプシン(赤オプシン)とMオプシン(緑オプシン)で、この2つの遺伝子はX染色体上に並んで存在する。SオプシンはX染色体上になく、分岐したのはもっと古いと思われる(塩基配列上も、最も違いが大きい)。現在、原猿や有胎盤類は、SオプシンとL/Mオプシンの2色性色覚のものが多いらしく、真猿類の中だと、新世界ザルのオスは(一部を除いて)2色型色覚で、旧世界ザル、類人猿、ヒトは3色型色覚を持つ。新世界ザルは、3〜4000年前に真猿類の一部がユーラシア大陸から南米に移住して分岐し、アフリカに残った真猿類の一種で、L/Mオプシンの遺伝子重複が起き、この種は、旧世界ザルと類人猿の共通祖先になったと考えられている。新世界ザルの祖先が、海を渡った方法は謎だけど、当時は今より大陸間は近かったと考えられている。

と、よく書いてあるのだけど、ヒトの場合、X染色体上にオプシンを3つ持つ人がいるらしい。
OPN1LW
https://www.ncbi.nlm.nih.gov/nuccore/NM_020061
OPN1MW
https://www.ncbi.nlm.nih.gov/nuccore/NM_000513
OPN1MW2
https://www.ncbi.nlm.nih.gov/nuccore/NM_001048181
OPN1MW3
https://www.ncbi.nlm.nih.gov/nuccore/NM_001330067

OPN1MW2とOPN1MW3は、目で見ると同じもののようなので、ヒトのX染色体上には、OPN1LWとOPN1MW、OPN1MW2の3つのオプシンがある。OPN1MW2は、ある人とない人がいるらしい(?)から、かなり最近になって出現した遺伝子と考えられる(数千年とか数万年とか)。OPN1MWとOPN1MW2は、アミノ酸配列が一致してるが、一部の人は、これによって4色型色覚となるようなので、多分、なんか多型があるのだろう。

チンパンジーには、(私が探した限り)2つしかないっぽい。ヒトでもチンパンジーでもXq28領域に遺伝子がある。
OPN1LW領域(チンパンジー)
https://www.ncbi.nlm.nih.gov/nuccore/CU467112.1?report=fasta&from=144100&to=156324
OPN1MW領域(チンパンジー)
https://www.ncbi.nlm.nih.gov/nuccore/CU467112.1?report=fasta&from=186659&to=199065
ヒトXq28領域
https://www.ncbi.nlm.nih.gov/nuccore/AC092402.2

どーでもいいけど、
チンパンジーX染色体
https://www.ncbi.nlm.nih.gov/nuccore/1348310337
ゴリラX染色体
https://www.ncbi.nlm.nih.gov/nuccore/NC_044625.1
を見ても、OPN1LW/OPN1MWに相同な遺伝子が一個しか見つからないので、アセンブルをミスってるような気がする。

ヒトのOPN1LWとOPN1MWは、共に364aaで、アミノ酸13個の違いしかない。通常、オプシンは、一つの視細胞に一種類しか発現してないらしいけど、4色型色覚の人で、そのへんの発現制御がどうなってるのかは謎。

OPN1MWとOPN1LWの場合は、上流に共通のエンハンサー領域LCRというのがあって、一方の遺伝子のプロモーターを(ランダムに)選択し、一旦行った選択は変化しないらしい(どうやって安定化するのかは分かってないらしい)。ゼブラフィッシュでは、RH2-LCRという領域が4つの下流にあるRH2遺伝子の発現制御に関与しているそうなので、ヒトでも、OPN1MW2が追加されたところで、同種の機構が働いているのかもしれないけど不明要素が多い

OPN1LW,OPN1MW,OPN1MW2とOPN1SWは、エクソンの数が合致していて、個々のエクソンのサイズもよく似ているから、これらの遺伝子が単一の遺伝子に由来するというのは尤もらしい。ヒトのオプシンファミリー遺伝子には、他にも、OPN2/RHOやメラノプシン(OPN4)などの複数の遺伝子があるが、これらではエクソン数が違っている。

OPN1以外のオプシンファミリー遺伝子も、元を辿れば、単一の遺伝子に起源があると信じられている。更に、オプシンは、GPCR(Gタンパク共役型受容体)で、最初のオプシンは、他のGPCR遺伝子から遺伝子重複と変異による新機能獲得(neofunctionalization)で生成された遺伝子かもしれない。ヒトには1000近いGPCR遺伝子があるが、277個のGPCRの系統樹を作成した論文(PMID 12429062)によれば、オプシンファミリーに最も近いのはメラトニン受容体のようである。カイメンには、メラトニン受容体もオプシン遺伝子もないらしいが、また別の論文(PMID23112152)には、平板動物はオプシンを持つと書いてある。

原核生物や古細菌も、バクテリオロドプシンやSchizorhodopsin、Heliorhodopsinと呼ばれる遺伝子を持つが、これらはGPCRではないようだ。レチナールと結合し発色団としている点は共通し、また7回膜貫通型受容体ではあるらしいから、全てのGPCRの起源という可能性はなくもない(?)(7回膜貫通型受容体という用語は、GPCRとほぼ同義で使われる)


他の例として、ホタルのルシフェラーゼ遺伝子は、遺伝子重複とneofunctionalizationで出現し、ついでに、その後、更に遺伝子重複が起きたと考えられるらしい。
cf)発光生物学　2つの最新トピックス
https://doi.org/10.1271/kagakutoseibutsu.57.409


単一遺伝子の重複だけでなく、ゲノム単位の重複も起きると考えられている。Hox gene clusterは、ショウジョウバエでは、8遺伝子からなる単一のクラスターだけど、多くの脊椎動物では、4つのクラスターがあり、それぞれのクラスターは10個前後の遺伝子からなる(進化の中で、一部が脱落したり増えたりしてるけど、並び順は保存されてる)。ヒトやマウスではHoxa,Hoxb,Hoxc,Hoxdと名前がついてて、異なる染色体上(ヒトではそれぞれ7,17,12,2番染色体)にあるけど、これはゲノム単位の重複が2回起きて生じたという仮説がある。クラスター単位で重複して転座しただけという可能性もあるけど、ParaHoxクラスターとかも4つあるようだ。

他にも、哺乳類のFOXP遺伝子は4つあるけど、ハエには一つしかない。2回ゲノム重複があったと考えると、丁度いいケースが結構ある。2010年代後半までの報告では、脊椎動物に近縁のナメクジウオやホヤが単一のHoxクラスターしか持たない一方、無顎類の生き残りである円口類のHox gene clusterやParaHox gene clusterは複数あるものの、どっちも2つとか、両方4つずつというわけではなく、ゲノム重複の起きた時期は、まだ完全には確定されていない。
cf) Chromosome evolution at the origin of the ancestral vertebrate genome
https://www.biorxiv.org/content/10.1101/253104v1
cf) Lampreys, the jawless vertebrates, contain only two ParaHox gene clusters
https://www.pnas.org/content/114/34/9146

ゲノム重複は魚類や両生類でも起こると考えられ、有性生殖を行うフナ類や金魚は、52対104本の染色体を持ち、1400万年前にゲノム重複があったと2019年に報告された(他の事例は、アフリカツメガエル)。また、多くの真骨魚類が、最低8つのHox gene clusterを持つらしく、数億年前に、初期の真骨魚類でゲノム重複が起きた可能性も考えられている。節足動物では、一部のクモやサソリに2つのHox gene clusterが見られるようである。

哺乳類だと、仮に4倍体が生じても、子孫を残せる可能性は絶望的と思われてるけど、21世紀になって、哺乳類でも、現在知られる中で最大の染色体数を誇るメンドサビスカーチャネズミ(学名:Tympanoctomys barrerae)及び近縁のゴールデンビスカーチャネズミ(学名:Pipanacoctomys aureus)というマイナーな齧歯類が、4倍体の可能性が議論されている。
Discovery of tetraploidy in a mammal
https://www.nature.com/articles/43815

Molecular cytogenetics discards polyploidy in mammals
https://doi.org/10.1016/j.ygeno.2004.12.004

Evolution of the Largest Mammalian Genome
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5569995/




ヒトの持つ遺伝子の多くが、重複遺伝子として生まれて進化してきたことは確からしいが、全部の遺伝子がそうだと証明されてるわけではない。21世紀になってから、既存遺伝子を雛形としない遺伝子(de novo gene)の生成があるという議論がされるようになった。orphan geneという言葉もあるけど、orphan geneは、遺伝子の水平伝播で獲得した遺伝子でも、そう見えることがあるのに対して、de novo geneは、DNAのnon-coding sequenceから生成されたものと定義される。

de novo gene birthという名前が付いてるけど、ヒト遺伝子の中で、de novo geneだと特定されてるものは今の所ないし、de novo gene birthのメカニズムもよく分かってない。

de novo geneに関する初期の論文として、2006〜7年頃の以下の論文が挙げられる。

Evidence for de novo evolution of testis-expressed genes in the Drosophila yakuba/Drosophila erecta clade
https://doi.org/10.1534/genetics.106.069245

論文中に、We refer to this class of orphan genes as “de novo,” to suggest the possibility that they may derive from ancestrally noncoding sequence.と書かれていて、de novoという接頭辞を使おうと提案したのは彼らっぽい。

彼らは、分子生物学でよく使われるモデル生物Drosophila melanogasterの近縁種Drosophilia yakubaの精巣で発現しているmRNAを解析し、D. melanogaster, D. erecta, D. ananassaeなどの近縁種に相同遺伝子の存在が確認できないD. yakuba固有の遺伝子を探索した。これらの種は、Drosophila melanogaster species subgroupもしくは、Drosophila melanogaster species groupに分類されている。species subgroupとかspecies groupはハエ業界以外ではあんまり使われてないっぽいけれど、属と種の中間の分類階級らしく、とりあえず、Drosophilia melanogasterに最も近縁な種たちを調べたということ。

Wikipediaによると、Drosophilia yakubaに最も近縁なのが、D. santomeaで、次が、D. teissieriとなっている。
Drosophila melanogaster species subgroup
https://en.wikipedia.org/wiki/Drosophila_melanogaster_species_subgroup

こうして見つかった遺伝子は、2〜4個の比較的少数のエクソンからなると予測されている。論文では、これらの遺伝子が、本当に、non-coding sequenceから出現したものであることが証明されてるわけではないし、タンパク質に翻訳されてない可能性もあるが、これらの遺伝子がタンパク質をコードしているならば、そのサイズは、100aa前後で、かなり小さく、突然、巨大な遺伝子が出現する可能性は低いという直感とは整合的である。


2012年の以下の論文では、酵母のORFの内、遺伝子としてアノテーションされてないものを調べている。
Proto-genes and de novo gene birth
https://www.nature.com/articles/nature11184

ゲノム配列からのORFのアノテーションでは、最初に、300塩基を閾値として取捨選択されていて、約6000の遺伝子に対して、約261000のアノテーションされてないORFがあったそうである(過去形だけど、多分、今も状況はそんなに変ってない)(個別に機能が同定された遺伝子の中には、閾値より小さいものもあると思われる)。論文では、長さが30塩基以上のアノテーションされてないORFについてタンパク質に翻訳されてるかどうかを調べた所、108000のORFの内、1139個について翻訳されている証拠を得たとしている("We found that 1,139 of 108,000 ORFs0 show such evidence of translation")。

割合としては1%だけど、そもそもプロモーターがなければ転写もされないので、少ないとも言い切れない(どれくらい転写されてるのかは分からないが)。実験方法としては、リボソームプロファイリングと呼ばれてる方法で、2009年に提案され、次世代シーケンサーを使うので、20世紀には出来なかった実験。但し、翻訳されているというだけで、それが何かの機能を果たしているかは分からない。

"weakly conserved annotated ORF"とアノテーションなしのORFでは、翻訳が確認されるものが1256個、300塩基以下の長さのものは1133個ある(Fig4a及び、aupplementary information fig8 middle)らしいので、100アミノ酸以下の割合が、かなり高い。


2018年の以下の論文は、ちょっと毛色が違って、以前から知られていた遺伝子がde novo geneであると主張している。
De Novo Gene Evolution of Antifreeze Glycoproteins in Codfishes Revealed by Whole Genome Sequence Data
https://dx.doi.org/10.1093%2Fmolbev%2Fmsx311

不凍タンパク質と呼ばれる凍結防止の機能を持つ遺伝子は、様々な生物で見つかっていて、独立に何度も獲得されたと考えられている。論文では、タイセイヨウダラ(Gadus morhua)やコダラ(Melanogrammus aeglefinus)が持つ不凍タンパク質AFGPが、1300〜1800万年前にde novo gene birthで獲得したと推論している。2019年の論文

Molecular mechanism and history of non-sense to sense evolution of antifreeze glycoprotein gene in northern gadids
https://doi.org/10.1073/pnas.1817138116

も、同様の結論だけど、こっちの方が読みやすい(結論だけ見たければ,Fig4)。Fig2(Boreogadus saidaという種のafgp遺伝子)を見ると、3つのエクソンからなり、冒頭2つのエクソンはシグナルペプチドで、3つ目の最も大きな本体は、TAAが繰り返し並ぶ特異な配列をしている。afgp遺伝子は、そこそこ大きいけど、例外的ケースじゃないかと思う。


2019年の以下の論文は、現存する自然界の遺伝子の話ではなく、ランダムに生成した人工的なsmall ORF(長さは10~50aa)によって、大腸菌が抗生物質耐性を獲得しうるという実験。
De Novo Emergence of Peptides That Confer Antibiotic Resistance
https://doi.org/10.1128/mbio.00837-19
10の9乗に満たない種類の配列から、3つの候補遺伝子が得られたそうだ。あくまで、抗生物質耐性の有無でスクリーニングしているだけなので、大腸菌にとって有用な機能を果たす遺伝子は、もっと多く含まれている可能性もある。

ヒトに、de novo geneが、どれくらいあるかは分からないが、現在のところ、小さい遺伝子は見過ごされてる可能性も高いので、意外と多い可能性もある。酵母の例を見ると、既知の遺伝子6000に対して、アノテーションされてないが翻訳されるORFが1000以上あり、これらの内、どれくらいが何らかの機能を果たしているのか分からないが、全遺伝子の一割くらい、de novo geneという可能性もある。



ドメインより小さい構造単位が存在するかどうかについて、以下のような論文を書いてる人たちがいる。

Centripetal modules and ancient introns
https://doi.org/10.1016/S0378-1119(99)00292-9

タンパク質の3次元構造に基づいて、各アミノ酸残基に対して定義される"F-value"(とは呼んでないが、勝手に命名した。周辺アミノ酸のα炭素との距離の二乗の平均)の極小点を境界として、centripetal moduleというものを定義して、moduleの境界は、フェーズ0のエクソン・イントロン境界(エクソンの始まりや終わりが、コドンの途中になってないような境界)と、よく相関すると主張する論文。

centripetal moduleの概念は、1987年に提案されたものらしいが、この論文では、module境界を計算するプログラムを書いて、統計的に相関があるのを確認したということらしい。プログラムの詳細は説明されてないし、公開されてもいない。

彼らは、ドメインについては特に何も言ってないけど、centripetal moduleは、直感的には、空間的にコンパクトにまとまっている単位に分割しようという発想で、ドメインも、しばしばコンパクトにまとまっているとされているので、ドメイン境界を含むものになるだろうと予想される。


色々疑問はあるけど、例えば、ヒトのピルビン酸キナーゼM2(PKM2)で、エクソン・イントロン境界とドメイン境界に、何か関わりがあるか見てみる。

PKM2は、グルコースを2つのピルビン酸に分解する解糖系に於いて、ホスホエノールピルビン酸からADPにリン酸を転移し、ピルビン酸を生成する過程を触媒する酵素である。ヒトは、複数のピルビン酸キナーゼを持つが、その内の一つ。

PKM2(PKM isoform a) domains
http://atlasgeneticsoncology.org/Genes/GC_PKM.html

によれば、PKM2は、4つのドメインからなり、A-domainは、配列上は2つに分断されていて、N-domain (aa 1-43), the A-domain (aa 44-116 and 219-389), the B-domain (aa 117-218) and the C-domain (aa 390-531)らしい。

pkm2遺伝子のエクソンは、
Homo sapiens pyruvate kinase M1/2 (PKM), transcript variant 1, mRNA
https://www.ncbi.nlm.nih.gov/nuccore/NM_002654.6
によれば、11個のエクソンからなるようで、その境界は、タンパク質のドメイン境界と、特に関係してそうには見えない。


PKM2は、解糖系に関わる酵素なので、その進化的起源は、物凄く古い可能性がある。

ヒトPKM2の構造は、
Functional cross-talk between allosteric effects of activating and inhibiting ligands underlies PKM2 regulation
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6636998/
のFig1(A)で見られ、

Evolutionary plasticity in the allosteric regulator-binding site of pyruvate kinase isoform PykA from Pseudomonas aeruginosa
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6802521/
のFig3Bに緑膿菌のピルビン酸キナーゼPykAの構造があるが、目で見る限り、ヒトPKM2に似ている。緑膿菌PykAのA-domainには、高度に保存された配列MVARGDLGVEがあり、触媒活性に重要な役割を果たすと書いてある。ヒトのPKM2のA-domainには、似た配列MVARGDLGIEが含まれている(触媒活性に関与するのかは知らない)


それはそれとして、PKM2のcentripetal moduleを調べることにする。構造データとして、
https://www.rcsb.org/structure/1T5A
を使用して、計算した。横軸は、残基で、縦軸は、F-value。kは、周辺アミノ酸何残基までを計算に入れるか決めるウインドウ幅(の半分。論文では、k=40~90がいいとか書いてるので、同じような範囲で計算した)。
f:id:m-a-o:20201213154649p:plain
は、平滑化なしの値をプロットしたもので、既知のドメイン境界も、重ねて書いてある。

PKM2のドメインは、目で見ても、何となく分かるので当然のことではあるが、ドメイン境界が極小点になってそうだとは分かる。この例では、極小点の中でも、ドメイン境界になっているものは、F-valueがウインドウ幅kの選び方に依存してないようである。いつでも、こんな風にドメイン境界が採れるという保証はない。

生データだと、極小点が多すぎるので、8残基移動平均をとって平滑化して、極小点を抽出してみた。この例だと、あまり、エクソン・イントロン境界と相関してそうに見えない。"+"が付いてるのは、フェーズ0で、"-"が付いてるのは、フェーズ0じゃないもの。
f:id:m-a-o:20201213154709p:plain
平滑化したグラフに、ドメイン境界を重ねてみる
f:id:m-a-o:20201213154712p:plain

もうちょっとencourageされる結果だったら、他の例も見てみようと思ったけど、大変微妙な結果。今の所、10〜20残基程度の短い構成単位が存在するという証拠はなさそうに思える。


描画に使用したコード

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal

def readcif(filepath):
   ret = []
   nprev = None
   for line in open(filepath):
       line = line.strip()
       if len(line)<5:continue
       if line[:5]!="ATOM ":continue
       ls = line.strip().split()
       if len(ls)<10:continue
       if ls[0]!="ATOM":continue
       if ls[3]!="CA":continue  #-- not alpha carbon
       if ls[4]!="." and ls[4]!="A":continue
       if ls[6]!="." and ls[6]!="A":continue
       aa = ls[5]
       n = int(ls[8])
       x,y,z=float(ls[10]),float(ls[11]),float(ls[12])
       assert(nprev==None or nprev+1==n),(nprev,n)
       ret.append( (aa,x,y,z) )
       nprev=n
   return ret


def getseq(r):
   dic = {'Ala': 'A', 'Asx': 'B', 'Cys': 'C', 'Asp': 'D', 'Glu': 'E', 'Phe': 'F', 'Gly': 'G', 'His': 'H', 'Ile': 'I', 'Lys': 'K', 'Leu': 'L', 'Met': 'M', 'Asn': 'N', 'Pro': 'P', 'Gln': 'Q', 'Arg': 'R', 'Ser': 'S', 'Thr': 'T', 'Sec': 'U', 'Val': 'V', 'Trp': 'W', 'Xaa': 'X', 'Tyr': 'Y', 'Glx': 'Z'}
   def normalize(s):
      s= s.lower()
      return s[0].upper()+s[1:]
   return "".join([dic[normalize(x[0])] for x in r])



def getF(r , k=50):
   ret = []
   for i,(aa,x,y,z) in enumerate(r):
       j1 = max(0 , i-k)
       j2 = min(len(r)-1 , i+k) + 1
       assert(j2>j1),(j1,j2)
       Fval = sum([(x-x2)**2+(y-y2)**2+(z-z2)**2 for(_,x2,y2,z2) in r[j1:j2]])/(j2-j1)
       ret.append(Fval)
   return ret


def phase(n):
  if n%3==0:return "+"
  else:return "-"


if __name__=="__main__":
  r = readcif("1t5a.cif")
  Nlen = 531  #-- PKM2 has 531aa
  offset = Nlen - len(r) + 1
  xtick = "module"
  smoothing = True
  for k,c in [(40,'r'),(50,'brown'),(60,'grey'),(70,'g'),(80,'b')]:
    d = getF(r,k=k)
    x = np.linspace(offset ,Nlen, len(d))
    y = np.array(d)
    if not smoothing:
      plt.plot(x,y,c,label='k={0}'.format(k))
    else:
      num=8  #移動平均の個数
      b=np.ones(num)/num
      y2=np.convolve(y, b, mode='same')#移動平均
      plt.plot(x,y2,c,label='k={0}'.format(k))
      minpos = signal.argrelmin(y2, order=3)
      plt.plot(x[minpos],y2[minpos],c, marker='o',linestyle='None')
    if xtick=="exon":
       exon_ticks=[((c-89)/3,str(n+3)+phase(c-89)) for (n,c) in enumerate([243,335,467,654,925,1076,1229,1396,1578])]
       plt.xticks([c[0] for c in exon_ticks],[c[1] for c in exon_ticks])
    elif xtick=="module":
       plt.xticks([44,117,219,390])
    plt.yticks([])
    plt.grid(True)
    plt.xlim(0,Nlen)
  plt.legend()
  plt.show()

対数共起頻度のSVD

NLP

対数共起頻度を用いた四項類推：word2vecとPMI との比較
https://doi.org/10.11517/pjsai.JSAI2020.0_4Rin177

という報告を読んだ。四項類推は、queen-woman+man ≒ kingみたいなタイプの類推問題。word2vecは、この手の問題によく答えられるということになっている。

四項類推の精度については、通常、最も類似度の高い単語を1〜数個取ってきて、妥当性を確認するという方法が取られ、類似度の値そのものは見てないことが多いけど、ちょっと類似度の値を見てみる。
Pretrained Embeddings
https://wikipedia2vec.github.io/wikipedia2vec/pretrained/

で公開されているword2vecで学習したenwiki_20180420_300dに含まれるデータでは、kingとqueenのコサイン類似度は、元々、0.63ある。kingもqueenも、最上級の支配者を表すという点では、共通していて同じような文脈で出現することが多いと考えられるので、それ自体は問題ない。しかし、king-man+womanとqueenのコサイン類似度は、0.65になるが、単にkingとqueenが、最初から類似度の高い単語であることが効いてるだけのように思える。類似度が上がってるように見えるのも多分偶然の産物で、queen-woman+manとkingのコサイン類似度は、0.60になって、むしろ小さくなる。

しかし、四項類推について全く意味のある情報を含んでいないとも言えない。(king,queen),(man,woman),(father,mother),(god,goddess),(boy,girl),(boys,girls),(son,daughter),(prince,princess),(waiter,waitress)という9組について、差のベクトルの平均を取って、(actor,actress)の差ベクトルとのコサイン類似度を計算すると、0.82という高い値になる。個別の差ベクトルとのコサイン類似度は、0.43〜0.69までバラバラであるが、平均化することで類似度が上がるので、"性差"を表現するベクトルが存在している可能性もある。

actorとactressのコサイン類似度は、元々、0.78ほどあり高いのだけど、actor-"性差"ベクトルとactressのコサイン類似度は、0.92とかなので、0.4~0.7くらいの胡散臭い数値を、うろうろしてるのとは違うようにも思える。個別のペアに関しては、性差以外の部分でも、用法の差があるのかもしれない(例えば、日本語で、"ボーイ"を給仕的な意味に使うことがあり、英語圏でも、この用例は存在するっぽいが、girlで、これに相当する使い方は、多分しないっぽいとか、そういった感じのこと)が、平均化することで、こうした微細な違いが消失するのかもしれない。

gloveの場合も、githubレポジトリのREADMEにリンクがあるglove.42B.300dで実験すると、同じようなことが起きる。例えば、上と同様にして"性差"ベクトルを計算して、actor-"性差"ベクトルとactressのコサイン類似度を見ると、0.89程度になっている。he-"性差"ベクトルとsheのコサイン類似度は、0.93以上ある。
stanfordnlp/GloVe
https://github.com/stanfordnlp/GloVe

gloveに同梱されてるdemo.shでは、小さなコーパスtext8をダウンロードして学習するようになっているけど、text8で学習したベクトルを使うと、同じことをやっても成績は悪いので、大きなコーパスを使うことは重要かもしれない。text8は、Wikipediaを元に作成したコーパスから100MBを抜粋したデータで、約1700万tokensを含む。glove.42B.300dは、420億トークンからなるコーパスで学習しているらしい。Wikipedia+αから作成したコーパス(60億トークン)で学習したglove.6B.300dでも、学習精度は、glove.42B.300dと大差ないように思える。

これらのコーパスは、通常の本などと比べて桁違いに大きく、例えば、King James版英訳聖書は、延べ80万tokens弱(unique wordsは約5400単語)らしい。
単語の分散表現を使った教師なし単語翻訳
https://m-a-o.hatenablog.com/entry/2019/04/01/221230
の末尾で、King James版英訳聖書で学習した分散表現がイマイチっぽいと書いたけど、本1〜2冊程度で学習できる分散表現は、(少なくとも、現在のアルゴリズムでは)やはりカスなんじゃないかと思う


で、冒頭の報告は、対数共起頻度の(trunated) SVDを使っても、四項類推問題で、それなりの成績を収めるということらしい。多くの点で、以下の2015年の論文を参照している。
Improving Distributional Similarity with Lessons Learned from Word Embeddings
https://www.aclweb.org/anthology/Q15-1016/

論文には、(LSAとかで古くから行われてたような)old-styleの"count-based"な手法でも、prediction-basedな手法(≒neuralなんとか)に近い性能を実現することができる的なことが書いてある。対数共起頻度も、count-based modelではあるけど、以下の論文では実験されてないので、それをやったのが冒頭の報告ということ。

対数共起頻度をSVDするだけというのは、理屈上は簡単なので自分でも試してみることにした。上にも書いた通り、コーパスサイズが小さいとゴミになる可能性があるので、報告と同じく、英語Wikipedia全文を利用することにする(少なくとも、word2vecやgloveは、これで、そこそこの性能が出てるらしいし)。報告では、メモリの都合が〜と言い訳しつつ、なんか回りくどいことをしているが、メモリを潤沢に利用できる環境を用意して直接計算した。

Wikipediaの前処理には、fasttextのレポジトリに含まれるwikifil.plを使用して、以下のように、コーパスを作成した。

wget https://raw.githubusercontent.com/facebookresearch/fastText/master/wikifil.pl
curl -L -O http://dumps.wikimedia.org/enwiki/latest/enwiki-latest-pages-articles.xml.bz2
bzcat enwiki-latest-pages-articles.xml.bz2 | perl wikifil.pl > enwiki_all.txt

次に、共起頻度の計算には、gloveで使われているcooccurを流用。デフォルトの設定では、distance-weightingという距離による重み付けが入ってたりするけど、オプションで切ることができる。
GloVe
https://github.com/stanfordnlp/GloVe

gloveのレポジトリに含まれてるdemo.shを改変して、以下のようなシェルスクリプトを使用した。メモリの設定は64GBとかにしてるので、環境に合わせて変更しないとメモリを使い切るかもしれない。shuffleとgloveは、やらなくてもいいけど、ついで。

#!/bin/bash
set -e

make

CORPUS=enwiki_all.txt
VOCAB_FILE=vocab.txt
COOCCURRENCE_FILE=cooccurrence.bin
COOCCURRENCE_SHUF_FILE=cooccurrence.shuf.bin
BUILDDIR=build
SAVE_FILE=glove.raw.300d
VERBOSE=2
MEMORY=64.0
VOCAB_MIN_COUNT=100
VECTOR_SIZE=300
MAX_ITER=15
WINDOW_SIZE=10
BINARY=2
NUM_THREADS=24
X_MAX=10

echo "$ $BUILDDIR/vocab_count -min-count $VOCAB_MIN_COUNT -verbose $VERBOSE < $CORPUS > $VOCAB_FILE"
$BUILDDIR/vocab_count -min-count $VOCAB_MIN_COUNT -verbose $VERBOSE < $CORPUS > $VOCAB_FILE
echo "$ $BUILDDIR/cooccur -memory $MEMORY -vocab-file $VOCAB_FILE -verbose $VERBOSE -window-size $WINDOW_SIZE -symmetric 1 -distance-weighting 0 < $CORPUS > $COOCCURRENCE_FILE"
$BUILDDIR/cooccur -memory $MEMORY -vocab-file $VOCAB_FILE -verbose $VERBOSE -window-size $WINDOW_SIZE -symmetric 1 -distance-weighting 0 < $CORPUS > $COOCCURRENCE_FILE
echo "$ $BUILDDIR/shuffle -memory $MEMORY -verbose $VERBOSE < $COOCCURRENCE_FILE > $COOCCURRENCE_SHUF_FILE"
$BUILDDIR/shuffle -memory $MEMORY -verbose $VERBOSE < $COOCCURRENCE_FILE > $COOCCURRENCE_SHUF_FILE
echo "$ $BUILDDIR/glove -save-file $SAVE_FILE -threads $NUM_THREADS -input-file $COOCCURRENCE_SHUF_FILE -x-max $X_MAX -iter $MAX_ITER -vector-size $VECTOR_SIZE -binary $BINARY -vocab-file $VOCAB_FILE -verbose $VERBOSE"
$BUILDDIR/glove -save-file $SAVE_FILE -threads $NUM_THREADS -input-file $COOCCURRENCE_SHUF_FILE -x-max $X_MAX -iter $MAX_ITER -vector-size $VECTOR_SIZE -binary $BINARY -vocab-file $VOCAB_FILE -verbose $VERBOSE

無事計算を終えると、cooccurrence.binに共起頻度が入っている。(symmetricオプションを有効にしたため)対角成分が、2倍になってたりするけど、どうせ対数を取るので、大きな問題はないと考えて、そのままにしてある

私がやったところ、共起頻度行列の次元は40万弱で、非ゼロ要素の数は、17.8億弱になった。非ゼロ要素の半分近くが1なので、対数を取ると0になるけど、それでも、疎行列を作るだけで10GB以上のメモリが要る(cooccurrence.binのサイズは28GB強)。疎行列のSVDには、RedSVDを使った(RedSVDの実装が正しいのかチェックしてないけど信じることにする)。使用した実装が、Eigenに依存してたので、Eigenも必要。
Eigen
http://eigen.tuxfamily.org/index.php?title=Main_Page

header-only version of RedSVD
https://github.com/ntessore/redsvd-h

実際のプログラムは、以下の通り。UとSを別々に出力してるけど、冒頭の報告では、UとSの平方根を掛けたもの(rootSVDと呼んでる)が一番成績がよいらしい。これは、後から計算する。U単体を、word vectorと見たものは、uSVDと呼んでいる。

#include <cmath>
#include <vector>
#include <iostream>
#include <Eigen/Sparse>
#include <RedSVD/RedSVD.h>

typedef struct cooccur_rec {
    int word1;
    int word2;
    double val;
} CREC;

int main(){
    std::vector< Eigen::Triplet<float> > triplets;
    FILE *fin = fopen("cooccurrence.bin","rb");
    fseek(fin , 0 , SEEK_END);
    size_t sz = (size_t)ftell( fin );
    triplets.reserve( sz / sizeof(CREC) );
    fseek(fin , 0 , SEEK_SET);
    int D = 0;
    while(!feof(fin)){
        CREC rec;
        fread(&rec, sizeof(CREC), 1, fin);
        if ( D < rec.word1 ) D = rec.word1;
        if ( D < rec.word2 ) D = rec.word2;
        if( rec.val > 1.0) {
            triplets.push_back( Eigen::Triplet<float>(rec.word1-1 , rec.word2-1 , logf((float)rec.val)) );
        }
    }
    fclose(fin);
    
    Eigen::SparseMatrix<float> M;
    M.resize(D,D);
    M.setFromTriplets(triplets.begin() , triplets.end());
    triplets.resize(0);
    triplets.shrink_to_fit();

    RedSVD::RedSVD<Eigen::SparseMatrix<float>> svd(M,300);
    auto U = svd.matrixU();
    auto S = svd.singularValues();
    std::cout << "S:" << std::endl << S << std::endl;
    std::cout << "U:" << std::endl << U << std::endl;
    return 0;
}

計算は、gloveより大分早かった(gloveは複数スレッド使って頑張ってるらしいのに数時間、rootSVDやuSVDは1スレッドで10分くらい)。後は、出力を適当に整形すれば、分散表現が得られる。で、四項類推とか以前に、そもそも、(king,queen),(man,woman),(father,mother),(god,goddess),(boy,girl),(boys,girls),(son,daughter),(prince,princess),(waiter,waitress)という9組で、単語のコサイン類似度がどうなってるのか見ると、以下のようになった。

名詞1 名詞2 rootSVD uSVD word2vec glove.42B.300d glove.text8
king queen 0.69 0.40 0.63 0.76 0.52
man woman 0.70 0.30 0.69 0.80 0.57
father mother 0.89 0.69 0.81 0.82 0.65
god goddess 0.64 0.45 0.61 0.47 0.42
boy girl 0.85 0.48 0.77 0.83 0.53
boys girls 0.84 0.65 0.86 0.86 0.54
son daughter 0.86 0.67 0.84 0.82 0.64
prince princess 0.79 0.64 0.68 0.66 0.50
waiter waitress 0.85 0.65 0.72 0.78 0.11

word2vecは、上の方で触れたenwiki_20180420_300dを使った時の数値。glove.text8は、gloveに同梱されてるdemo.shの学習結果で、上に書いたようにコーパスが小さい。雰囲気としては、uSVDやglove.text8は、他よりコサイン類似度が小さい傾向にあるように見える。

rootSVDでは、king-man+womanとqueenのコサイン類似度は、0.64で下がる。king-male+femaleとqueenのコサイン類似度も、0.67で、元々のkingとqueenの類似度が効いてるだけのように見える。性差ベクトルを計算して、king-"性差"ベクトルとqueenのコサイン類似度は、0.75になり改善はする。

網羅的なテストはしてないけど、rootSVDが、四項類推課題に対して、そこそこの性能を持つという報告は正しそうに見える。ただ、それは、元々、類似度の高い単語を取ってるだけでないかという疑惑はある。報告を見る限り、skip-gramよりは多少精度が悪そうなのだけど、共起頻度のみを使う場合、共起単語の出現位置に関する情報(the kingとking theという単語列を区別できない)は完全に捨てられていて、入力の時点で、情報が少ないので、当然という気はする。SVDベースの方法で、このような情報を取り込むには、前方共起頻度と後方共起頻度を別にカウントするとか、もっと詳細には、窓幅4以下での共起と窓幅8以下での共起を分けるとか、複数のcontextを同時に使う方法が思いつくけど、実際有効かどうかは不明

NVIDIA Thrustで分子動力学(1)単原子分子集団のMD

MDは時々使うけど、自分で実装したことがなかったので、一回くらいは、やっておくべき気がした。

実装
https://gist.github.com/vertexoperator/1d2b4c6b71138d59484665cb524143ac

最初なので、相互作用はLennard-Jonesポテンシャルのみで、周期境界条件、NVE(粒子数、体積、総エネルギー)一定という簡単なやつ。粒子の軌跡を追跡するだけでは詰まらないので、任意時点での温度と圧力を計算できるようにしてある。

普通に、全粒子間のLJ相互作用を全部計算すると、粒子数Nに対して、O(N^2)の計算量が必要になるので、分子動力学では、適当なカットオフ距離を決めて、カットオフ距離内の粒子との相互作用のみ計算することが多い(バッキンガムポテンシャルを使う場合も同様の方法が使えると思うけど、クーロン相互作用を計算する場合は、LJ相互作用より減衰しにくいので、ももう少し工夫がいる)。この場合、計算全領域をカットオフ距離幅で分割することで、計算量をO(N)に減らせる。

実装に難しいところはないけど、生のCUDAカーネルを直接書かないという目標は、一箇所だけ、いい方法が思いつかなくて達成できなかった(sortされた整数の配列がある時に、各整数の開始indexを計算する処理で、scan_by_key→unique_by_key→scan→gatherとかで出来るけど、これだけのことに、4つも関数を呼びたくない)。力の計算も、なんか、thrustの正しい使い方じゃない気もするけど、一応正しく、それなりの性能で動いてる。

時間積分は、velocity Verletだけど、今の所、知られてる時間積分では、どの方法であっても、長時間計算すると、エネルギーが発散するという問題は避けられないらしい。エネルギーを保存する時間積分法が欲しいところではある。

カットオフ距離2nm、タイムステップ1fsで、うちのマシンで適当に性能を測定すると、
倍精度、N=100万:1ステップ20msくらい
倍精度、N=1000万:1ステップ140msくらい
単精度、N=100万:1ステップ10msくらい
単精度、N=1000万:1ステップ70msくらい
という感じ。計算は、常温・約8気圧で、カットオフ距離は2nmと、やや大きめにしたものの、単精度と倍精度で、性能比が、ほぼ2：1なことから、残念な感じが察せられる。

力の計算の演算数を数えると、一粒子あたり
・3個の除算（最適化すれば、3個の乗算)と3個のfloor
・8〜11個の加減乗算が近傍グリッドにいる粒子の個数回
・(25個の加減乗算+1個の除算)がカットオフ距離内にいる粒子の個数回
となる。

100万粒子の計算では、系の一辺の長さを、約167nmに設定しているので、粒子が一様に分布していれば、近傍グリッドにいる粒子の個数は、平均46個。カットオフ距離内にいる粒子は(自分自身を除いて)6個強になる。使ったGPUの倍精度演算性能は420Gflopsと公表されていて、100万粒子に対して、420Gflops出れば、これくらいの演算をやるのに、2msあれば十分というくらい。

近傍グリッドは27個あるので、平均すると、1グリッドあたり、1〜2個の粒子しか入ってない。より計算需要があるだろう超高圧の気体や液体になれば、粒子の密度は、数千倍にもなって、一つのグリッドにいる粒子数が増大するので、多分、性能は向上する。希薄気体であっても、全ペアの計算をするよりは、ずっとマシ(100万粒子に対して、全ペアの計算をすれば、10Tflopsの演算性能があっても、1ステップ1秒以上は必要になるはず)なので、こんなんでもいいだろう。


一番はまったのは、初期位置を、単純に一様分布で生成すると、かなり近い位置に粒子が落ちることが結構あって、そうすると、運動エネルギーが普通でも、総エネルギーが異常に大きくなって困る。こういう希薄な気体でも、こういう問題が起きるのは、意外だった。

もう一つ、本来は初期化でやるべき気がすることとして、系の運動量及び角運動量を0にすることがあるが、やっていない。(初期位置を決めれば)運動量と角運動量が0という条件は、速度に関する6個の線形な拘束条件になるので、この拘束条件が6x3N行列Cで書けるとすると、I-C^{T}(CC^{T})^{-1}Cを掛ければ、運動量と角運動量が0の部分空間に落とせる(射影すると、速度分布は、Maxwell-Boltzmann分布からは、ずれると思われる)。LJポテンシャルは方向に依存しない二体相互作用なので、運動量と角運動量は(原理的には)保存するはずである。


NVE条件で計算すると、温度や圧力は別に保存量じゃないので一定値にはならないが、系が平衡状態にあれば、(相転移点の近傍にいない限りは)ある平均値の周りで変動するものと思われる。日常的な感覚に基づく推測でも、外部とエネルギーのやり取りがない孤立系で、突然、温度や圧力が不安定化することは、あんまりないと思われる。けれど、所与の温度や圧力が、どのようなNVE条件で成立するか事前に知ることは、通常、そんなに簡単じゃない。普通、計算したい物理量は、特定の温度や圧力条件下における値なので、これは、少し不便。

なので、分子動力学でも、温度や圧力を制御する方法が色々と考えられている。これらの方法の多くは、1980年台に提案されたものらしい。Nose-Hoover thermostat, Andersen barostat, Berendsen thermostatなどがよく使われるが、実装は、別に難しくない。

何らかのthermostatを使ったNVT計算を、統計力学に於けるカノニカル・アンサンブルの実現と見なすようなことは理論的根拠がないので、これらのbarostat/thermostatは、所与の温度、圧力に対応する密度とエネルギーを発見するために使われる(この値が間違ってないことは理論的に保証されないと思うので、NVE条件による"production run"で、温度や圧力はチェックすべきである)。

生体高分子の計算をやる場合は、NVE一定条件でやると、タンパク質のフォールディングによって温度が変わることがある(十分多くの水分子を用意すれば、系の温度変化は極めて小さくなるだろうが、計算量が増える)ようなので、production runを、NVE条件でやらずに、上に書いたようなthermostat/barostatを使っているのを見ることがあるかもしれない。


実用的な分子動力学に必要な一般化として、多原子分子の計算もある。(古典)分子動力学では、電子の自由度を無視しているので、多原子分子の扱いは、色々無理があるけど、電子の自由度を考慮するには、量子論的な計算が必要になって計算量が爆発的に増えるので、古典分子動力学は、色んな業界で使われている。

電子の自由度を無視した結果として、分子内の相互作用は、二体相互作用の形ではなく、多体相互作用の形でモデル化されることも、しばしばある。こういう相互作用を許すと、運動量や角運動量の保存則が壊れる気がする(そもそも、原子核だけでは、これらの量は保存しないので問題ないとも言える)けど、計算結果は、それなりに有用なものとして扱われる。分子内の相互作用の計算は、実装としては、そんなに難しいものでもない。

分子間の相互作用は、Lennard-Jonesポテンシャルに加えて、Coulomb力が考慮されることが多い。Coulomb力は、LJポテンシャルより、減衰しにくいので、LJポテンシャルのように単純に小さなカットオフを使って遠くの粒子からの影響を切り捨てるわけにはいかないし、それどころか、Coulombポテンシャルは周期境界条件を考慮して、周期対称性を持つ形にしなければならない。つまり、1/r^2ではなく、\sum_{\mathbf{k} \in L} 1/(\mathbf{r} + \mathbf{k})^2の形を取る(Lは周期格子)

Coulombポテンシャルを考慮した計算も、O(N^2)で良ければ、定義どおりに実装してみることは、特に難しくない。計算量を減らすアルゴリズムとしては、particle Mesh Ewald法というものが、よく使われて、O(N log N)になる。最近は、fast multipole method(FMM)を使った実装が増えてきていて、こっちは、O(N)になる。

cf)FMMについては、2020年に実装しましたというだけの論文が出てる程度には枯れてない話題
A GPU-Accelerated Fast Multipole Method for GROMACS: Performance and Accuracy
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7660746/

分子内と分子間の相互作用をモデル化する方法は多分色々提案されているが、特に、水分子については、専用のモデルがいくつも作られている。有名なものとしては、SPC,SPC/E, TIP4P,TIP5Pなどがある。Wikipediaで
水モデル
https://ja.wikipedia.org/wiki/%E6%B0%B4%E3%83%A2%E3%83%87%E3%83%AB
を読めば、大体、どういうものかはわかる。

次は、とりあえず、水の分子動力学計算を実装して、いくつかの物理量を計算して、水モデルによる違いでも見ることにする。

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トン相当と結論付けたのかもしれない。

Tales of Tetranitrogenia

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

f:id:m-a-o:20200531182824p:plainf:id:m-a-o:20200531182850p:plain

周期表で、窒素の下にあるリン元素の場合は、リン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スペクトルしかないっぽいので若干微妙だけど、合成の方法は、フラーレンとかを最初に合成したのと同じような方法っぽい(?)

コンピュータ以前の数値計算(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の結果のトレースとは趣旨が外れるので、それも省略