Previous ToC Next

11. ハミルトニアン分割

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

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

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

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

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

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

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

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

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

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

学生C:えーと、、、

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

11.1. 従来の考え方

赤木 : まず、「時間刻みを可変にする」っていうのは、どっかで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:と思いますが、じゃあどうするんですか?

赤木 :1960年代以降の基本的な考え方は、「独立時間刻み」というものね。 要するに、それぞれの星に、「自分の時間」と「自分の時間刻み」をもたせる の。図 20 みたいな感じ。

Figure 20: 独立時間刻みの考え方

これ、それぞれの星が、今の自分の時刻と次の時刻というのをもってて、 次の時刻が一番小さい星を積分するのね、積分するには、少なくとも新しい時 刻で他の星からの重力を計算する必要があるわけで、それは実際に他の星をそ こまで積分はできないから場所を「予測」するの。で、他の星の予測した位置 から自分の予測した位置への重力を計算して、それ使って自分の位置・速度は 計算しなおして、で、新しい「次の時刻」を決めるわけ。そうしたら、 1個積分すると次の時刻が一番小さい星が変わるから、それを次々に積分して いくのね。

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

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

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

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

学生C:はあ。

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

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

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

学生C:はあ。

赤木 : あともうちょっと独立時間刻みの話ね。最近のスパコンってみんな並 列計算機だから、こういう、1個星選んでそれを積分、では効率でないの。 これは1980年代くらいから問題になりはじめて、なのでブロックステップとい うのが発明されたのね。これは 21 にある通りで、時間刻みを 2のべき乗に合わせておくと、 時間や時間刻みの値が同じになる星がでてくるから、それらは並列に計算でき る、というわけ。

Figure 21: ブロックステップの考え方

学生C:そりゃ、元々よりましそうですが、どれくらい上手くいくんですか?

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

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

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

11.2. 新しい考え方

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

学生C:それ、大丈夫なんですか?

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

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

という形で書いたけど、もともと2つの公式を組合せてた、ということを思い 出すと、

とも書けるわけね。これって、

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

となるじゃない?

学生C:そうですが、だからどうと?

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

と、加速度が2つの関数の和だとするわね。で、 f は計算が大変だけどあんま り速くは変わらない、 g はそれほど大変じゃないけど f より急に変わる、あ るいはずっと大きくて正確に計算する必要がある、という場合を考ると、

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

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

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

例えば、太陽系の惑星とかだとこれ上手くいくのね。 f を惑星間の重力相互 作用、 g を太陽重力とすると、惑星は沢山あるから f の計算は大変で、gは 数値計算といっても解析解もあるからそっちつかうこともできるし、高精度の 公式で数値積分してもいいわけ。

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

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

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

つまり、重力加速度は

として、適当な関数 をもってきて、

とするわけ。で、前に書いた、リープフロッグで f を積分して、但し途中で は g を積分するのね。 が、 $x$ が小さい、つまり太陽に近いところで は0になるようにすれば、太陽に近いところでは g の、なんか正確な方法での 積分、 遠いところでは逆に f だけの、普通のリープフロッグでの積分になって、リー プフロッグでは太陽に近くなると上手くいかない、という問題が回避されるわ け。ここでは g は0だから、計算しなくていいのね。

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

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

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

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

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

11.3. 遠距離力の高速計算1 FFT と球面調和関数展開法

赤木 :そうなんだけど、もうちょっと上手く計算できないか、という話ね。 伝統的にというか、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次元極座標を で書いた時に、 の三角関数と のルジャンドル陪関数っていう直交多項 式系の積で書けるんだけど、ここはまあ具体的な形はよくて要するにそういう、 円周上のフーリエ級数展開を球面上に拡張したものと思って。 これを多重極展開っていうの。

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

11.4. 遠距離力の高速計算2 ツリー法

学生C:もうちょっと新しい方法とかはないんですか?

赤木 :もちろんあって、ちょっとその話ね。今割合よく使われる方法の1つは ツリー法というものね。

Figure 22: ツリー法の考え方。2次元の場合。

3次元面倒なので2次元で説明するわね。図 22 みたいに、なにか粒 子分布があったとして、それ全体がはいる正方形を、再帰的に4個に分けるの を、それぞれ正方形に粒子が1または0個になるま繰り返してできる木構造を 4分木、3次元で立方体だと8分木というのね。これ使うと、例えば右下の粒子 から、左上の、全体の 1/4 の面積の正方形をみると、割合離れてるから、 そうすると多重極展開で重力評価してもあんまり誤差大きくないかもしれないわ け。つまり、この木構造の各節点に対応する粒子分布について全部多重極展開 を求めておけば、それらをを上手く使って各粒子への力を高速に求められる、 というのがツリー法ね。

上手く使う、というのが実は簡単で、ツリーを上から辿っていって、精度的に 大丈夫ならそこで多重極展開を評価する、というだけね。まあ、並列計算機と か高速化とかいうと無限に色々なことがあるんだけど、それはちょっとまた別の話。

1980年代に、よく似た考え方で少なくとも2つ論文がでてるんだけど、8分木を 使う、というのは Barnes と Hut のアイディアで、こちらがプログラム書く のが割合楽でツリー作るのも高速ということがあって圧倒的に一番使われてる わね。

ツリー法のいいところは、高速フーリエ変換使う方法みたいに格子を切るわけ ではなくて、粒子分布に合わせて木構造作るので、密度が高いところと低いと ころがあっても大丈夫で、計算時間もあんまりのびないということね。

というわけで、この、ツリー法と、 みたいな関数使って重力相互作 用を分けて、短い時間刻み必要な粒子というか粒子間相互作用を別途計算する、 というのが、まあ一応現状で最先端の計算法ということになるわね。 これは作者のところでドクターとった人の博士論文。

学生C:えー、それ、本当に大丈夫な方法なんですか?

11.5. ハミルトニアン分割の条件

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

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

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

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

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

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

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

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

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

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

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

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

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

わけ。で、多項式でなんかできそうだけど、多項式だとその次数まで微分した ら必ず境界のところで不連続になるでしょ?例えば、 で0、 としたら、 n-1次導関数までは で0だけど、 n次導関数は0じゃないわけ。

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

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

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

な関数というのがあって、 bump function (隆起関数) っていうの、って私も 昨日まで知らなかったんだけど。だから、これの不定積分をつくると、上の 欲しい性質満たして無限回微分可能にはなるのね。

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

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

学生C:解析的ってなんでしたっけ?

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

学生C:はあ。

赤木 :例えば、

はそういうものなんだって。図
23 みたいな感じ。

Figure 23: 隆起関数の例(ウィキペディアの隆起関数の項から)

学生C:よくこんな関数考えつくものですね、、、

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

とかね。これを右端で1になるように規格化すればいいわけ。

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

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

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

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

で、色々やってみる、ということでどうかしら?

学生C:はあ、、、なんか の積分とかでてきてますけどこれ どうするんですか?

赤木 :まあその辺も考えてみて。もちろん、がんばって紙と鉛筆で計算してくれてもいいわよ。

学生C:えー、ちょっと考えます。

11.6. 多項式クラス、多項式の乗算と積分

以下学生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
と。 T.zero は、この多項式クラスの要素が何か分からないけど、 zero をクラスメソッドとしてもってるとしてそれはゼロにあたるものである という約束で使いましょうということ。整数とか浮動小数点数以外に、 有理数でもこれで。あと、理屈としてはこれで多項式の多項式もできるから、 多変数にも使えるかな?例えば2次と3次の多項式の積は5次だけど、 係数の数としては3、4 で6になるから、まずゼロで埋まった サイズ self.size+a.size-1 の配列を作って、あとはそこに加算していくだけ。 to_poly は Array にメソッド追加して多項式に変換するやつ。これは MathVector を真似して作ってみたけどなんか一応これで動くんだけどもうちょっ といい方法がありそう。

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

  def integrate
    ([T.zero] +  self.map_with_index{|x,k| x/(k+1)}).to_poly
  end
なので、これだけ。あとべき乗も作ると、単に n 個の間の掛け算するだけで

  def ^(n : Int32)
    newp=self
    (n-1).times{newp *= self}
    newp
  end
かな。これももうちょっと賢く書ける気が。あと、多項式に実際に x の値い れて評価する関数がいるから、これを evaluate で、これホーナーの公式だっ たかな、 ax^2+bx+c = (ax+b)*x+c みたいな次数が高いほうから順番にするの を使うと、、、配列の順番ひっくり返す必要があるのか。 reverse って、、、 あ、あるから、これと reduce で

  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)})
必要なのは、 の積分で、あとこれを で 0から1に変わるように入力の値を細工し てから計算するのを前にやった を使って作るのかな。

 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
で、最初の行で多項式作ってそれ積分して、次ではそれの x=1 での値をだし てあとでスケールするのにつかう、3行目は x のほうのスケール用の係数を 計算しておく。で、あとは定義通りに範囲内なら poly.evaluate で、左側は0、 右側は1と。

実行してみると、

 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
グラフは
24 と。これで大丈夫かな?

Figure 24: 5次で作ってみたスイッチ関数のグラフ

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

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

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

でよかったはずで、さっきの話だと「位置を更新する」のところがルンゲクッ タで軌道積分するになればいいのか。だから、

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

ライブラリで作ったリープフロッグが2つあるけど、よくわかんないから まず x, v があるほうでやるか。元々が

  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
みたいにすればいいのかな。 q はルンゲクッタのほうの精度パラメータとし て。f, g がそれぞれリープフロッグ部分とルンゲクッタ部分の加速度。 これあと、gが0の時に無駄な計算しないで、というのをいれないともったいな いけど、それは後回しにすることで。あれ、これだとルンゲクッタに入る前の 速度の半分更新がないな。えーと。

  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
で、エネルギー誤差が図
25 です。

Figure 25: ハミルトニアン分割での積分結果

赤木 :離心率は、、、 0.99999 ってこと?それでちゃんと積分できたなら割 とあってるかもね。プログラムはどんなの?

学生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:
です。元々のケプラー軌道の積分のプログラムから結構あんまり変わってない ので、簡単に説明します。プログラム本体は154行目からで、 158行目からの

  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)))}
は、スイッチ関数と、それ使った、ハミルトニアン分割の f と g を定義してます。 create_switch_function は、多項式扱うクラス自分で作って、関数定義しま した。次数と、内側、外側の境界を与えると、それぞれで 0, 1 に、与えられ た次数まで微分可能につなげる関数を生成します。

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

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

  x, v = integrator.call(x,v,h)
だったのが

  x, v = hybrid(x, v, options.h, fout, fin, options.eta)
になってます。今のところ、積分公式はリープフロッグ+4時間ルンゲクッタだ けなので、そのへんはオプションないですが、リープフロッグの時間刻み、外 側の関数、内側の関数、内側の積分の時間刻みパラメータの順番です。 この hybrid の本体は

  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
です。システムじゃなくて x, v, f をもらうほうのリープフロッグが

  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

これあんまり美しくないんですが、ライブラリで定義したルンゲクッタを そのまま使うためとりあえずこうしてます。ライブラリで定義したルンゲクッ タはケプラー問題なら x と v を一緒に1つのベクトルと思う関数で、 これには MathVector クラスを渡す必要がありますが、ここまでは リープフロッグで x, v を別の Vector3 クラスでもってます。なので、加速 度を計算する関数とかもその辺変える必要がありました。それが

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

で、ここでは、 MathVector として、 Vector3 を2つ要素に持つようにしてま す。最初は6要素のベクトルを作って書いてたんですが、こっちでいけるかな と思ってやったらコンパイルできたので。多分大丈夫かなあと。

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

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

      dt = sqrt((x*x)/(v*v))*q
でだして、 t+dt が h を超えないように

      dt = h-t if t+dt > h
で小さくします。

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

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

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

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

11.8. 課題

  1. 関数 44 が両端で n 回導関数まで0であることを証明せよ。
  2. 4次ルンゲクッタの場合に、スイッチ関数の次数を変えるとエネルギー誤差の 振舞いがどう変わるか調べよ
  3. 8次程度のルンゲクッタを実装して、4次の場合と同様にエネルギー誤差と スイッチ関数の次数の関係を調べよ。

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