赤木 ： 今回は多体問題ということで。まあその、本当に大規模な数値計算す るなら最近だと全部自分でプログラム書くってだいぶ無理で、並列化のところ にパッケージ使うとか、初めから全部人が書いたのをもらってきて使うか なんだけど、基本的な仕掛けは自分でわかってて欲しいのね。

学生C：そういうものですか？

赤木 ：そういうものよ。ここではあんまり科学的な意味とかはを追求しないこ とにして、初期条件としては沢山星がある、例えば半径1の球の中の一様分布 とかにしてみて。

学生C： 質量と速度はどうしますか？あと1ってどういう単位ですか？光年とかですか？

赤木 ：ケプラー問題では太陽質量を1、重力定数 G も1にしたじゃない？ そういう感じで、星団の質量を 1、重力定数 G も1としてみて。 星の質量は、まずは N 個ならみんな等しく ね。

学生C：ケプラー問題だと、長さの単位は天文単位で、 が1年という ことだったじゃないですか？

赤木 ：あら、別にそんなことはないわよ。それだと地球の軌道しかし計算で きないみたいじゃない？

学生C：あれ、そうじゃないんですか？木星なら軌道半径を 5.2 にしないと、周期が12年とかにならないですよね？

赤木 ：あー、これ割合説明がというか理解が難しい？1が1天文単位だと思う と地球の軌道だけど、1が 5.2 天文単位だと思うと、半径1の軌道が木星の軌 道になるの。

学生C：？

赤木 ：うーん、これ意外に説明難しいわね。運動方程式は

今、物理学では単位といえば SI単位のことで、それ以外駄目みたいな話になっ てるから、単位系を勝手にとるってあんまり馴染みがないかもしれないけど、 天文学では系の大きさにみあった単位使わないとすごく計算間違いしやす くなるし、計算機でも上手く計算できなかったりするから、単位系とりなおす 話は重要なの。

SI だと、Gの値は、1 kg のものが1m 離れたところにある時にそれから受ける 重力加速度を を単位にして測ったものでしょ？そうすると、 例えば を単位としたらどうなるかしら？加速度の単位も になるわけね。

学生C： のところにあるから、 1m の時より に小さくなりま すね。で、加速度が1というのも を1と思うわけですね。だか らもうひとつ X がはいって G はに小さくなります。

赤木 ：多分あってるんじゃないかしら。じゃあ時間の単位のほうを、 1秒でなくて秒にしたら？

学生C： 加速度の単位が 1秒あたり、「1秒あたり動く距離」の変化だったの が、秒あたり、「秒あたり動く距離」の変化に変わるわけですね。例えば なら、1m/s は 2m/(2s) だから、 は秒後に で、それはで、だから で大きくな ればいいわけですね。

赤木 ：そう。そうすると、とを上手くとれば が変わらない でしょ？

学生C： ということですね。だから ですね。

赤木 ：そう。これケプラーの第三法則と同じ形でしょ？

学生C：第三法則ってなんでしたっけ？

赤木 ： 軌道周期が軌道長半径の 3/2乗ってやつ。

学生C：確かに同じですが、でも？

赤木 ：つまり、天文単位を1、1年をにとる代わりに、 木星の軌道長半径を1、木曜の軌道周期をにする単位系でもいいし、 他の惑星でもなんでもいいわけ。

星団でも、例えば1パーセクの中に合計で太陽質量の100倍があると思って、それを で計算するなら、時間の単位が決まるわけ。

あ、質量の単位を 倍にしたらどうなるかしら？

学生C：加速度は 倍だから、 も 倍？

赤木 ：そうね。だから、地球の周りの人工衛星でも月でも、木星の衛星でも、 土星のリングでも、 として、あと長さの単位を決めると、それから 時間の単位が決まるわけ。で、それはもちろん、円運動の周期が になるように、つまり、その長さのところでの円運動 の角速度が1になるように、ということになるわけ。

これは太陽系でなくてお互いの重力で運動する系でも式変形の原理は同じだか ら、同じように として、あと長さの単位を決めると、それから 時間の単位が決まるわけ。

学生C：でも星団みたいなのだと星は円運動するわけじゃないし、ケプラー運 動でもないですよね？時間の意味ってなんですか？

赤木 ：まず、形式的にはとにかく、実際の何年というのに対応するわけ。あと、星 の運動は色々だけど、平均の速度とか、系の特徴的な半径とか、そういうのは定義できるわね。

学生C：うーん、そういうものですか。

赤木 ：まあ、実際に星動かしてみるともうちょっと気分わかるかも。位置は 半径1の球の中の乱数でやるとして、速度はなんかパラメータ与えてやっぱり その半径の球の中で一様とかしてみて。熱力学的にはボルツマン分布とかがもっ ともらしいんだけど、星団でボルツマン分布って重力で束縛されない(無限遠 までいっちゃう)粒子がでてきてちょっと取り扱いが厄介なの。

学生C： だからといって速度分布や粒子分布を適当に作ることになんか意味あるんですか？

赤木 ：うん、そうね、これ、わりと意味ないわけでもないの。これは要する に、重力多体系というものはどう進化するか、という話ね。まず、理論的には、 沢山の粒子でできてるわけだから、統計力学で記述されるはずで、古典統計だ からボルツマン分布になるんだけど、さっきの話でボルツマン分布だと粒子が 無限遠にいっちゃうからそうならないわけ。でも、1つ1つの星についてみると、 ランダムに速度変化して、速度というかエネルギーの空間の中で拡散していく のね。なので、そのうちに、初期条件がなんでも、同じような分布になるとい うか、同じような進化をするようになるの。

学生C：はあ。

赤木 ：まあ、ちょっと実際にやってみて。

学生C：了解です。

赤木 ：あ、あとね、ケプラー問題では、軌道が楕円軌道で変化しなかったけ ど、多体問題ではみんな勝手に動くわけで、2つの粒子がすごく近くを通るこ とがあるのね。そうすると、タイムステップいくら短くしておいても破綻する の。

この、近くを通るのを、英語では close encounter っていうの。ちょっと余 談だけど、君くらい若いと知らないかもだけど、スピルバーグの映画に 「未知との遭遇」というのがあって、英語の原題は Close encounters of the third kind なのね。この第三種というのは誰かが勝手に作った宇宙人だかと の遭遇の分類だけど、 close encounter 自体は普通につかう言葉なのね。

学生C：はい？えーと、、、

赤木 ：なので、今回、昔から使われている方法なんだけど、重力ポテンシャ ルの式を、 じゃなくて としてみて。 は定数で、パラメータとして入力できるようにして。

学生C：そうすると、軌道とか変わるわけで、なんか違うものを計算している ことにならないですか？

赤木 ：もちろん違うんだけど、統計的には同じ性質になると信じましょうみ たいな感じ。まあ、 いくつか変えてみて、同じような性質がで てきたら大丈夫とかね。ちなみにこの を重力 softening (ソフトニング)というの。

学生C：そんなものでしょうか？

赤木 ：あ、もちろん、もうちょっと色々もっともらしい議論はできるし、じゃ あ結果が を0にもっていった極限で収束するのかとかそういう 話はできるし、しないといけないわけよ。まあそういうのはそのうちにね。

学生C：わかりました。では実際に作ってみます。

(このあと無限に色々あったが全部省略)

学生C：あ、なんか動いてるような気がします。プログラム particle0.cr と いうのを作りました。入力パラメータはこんな感じです

  -n     --numner-of-particles          int        2 Number of particles
  -s               --softening        float      0.0 Size of softening
  -e            --eccentricity        float      0.0 Eccentricity. Used only when n=2
  -d               --step-size        float     0.01 Size of timestep
  -T                --end-time        float      1.0 Time to stop integration
  -w             --window-size        float      1.5 Window size for plotting
  -e                     --ecc        float      0.0 Initial eccentricity of the orbit
  -v          --velocity-scale        float      0.5 Scaling factor for the initial velocity

で、アニメーションがでるんですが、2体だとそれっぽく回って、3とかそれ以 上だとなんか動いてました。

赤木 ：んーと、計算正しいかどうかってなんかわからない？まあその前にプ ログラム見せて。

学生C：はい。

    1:require "grlib"
    2:require "clop"
    3:include Math
    4:include GR
    5:
    6:optionstr= <<-END
    7:  Description: Test integrator  driver for Kepler problem
    8:  Long description:
    9:    Test integrator driver for Kepker problem
   10:    (c) 2020, Jun Makino
   11:
   12:  Short name:           -n
   13:  Long name:  --numner-of-particles
   14:  Value type:int
   15:  Default value: 2
   16:  Variable name: n
   17:  Description:Number of particles
   18:  Long description:     Number of particles
   19:
   20:  Short name:           -s
   21:  Long name:  --softening
   22:  Value type:float
   23:  Default value: 0.0
   24:  Variable name: eps
   25:  Description:Size of softening
   26:  Long description:     Size of softening
   27:
   28:  Short name:           -e
   29:  Long name:  --eccentricity
   30:  Value type:float
   31:  Default value: 0.0
   32:  Variable name: ecc
   33:  Description:Eccentricity. Used only when n=2
   34:  Long description:     Eccentricity. Used only when n=2
   35:
   36:  Short name:           -d
   37:  Long name:  --step-size
   38:  Value type:float
   39:  Default value: 0.01
   40:  Variable name: h
   41:  Description:Size of timestep
   42:  Long description:     Size of timestep
   43:
   44:  Short name: -T
   45:  Long name:  --end-time
   46:  Value type:int
   47:  Default value:1.0
   48:  Value type:float
   49:  Variable name:tend
   50:  Description:Time to stop integration
   51:  Long description:     Time to stop integration
   52:
   53:  Short name:-w
   54:  Long name:  --window-size
   55:  Value type:  float
   56:  Variable name:wsize
   57:  Default value:1.5
   58:  Description:Window size for plotting
   59:  Long description:
   60:    Window size for plotting orbit. Window is [-wsize, wsize] for both of
   61:    x and y coordinates
   62:
   63:  Short name:-e
   64:  Long name:--ecc
   65:  Value type:float
   66:  Default value:0.0
   67:  Variable name:ecc
   68:  Description:Initial eccentricity of the orbit
   69:  Long description:     Initial eccentricity of the orbit
   70:
   71:
   72:  Short name:-v
   73:  Long name:--velocity-scale
   74:  Value type:float
   75:  Default value:0.5
   76:  Variable name:vscale
   77:  Description:Scaling factor for the initial velocity
   78:  Long description:
   79:    Scaling factor for the initial velocity.
   80:    positions and velocities are set from random vectors within
   81:    spheres of radius 1 and vscale.
   82:
   83:END
   84:
   85:clop_init(__LINE__, __FILE__, __DIR__, "optionstr")
   86:options=CLOP.new(optionstr,ARGV)
   87:
   88:require "./vector3.cr"
   89:class Particle
   90:  property :m, :x, :v, :acc, :phi
   91:  def initialize(m : Float64=0, x : Vector3=Vector3.new(0,0,0),
   92:                 v : Vector3=Vector3.new(0,0,0))
   93:    @m = m; @x = x;  @v = v
   94:    @acc = Vector3.new(0,0,0); @phi=0.0
   95:  end
   96:  def self.random(m, rx, rv)
   97:    Particle.new(m.to_f, randomvector(rx), randomvector(rv))
   98:  end
   99:  def calc_gravity(p,eps2)
  100:    dr = @x - p.x
  101:    r2inv = 1.0/(dr*dr+eps2)
  102:    rinv = sqrt(r2inv)
  103:    r3inv = r2inv*rinv
  104:    @phi -= p.m*rinv
  105:    p.phi -= @m*rinv
  106:    @acc -= p.m*r3inv*dr
  107:    p.acc += @m*r3inv*dr
  108:  end
  109:    
  110:end
  111:
  112:class ParticleSystem
  113:  property :particles
  114:  def initialize()
  115:    @particles = Array(Particle).new
  116:  end
  117:  def +(p : Particle)
  118:    @particles.push p
  119:    self
  120:  end
  121:  def self.random(n, rx, rv)
  122:    ps=ParticleSystem.new
  123:    m = 1.0/n
  124:    n.times{ ps += Particle.random(m,rx,rv)}
  125:    ps.adjust_cm
  126:    ps
  127:  end
  128:  def self.twobody(x,v)
  129:    ps=ParticleSystem.new
  130:    n=2
  131:    m = 1.0/n
  132:    ps += Particle.new(m,x,v)
  133:    ps += Particle.new(m,-x,-v)
  134:    ps
  135:  end
  136:  def clear_gravity
  137:    @particles.each{|p|p.acc=Vector3.new; p.phi=0.0}
  138:  end
  139:  def calc_gravity(eps2)
  140:    @particles.each_with_index{|p,i|
  141:      ((i+1)..(particles.size-1)).each{|j|
  142:        p.calc_gravity(@particles[j], eps2)
  143:      }
  144:    }
  145:  end
  146:  def leapfrog(h,eps2)
  147:    @particles.each{|p|p.x += p.v*h + p.acc*(h*h*0.5)}
  148:    @particles.each{|p|p.v += p.acc*(h*0.5)}
  149:    clear_gravity
  150:    calc_gravity(eps2)
  151:    @particles.each{|p|p.v += p.acc*(h*0.5)}
  152:    self
  153:  end
  154:  def calc_cm
  155:    sumx=Vector3.new
  156:    sumv=Vector3.new
  157:    summ=0.0
  158:    @particles.each{|p|
  159:      sumx += p.m*p.x
  160:      sumv += p.m*p.v
  161:      summ +=p.m
  162:    }
  163:    {sumx*(1/summ),sumv*(1/summ)}
  164:  end
  165:  def adjust_cm
  166:    cmx,cmv = calc_cm
  167:    @particles.each{|p|
  168:      p.x -= cmx
  169:      p.v -= cmv
  170:    }
  171:  end
  172:end
  173:  def randomvector(r)
  174:    sqsum= r*r*2
  175:    v=Vector3.new
  176:    while sqsum > r*r
  177:      v = Array.new(3){rand(r)*2-r}.to_v
  178:      sqsum = v*v
  179:    end
  180:    v
  181:  end
  182:
  183:
  184:ENV["GKS_DOUBLE_BUF"]= "true" 
  185:
  186:#ps = ParticleSystem.random(2,1.0,0.4)
  187:if options.n == 2
  188:  ps = ParticleSystem.twobody([0.5*(1+options.ecc),0.0,0.0].to_v,
  189:                              [0.0,0.5*sqrt((1-options.ecc)/(1+options.ecc)),0.0].to_v)
  190:else
  191:  ps =ParticleSystem.random(options.n,1.0,options.vscale)
  192:end
  193:eps2= options.eps*options.eps
  194:ps.calc_gravity(eps2)
  195:wsize=options.wsize
  196:setwindow(-wsize, wsize,-wsize, wsize)
  197:
  198:time=0
  199:while time < options.tend+options.h/2
  200:  ps = ps.leapfrog(options.h, eps2)
  201:  time += options.h
  202:  clearws() 
  203:  box
  204:  setcharheight(0.05)
  205:  mathtex(0.5, 0.06, "x")
  206:  mathtex(0.06, 0.5, "y")
  207:  text(0.6,0.91,"t="+sprintf("%.3f",time))
  208:  setmarkertype(4)
  209:  setmarkersize(1)
  210:  polymarker(ps.particles.map{|p| p.x[0]}, ps.particles.map{|p| p.x[1]})
  211:  updatews() 
  212:end
  213:c=gets 
  214:

赤木 ：お願い。

学生C：最初の4行はグラフィクスとかコマンドラインとか数学ライブラリ使い ますですね。6から83行目まではオプション定義です。上の -n から -v まで。 88行は前に作ってもらった3次元ベクトルを使います。

89行から粒子クラスです。こういうの作るのがいいかどうかよくわかってない ですが、質量 m、 位置 x、速度 v、加速度 acc、ポテンシャル phi がメンバー 変数としてあって、それら全部普通に p.phi みたいなふうにアクセスできる としてます。

initialize は粒子作る時に Particle.new するわけですが、それから呼ばれ る関数でした。なにも引数なければ全部0に、ということでデフォルト値と、 あとここでは引数に型を与えます。

その次の self.random は、クラスメソッドというやつですね。 Particle.random(m, rx, rv) で、質量 m、位置と速度はそれぞれ半径 rx, rv の球の中の一様乱数です。

中では randomvector っていう関数で球の中の乱数生成します。

赤木 ：あら、その randomvector ってどこにあるの？

学生C：173行からですね。これは、球に外接する立方体の中の乱数は、3成分 独立に作ればいいのでそうして、半径1の外側ならその乱数使わないで、球の 中に入るまで繰り返す、という方法です。わりとこれ標準的みたいです。

赤木 ：そうね。 rejection method ね。乱数の半分くらい捨てるわけだから 計算量的にはもったいないけど、まあ、ちゃんと書けてる可能性が高いから。 そういえば、これ Crystal の面白いとこで、関数定義ってどこにあってもい いのね。

これが普通のコンパイル言語だと、よっぽど古いのでない限り、あるところで 使う関数はその前で定義がでてこないといけないのね。で、その定義って、 関数自体の型と引数の型なの。古い Fortran (Fortran77 とかその前)だと、 使うところで関数自体の型は宣言できるけど、引数は特に宣言する方法とか ないから、呼ぶところと定義するところでちゃんとつじつまがあってるかどう かはプログラム書く人任せで、もちろん間違えるわけ。で、 Algol 系の言語 のほとんどでは、関数定義と使うところで矛盾がないかどうかチェックするよ うになってるの。でも、そうすると、コンパイラが、ある関数が呼ばれるところを処 理するにはあらかじめその関数がどういうものか知らないといけないわけ。

単純には、これ、関数の定義が関数使うところより前にあれば いいわけで、Pascal とかは基本的にそうなってるのね。でも、 例えば2つの関数がお互いを自分の中で呼びたいとかあるとこれでは 不足だから、そういう時のために関数と引数の型だけ宣言できるような 文法があるの。で、最近の言語、特に C や C++ では、関数の引数と型の 宣言のことをプロトタイプ宣言といって、それが書いてあるファイルを、その 関数使うところの前で読み込む、というふうにするわけ。 C だと

  #include <stdio.h>

 #include <iostream>

学生C： Crystal でも require でなんかインクルードしてるからおんなじじゃないですか？

赤木 ： あ、でも、それ、必ずしも使うところの前にないといけないわけでは ないの。普通

  require "./vector3.cr"
  pp! Vector3.new(0)

  pp! Vector3.new(0)
  require "./vector3.cr"

C++ なんかだと、そういう、コンパイル時に呼び出し側から型がつたわってく る関数はテンプレート関数というものになって、「これはテンプレート関数で す」という形の宣言が必要ですごく複雑になって、それをもうちょっと簡単に するための新しい文法ができたりして一層わけがわからなくなんるんだけど、 Crystal だとその辺コンパイラが何とかするように頑張って作ってあって、見 かけ上実行時に型がわたってくる Ruby と概ね同じなんだけど、コンパイル言 語である、ということからこんなふうに Ruby でもできない気持ちが悪いこと もできるのね。

学生C：(あんまり良くわかってない)はあ、、、

赤木 ：と、ちょっと余談よねこれ。次の関数は？

学生C：次は calc_gravity で、2粒子間の重力を計算します。式の通りであん まりいうことないですが、自分と、もうひとつの同じ型の粒子で相互重力計算 して、両方の acc と phi にそれぞれが受ける加速度とポテンシャルを加算し ていきます。粒子クラスはこれだけです。

赤木 ：その次は？

学生C： 112行目からは ParticleSystem です。これ基本的には単に粒子の配 列をもっているクラスで、その配列に particles という名前をつけてます。 initialize は、なので、particles を長さ0のParticleクラスの配列、として ます。

その次の self.random と self.twobody は、それぞれ乱数でn個の粒子を作る のと、指定した位置、速度(一方は符号ひっくり返して)で2粒子、という関数 です。2粒子のほうは、このプログラムでちゃんとケプラー問題できるかどう かのチェック用です。

で、次の clear_gravity は、関数でなくてもいいんですが、各粒子の acc, phi をゼロにして、 calc_gravity は粒子間重力の計算です。

赤木 ： その次は？

学生C： はい、これがなんかちょっと嫌なんですが、この ParticleSystem ク ラス用のリープフロッグ積分関数になってます。

赤木 ：嫌ってのは？

学生C：時間積分法としては同じなのに、折角作った integratorlib.cr の中 の leapflog とはちがうものをもう一度書いてるわけで、なんか無駄というか 格好悪くないですか？

赤木 ：うーん、そうねえ、これ確かにちょっと難しいところね。 intergratorlib.cr の leapfrog は

  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(h,eps2)
    @particles.each{|p|p.x += p.v*h + p.acc*(h*h*0.5)}
    @particles.each{|p|p.v += p.acc*(h*0.5)}
    clear_gravity
    calc_gravity(eps2)
    @particles.each{|p|p.v += p.acc*(h*0.5)}
    self
  end

数学的というか、古典力学の定式化のほうからみると、シンプレクティック公 式だから、ハミルトニアンに対して定義されるわけで、ハミルトニアンって じゃない？ が一般化運動量で が一般化座標で。 そう思うと、intergratorlib.cr の leapfrog は割とそれに近い形で、 f が正準方程式の になってるわけね。 ParticleSystem のは、粒子、つまり Particle クラスのデータ構造が、 ハミルトニアンの型式とはあってないのね。

じゃあ、粒子クラスなんてのは止めてしまって、x とか v を 3N 次元のベク トルにしちゃえばいいか、っていうと、それもなんだかなあというところがあっ て。だって、「1つの粒子」ってやっぱり実体としてある気がするわけで、 それがデータ構造として存在しないのはなんか気持ちが悪いし、 初期条件を作るとか結果を解析するとかの時にはやっぱり粒子を単位に扱いた いじゃない。

学生C：はい、でも、それと今のシンプレクティック公式のライブラリは上手 くあわないですよね？

赤木 ：そうねえ、、、いろいろ考え方があると思うんだけど、一つは、 ParticleSystem みたいな物理系をあらわすクラスのほうに、 リープフロッグ公式を作るのに必要なだけのメソッドを用意する、ということ ね。つまり、リープフロッグで必要なのは

• 加速度を計算する
• 加速度と時間刻みを使って速度をアップデートする
• 速度と時間刻みを使って位置座標をアップデートする

 def leapfrog(s, h)
   s.inc_vel(h*0.5)
   s.inc_pos(h)
   s.calc_accel
   s.inc_vel(h*0.5)
 end

学生C： じゃあやってみます。

(次の日くらい)

一応できて動いてると思いますが、、、

赤木 ：あら、速いわね。プログラムは？

    1:require "grlib"
    2:require "clop"
    3:require "./integratorlib.cr"
    4:include Math
    5:include GR
    6:
    7:optionstr= <<-END
    8:  Description: Test integrator  driver for Kepler problem
    9:  Long description:
   10:    Test integrator driver for Kepker problem
   11:    (c) 2020, Jun Makino
   12:
   13:  Short name:           -n
   14:  Long name:  --numner-of-particles
   15:  Value type:int
   16:  Default value: 2
   17:  Variable name: n
   18:  Description:Number of particles
   19:  Long description:     Number of particles
   20:
   21:  Short name:           -s
   22:  Long name:  --softening
   23:  Value type:float
   24:  Default value: 0.0
   25:  Variable name: eps
   26:  Description:Size of softening
   27:  Long description:     Size of softening
   28:
   29:  Short name:           -e
   30:  Long name:  --eccentricity
   31:  Value type:float
   32:  Default value: 0.0
   33:  Variable name: ecc
   34:  Description:Eccentricity. Used only when n=2
   35:  Long description:     Eccentricity. Used only when n=2
   36:
   37:  Short name:           -d
   38:  Long name:  --step-size
   39:  Value type:float
   40:  Default value: 0.01
   41:  Variable name: h
   42:  Description:Size of timestep
   43:  Long description:     Size of timestep
   44:
   45:  Short name: -T
   46:  Long name:  --end-time
   47:  Value type:int
   48:  Default value:1.0
   49:  Value type:float
   50:  Variable name:tend
   51:  Description:Time to stop integration
   52:  Long description:     Time to stop integration
   53:
   54:  Short name:-w
   55:  Long name:  --window-size
   56:  Value type:  float
   57:  Variable name:wsize
   58:  Default value:1.5
   59:  Description:Window size for plotting
   60:  Long description:
   61:    Window size for plotting orbit. Window is [-wsize, wsize] for both of
   62:    x and y coordinates
   63:
   64:  Short name:-e
   65:  Long name:--ecc
   66:  Value type:float
   67:  Default value:0.0
   68:  Variable name:ecc
   69:  Description:Initial eccentricity of the orbit
   70:  Long description:     Initial eccentricity of the orbit
   71:
   72:
   73:  Short name:-v
   74:  Long name:--velocity-scale
   75:  Value type:float
   76:  Default value:0.5
   77:  Variable name:vscale
   78:  Description:Scaling factor for the initial velocity
   79:  Long description:
   80:    Scaling factor for the initial velocity.
   81:    positions and velocities are set from random vectors within
   82:    spheres of radius 1 and vscale.
   83:
   84:END
   85:
   86:clop_init(__LINE__, __FILE__, __DIR__, "optionstr")
   87:options=CLOP.new(optionstr,ARGV)
   88:
   89:require "./vector3.cr"
   90:class Particle
   91:  property :m, :x, :v, :acc, :phi
   92:  def initialize(m : Float64=0, x : Vector3=Vector3.new(0,0,0),
   93:                 v : Vector3=Vector3.new(0,0,0))
   94:    @m = m; @x = x;  @v = v
   95:    @acc = Vector3.new(0,0,0); @phi=0.0
   96:  end
   97:  def self.random(m, rx, rv)
   98:    Particle.new(m.to_f, randomvector(rx), randomvector(rv))
   99:  end
  100:  def calc_gravity(p,eps2)
  101:    dr = @x - p.x
  102:    r2inv = 1.0/(dr*dr+eps2)
  103:    rinv = sqrt(r2inv)
  104:    r3inv = r2inv*rinv
  105:    @phi -= p.m*rinv
  106:    p.phi -= @m*rinv
  107:    @acc -= p.m*r3inv*dr
  108:    p.acc += @m*r3inv*dr
  109:  end
  110:end
  111:
  112:class ParticleSystem
  113:  property :particles, :eps2
  114:  def initialize(eps2 : Float64 = 0.0 )
  115:    @particles = Array(Particle).new
  116:    @eps2 = eps2
  117:  end
  118:  def +(p : Particle)
  119:    @particles.push p
  120:    self
  121:  end
  122:  def self.random(n, rx, rv)
  123:    ps=ParticleSystem.new
  124:    m = 1.0/n
  125:    n.times{ ps += Particle.random(m,rx,rv)}
  126:    ps.adjust_cm
  127:    ps
  128:  end
  129:  def self.twobody(x,v)
  130:    ps=ParticleSystem.new
  131:    n=2
  132:    m = 1.0/n
  133:    ps += Particle.new(m,x,v)
  134:    ps += Particle.new(m,-x,-v)
  135:    ps
  136:  end
  137:  def calc_accel
  138:    @particles.each{|p|p.acc=Vector3.new; p.phi=0.0}
  139:    @particles.each_with_index{|p,i|
  140:      ((i+1)..(particles.size-1)).each{|j|
  141:        p.calc_gravity(@particles[j], @eps2)
  142:      }
  143:    }
  144:  end
  145:  def inc_vel(h)
  146:    @particles.each{|p|p.v += p.acc*(h)}
  147:  end
  148:  def inc_pos(h)
  149:    @particles.each{|p|p.x += p.v*h}
  150:  end
  151:  def calc_cm
  152:    sumx=Vector3.new
  153:    sumv=Vector3.new
  154:    summ=0.0
  155:    @particles.each{|p|
  156:      sumx += p.m*p.x
  157:      sumv += p.m*p.v
  158:      summ +=p.m
  159:    }
  160:    {sumx*(1/summ),sumv*(1/summ)}
  161:  end
  162:  def adjust_cm
  163:    cmx,cmv = calc_cm
  164:    @particles.each{|p|
  165:      p.x -= cmx
  166:      p.v -= cmv
  167:    }
  168:  end
  169:end
  170:  def randomvector(r)
  171:    sqsum= r*r*2
  172:    v=Vector3.new
  173:    while sqsum > r*r
  174:      v = Array.new(3){rand(r)*2-r}.to_v
  175:      sqsum = v*v
  176:    end
  177:    v
  178:  end
  179:
  180:
  181:ENV["GKS_DOUBLE_BUF"]= "true" 
  182:
  183:if options.n == 2
  184:  ps = ParticleSystem.twobody([0.5*(1+options.ecc),0.0,0.0].to_v,
  185:                              [0.0,0.5*sqrt((1-options.ecc)/(1+options.ecc)),0.0].to_v)
  186:else
  187:  ps =ParticleSystem.random(options.n,1.0,options.vscale)
  188:end
  189:ps.eps2= options.eps*options.eps
  190:ps.calc_accel
  191:wsize=options.wsize
  192:setwindow(-wsize, wsize,-wsize, wsize)
  193:
  194:time=0
  195:while time < options.tend-options.h/2
  196:#  SymplecticIntegrators.leapfrog(ps,options.h)
  197:  SymplecticIntegrators.yoshida4(ps,options.h)
  198:  time += options.h
  199:  clearws() 
  200:  box
  201:  setcharheight(0.05)
  202:  mathtex(0.5, 0.06, "x")
  203:  mathtex(0.06, 0.5, "y")
  204:  text(0.6,0.91,"t="+sprintf("%.3f",time))
  205:  setmarkertype(4)
  206:  setmarkersize(1)
  207:  polymarker(ps.particles.map{|p| p.x[0]}, ps.particles.map{|p| p.x[1]})
  208:  updatews() 
  209:end
  210:c=gets 
  211:

  gravity> diff particle0.cr particle1.cr
  2a3
  > require "./integratorlib.cr"
  109d109
  <     
  113,114c113,114
  <   property :particles
  <   def initialize()
  ---
  >   property :particles, :eps2
  >   def initialize(eps2 : Float64 = 0.0 )
  115a116
  >     @eps2 = eps2
  136c137
  <   def clear_gravity
  ---
  >   def calc_accel
  138,139d138
  <   end
  <   def calc_gravity(eps2)
  142c141
  <         p.calc_gravity(@particles[j], eps2)
  ---
  >         p.calc_gravity(@particles[j], @eps2)
  146,152c145,149
  <   def leapfrog(h,eps2)
  <     @particles.each{|p|p.x += p.v*h + p.acc*(h*h*0.5)}
  <     @particles.each{|p|p.v += p.acc*(h*0.5)}
  <     clear_gravity
  <     calc_gravity(eps2)
  <     @particles.each{|p|p.v += p.acc*(h*0.5)}
  <     self
  ---
  >   def inc_vel(h)
  >     @particles.each{|p|p.v += p.acc*(h)}
  >   end
  >   def inc_pos(h)
  >     @particles.each{|p|p.x += p.v*h}
  186d182
  < #ps = ParticleSystem.random(2,1.0,0.4)
  193,194c189,190
  < eps2= options.eps*options.eps
  < ps.calc_gravity(eps2)
  ---
  > ps.eps2= options.eps*options.eps
  > ps.calc_accel
  199,200c195,197
  < while time < options.tend+options.h/2
  <   ps = ps.leapfrog(options.h, eps2)
  ---
  > while time < options.tend-options.h/2
  > #  SymplecticIntegrators.leapfrog(ps,options.h)
  >   SymplecticIntegrators.yoshida4(ps,options.h)

  module SymplecticIntegrators
    extend self
    def leapfrog(s, h)
      s.inc_vel(h*0.5)
      s.inc_pos(h)
      s.calc_accel
      s.inc_vel(h*0.5)
    end
    D1 = 1.0 / (2-exp(log(2.0)/3))
    D2 = 1 - 2*D1

   def yoshida4(s,h)
      leapfrog(s,h*D1)
      leapfrog(s,h*D2)
      leapfrog(s,h*D1)
    end
    
    def  yoshida6(x,v,h,f)
      d = {0.784513610477560, 0.235573213359357,
           -1.17767998417887, 1.31518632068391};
      4.times{|i|leapfrog(s,h*d[i])}
      3.times{|i|leapfrog(s,h*d[2-i])}
    end
  end

赤木 ：動くのはいいとしてして、これちゃんと正しく計算できてる？

学生C：と思いますが、、、

赤木 ：こういうのは思いますでは駄目で、根拠がいるのよ。だって、 シミュ レーションしてこう答がでました、っていっても、それだけではそれが本当か どうかわからないでしょ？

学生C：でも、元々解析解がなくて答がわからないからシミュレーションする んだから、その答が正しいかどうかってわからないんではないですか?

赤木 ：それはそうなんだけど、それでも、全然間違ってればわかるチェック、 というのはあるの。例えば、古典力学系だから保存量はあるでしょ？

学生C：3次元なので、エネルギー、全体の運動量、全体の角運動量ですね？

赤木 ：そう。だから、その辺のチェックは最低限して欲しいの。あ、あと、 恒星系だと、「重力ポテンシャルエネルギーと運動エネルギーの比」が大事な 量なの。これは「ビリアル比」っていって、自己重力系の定常状態だと変化し ない、というのが証明できるの。

学生C：すみません、定常状態ってなんですか？

赤木 ：分布関数が、、、あ、えーと、つまり、統計力学と同じで、星の数が 無限に大きな極限を考えると、恒星系って、位置と速度の6次元空間の中での 密度分布とみなせるでしょ？

学生C：でしょって言われても、、、統計力学なら原子はみんな同じですが、 星はひとつひとつ違うじゃないですか?統計力学って成り立たつんですか?

赤木 ：今は星は重力相互作用する質点だから、違うのは質量だけじゃない？ そうすると、星の数が無限に多いなら、分布関数が質量にも依存するとして 統計力学の式をたてればいいわけ。ただ、まずここで問題なのは、 統計力学的な平衡、熱平衡ね、それじゃなくて、力学平衡というものなの。

学生C：どう違うんですか?

赤木 ：統計力学的平衡はエントロピー最大なんだけど、力学平衡は単に分布 関数が定常というだけね。普通の気体だと、力学平衡は圧力がつりあってる状 態、静水圧平衡にあたるものね。静水圧平衡だと温度分布はあってもいいわけ。

学生C： でも星は動くわけじゃないですか?

赤木 ：そう。だから、割合その辺違うのね。気体分子運動論から流体力学を 導く時には、分子同士は頻繁に衝突するから、局所的には統計力学的に平衡だ として、局所的なマクロな物理量を導くんだけど、星同士は滅多にぶつからな いから局所的にも大局的にも熱平衡にならないの。

1. 粒子数 100 程度で、時刻 10 程度まで積分した時のエネルギーの誤差の絶 対値の最大値を、タイムステップの関数としてプロットしてみよ。ソフト ニングは 0.01 で固定でよいが、変化させたらどのように変わるかも 余力があれば調べてみよ。

2. 全運動量、全角運動量について、保存精度を調べよ。

3. 4次シンプレクティック公式も選択できるようにプログラムを拡張し、 課題 1, 2 についてどうなるか調べよ。

10.6. まとめ

1. シンプレクティック積分公式を、ハミルトン力学系を表現するクラスに対 して一般的に適用できる形のライブラリとして作成した。
2. 重力多体問題を例に、上の公式を適用して系の性質を学んだ。
3. 重力多体問題では、ソフトニングを導入してポテンシャルの発散をふせいだ。

10.7. 参考資料

http://jun-makino.sakura.ne.jp/kougi/stellar_dynamics/index.html|

Previous

ToC

Next