Rubyで数値計算2

運動方程式なんぞ、解いてみますか。

require 'matrix'

class RungeKutta
  def initialize(dy)
    @dy = dy
  end
  
  def step(x, y, h)
    k1 = h * @dy.call(x, y)
    k2 = h * @dy.call(x + h / 2, y + k1 / 2)
    k3 = h * @dy.call(x + h / 2, y + k2 / 2)
    k4 = h * @dy.call(x + h, y + k3)
    y + k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6
  end
end

# force
f = proc{|t|
  if t < 0.2 or (t > 5.0 and t <= 5.1)
  	-1.0
  else
  	0.0
  end
}

# y = [[v], [x]]
# d^2x/dt^2 + dx/dt + 3x + f(t) = 0
# -> dv/dt = d^2x/dt^2 = - v - 3x - f(t) 
#    dx/dt = v
dy = proc{|t, y|
  Matrix.rows([[- y[0,0] - 3 * y[1,0] - f.call(t)], [y[0,0]]])
}

t = 0.0
y = Matrix.rows([[0.0], [0.0]])
dt = 0.1

rk = RungeKutta.new(dy)

for i in 1..100
  y = rk.step(t, y, dt)
  t += dt
  printf("%20.10f, %20.10f, %20.10f", t, y[0,0], y[1,0])
  puts
end

Motion eq.