微分方程式の数値解析 #12

滑らかなものばかりだったので滑らかさのないものを扱ってみる。
両側を壁で囲まれた一次元の質点の運動について、
壁で瞬間的に衝突、反射するとすれば運動方程式からの数値解析が難しくなるので、
壁がばねで支えられているとして、その復元力をもって有限時間で完全弾性衝突することにする。
ばねの自然長を中心とし両側の壁の間が不感帯であるようなばね-質点運動と考えればよい。
位置±5にばね仕込みの壁を配置し、質点は時刻0で位置0、速度1の初期値を持つとする。
壁のばね定数と質点質量の比として10, 10000の二種類について見てみる。

program ode
  use runge_kutta_4
  implicit none
  double precision, parameter :: t0 = 0           ! initial time
  double precision, parameter :: x0(2) = (/0, 1/) ! initial value at t0
  double precision :: t
  double precision :: x(size(x0))
  double precision, parameter :: udiv = 1         ! division unit
  integer, parameter :: ndiv = 1000               ! number of divisions per udiv
  double precision, parameter :: h = udiv / ndiv  ! delta-t
  integer, parameter :: niter = 50 / udiv * ndiv  ! number of iterations
  integer, parameter :: pdiv = ndiv / 100         ! print interval
  integer :: n
  type (func_type), dimension(size(x0)) :: fns
  fns(1)%f => f1
  fns(2)%f => f2
  call set_func(fns)
  x = x0
  do n = 0, niter
    t = t0 + n * h
    if (mod(n, pdiv) == 0) print '(4F13.8)', t, x(1), x(2), f2(t, x) ! time, location, velocity, force/mass
    x = next(t, x, h)
  end do
contains
  double precision function f1(t, x)
    double precision, intent(in) :: t
    double precision, dimension(:), intent(in) :: x
    f1 = x(2)
  end function f1
  double precision function f2(t, x)
    double precision, intent(in) :: t
    double precision, dimension(:), intent(in) :: x
    double precision, parameter :: kpm = 1d4 ! spring constant / mass
    if (x(1) .gt. 5) then
      f2 = -kpm * (x(1) - 5)
    else if (x(1) .lt. -5) then
      f2 = -kpm * (x(1) + 5)
    else
      f2 = 0
    end if
  end function f2
end program ode

質点の位相空間での軌跡を示すと、

kpm = 1d4で図的には瞬間的に壁で反射しているようにみえる。
質点に働く力の時間変化を示すと、


単位質量当たりの力が0ラインとの間に作る面積が壁との衝突による単位質量当たりの力積ということになる。
kpm = 1d4ではかなり短時間だけ力が働いており、まあ撃力と呼んでよいだろう。
最後にkpm = 1d4での質点位置の時間変化を示す。

壁間距離10を速さ1で運動すれば周期20になるはずだが、
時刻50で位置0に戻り切れていないことが分かる。
周期は約20.07であり、ばね定数-質量比10000で実効壁間距離が0.35%ほど長くなるとみなせる。