スリンキーの自由落下（４）

スリンキーの自由落下（３）からの続き

これまで，Mathematicaでの計算を示してきた。黒木さんのようにJuliaをスイスイ使えるようになりたいので，微分方程式をJuliaで数値的に解いてみよう。これがまた，なかかな適当な初心者向けのテキストがないのだ。まあ，オリジナルのチュートリアルとかそのビデオ版をみて勉強すればよいのだが，なまけものには日本語の簡単な入門書が欠かせない。

とりあえず，DifferentialEquationsパッケージの連立常微分方程式を使うための，最小限の手続きを見様見まねであてはめてみると次のようになった。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
using DifferentialEquations
using ParameterizedFunctions
using Plots; gr()

sky = @ode_def slinky begin
  dx1 = v1
  dx2 = v2
  dv1 = -g + k/m1*(x2-x1+L)
  dv2 = -g - k/m2*(x2-x1+L)
  end g k m1 m2 L
L=1
m2=1
k=5
g=10
u0 = [0,-(L+m2*g/k),0,0]
p = (10.0,5.0,1.0,1.0,1.0)
tspan = (0.0,1.0)
prob = ODEProblem(sky,u0,tspan,p)
sol = solve(prob)
plot(sol,vars=(0,1))
plot!(sol,vars=(0,2))
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -