1 万有引力による惑星の軌道

　太陽と惑星の間には、万有引力が働いています。

地球の重さを \(m_e\) 、太陽の重さを \(m_s\)、距離を \(r\)、万有引力定数を \(G=6.67259\times 10^{-11}~［\mathrm{N\cdot m^2/kg^2}］\)として
　　　\(F=\dfrac{Gm_em_s}{r^2}\)
地球の位置を \((x,y)\)、太陽の位置を \((0,0)\)とすると、　

　　　\(\overrightarrow{F}=-\dfrac{Gm_em_s}{r^2}(x,y)/\sqrt{x^2+y^2}\) 引力は、成分ごとに計算すればよく、
　　　\(F_x=-\dfrac{Gm_em_s}{(x^2+y^2)^{\frac{3}{2}}}x~~,~~~F_y=-\dfrac{Gm_em_s}{(x^2+y^2)^{\frac{3}{2}}}y\)
\(F_x=m_ea_x,~~F_y=m_ea_y\)　なので、
　　　　\(a_x=-\dfrac{Gm_s}{(x^2+y^2)^{\frac{3}{2}}}x~~,~~~a_y=-\dfrac{Gm_s}{(x^2+y^2)^{\frac{3}{2}}}y\)
地球の軌道は、２階微分方程式
　　　\(\dfrac{d^2x}{dt^2}=-\dfrac{Gm_s}{(x^2+y^2)^{\frac{3}{2}}}x(t)~~,~~~\dfrac{d^2y}{dt^2}=-\dfrac{Gm_s}{(x^2+y^2)^{\frac{3}{2}}}y(t)\)　　で表現されます。

数値計算法を用いて、地球の軌道を算出しましょう。
太陽と地球の距離を、149600000000メートルとします。

天文学では、\(1.496\times 10^{11}\)メートル、コンピュータでの表記は 1.496E+11 となります。

G=6.67259*10^{-11};me=5.974*10^{24};ms=1.989*10^{30}
x=1.496*10^{11};y=0;vx=0;vy=2.3*10^{4}
dt=86400
xlist=c();ylist=c()
ax <- function(x, y) {
  ax <- -(G*ms*x)/(x^2+y^2)^{3/2}            # x方向の加速度を求める
  return(ax)
}
ay <- function(x, y) {
  ay <- -(G*ms*y)/(x^2+y^2)^{3/2}              # y方向の加速度を求める
  return(ay)
}

for (i in 1:200){
  xlist = c(xlist, x)
  ylist = c(ylist, y)
  vx = vx+ax(x,y)*dt
  vy = vy+ay(x,y)*dt
  x = x+vx*dt
  y = y+vy*dt
}  
df = data.frame(x=xlist,y=ylist)
head(df)
x y
149600000000 0
149555731631 1987200000
149467180410 3973811598
149334329088 5959245411
149157157655 7942910612
148935643353 9924214469 
plot(df$x,df$y,col="blue",asp=1)

G=0.001;me=1;ms=10000
x=10;y=0;vx=0;vy=1
xlist=c();ylist=c()
ax <- function(x, y) {
  ax <- -(G*ms*x)/(x^2+y^2)^{3/2}            # x方向の加速度を求める
  return(ax)
}
ay <- function(x, y) {
  ay <- -(G*ms*y)/(x^2+y^2)^{3/2}              # y方向の加速度を求める
  return(ay)
}
for (i in 1:60){
  xlist = c(xlist, x)
  ylist = c(ylist, y)
  vx = vx+ax(x,y)
  vy = vy+ay(x,y)
  x = x+vx
  y = y+vy
}  
df = data.frame(x=xlist,y=ylist)
head(df)
x y
10.000000 0.000000
9.900000 1.000000
9.699511 1.989850
9.399105 2.959201
9.000465 3.897625
8.506435 4.794741 
plot(df$x,df$y,asp=1)