微分方程式の数値解析 #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%ほど長くなるとみなせる。