ここで訂正、減衰振動の速度の理論解は、筆者の計算ミスで正解ではありません。変位の理論解は正解です。

微分の計算を間違えました。

減衰振動の速度の理論解は、読者が挑戦して下さい。

読者の楽しみを残しました。

Excel macro  sub  srk-2()

       Sub srk-2()

       [H1].Select

       For i = 1 To 1000

       ActiveCell.Offset(i, 0).Characters.Text = "             "

       Next i

       tend = Range("B5")

       step = Range("B6")

       m = Range("B7")

       k = Range("B8")

       zeta = Range("B9")

       x0 = Range("B10")

       v0 = Range("B11")

       c = 2# * zeta * Sqr(m * k)

       omega = Sqr(k / m)

       root_1_zeta = Sqr(1# - zeta * zeta)

       [F5].Select

       dt = tend / step

       ActiveCell.Offset(2, 0) = dt

       ActiveCell.Offset(3, 0) = omega * root_1_zeta / 2# / 3.14159265358979

       [O1].Select

       t = 0#

       x = x0

       v = v0

       For i = 1 To step

       Call rk2(dt, t, x, v, xnew, vnew, m, k, c)

       ActiveCell.Offset(i, 0) = t

       ActiveCell.Offset(i, 1) = x

       ActiveCell.Offset(i, 2) = v

       ActiveCell.Offset(i, 3) = Exp(-zeta * omega * t) * (x0 * Cos(omega * t * root_1_zeta) + (omega * zeta * x0 + v0) / (omega * root_1_zeta) * Sin(omega * t * root_1_zeta))

       a1 = zeta * (1# - omega) * x0 + v0

       b1 = -x0 * (zeta * zeta / Sqr(1# - zeta * zeta) + omega * Sqr(1# - zeta * zera)) + v0 * zeta / Sqr(1# - zeta * zeta)

       ActiveCell.Offset(i, 4) = Exp(-zeta * omega * t) * (a1 * Cos(omega * Sqr(1# - zeta * zeta) * t) + (b1 * Sin(omega * Sqr(1# - zeta * zeta) * t)))

       t = t + dt

       x = xnew

       v = vnew

       Next i

       End Sub

      Sub rk2(dt, t, x, v, xnew, vnew, m, k, c)

      k1 = dt * f(t, x, v)

      l1 = dt * g(t, x, v, m, k, c)

      k2 = dt * f(t + dt / 2, x + k1 / 2, v + l1 / 2)

      l2 = dt * g(t + dt / 2, x + k1 / 2, v + l1 / 2, m, k, c)

      k3 = dt * f(t + dt / 2, x + k2 / 2, v + l2 / 2)

      l3 = dt * g(t + dt / 2, x + k2 / 2, v + l2 / 2, m, k, c)

      k4 = dt * f(t + dt, x + k3, v + l3)

      l4 = dt * g(t + dt, x + k3, v + l3, m, k, c)

      xnew = x + (k1 + 2 * (k2 + k3) + k4) / 6

      vnew = v + (l1 + 2 * (l2 + l3) + l4) / 6

      End Sub

      Function f(t, x, v)

      f = v

      End Function

      Function g(t, x, v, m, k, c)

      g = k / m * (-x) + c / m * (-v)

      End Function