julia(1.4.1)~二階微分方程式（バネマスダンパ系）の数値積分~

概要

バネマスダンパ系の運動を数値計算により解く

環境

• julia v(1.4.1)
  • "DifferentialEquations" => v"6.13.0"
  • "Plots" => v"0.29.9"

問題設定

質量m，ばね定数k，粘性cのバネマスダンパ系を考える．

運動方程式は，となる．

この式を，以下のように式変形し，数値積分できる形にする．

とおくと，以下のように表せる． なお，はの時間微分である．

数値積分の実装

使用するモジュールを指定する． 今回，結果のプロットにはPlots.jlを使用する．これは最後にアニメーションで記載する都合．

using DifferentialEquations
using Plots

次に，問題を定義する．

#問題の定義
function SMD_ODE(du,u,p,t)
    m = p[1];
    k = p[2];
    c = p[3];
    du[1] = -(k/m)*u[2] - (c/m)*u[1]
    du[2] = u[1]
end

ここで，pは外部から与えるパラメータ．わかりやすいように一度別の変数に展開している．

次に，シミュレーションを行うパラメータを指定する．

#シミュレーションパラメータ
m = 1.0
k = 2.0
c = 0.5
#初期値の定義
u0 = [0.0,1.0] #[dx0,x0]
param = [m,k,c]
#積分区間の指定
tspan = (0.0,20)
#積分問題の定義
prob = ODEProblem(SMD_ODE,u0,tspan,param)

ここでODEproblemに与える引数の順番は，関数SMD_ODEと異なることに注意．

そして数値積分を実行する．今回はsolverにVern7を使用する．

sol = solve(prob,Vern7(),abstol = 1e-12,reltol = 1e-12,saveat = 0.1);

これは以前の記事を参照して決めた．

Julia(1.4.1)~常微分方程式の計算~ - Tagebuch von Spargel

julia(1.4.1)~数値積分の精度確認~ - Tagebuch von Spargel

結果をプロットする．

x = [u[2] for u in sol.u]
plot(sol.t,x)
ylabel!("Position")
xlabel!("Time [s]")
ylims!((-1,1))

結果を取り出すのにリスト内包表記を使用している． 単純な参照をやる方法がわからなかったため，おとなしくx = [u[2] for u in sol.u]とするのがよいと思う．

位置の変化は以下のようになる． きちんと計算されていることが確認できる．

アニメーション（おまけ）

最後にアニメーションで視覚化してみる．ばねは適当に書いたものでバネっぽい以外の意味はない．

今記事ではアニメーションについて詳細に取り扱わないが，これをじっこうするためにPlots.jlを使用する．

using Plots
anim = Animation()
function spring_pos(x)
    sx = range(-2,x-0.1,length=14)
    sy = [0,0.1,-0.1,0.1,-0.1,0.1,-0.1,0.1,-0.1,0.1,-0.1,0.1,-0.1,0]
    return sx,sy
end
for i in 1:100
    spring_position = spring_pos(sol.u[i][2])
    plt = plot(spring_position[1],spring_position[2],color = "black")
    plot!([sol.u[i][2]],[0],marker = :circ,markersize=20)
    ylims!((-1,1))
    xlims!((-2,1.5))
    plot!(legend=false)
    frame(anim,plt)
end

gif(anim,"test.gif",fps = 10)

アニメーションを使わないならPyPlot.jlのほうが使いやすいような・・・．