赤木 ： 前回、多体問題をやったわけだけど、重力ソフトニングっていうのをいれたじゃない？

学生C：はい。そうしないと上手くいかないからとかで。あんまりそれでいいのかどうかよくわかってないですが。

赤木 ： まあその、簡単にいっちゃうと、いい時もあるけど駄目な時もあるの ね。例えば、星の数がすごく少ない、3とか10とかだと、なんか全然違いそう でしょ？2つの星が近づいて、お互いの重力で軌道変えるんだけど、 そのたびに曲がりかたが違うからちょっと時間たつと全然変わっちゃうわね。 あと、例えば3体でやってみるとわかるけど、連星ができることがあるでしょ？ 3つがぐるぐる動いているうちに、1つがはねとばされて、あとの2つが結合状 態になるのね。これ、最終的にどうなるかはもちろんソフトニングとかいれる と全然変わるから。

逆にいうと、星とダークマターだけの楕円銀河とかで、星の数もすごく多い時 には、星は銀河全体の重力うけて運動するから、近くの星の影響をそんなにう けなくて、ソフトニングがあっても軌道そんなにすぐにはかわらないわけ。

学生C：そうすると、銀河とかならこれでいい、ということですか?

赤木 ： 楕円銀河でガスがなければそうなんだけど、我々の銀河系みたいな、ガ スがあって円盤がある銀河だとそうでもないのね。というのは、今でも温度が 低いガスが重力で集まって、そこで星ができてるわけ。ガスが集まってるのを 分子雲といって、星ができているところを星形成領域というんだけど、 今の銀河系dも数百個とかの星ができてつつある星形成領域とかあるの。

そうすると、新しくできた星は、それらが自己重力で集まった星団になってる のね。これは星の数少ないから、銀河全体は沢山星があっても星団で生まれた 星がどうなるかを計算しようと思ったらソフトニングとか使えないわけ。

学生C：じゃあどうするんですか?

赤木 ： というのが今日の話ね。こういう、星が集まった星団とか銀河とかを 計算シミュレーションしたい、というのはわりと昔からやられていて、 1960年くらいにもう最初の論文があるのね、で、色々発展があったんだけど、 基本的な考え方は3つなの。

1. 時間刻みを可変にする
2. 時間刻みを星毎ににバラバラに可変にする
3. 二つの星の間の重力ポテンシャル自体を、複数の成分に分けてそれぞれバ ラバラに時間積分する。

赤木 ： 今回やってみて欲しいのは最終的には3つめなんだけど、基本的な考 え方とししては最初の2つもわかっておく必要があるから、簡単にに説明するわね。

赤木 ： まず、「時間刻みを可変にする」っていうのは、どっかで2つの星が 近づいて、今のタイムステップで上手く計算できなくなりそうなら、タイムステップ 短くすればいい、というものね。常微分方程式の数値解法の研究としては、 どういうふうに時間刻みを変えればいいか、というのはすごく大事なことで、 そういうのがはいった実用的なアルゴリズムや数値計算パッケージは一杯ある の。

例えばPython の数値計算用のパッケージにscipy ってのがあって、この中に odeint という常微分方程式の初期値問題用のライブラリがあるのね。 これは ODEPACK っていう、1980年代に Fortran77 で書かれたライブラリを中 で使ってて、これは普通の問題、 non-stiff っていうんだけど、には時間刻 み可変の ABM公式、 stiff な問題には BDF を使うのね、 ABM 公式のうち Aの部分は19世紀に、、、

学生C：あの、すみません、日本語で話を、、、

赤木 ： あー、ごめんなさい。例えば星がN個ある星団で、みんなが勝手に 動いているとして、半径1、全質量1、重力定数1で、速度も1ぐらいだとしましょ う。隣の星までの距離はオーダーとしては だから、普通の星はそれ くらいの時間刻みで積分すればいいわけね。でも、実際やってみたら起こった みたいに、時々2つの星が重力でどんどん近づいて、また離れていく「近接遭 遇」が起こるから、そういうのをちゃんとつかまえるように時間刻み短くすれ ばいいだろうというわけ。

例えばルンゲクッタ公式だと、1ステップ進めるたびに時間刻みかえても大丈 夫だから、そういうことができるように工夫された公式があるの。 つまり、積分誤差を推定できるようにして、その推定された誤差がこちらで 設定した値を超えないようにするわけ。

一番簡単には、まずはある時間刻みで1ステップ積分して、また、元の値から、 時間刻み半分にして2ステップ積分したら、こっちのほうが正確なはずでしょ？ だから、この2つ比べると誤差が大体わかるから、その誤差があまり大きくな らないように、次のステップを調整するわけ。

学生C：なるほど。それでよさそうな気がしますが、駄目なんですか?

赤木 ： これは、星団とか銀河だと駄目なの。理由は2つあって、一つは、星 の数が増えると、平均の距離は でも、「一番近いペア同士」の 距離はもっと小さいのね。これは確率の問題だけど、 で小さくなるはず。だから、沢山の星の中でたまたま2つが近 いだけで、全部の星を無駄に短いステップで積分することになるわけね。

学生C：でも、数値積分ってそういうものじゃないんですか?

赤木 ： そうじゃない、っていう話をこれからするからもうちょっとまって。 もうひとつの問題は、銀河系でも球状星団でもそうなんだけど、密度構造があ るのね。つまり、中心のほうが星の数密度が高いの。銀河だと、銀河中心 にブラックホールがあるけど、それ別にしても太陽から内側くらいだと大体半 径の2乗に反比例して星というか質量密度が上がるのね。

学生C：2乗に反比例だと、、、質量自体だと半径に比例ですね？

赤木 ： そう。そうすると、例えば円軌道で回ってる星の速度は半径でどう変わるかしら？

学生C：えーと、重力は質量に比例で距離の2乗に反比例だから、半径に反比例 ですね。速度を v, 半径 r として遠心力は でこれが に 比例と、速度一定ですか?

赤木 ： そう、我々の銀河系は大体そうなってて、半径どこでも回転速度は 240km/s なの。この値自体は銀河の観測が精密になると段々変わるんだけど、一定、 というのはあんまり変わらないはず。

学生C：えーと、で、なんの話でしたっけ？

赤木 ： つまり、速度が同じで星の密度が高いと、タイムステップ短くしないといけ ないでしょ？なので、星団でも銀河でも、中心に近い星ほど短い刻みが必要だ から、全体を中心で必要な短い時間刻みで積分しないといけないことになるの ね。

我々の銀河系だと、中心に太陽質量の400万倍の質量の超巨大、天文学者の用語 で英語で supermassive っていうんだけど、まあだから超大質量かな、ブラッ クホールがあって、見える星で、15年周期でそのブラックホールの周りを回っ てるのがあるのね。ちょっと余談だけど 2020年4月16日にプレスリリースがあっ て、その星の近点移動が観測できて、一般相対性理論の予言と10%くらいの精 度であったとのことで、移動自体は1周期で12分だったかな、それくらいです ごく大きいんだけど、それが地上の望遠鏡で観測できたのもすごいわよね。

と、それはともかく、15年で1周期の星を積分する時間刻みで太陽みたいな2億 年周期で銀河内運動する星の軌道積分をしたら計算おわらないでしょ？

学生C：と思いますが、じゃあどうするんですか?

学生C：えーと、そんなの本当に上手くいくんですか?大体その予測ってどうするんですか?

赤木 ： 1960年代から90年代くらいまで使われてたのは、3ステップ前くらい までの加速度をとっておいて、それぞの時刻ではその値になるようなの時間 の3次多項式を作って、それで予測する、という方法ね。 これはさっきの君のいうところの日本語じゃないところの、ABM公式の Aのと ころにあたる方法なの。

学生C：補間多項式っていうやつですね。聞いたことあります。

赤木 ： そう。細かいことをいうと、加速度そのものじゃなくて、 divided difference っていう、なんかややこしい量を計算してとってもおくと、補間 計算が簡単になるのでそっち使うの。考え方としてはニュートン補間とかニュー トン・コーツ補間ね。

学生C：はあ。

赤木 ： で、90年代くらいまでは、といったのは、その辺からもうちょっと簡 単な方法をみんな使うようになったから。これでは、重力加速度自体の時間微 分を2階微分まで直接計算して、それでテイラー展開を作ってそれをそのまま 計算する方法。常微分方程式の数値解法としては、そういう、導関数の微分を 使う方法はもちろん昔から知られていたんだけど、重力多体問題に実際に 実用的な形で使ったのは作者氏と、元々の独立時間刻みを発明した Aarseth 先生の共著の論文ね。

学生C：へえ。偉い人なんですね。

赤木 ： うーん、どうかしら？あと、2回導関数までだと4次精度なんだけど、 6次と8次は似鳥さんっていう人が論文にしてて、あとそれ以上も論文にはして ないけどセミナーで発表した資料はあって、それ引用して論文かいてた人がい たわ。

学生C：はあ。

赤木 ： まあ理論的には、さっきの話で一番短いステップが くらいで他が平均的に なら、 一番短いステップで進めていった時に ステップで N 粒子だから、平均的には 粒子よね。結構沢山。

でも、中心密度が高いとかになるとあんまり上手くいかないのね。

学生C：そうですよね。もっとうまい方法はないんですか？

赤木 ： その辺をこれからね。ただ、今のところ世界的にも、星団でも銀河で も、この独立時間刻みが基本的に使われている方法なの。いわゆる state of the art というのね。だから、これから話をするのは、もちろん論文とかあっ て 上手くいくのもわかってるけど、世界で広く使われている、というわけではま だないのね。

学生C：それ、大丈夫なんですか？

赤木 ： まあ、その辺ちょっと調べるのもやってみようと思ってるわけ。 考え方は単純なのよ。といってもいきなり重力多体問題だと説明もちょっとや やこしいので、まず2体問題を例にして考えることにするわね。

話をリープフロッグ公式に戻すと、

1. 速度を今の位置の加速度で時間刻みの半分だけ進める
2. 位置を上の速度で時間刻み分進める。
3. 速度を今の位置の加速度で時間刻みの半分だけ進める

学生C：そうですが、だからどうと？

赤木 ： ここで、微分方程式が、

学生C：えーと、でも、2体問題では重力相互作用って1つの関数じゃないです か？この話どう関係するんですか？

赤木 ： そうなんだけど、もうちょっと待って。とりあえず、2つだとすると、 こういう手順もありえるわけ。

1. 速度を今の位置の f で時間刻みの半分だけ進める
2. gだけ使って次の時刻まで数値積分する
3. 速度を今の位置の加速度で時間刻みの半分だけ進める

学生C：それはなんとなく上手くいきそうですが。惑星間相互作用は太陽から の重力に比べてずっと小さいわけですよね？

赤木 ：そう。この方法は MIT の Jack Wisdom と Matthew Holman が 1991 年 の Symplectic maps for the N-body problem っていう論文で提案して、すごく広く使われるようになったの。

でね、それとちょっと似たことを、でも太陽重力と惑星間みたいな違うものの 和ではなくて、ケプラー問題みたいな、1つなんだけど、例えば離心率がすご く大きくて、太陽に近い時と遠い時で時間ステップ変えたい、みたいな時に 適用したのが Skeel と Biesiadecki の 1994年の、 Symplectic Integration with Variable Stepsize っていう論文なのね。どういう 考え方かというと、太陽からの重力を無理矢理2つとかそれ以上の成分に分けるのね。

つまり、重力加速度は

学生C：なんか無駄に面倒じゃないですか?独立時間刻みとか普通の可変時間刻みでは駄目なんですか?

赤木 ：まあだから、少なくとも重力多体系ではこれすごく得なの。なぜかと いうと、どこかで2つの粒子がすごく近づいいたとして、その2つの軌道とその 2つの間の重力だけ細かく計算されるのね、他の粒子は全く影響うけなくてす むわけ。だから、独立時間刻みだと、動かす粒子が2個しかなくても他の全部 の粒子からの力計算してたけど、こっちだと2個の間だけでいいわけ。

もうちょっというと、あちこちで同時とかちょっと違う時間に近接遭遇がおき ても、それぞれちょっとづつ計算時間増えるだけだから、あんまり全体に影響 しないのね。

これは、並列計算機使うとか、あと重力の計算法工夫して減らすとかの時にす ごくメリット大きいのね。ちょっと減らすほうの話をすると、前に作ってもらっ た重力多体問題のプログラムだと、運動方程式の通りに1つの粒子は他の 全粒子からの重力を受けるとしたじゃない？

学生C：はい。だって、実際にそうですよね？

赤木 ：そうなんだけど、もうちょっと上手く計算できないか、という話ね。 伝統的にというか、1970年代くらいに考えられた方法は2つあって、一つは 高速フーリエ変換を使う方法ね。

物理数学で偏微分方程式はやった？

学生C：ええ、熱伝導とか波動方程式とか。

赤木 ：ポアソン方程式ってやった？

学生C：あんまり記憶がないんですが、、、

赤木 ：ラプラス方程式は？

学生C：そっちは記憶があるような気がします。

赤木 ：あらら、電磁気でマクスウェルの方程式は？

学生C：えーと、やったはずですね。

赤木 ：それに、電場の方程式ってのがあるでしょ？

学生C：そういわれるとそんな気がしてきました。

赤木 ：これで、 は電荷密度で、 は静電ポテンシャルで、 は電場で、そこに電荷 があると の力を受ける わけね、これは、電荷が2つあると距離の二乗に反比例する力をお互いに及ぼ す、というのと同じことになってるのはわかるわね？

学生C： あんまりわかってないですが、、、えーと、そうなるんですか?

赤木 ：そうなるの。で、重力相互作用と静電相互作用って、質量と電荷とで 力の種類は違うけど式の形は同じでしょ？だから重力場と重力ポテンシャルも、 静電場の方程式と同じように質量密度に対するポアソン方程式を解けば求めら れるのね。なので、ポアソン方程式を解くのにフーリエ変換を使う、というの がこの方法の考え方で、具体的には、質量密度を、空間を格子に切って、各点 での値での値と思うことにして、それを数値的にフーリエ級数展開というかフー リエ変換して、フーリエ級数のほうで係数掛けるとポアソン方程式解けるで しょ？それで逆フーリエ変換してポテンシャルにするわけ。力も逆フーリエ変 換でもいいし、精度ちょっと落ちるけどポテンシャルを数値微分でもまあ求め るわね。

高速フーリエ変換は、格子点の数を として計算量が ですむから、まともに粒子間相互作用を計算すると になるのに比べてずっと速く計算できるの。

学生C：えーと、でも、フーリエ変換だと周期境界になってなんか変じゃないですか?

赤木 ：そこはなんとかする方法があるの。私よく知らないけど。

学生C：そうなんですか。あと、これ、2つの粒子が近づいたとかあったらどう するんですか?

赤木 ：基本的にはどうにもならなくて、初めてから粒子の質量は空間的に広 がっている、と思うのが1つの方法で、上にでてきた みたいなのを使っ て相互作用2つにわけて、近くは直接計算、遠くだけフーリエ変換で、という のがもうひとつの方法ね。これは particle-particle-particle-tree ()法っていわれていて 1980年代に提案されたの。

学生C： よさげな方法にみえますが、どうなんでしょう？

赤木 ： いいんだけど、天文のシミュレーションだと、中心の密度がどんどん 上がるとか、銀河だとずーっと外側まで広がってるとかあって、あんまり上手 くいかないの。

学生C：あ、あともうひとつあるっていってませんでしたっけ？

赤木 ：あるけど、これも同じような問題があるので簡単にね。物理数学で球面調和関数 展開ってやった？量子力学とかでは？

学生C：やったのかもしれないですがえーと、、、

赤木 ：これもやっぱりラプラス方程式とかポアソン方程式と関係するんだけ ど、ラプラス方程式って境界条件与えると解決まるでしょ？2次元で、円板領 域だと、境界は円で、円周上での周期関数が境界条件で、その固有関数を求め ると みたいなのがでてくるでしょ？

学生C：そうだったような気がします。

赤木 ：で、円板の外側領域だと で書けて、こ れを外部解っていうんだけど、これ、円板の中の任意の粒子分布について、そ れが作る電場とか重力場の角度方向のフーリエ変換をすれば、無限遠方までの 解が求められることになるでしょ？

学生C：なんか騙されてる気がしますが、そうなんですか?

赤木 ：数学としてはそうなるでしょ？でね、2次元で円板だとフーリエ変換で いいけど、3次元だと球面になるじゃない？なので、球面上にフーリエ変換を 拡張したものが必要で、それが球面調和関数なのね。具体的な形は、 3次元極座標を で書いた時に、 の三角関数と のルジャンドル陪関数っていう直交多項 式系の積で書けるんだけど、ここはまあ具体的な形はよくて要するにそういう、 円周上のフーリエ級数展開を球面上に拡張したものと思って。 これを多重極展開っていうの。

そうすると、ある半径での重力ポテンシャルは、その内側の粒子が作る重力場 の外部解と、その外側の粒子が作るポテンシャルの内部解の合計になるから、 それを順番に、あるいはなんか上手く並列化して計算すればいいわけ。 これも、じゃあ展開項数どこまでとるの問題というのがあって、楕円銀河みた いな丸いのならこれでいいけど、、、みたいな話で最近あんまり使われてない かも。

赤木 ：まあ、どうなんだろう？ということで、色々研究されてるわけ。 天文での応用としては、2つの粒子が近づくというのはそこそこ稀だから、 そこは短い刻みとか高精度の公式とかでちゃんとやって、特に2つの粒子がす ごく近づくとかあるから可変時間刻みの高精度公式を使いたいのね。一方、 それ以外の遠くの粒子からの重力は、一定刻みで、まあ簡単な話としてはリー プフロッグですませたいわけ。そうすると、ひとつ問題なのは、 近接遭遇に対して使う公式と の関係なのね。

学生C：えーと、どういうことでしょうか?

赤木 ：高精度の公式って、軌道や加速度が少なくともその精度の次数くらい の回数微分可能じゃないと上手くいかない「かもしれない」のね。

学生C：何回も微分可能ってどういうことでしたっけ？

赤木 ：ちょっと確認だけど、微分可能ってどういうこと？

学生C：1変数でいいですよね？各点で微分可能であればいいんじゃないですか?

赤木 ：それだと上手くいかない変な関数沢山つくれるから、普通微分可能っ ていう時には、微分した関数が連続関数になることをさすのね。連続的微分可 能とも。

連続関数の定義はいいわよね？

学生C：こっちは各点でいいですよね？その点での値に、その点でに近傍の値 が収束すること、つまり、任意の に対して、 で になることです。

赤木 ：はい、よくできました。なので、n回連続的微分可能っていうのは、 1回からn回導関数までが全て連続関数ということね。級って書くこと が数学では多いかも。

なので、数値積分法がn次精度だったとしても、積分する関数のほうが じゃなくてもっと微分できる回数少ないと上手くいかないかもしれな いじゃない？

学生C：それなら、無限回微分可能ならいいんじゃないですか？なんか とかいうのがなかったでしたっけ？

赤木 ：あ、それは、「良い質問」というやつね。どういう関数が として欲しいかというと に対して

• で0
• で1
• その間を連続かつ微分可能につなぐ

逆に、誤差関数、って知ってるわね？正規分布を積分して累積分布関数にした やつ、これは無限回微分可能で左側で0で右側で1になるからいい感じなんだけ ど、数学的に厳密には無限遠まで0じゃないから、

学生C：それ、多項式にするしかないっていってます？

赤木 ：概ねそうなんだけど、一応この辺数学的に厳密には、そうじゃなくっ て上の性質をもって無限回連続的微分可能な関数って存在するのね。 数学の話としては、

• で正、その区間での積分が有限値
• その両側で0
• 実数全体で無限回連続的微分可能

学生C：じゃあそれ使えるんじゃないですか？

赤木 ：と一瞬私も思ったんだけど、よく解説読むと隆起関数は必ず解析的じゃないって書い てあってね、、、

学生C：解析的ってなんでしたっけ？

赤木 ：テイラー展開が、少なくとも近傍で、収束することが、ある関数があ る点で解析的であることの定義。隆起関数は 1 と -1 以外のどこでも解析的 なんだけど、この2点では解析的じゃないことが数学的に証明できるんだって。

学生C：はあ。

赤木 ：例えば、

赤木 ：そうねえ。まあ、世界には色々な人がいるのよ。なんか話がそれたけ ど、なんだっけ？えーと、だから、無限回は無理として有限回ならどうかとい うことね。隆起関数と同じように、 で正で微分可能で両端で n回導関数まで0、みたいなのを考えればいいから、多分一番簡単な例は

学生C：えーと、これでいいってのはわかる気がしますが、他にもっといいの とかあったりしないんですか？

赤木 ：それはまあ、良いとは何かということね。まず、「条件を満たす、最 低次の多項式」が、まあ、計算が楽という観点では一番良いから、そういうも のを、としてみると、両端で n回導関数まで0、関数の値は 0, 1 というこ とから条件は 2n+2 個あるから、一般には係数が 2n+1 個ある、2n+1 次多項式が いるでしょ？で、上のは 2n次多項式を1回積分して2n+1 次多項式になってる から、最低次の条件を満たしてて、次数の意味ではよい、ということにはなる わね。

学生C：次数がもっと高いののほうがなんかいいことがあるとかそういうのは どうなんでしょう？

赤木 ：その辺あんまりわかってないから、調べてみましょうみたいなね。 というわけで、今日の話としては、

以下学生C独白。

そんなこといっても、今なら Mathematica とか、Web の Wolfram Alpha でちょ いちょいと。Alpha にいれてテキストでだすと、

  \int x^4 (1 - x)^4 dx

  integral x^4 (1 - x)^4 dx = x^9/9 - x^8/2 + (6 x^7)/7 - (2 x^6)/3 + x^5/5 + 定数

(ちょっと捜すが Ruby のは一杯見つかったが Crystal のは発見できなかったっぽく)

多項式って、要するに x のべき乗の係数の配列だから、MathVector みたいに Array から適当につくればいいか、必要なのは、、、掛け算と積分だけかな。 一応加減算と微分もいれておくと、、、

    1:#
    2:# polynomial handling library 
    3:# Copyright 2020- Jun Makino
    4:#
    5:class Polynomial(T) < Array(T)
    6:  def Polynomial.zero
    7:    [T.zero].to_poly
    8:  end
    9:  def extended(i)
   10:    (i < self.size) ? self[i] : T.zero
   11:  end
   12:  def +(a)
   13:    newsize = {self.size, a.size}.max
   14:    newp=(0..(newsize-1)).map{|k|  self.extended(k)+ a.extended(k)}
   15:    while newp[(newp.size) -1] == T.zero
   16:      newp.pop
   17:    end
   18:    newp.to_poly
   19:  end
   20:  def -() self.map{|x| -x}.to_poly  end
   21:  def -(a) self + (-a) end
   22:  def +() self end
   23:  def *(a)
   24:    newp = Array.new(self.size+a.size-1){T.zero}
   25:    self.size.times{|i|
   26:      a.size.times{|j|
   27:        newp[i+j] += self[i]* a[j]
   28:      }
   29:    }
   30:    newp.to_poly
   31:  end
   32:  def ^(n : Int)
   33:    newp=self
   34:    (n-1).times{newp *= self}
   35:    newp
   36:  end
   37:  def differentiate
   38:    newp= self.map_with_index{|x,k| x*k}
   39:    newp.shift
   40:    newp.to_poly
   41:  end
   42:  def integrate
   43:    ([T.zero] +  self.map_with_index{|x,k| x/(k+1)}).to_poly
   44:  end
   45:  def evaluate(x)
   46:    self.reverse.reduce{|p,a| p = p*x+a}
   47:  end
   48:end
   49:class Array(T)
   50:  def to_poly
   51:    x=Polynomial(T).new
   52:    self.each{|v| x.push v}
   53:    x
   54:  end
   55:end

  def *(a)
    newp = Array.new(self.size+a.size-1){T.zero}
    self.size.times{|i|
      a.size.times{|j|
        newp[i+j] += self[i]* a[j]
      }
    }
    newp.to_poly
  end

積分は多項式作るから不定積分で、定数は0にするとして、 今の多項式の各係数は添字+1次の項になるからそれで割ると。

  def integrate
    ([T.zero] +  self.map_with_index{|x,k| x/(k+1)}).to_poly
  end

  def ^(n : Int32)
    newp=self
    (n-1).times{newp *= self}
    newp
  end

  def evaluate(x)
    self.reverse.reduce{|p,a| p = p*x+a}
  end

    1:#
    2:# test_polynomial.cr
    3:#
    4:require "./polynomial.cr"
    5:require "grlib"
    6:include GR
    7:
    8:poly=[1,1].to_poly
    9:pp! poly
   10:pp! poly+poly
   11:pp! poly - [1,1,1].to_poly
   12:pp! poly*poly*poly
   13:pp! ([1,1].to_poly)^10
   14:pp! ([0.1,0.1].to_poly)^10
   15:pp! ([1,1].to_poly)^5
   16:pp! (([1,1].to_poly)^5).differentiate
   17:pp! (([1.0,1.0].to_poly)^5).integrate
   18:base= [1.0,-1.0].to_poly*[0.0,1.0].to_poly
   19:pp! base
   20:pp! (base^3).integrate
   21:pp! [1.0,2.0,1.0].to_poly.evaluate(1.0)
   22:pp! [1.0,2.0,1.0].to_poly.evaluate(2.0)
   23:pp! [1.0,2.0,1.0].to_poly.evaluate(0.5)
   24:a = [1.0,2.0,1.0].to_poly
   25:b = [a,a].to_poly
   26:pp! a
   27:pp! b
   28:
   29:def create_switch_function(order, rin : Float64, rout : Float64)
   30:  poly= (([1.0,-1.0].to_poly*[0.0,1.0].to_poly)^order).integrate
   31:  normalization = 1.0/poly.evaluate(1.0)
   32:  scale = 1.0/(rout-rin)
   33:  -> (x : Float64){
   34:    if x > rout
   35:      1.0
   36:    elsif x < rin
   37:      0.0
   38:    else
   39:      poly.evaluate((x-rin)*scale)*normalization
   40:    end
   41:  }
   42:end
   43:
   44:f5 = create_switch_function(5,0.5,1)
   45:11.times{|i| pp! f5.call(i*0.1)}
   46:
   47:setwindow(-0.2,1.2,-0.2,1.2)
   48:box
   49:setcharheight(0.05)
   50:mathtex(0.5, 0.06, "x")
   51:mathtex(0.03, 0.5, "K(x)")
   52:x = Array.new(100){|i| 1.4*i/100 -0.2}
   53:polyline(x, x.map{|v| f5.call(v)})

 def create_switch_function(order, rin : Float64, rout : Float64)
   poly= (([1.0,-1.0].to_poly*[0.0,1.0].to_poly)^order).integrate
   normalization = 1.0/poly.evaluate(1.0)
   scale = 1.0/(rout-rin)
   -> (x : Float64){
     if x > rout
       1.0
     elsif x < rin
       0.0
     else
       poly.evaluate((x-rin)*scale)*normalization
     end
   }
 end

実行してみると、

 gravity> crystal test_polynomial.cr
 poly # => [1, 1]
 poly + poly # => [2, 2]
 poly - [1, 1, 1].to_poly # => [0, 0, -1]
 (poly * poly) * poly # => [1, 3, 3, 1]
 ([1, 1].to_poly) ^ 10 # => [1, 10, 45, 120, 210, 252, 210, 120, 45, 10, 1]
 ([0.1, 0.1].to_poly) ^ 10 # => [1.0000000000000006e-10,
  1.0000000000000009e-9,
  4.500000000000003e-9,
  1.200000000000001e-8,
  2.1000000000000016e-8,
  2.520000000000002e-8,
  2.1000000000000016e-8,
  1.200000000000001e-8,
  4.500000000000003e-9,
  1.0000000000000009e-9,
  1.0000000000000006e-10]
 ([1, 1].to_poly) ^ 5 # => [1, 5, 10, 10, 5, 1]
 (([1, 1].to_poly) ^ 5).differentiate # => [5, 20, 30, 20, 5]
 (([1.0, 1.0].to_poly) ^ 5).integrate # => [0.0, 1.0, 2.5, 3.3333333333333335, 2.5, 1.0, 0.16666666666666666]
 base # => [0.0, 1.0, -1.0]
 (base ^ 3).integrate # => [0.0, 0.0, 0.0, 0.0, 0.25, -0.6, 0.5, -0.14285714285714285]
 [1.0, 2.0, 1.0].to_poly.evaluate(1.0) # => 4.0
 [1.0, 2.0, 1.0].to_poly.evaluate(2.0) # => 9.0
 [1.0, 2.0, 1.0].to_poly.evaluate(0.5) # => 2.25
 a # => [1.0, 2.0, 1.0]
 b # => [[1.0, 2.0, 1.0], [1.0, 2.0, 1.0]]
 f5.call(i * 0.1) # => 0.0
 f5.call(i * 0.1) # => 0.0
 f5.call(i * 0.1) # => 0.0
 f5.call(i * 0.1) # => 0.0
 f5.call(i * 0.1) # => 0.0
 f5.call(i * 0.1) # => 0.0
 f5.call(i * 0.1) # => 0.011654205440003432
 f5.call(i * 0.1) # => 0.24650186752007183
 f5.call(i * 0.1) # => 0.7534981324802157
 f5.call(i * 0.1) # => 0.9883457945601857
 f5.call(i * 0.1) # => 1.0

11.7. ハミルトニアン分割によるケプラー問題数値解

じゃあこれいれてケプラー問題か。どうするんだこれ？単純に考えると、リー プフロッグは

1. 加速度計算する
2. 速度を半分更新する
3. 位置を更新する
4. 加速度計算する
5. 速度を半分更新する

1. 加速度計算する
2. 速度を半分更新する
3. ルンゲクッタで軌道積分
4. 加速度計算する
5. 速度を半分更新する

  def leapfrog(x,v,h,f)
    f0 = f.call(x)
    x+= v*h + f0*(h*h/2)
    f1 = f.call(x)
    v+= (f0+f1)*(h/2)
    {x,v}
  end

  def hybrid(x,v,h,f,g,q)
    f0 = f.call(x)
    variable_step_rk4(x,v,h,g,q)
    f1 = f.call(x)
    v+= (f0+f1)*(h/2)
    {x,v}
  end

  def hybrid(x,v,h,f,g,q)
    f0 = f.call(x)
    v+= f0*(h/2)
    variable_step_rk4(x,v,h,g,q)
    f1 = f.call(x)
    v+= f1*(h/2)
    {x,v}
  end

(数日後)

学生C：なんかできたみたいです。

赤木 ：え、本当？

学生C：本当って、、とりあえずエネルギーのグラフだけ。

  splitintegrator -g -w 0.0001 -e 0.99999 -q 0.005 -s 0.01 -E -o 10 -p 4

学生C：

     1:#
     2:# splitintegrator.cr
     3:# Copyright 2020 J. Makino
     4:#
     5:
     6:require "grlib"
     7:require "./integratorlib.cr"
     8:require "./vector3.cr"
     9:require "./mathvector.cr"
    10:require "./polynomial.cr"
    11:require "clop"
    12:include Math
    13:include GR
    14:
    15:optionstr= <<-END
    16:  Description: Test integrator for Kepker problem with split integrator
    17:  Long description:
    18:    Test integrator for Kepker problem with split integrator
    19:    (c) 2020, Jun Makino
    20:
    21:  Short name:           -s
    22:  Long name:  --soft-step
    23:  Value type:float
    24:  Default value: 0.1
    25:  Variable name: h
    26:  Description:timestep for leapfrog part
    27:  Long description:     timestep for leapfrog part
    28:
    29:  Short name: -o
    30:  Long name:  --norbits
    31:  Value type:int
    32:  Default value:1
    33:  Variable name:norb
    34:  Description:Number of orbits
    35:  Long description:     Number of orbits
    36:
    37:  Short name:-w
    38:  Long name:  --window-size
    39:  Value type:  float
    40:  Variable name:wsize
    41:  Default value:1
    42:  Description:Window size for plotting
    43:  Long description:
    44:    Window size for plotting orbit. Window is [-wsize, wsize] for both of
    45:    x and y coordinates
    46:
    47:  Short name:-e
    48:  Long name:--ecc
    49:  Value type:float
    50:  Default value:0.0
    51:  Variable name:ecc
    52:  Description:Initial eccentricity of the orbit
    53:  Long description:     Initial eccentricity of the orbit
    54:
    55:  Short name:-g
    56:  Long name:--graphic-output
    57:  Value type:        bool
    58:  Variable name:gout
    59:  Description:
    60:    whether or not create graphic output (default:no)
    61:  Long description:
    62:    whether or not create graphic output (default:no)
    63:
    64:
    65:  Short name:-E
    66:  Long name:--plot-energy-error
    67:  Value type:        bool
    68:  Variable name:eplot
    69:  Description:
    70:    plot de instead of orbit
    71:  Long description:
    72:    plot de instead of orbit
    73:
    74:  Short name:-r
    75:  Long name:--outer-cutoff-radius
    76:  Value type:        float
    77:  Variable name:rout
    78:  Default value:0.5
    79:  Description:          outer cutoff radius for splitting
    80:  Long description:     outer cutoff radius for splitting
    81:
    82:
    83:  Short name:-p
    84:  Long name:--cutoff-order
    85:  Value type:        int
    86:  Variable name:cutoff_order
    87:  Default value:4
    88:  Description:          Smoothness order for the cutoff function
    89:  Long description:     Smoothness order for the cutoff function
    90:
    91:  Short name:-q
    92:  Long name:--accuracy-parameter
    93:  Value type:        float
    94:  Variable name:eta
    95:  Default value:0.1
    96:  Description:          Parameter for timestep criterion for hard part
    97:  Long description:     Parameter for timestep criterion for hard part
    98:END
    99:
   100:clop_init(__LINE__, __FILE__, __DIR__, "optionstr")
   101:options=CLOP.new(optionstr,ARGV)
   102:
   103:def variable_step_rk4(x,v,h,g,q)
   104:  t=0.0
   105:  ff =  ->(x : MathVector(Vector3), t : Float64){ [x[1], g.call(x[0])].to_mathv} 
   106:  while t < h
   107:    dt = sqrt((x*x)/(v*v))*q
   108:    dt = h-t if t+dt > h
   109:    xv = [x,v].to_mathv
   110:    xv,t = Integrators.rk4(xv,t,dt,ff)    
   111:    x,v= xv
   112:  end
   113:  {x,v}
   114:end
   115:
   116:def hybrid(x,v,h,f,g,q)
   117:  f0 = f.call(x)
   118:  v+= f0*(h/2)
   119:  x,v = variable_step_rk4(x,v,h,g,q)
   120:  f1 = f.call(x)
   121:  v+= f1*(h/2)
   122:  {x,v}
   123:end
   124:
   125:def kepler_acceleration(x,m)
   126:  r2 = x*x
   127:  r=sqrt(r2)
   128:  mr3inv = m/(r*r2)
   129:  -x*mr3inv
   130:end
   131:def energy(x,v,m)
   132:  m*(-1.0/sqrt(x*x)+v*v/2)
   133:end
   134:
   135:def create_switch_function(order, rin : Float64, rout : Float64)
   136:  poly= (([1.0,-1.0].to_poly*[0.0,1.0].to_poly)^order).integrate
   137:  normalization = 1.0/poly.evaluate(1.0)
   138:  scale = 1.0/(rout-rin)
   139:  -> (x : Float64){
   140:    if x > rout
   141:      1.0
   142:    elsif x < rin
   143:      0.0
   144:    else
   145:      poly.evaluate((x-rin)*scale)*normalization
   146:    end
   147:  }
   148:end
   149:
   150:m=1.0
   151:switch= create_switch_function(options.cutoff_order,
   152:                               options.rout*0.5,
   153:                               options.rout)
   154:fout =  -> (xx : Vector3){kepler_acceleration(xx,m)*switch.call(sqrt(xx*xx))}
   155:fin=  -> (xx : Vector3){kepler_acceleration(xx,m)*(1-switch.call(sqrt(xx*xx)))}
   156:t=0.0
   157:x= Vector3.new(1.0+options.ecc,0.0,0.0)
   158:v= Vector3.new(0.0,sqrt((1-options.ecc)/(1+options.ecc)),0.0)
   159:if options.gout
   160:  wsize=options.wsize
   161:  if options.eplot
   162:    setwindow(0.0, options.norb*PI*2,-wsize, wsize)
   163:  else
   164:    setwindow(-wsize, wsize,-wsize, wsize)
   165:  end
   166:  box
   167:  setcharheight(0.05)
   168:  if options.eplot
   169:    mathtex(0.5, 0.06, "t")
   170:  mathtex(0.06, 0.5, "|de|")
   171:  else
   172:    mathtex(0.5, 0.06, "x")
   173:    mathtex(0.06, 0.5, "y")
   174:  end
   175:end
   176:e0 = energy(x,v,m)
   177:emax = 0.0
   178:de=0.0
   179:while t < options.norb*PI*2 - options.h
   180:  dep=de
   181:  xp=x
   182:  tp=t
   183:  x, v = hybrid(x, v, options.h, fout, fin, options.eta)
   184:  t+=options.h
   185:  de =energy(x,v,m)-e0
   186:  emax = {de.abs, emax}.max
   187:  if options.gout
   188:    if options.eplot
   189:      polyline([tp, t], [dep, de])
   190:    else
   191:      polyline([xp[0], x[0]], [xp[1], x[1]])
   192:    end
   193:  end
   194:  pp! x, v, energy(x,v,m)
   195:end
   196:p! -emax/e0
   197:c=gets if options.gout
   198:

  switch= create_switch_function(options.cutoff_order,
                                 options.rout*0.5,
                                 options.rout)
  fout =  -> (xx : Vector3){kepler_acceleration(xx,m)*switch.call(sqrt(xx*xx))}
  fin=  -> (xx : Vector3){kepler_acceleration(xx,m)*(1-switch.call(sqrt(xx*xx)))}

fout, fin は外側、内側の関数で、もともとの相互作用がケプラーなので距離 の2乗に反比例の力、これは前に使ったのと同じですね。それにスイッチ関数 と、1からスイッチ関数をひいたものをそれぞれ掛けてます。

そのあとは、エネルギー誤差を書けるようにしたのの他は前のケプラー問題の 数値積分とほとんど同じなのですが、数値積分するとこの本体が、前のは

  x, v = integrator.call(x,v,h)

  x, v = hybrid(x, v, options.h, fout, fin, options.eta)

  def hybrid(x,v,h,f,g,q)
    f0 = f.call(x)
    v+= f0*(h/2)
    x,v = variable_step_rk4(x,v,h,g,q)
    f1 = f.call(x)
    v+= f1*(h/2)
    {x,v}
  end

  def leapfrog(x,v,h,f)
    f0 = f.call(x)
    x+= v*h + f0*(h*h/2)
    f1 = f.call(x)
    v+= (f0+f1)*(h/2)
    {x,v}
  end

  def leapfrog(x,v,h,f)
    f0 = f.call(x)
    v+= f0*(h/2)
    x+= v*h
    f1 = f.call(x)
    v+= f1*(h/2)
    {x,v}
  end

    x+= v*h

    x,v = variable_step_rk4(x,v,h,g,q)

赤木 ：なるほど。もっともらしいわね。

学生C：variable_step_rk4 は gのほうを時間刻み可変のルンゲクッタで時間 h まで積分するものです。これはこんな感じです。

  def variable_step_rk4(x,v,h,g,q)
    t=0.0
    ff =  ->(x : MathVector(Vector3), t : Float64){ [x[1], g.call(x[0])].to_mathv} 
    while t < h
      dt = sqrt((x*x)/(v*v))*q
      dt = h-t if t+dt > h
      xv = [x,v].to_mathv
      xv,t = Integrators.rk4(xv,t,dt,ff)    
      x,v= xv
    end
    {x,v}
  end

    ff =  ->(x : MathVector(Vector3), t : Float64){ [x[1], g.call(x[0])].to_mathv} 

赤木 ：まあいいんじゃないかしら？速度気にするとまた違う書き方がいるかもだけど。

学生C： はい。で、 while の中で h だけ積分します。時間刻みを

      dt = sqrt((x*x)/(v*v))*q

      dt = h-t if t+dt > h

赤木 ：これ、このあとで t+dt が永久に h にならないとかは起こらない？

学生C： 引き算で、 h のほうが大きいと大丈夫じゃないですか？

赤木 ：そうね、IEEE-754 の規定からは大丈夫なはずなんだけど、何故そうい えるかは調べといてね。

学生C：読者がですね？で、あとは、 x,v を MathVector にして rk4 に渡し てまた Vector3 2つに戻すだけです。

11.9. まとめ

1. ハミルトニアン分割によって、重力多体問題を高速かつ高精度に計算する 方法を、まずケプラー問題で振舞いを調べた。
2. この目的のため、多項式の乗算・不定積分を行なうライブラリを作成した。

11.10. 参考資料

Skeel, R. D. and Biesiadecki, J. J., "Symplectic Integration with Variable Stepsize", Annals of Numer. Math. 1994, 1, 191-198. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.28.9532&rep=rep1&type=pdf

Hernandez, David M. Improving the accuracy of simulated chaotic N-body orbits using smoothness, Monthly Notices of the Royal Astronomical Society, Volume 490, Issue 3, p.4175-4182 https://academic.oup.com/mnras/article-abstract/490/3/4175/5575201?redirectedFrom=fulltext

Previous

ToC

Next