Previous ToC Next

3. 「とんでる力学」8 月分

3.1. プログラムと解説

赤木: というわけで、6月に続いてプログラムとかの解説ね。

学生S: えーと、そうですね。今回結構大変なんですよね。


   august.rb : numerical integration for three bodies

  require 'vector.rb'
  include Math
  class Body
    attr_accessor :mass, :pos, :vel, :acc, :pot, :withring
    def initialize
      @pos=[]
      @vel=[]
      @withring=nil
    end
    def froms(s)
      a=s.split.collect{|x| x.to_f}
      p a
      n=(a.size-1)/2
      @mass = a[0]
      @pos  = a[1..n].to_v
      @vel  = a[(n+1)..(2*n)].to_v
      p self
    end

    def convertfromelements(a,e)
      @pos[0]=a*(1+e)
      @pos[1]=0
      @vel[0]=0
      @vel[1]= sqrt(2*(1/@pos[0]-0.5/a))
      @pos=@pos.to_v
      @vel=@vel.to_v
    end
    def fromvectortocontactelements(mass,pos,vel)
      raise "contactelements works only for 2D vector" if pos.size != 2
      r=sqrt(pos*pos)
      be = -mass/r+0.5*(vel*vel)
      # be = -mass/2a, a = -mass/2be
      a = -mass/(be*2)
      vcirc = sqrt(mass/a)
      jmax = a*vcirc
      j = pos[0]*vel[1]-pos[1]*vel[0]
      vbyj = vel*j
      vbyh=[vbyj[1],-vbyj[0]].to_v
      e = vbyh-pos*(1/r)
      eabs=sqrt(e*e)
      argp = atan2(-e[1],-e[0])
      [e,eabs,a,argp]
    end

    def contactelements(ref)
      fromvectortocontactelements(@mass+ref.mass,
                                  @pos-ref.pos,@vel-ref.vel)
    end
      #
   relation between J and e

   assume a=1,j=1
   ra = (1+e)
   E= -1/r+v^2/2
   va = sqrt(1/ra - 1/2a)
    def accumulate_force(otherbody)
      dr= @pos-otherbody.pos
      r2 = dr*dr
      r = sqrt(r2)
      r3 = r2*r
      mfact=@mass/r3
      @acc += dr*(-otherbody.mass/r3)
      otherbody.acc += dr*(@mass/r3)
      @pot -= otherbody.mass/r
      otherbody.pot -= @mass/r
    end

    def forcefromring
      dummyp = Body.new
      dummyp.mass=@ringmass/@ringndivide
      dummyp.acc=@acc
      dummyp.pot=@pot
      dt = 2*Math::PI/@ringndivide
      for i in 0..(@ringndivide-1)
        theta = dt*i
        dummyp.pos[0]=@ringradius*cos(theta)
        dummyp.pos[1]=@ringradius*sin(theta)
        accumulate_force(dummyp)
      end
    end
    def ringinit(radius,mass,ndivide)
      @ringradius,@ringmass,@ringndivide = radius,mass,ndivide
      @withring=1
    end

    def makenull
      null = ([0.0]*(@pos.size)).to_v
    end
  end
  class Nbody
    attr_accessor :ba,  :time

    def initialize
      @ba=[]
    end
    def calculate_force
      @ba.each{|x|
        x.pot=0; x.acc = x.makenull
        x.forcefromring if x.withring
      }
      for i in 0..(@ba.size-2)
        for j in (i+1)..(@ba.size-1)
          @ba[i].accumulate_force(@ba[j])
        end
      end
    end

    def leapfrog(first_call,dt)
      if first_call
        calculate_force
        first_call=nil
      end
      @ba.each{|p| p.vel += p.acc*0.5*dt}
      @ba.each{|p| p.pos += p.vel*dt}
      calculate_force
      @ba.each{|p| p.vel += p.acc*0.5*dt}
    end
    def yoshida6b(first_call,dt)
      d =   [0.784513610477560e0,
        0.235573213359357e0,
        -1.17767998417887e0,
        1.31518632068391e0]
      for i in 0..2 do leapfrog(first_call,dt*d[i]) end
      leapfrog(first_call,dt*d[3])
      for i in 0..2 do leapfrog(first_call,dt*d[2-i]) end
    end

    def calculate_energy
      ke=pe=0
      @ba.each{|p|
        ke += 0.5*p.mass*(p.vel*p.vel)
        pe += p.pot*p.mass
      }
      pe *= 0.5
      [ke+pe,ke,pe]
    end
    def cmadjust
      cmp=@ba[0].makenull
      cmv=@ba[0].makenull
      cmm=0
      @ba.each{|p|
        mass = p.mass
        cmp += p.pos*mass
        cmv += p.vel*mass
        cmm += mass
      }
      cmp *= 1/cmm
      cmv *= 1/cmm
      @ba.each{|p|
        p.pos -= cmp
        p.vel -= cmv
      }
    end

    def plot
      @ba.each{|p|
        relocate(*p.pos)
        dot
      }
    end
    def run(dt,tend,interval,plotendtime=0)
      steps = ((tend+dt*0.5)/dt).to_i
      calculate_force
      e0= calculate_energy
      print "Initial energy = #{e0[0]}\n"
      i=0
      @time=0
      first_call = 1
      steps.times{
        yoshida6b(first_call,dt)
        @time += dt
        i+= 1
        if i==interval or @time < plotendtime
          i=0
          yield
        end
      }
      e= calculate_energy
      print "Final energy = #{e[0]}, de = #{e0[0]-e[0]}\n"
    end
  end

  the following is a test config for circular hierarchical binary
  if nil
    b.ba= [Body.new,Body.new]
    b.ba[0].froms("0.5 1 0.0 0   1")
    b.ba[1].froms("0.5 0 0   0   0")
    p b
    b.cmadjust
    b.ba += [Body.new]
    #
    # M=2, R=5, a= 2/25, v^2/=a v=sqrt(2/5)
    # pe = 2/5 ke=1/5
    b.ba[2].froms("1 5 0   0   0")
    b.ba[2].vel[1]=sqrt(2.0/5.0)
    b.cmadjust
  end
  def marssetup(jupiterascale=1,marsecc=0.093,jupiterecc=0)
    b = Nbody.new
    b.time=0
    b.ba= [Body.new,Body.new,Body.new]
    b.ba[0].froms("1 0.0 0.0  0  0")
    b.ba[1].convertfromelements(5.203/1.524*jupiterascale,jupiterecc)
    b.ba[1].mass=0.001
    b.ba[2].convertfromelements(1,marsecc)
    b.ba[2].mass=0.00000
    p b
    b.cmadjust
    b
  end

  class Nbody
    def connectorplot(dt,tconnect)
      if @time <= dt or @time > tconnect
        relocate(*ba[2].pos)
        dot
      else
        draw(*ba[2].pos)
      end
    end
  end
  unless File.exist?("august-fig1b.ps")
    printer august-fig1b.ps/vcps
    term
    endtime=10000
    endtime=100
    tpend=6.38
    b=marssetup
    square
    expand 1.5
    viewport 0 0.45 0.55 1
    limit -1.2 1.2  -1.2 1.2
    box BCTSV  BCNTSV
    ylabel y
    expand 0.1
    dt = 0.07
    b.run(dt,endtime,50,tpend){
      b.connectorplot(dt,tpend)
    }

    viewport 1.2 2.2 0 1
    b=marssetup
    b.ba[1].mass*10
    expand 1.5
    limit -1.2 1.2  -1.2 1.2
    box BCTSV  BCTSV
    expand 0.1
    b.run(dt,endtime,50,tpend){
      b.connectorplot(dt,tpend)
    }
    viewport -1.2 -0.2 -1.2 -0.2
    b=marssetup(0.5)
    expand 1.5
    limit -1.2 1.2  -1.2 1.2
    box BCNTSV  BCNTSV
    xlabel x
    ylabel y
    expand 0.1
    b.run(dt,endtime,50,tpend){
      b.connectorplot(dt,tpend)
    }
    expand 1

    viewport 1.2 2.2 0 1
    b=marssetup(0.5)
    b.ba[1].mass*10
    expand 1.5
    limit -1.2 1.2  -1.2 1.2
    box BCNTSV  BCTSV
    xlabel x
    expand 0.1
    b.run(dt,endtime,50,tpend){
      b.connectorplot(dt,tpend)
    }
    expand 1
    pgend
    psfix
  end
  unless File.exist?("august-fig2.ps")
    printer august-fig2.ps/vcps
    b = marssetup
    square
    viewport 0.1 1.1 0.1 1.1
    limit 0 500 -0.005 0.015
    expand 1.2
    box BCNSTV  BCNSTV
    expand 2
    xlabel T
    reloc -100 0.007
    label \\gh
    a=b.ba[2].contactelements(b.ba[0])
    p b.time
    p a[3]
    reloc b.time a[3]
    b.run(0.07,500,5){
      a=b.ba[2].contactelements(b.ba[0])
      draw(b.time,a[3])
    }

    b = marssetup
    b.ba[1].mass*=0.5
    a=b.ba[2].contactelements(b.ba[0])
    reloc b.time a[3]
    b.run(0.07,500,5){
      a=b.ba[2].contactelements(b.ba[0])
      draw(b.time,a[3])
    }
    b = marssetup
    b.ba[1].mass*=0
    a=b.ba[2].contactelements(b.ba[0])
    reloc b.time a[3]
    b.run(0.07,500,5){
      a=b.ba[2].contactelements(b.ba[0])
      draw(b.time,a[3])
    }
    expand 1
    pgend
    psfix
  end

  unless File.exist?("august-fig3.ps")
    printer august-fig3.ps/vcps
    $plotlylimit = -6
    $ymin = exp($plotlylimit*log(10.0))
    $lymax =   $plotlylimit
    def logp(x,y)
      p x
      p y
      if y > 0
        ly = log10(y)
        if ly > $lymax
          draw x ly
          $lymax = ly
        end
      end
    end
    square
    viewport 0.1 1.1 0.1 1.1
    limit 0 500 $plotlylimit -0.5
    expand 1.2
    box BCNSTV  BCNSTVL
    expand 2
    xlabel T
    ylabel \\gh
    b = marssetup
    a=b.ba[2].contactelements(b.ba[0])
    reloc b.time $plotlylimit
    b.run(0.07,500,5){
      a=b.ba[2].contactelements(b.ba[0])
      logp(b.time, a[3])

    }
    $lymax =   $plotlylimit
    b = marssetup
    b=marssetup(0.5)
    a=b.ba[2].contactelements(b.ba[0])
    reloc b.time $plotlylimit
    b.run(0.07,500,5){
      a=b.ba[2].contactelements(b.ba[0])
      logp(b.time, a[3])
    }

    $lymax =   $plotlylimit
    b = marssetup
    b=marssetup(2)
    a=b.ba[2].contactelements(b.ba[0])
    reloc b.time $plotlylimit
    b.run(0.07,500,5){
      a=b.ba[2].contactelements(b.ba[0])
      logp(b.time, a[3])
    }
    $lymax =   $plotlylimit
    b = marssetup
    b=marssetup(4)
    a=b.ba[2].contactelements(b.ba[0])
    reloc b.time $plotlylimit
    b.run(0.05,500,5){
      a=b.ba[2].contactelements(b.ba[0])
      logp(b.time, a[3])
    }
  end

  unless File.exist?("august-fig4.ps")
    printer august-fig4.ps/vcps
    square
    viewport 0.1 1.1 0.1 1.1
    limit 0 500 -0.005 0.01
    expand 1.2
    box BCNSTV  BCNSTV
    expand 2
    xlabel T
    ylabel \\gh
    b = marssetup(1,0.1)
    a=b.ba[2].contactelements(b.ba[0])
    reloc b.time 0
    b.run(0.07,500,5){
      a=b.ba[2].contactelements(b.ba[0])
      draw b.time a[3]

    }
    b = marssetup(1,0.4)
    a=b.ba[2].contactelements(b.ba[0])
    reloc b.time 0
    b.run(0.03,500,5){
      a=b.ba[2].contactelements(b.ba[0])
      draw b.time a[3]

    }
  end
  unless File.exist?("august-fig5.ps")
    printer august-fig5.ps/vcps
    term
    square
    viewport 0.1 1.1 0.1 1.1
    tend = 500
    limit 0 tend -0.005 0.01
    expand 1.2
    box BCNSTV  BCNSTV
    expand 2
    xlabel T
    ylabel \\gh

    color 1
    b = marssetup
    reloc b.time 0
    b.run(0.07,tend,5){
      a=b.ba[2].contactelements(b.ba[0])
      draw b.time a[3]
    }

    b = marssetup
    reloc b.time 0
    b.ba[1].mass=0
    b.ba[2].ringinit(5.203/1.524,0.001,32)
    color 2
    b.run(0.07,tend,5){
      a=b.ba[2].contactelements(b.ba[0])
      draw b.time a[3]

    }
    color 1
  end

赤木: これちょっと長いどころの騒ぎじゃないわね。 500行近いわよ。

学生S: これ1つで8月号の図全部書くんだから、しょうがないと思いま す。頭から見ていってもいいんですが、どう使われるかをみないとわけがわか らないところがおおいと思いますから、実際に図を書くところから見ていき ませんか?

赤木: そうね。じゃあ、最初の図のところね。

学生S: はい、メインプログラムというか、トップレベルの、関数とか の中でなくて図1を書いてるのはここです。


 unless File.exist?("august-fig1b.ps")
   printer august-fig1b.ps/vcps
   term
   endtime=10000
   endtime=100
   tpend=6.38
   b=marssetup
   square
   expand 1.5
   viewport 0 0.45 0.55 1
   limit -1.2 1.2  -1.2 1.2
   box BCTSV  BCNTSV
   ylabel y
   expand 0.1
   dt = 0.07
   b.run(dt,endtime,50,tpend){
     b.connectorplot(dt,tpend)
   }

   viewport 1.2 2.2 0 1
   b=marssetup
   b.ba[1].mass*10
   expand 1.5
   limit -1.2 1.2  -1.2 1.2
   box BCTSV  BCTSV
   expand 0.1
   b.run(dt,endtime,50,tpend){
     b.connectorplot(dt,tpend)
   }
   viewport -1.2 -0.2 -1.2 -0.2
   b=marssetup(0.5)
   expand 1.5
   limit -1.2 1.2  -1.2 1.2
   box BCNTSV  BCNTSV
   xlabel x
   ylabel y
   expand 0.1
   b.run(dt,endtime,50,tpend){
     b.connectorplot(dt,tpend)
   }
   expand 1

   viewport 1.2 2.2 0 1
   b=marssetup(0.5)
   b.ba[1].mass*10
   expand 1.5
   limit -1.2 1.2  -1.2 1.2
   box BCNTSV  BCTSV
   xlabel x
   expand 0.1
   b.run(dt,endtime,50,tpend){
     b.connectorplot(dt,tpend)
   }
   expand 1
   pgend
   psfix
 end

赤木: で、これが何をしているの?

学生S: 要するに、"august-fig1b.ps" という ファイルに図1を書く、というものですね。

赤木: 1b なのはなんか深いわけがあるの?

学生S: えーと、すみません、だいぶ前の話なんで忘れました。という か、これ編集から注文が来て書き直したんで変わってるんだったと思います。

赤木: プログラムのほうが残してあれば、図とかのファイルの名前は 変えなくてもいいような気がするけど、まあいいわ。最初の

  unless File.exist?("august-fig1b.ps")

っていうのは、図があれば新しくは作らないというのね?

学生S: そうです。トップレベルで

  unless File.exist?("foo.ps")
     ...
  end

というのが一杯あって、それぞれ別の図を書くというふうになってます。図毎 にプログラムわけてもいいんですが、そうするとかえって複雑になるような気 がして。

赤木: まあ、原稿一つの図が一つのプログラムで全部かけるというの は便利は便利かもね。 次の

  printer august-fig1b.ps/vcps

は Rongo のコマンドになってるわけね?普通の Ruby/PGPLOT ではなくて。

学生S: です。

赤木: その後いくつか変数を初期化して、、、あ、

  b=marssetup

ってのは何?

学生S: これは時間積分の初期設定です。marssetup の実体は


def marssetup(jupiterascale=1,marsecc=0.093,jupiterecc=0)

  b = Nbody.new
  b.time=0
  b.ba= [Body.new,Body.new,Body.new]
  b.ba[0].froms("1 0.0 0.0  0  0")
  b.ba[1].convertfromelements(5.203/1.524*jupiterascale,jupiterecc)
  b.ba[1].mass=0.001
  b.ba[2].convertfromelements(1,marsecc)
  b.ba[2].mass=0.00000
  p b
  b.cmadjust
  b

end


で、Nbody というクラスの変数 b を作って、時刻 @time を 0 にして粒子3個 の質量、位置、速度を初期設定します。粒子は @ba という配列にはいってい て、まず Body.new で初期化してから、粒子0(太陽)は質量 1、位置、速度は 0 で、 木星、火星(1、2)は軌道長半径と離心率から位置、速度を計算します。計算す るのは


  def convertfromelements(a,e)
    @pos[0]=a*(1+e)
    @pos[1]=0
    @vel[0]=0
    @vel[1]= sqrt(2*(1/@pos[0]-0.5/a))
    @pos=@pos.to_v
    @vel=@vel.to_v
  end

で、最初は x 軸上で、遠日点であるとして位置、速度を計算するようになっ てます。

赤木: その後しばらくは画面の設定ね。それはいいとして

  b.run(dt,endtime,50,tpend){
    b.connectorplot(dt,tpend)
  }

だけど、これで {} の中は?なんかいかにもこれで絵を書いてるって感じだけ ど、、、

学生S: そうですね、これは Ruby のちょっと便利な機能を使ってみま した。 Nbody#run (Nbody というクラスの run というメソッドという意味) は


  def run(dt,tend,interval,plotendtime=0)
    steps = ((tend+dt*0.5)/dt).to_i
    calculate_force
    e0= calculate_energy
    print "Initial energy = #{e0[0]}\n"
    i=0
    @time=0
    first_call = 1
    steps.times{
      yoshida6b(first_call,dt)
      @time += dt
      i+= 1
      if i==interval or @time < plotendtime
        i=0
        yield
      end
    }
    e= calculate_energy
    print "Final energy = #{e[0]}, de = #{e0[0]-e[0]}\n"
  end

というものなんですが、ここで

        yield

というのがありますよね?プログラムのここにくると、 {} の中身のコードが 実行されるんです。Ruby ではこういうのをブロックとかコードブロックといいます。

赤木: えーと、関数をパラメータで渡すみたいなもの?

学生S: まあ、この場合はそうともいえますね。中でメソッド一つ呼ん でるだけですし。でも、ここでしか使わないものだといちいち関数程度しなく ても使えるから簡単でわかりやすいし。あと、これがちょっと面白いのは、 ブロックが、呼ばれた関数、この場合 run ですね、のコンテキストではなく て呼んだほう、この場合はトップレベルですが、のほうのコンテキストで実行 されるということです。つまり、呼んでほうのローカル変数とかは見えるけど 呼ばれたほうのは見えないんです。

赤木: それって不便っていうか、やりたいことと違うんではないの?

学生S: 実はこの場合にはちょっとそういうところもあるんですけど。 一般的にはこっちのほうが便利だと思います。Ruby のいろんな機能自体がこ れを使って実現されてますから。例えば each とかって使ったと思いますけど、 これは実はメソッドなわけでブロックが呼んだ側のコンテキストで実行されな いと each に渡すブロックに何も書けなくなってしまいます。

赤木: それはそうね。

学生S: で、先に、、、といっても、 run を呼んだ後はもう絵を書く 処理を終えてるだけですね。 run のほうの中身にいってみると、

まず

    steps = ((tend+dt*0.5)/dt).to_i

は繰り返しの回数計算ですね。で、ここではいわゆる吉田6次公式というシン プレクティック積分公式を使うのですが、そのためには最初に加速度を計算し ておく必要があるので

    calculate_force

で加速度を計算してます。これは後で説明します。次はチェック用ですが、

    e0= calculate_energy
    print "Initial energy = #{e0[0]}\n"

で系の全エネルギーを計算してます、後は繰り返しです。

    i=0
    @time=0
    first_call = 1
    steps.times{
      yoshida6b(first_call,dt)

これが時間積分ですね。

      @time += dt
      i+= 1
      if i==interval or @time < plotendtime

ここで、条件がややこしいのは、編集から最初の1周期は線をひいて欲しいと いう要請があったからです。それ以降は何回かに一回だけ点を打つわけです。 その判断は、問題の

        yield

で呼ばれた

    b.connectorplot(dt,tpend)

のほうがしているわけです。

赤木: なるほど。まあ、プログラム自体は割合素直ね。とすると後は

    calculate_force
    calculate_energy
    yoshida6b(first_call,dt)
    connectorplot(dt,tpend)

ね、順番にいってみて。

学生S: まず calculate_force ですね。


   def calculate_force
     @ba.each{|x|
       x.pot=0; x.acc = x.makenull
       x.forcefromring if x.withring
     }
     for i in 0..(@ba.size-2)
       for j in (i+1)..(@ba.size-1)
         @ba[i].accumulate_force(@ba[j])
       end
     end
   end

 うーん、見ての通りですが、、、

赤木: 最初は全粒子についてポテンシャル、加速度を 0 にしてるのね。 makenull ってのは?

学生S: 0 ベクトルを作るだけの関数です。その次に、後で木星をリン グで置き換えるのがあるのでその時にはリングからの力を計算します。これは こんなの


   def forcefromring
     dummyp = Body.new
     dummyp.mass=@ringmass/@ringndivide
     dummyp.acc=@acc
     dummyp.pot=@pot
     dt = 2*Math::PI/@ringndivide
     for i in 0..(@ringndivide-1)
       theta = dt*i
       dummyp.pos[0]=@ringradius*cos(theta)
       dummyp.pos[1]=@ringradius*sin(theta)
       accumulate_force(dummyp)
     end
   end

で、円環上にダミー粒子を並べて、そこからの力を計算するというのです。

赤木: これ、呼ばれる度にダミー粒子作り直すの?

学生S: そうです。無駄は無駄なんですが、そんなに時間かからなかっ たというか、書き直すのにかかりそうな時間より速く計算終わったので。

赤木: まあ、いいわ。 accumulate_force が実際に加速度を計算して るわけね。

学生S: そうです。


   def accumulate_force(otherbody)
     dr= @pos-otherbody.pos
     r2 = dr*dr
     r = sqrt(r2)
     r3 = r2*r
     mfact=@mass/r3
     @acc += dr*(-otherbody.mass/r3)
     otherbody.acc += dr*(@mass/r3)
     @pot -= otherbody.mass/r
     otherbody.pot -= @mass/r
   end

です。これはもう式の通りです。引数に相手の粒子をとって、そちらにも積算 します。ですから、 calculate_force では内側のループは i+1 から回すよう になってます。

calculate_energy は、


   def calculate_energy
     ke=pe=0
     @ba.each{|p|
       ke += 0.5*p.mass*(p.vel*p.vel)
       pe += p.pot*p.mass
     }
     pe *= 0.5
     [ke+pe,ke,pe]
   end

これも見ての通りですが、、、あ、最後の

    [ke+pe,ke,pe]

は、これらを要素にする配列を作っていて、このメソッドの返り値がこの配列 になります。こういうふうに気軽に配列とかを返せるのは便利だと思います。時間積分は


   def yoshida6b(first_call,dt)
     d =   [0.784513610477560e0,
       0.235573213359357e0,
       -1.17767998417887e0,
       1.31518632068391e0]
     for i in 0..2 do leapfrog(first_call,dt*d[i]) end
     leapfrog(first_call,dt*d[3])
     for i in 0..2 do leapfrog(first_call,dt*d[2-i]) end
   end

で、これは都合7回、時間刻みをいろいろ変えた leapfrog を呼ぶというもの ですね。leapfrog は


   def leapfrog(first_call,dt)
     if first_call
       calculate_force
       first_call=nil
     end
     @ba.each{|p| p.vel += p.acc*0.5*dt}
     @ba.each{|p| p.pos += p.vel*dt}
     calculate_force
     @ba.each{|p| p.vel += p.acc*0.5*dt}
   end

で、公式通りです。グラフを書く connectorplot は


   def connectorplot(dt,tconnect)
     if @time <= dt or @time > tconnect
       relocate(*ba[2].pos)
       dot
     else
       draw(*ba[2].pos)
     end
   end

で、指定した時刻の前なら線を引いて、後なら点を打つというだけです。

   relocate(*ba[2].pos)

で、

   *ba[2].pos

というのは見なれないと思いますが、これは、ここで配列が、なんというのか、 展開されて、要素が2つなので

   ba[2].pos[0], ba[2].pos[1]

と書くのと同じことになります。

赤木: んー、なくても困らないような機能ね。

学生S: 困ることはないかもしれないですけど、短く書けるのはいいと 思います。

赤木: まあ、それはそうね。他の図は?

学生S: 全部 run でブロックを渡して書くというのは今回同じなので、 ここまででいいことにしませんか?

赤木: ちょっと手抜きな気がするけど、、、

学生S: でも、ほら、これ8月号分で、9月号が今週でちゃいますよね? そっち書いて時間があったらこっちももうちょっと追加ということで、、、

赤木: そのほうがいいのは確かね。じゃあ、そうさせてもらいましょ。 あ、ルンゲ・レンツベクトルの話は?

学生S: え、でも、吉田さんの教科書とか見てもらったほうがいいと思 いますけど、、、

赤木: まね。リンクは?

学生S:これですね。

www.iwanami.co.jp/.BOOKS/00/5/0079560.html

Previous ToC Next