斜方投射の数値計算と描画(Ruby)

4次のルンゲ−クッタ法の練習としてやってみました。斜方投射(斜め上にものを投げるということです)じゃ全然威力を発揮しないのですが、とりあえず動くかどうかだけ確かめました。描画は gnuplot でやっています。
 

 
runge_kutta_sample1.rb

require 'numo/gnuplot'
require 'matrix'

def runge_kutta(f, dt, x, v, t, m)    #x, v は Vector, f[x, v, t] は Vector を返す
  k1 = v * dt
  l1 = f[x, v, t] * dt / m
  k2 = (v + l1 / 2) * dt
  l2 = f[x + k1 / 2, v + l1 / 2, t + dt / 2] * dt / m
  k3 = (v + l2 / 2) * dt
  l3 = f[x + k2 / 2, v + l2 / 2, t + dt / 2] * dt / m
  k4 = (v + l3) * dt
  l4 = f[x + k3, v + l3, t + dt] * dt / m
  
  nx = x + (k1 + 2 * k2 + 2 * k3 + k4) / 6
  nv = v + (l1 + 2 * l2 + 2 * l3 + l4) / 6
  
  [nx, nv]
end

g = 9.80665
m = 1.0
x = Vector[0.0, 0.0]
v = Vector[1.0, 1.0]
f = lambda {|x, v, t| Vector[0, -m * g]}
dt = 0.001

t = 0
ax, ay = [], []
while x[1] >= 0
  nx, nv = runge_kutta(f, dt, x, v, t, m)
  ax << x[0]
  ay << x[1]
  x, v = nx, nv
  t += dt
end

Numo.gnuplot do
  unset :key
  plot ax, ay, w: :l
end
gets    #終了待ち

初期値はいい加減に設定したのですが、グラフを見ると 5cm くらいの高さまでしか上がらないですね(笑)。初速はだいたい 1.5m/s くらいなのですが、これは時速 5km くらい(おおよそ人の歩く速さ)で、これが小さすぎるのですね。高校生で速い球を投げる人だと時速 100km 以上出るので、それくらいにしてみると(斜めに投げ上げても)高さ 20m くらいまで上がります。また、80m くらい先まで飛びます。なお、物理を知っている人には自明ですが、これらは球の質量には関係しません。重い球だろうが軽い球だろうが関係ないのですね。不思議ですか? このプログラムで実際に質量 m の値を変えてやれば、はっきりとわかります。
 
空気抵抗を入れてみます。はっきりとグラフに出るよう、非常に強く効かせています。力を次の関数で与えます。

f = lambda {|x, v, t| Vector[0, -m * g] - v * 10}

グラフはこんな風になります。なかなか飛ばない感じが出ているでしょうか。

 
RubyVector クラスはじつに手軽ですね。こういうことに使うにはすばらしいと思います。以上は 2次元でやっていますが、このまま 1次元でも 3次元でも動きます。
 

※参考
運動方程式の数値計算法(pdf)