B..3.3 単振動の方程式

shm1.jl

# shm1.jl --- simple harmonic motion
#

using DifferentialEquations
using Plots
using Printf

function shm!(y,x,p,t)
  omega2,=p
  y[1] = x[2]
  y[2] = - omega2 * x[1]
end

omega = 1; omega2 = omega * omega; p=(omega2,)
t0=0; T=20; tspan=(t0, T)
x0=1; v0=0; X0=[x0,v0]
prob = ODEProblem(shm!,X0,tspan,p)
sol = solve(prob, reltol=1e-8, abstol=1e-8);

plot(sol, xlim=(t0, T), ylim=(-1.2, 1.2),
     title="simple harmonic motion: dx/dt=y, dy/dt=-ω^2 x, x(0)=1, y(0)=0",
     xlabel="t", ylabel="x and y", size=(800,500))

savefig("shm1_julia.png")
savefig("shm1_julia.pdf")

図 18: 単振動の方程式: $x(t)$, $y(t)$ のグラフを描く

shm1.jl

# shm2.jl --- simple harmonic motion
#

using DifferentialEquations
using Plots
using Printf

function shm!(y,x,p,t)
  omega2,=p
  y[1] = x[2]
  y[2] = - omega2 * x[1]
end

function draworbits(omega)
  omega2 = omega * omega; p=(omega2,)
  t0=0; T=2*pi/omega; tspan=(t0, T)

  for i=0:30
    x0=0.1*i; v0=0; X0=[x0,v0]
    prob = ODEProblem(shm!,X0,tspan,p)
    sol = solve(prob, reltol=1e-4);
    if i==0
      plot(sol, vars=(1,2), xlim=(-2, 2), ylim=(-2, 2), aspect_ratio=1,
           title="simple harmonic motion: dx/dt=y, dy/dt=-ω^2 x",
           xlabel="t", ylabel="x and y", legend=:none)
    else
      global fig=plot!(sol, vars=(1,2), xlim=(-2,2), ylim=(-2,2),legend=:none)
    end
  end
  return fig
end

if abspath(PROGRAM_FILE) == @__FILE__
  print("こうやって実行すると時間がかかります。")
  draworbits(1.0)
  savefig("shm2_julia.png")
  savefig("shm2_julia.pdf")
else
  println("draworbits(omega)")
  println("savefig(\"なんとか.png\")")
end

図 19: 単振動の方程式: $x(t)$, $y(t)$ のグラフを描く