Sundialsを使ったJuliaでの微分方程式の数値計算 ref: http://qiita.com/so1_tsuda/items/12c945db5759d8229ec4

file0.txt

;; ESS - Emacs Speaks Statistiscs (for R and Julia)
(add-to-list 'load-path "ESSへのパス")
(load "ess-site")
;; julia 
(setq inferior-julia-program-name "juliaバイナリへのパス")

file1.txt

Pkg.update()
Pkg.add("DifferentialEquations") # Sundials.jlはここで一緒にインストールされる
Pkg.add("Plots")
Pkg.add("DataFrames") # optional
# ちゃんと読み込めるか確認（＆初回プリコンパイル）
using DifferentialEquations
using Sundials
using Plots

file10.txt

# Sundials.jl
  0.000547 seconds (1.25 k allocations: 33.453 KB)
  0.000481 seconds (1.25 k allocations: 33.453 KB)
  0.000477 seconds (1.25 k allocations: 33.453 KB)
  0.000584 seconds (1.25 k allocations: 33.453 KB)
# DifferentialEquations.jl
  0.001087 seconds (1.55 k allocations: 45.375 KB)
  0.000908 seconds (1.55 k allocations: 45.375 KB)
  0.000871 seconds (1.55 k allocations: 45.375 KB)
  0.001277 seconds (1.55 k allocations: 45.375 KB)
# ODE.jl
  0.001204 seconds (3.37 k allocations: 109.469 KB)
  0.001843 seconds (3.37 k allocations: 109.469 KB)
  0.001066 seconds (3.37 k allocations: 109.469 KB)
  0.001506 seconds (3.37 k allocations: 109.469 KB)

file11.txt

using ParameterizedFunctions
f = @ode_def VdP begin
    du[1] = u[2]
    du[2] = mu * (1.0 - u[1]^2.0) * u[2] - u[1]
end mu=>0.2
# デフォルトではmu=0.2
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, CVODE_BDF(), dt=0.01)
x = [u[1] for u in sol.u]
y = [u[2] for u in sol.u]
plot(x,y)
# パラメータmuを変更
f.mu = 5.0
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, CVODE_BDF(), dt=0.01)
x = [u[1] for u in sol.u]
y = [u[2] for u in sol.u]
plot!(x,y)
# mu=1.0のインスタンス g を生成
g = VdP(mu=1.0)
prob = ODEProblem(g,u0,tspan)
sol = solve(prob, CVODE_BDF(), dt=0.01) 
x = [u[1] for u in sol.u]
y = [u[2] for u in sol.u]
plot!(x,y)
savefig("vdp_parameterized.png")

file2.txt

Pkg.update()
Pkg.add("Sundials")

file3.txt

LoadError: failed process: Process(`powershell -Command '(new-object net.webclient).DownloadFile("https://cache.e.ip.saba.us/https://bintray.com/artifact/download/tkelman/generic/sundials-2.5.0.7z", "C:\cygwin64\home\user\.julia\v0.4\Sundials\deps\downloads\sundials-2.5.0.7z")'`, ProcessExited(1)) [1]
    while loading C:\cygwin64\home\user\.julia\v0.4\Sundials\deps\build.jl, in expression starting on line 37

file4.txt

Pkg.build("Sundials")
Pkg.test("Sundials")

file5.txt

\frac{dx}{dt} = y \\

\frac{dy}{dt} = -\mu (x^2 - 1)y - x

file9.txt

sol = solve(prob, CVODE_BDF(), dt=0.01);
for i in 0:3 
    @time sol = solve(prob, CVODE_BDF(), dt=0.01);
end
res = Sundials.cvode(vdp_sun, u0, t);
for i in 0:3 
    @time res = Sundials.cvode(vdp_sun, u0, t);
end
(t, pos) = ode45(vdp_ODE, u0, t);
for i in 0:3 
    @time (t, pos) = ode45(vdp_ODE, u0, t);
end

vdp_diffeq.jl

using DifferentialEquations
using Plots 
# 非線形パラメータ
mu = 0.2
# 関数定義
function vdp(t, u, du)
    # dx/dt
    du[1] = u[2]
    # dy/dt
    du[2] = mu * (1.0 - u[1]^2.0) * u[2] - u[1]
end
# パラメータ設定
u0 = [0.2; 0.2]
tspan = (0.0,10.0)
# メインの計算
prob = ODEProblem(vdp,u0,tspan)
sol = solve(prob, CVODE_BDF(), dt=0.01) # Sundialsのアルゴリズムを使う
# プロット
x = [u[1] for u in sol.u]
y = [u[2] for u in sol.u]
plot(x,y)
savefig("vdp_diffeq.png")

vdp_ode.jl

using ODE
mu = 0.2
function vdp_ODE(t, u)
    # dx/dt
    dx_dt = u[2]
    # dy/dt
    dy_dt = mu * (1.0 - u[1]^2.0) * u[2] - u[1]
    # return du
    [dx_dt; dy_dt]
end
u0 = [0.2; 0.2]
t = 0:0.01:10.0
(t, pos) = ode45(vdp_ODE, u0, t)
x = [u[1] for u in pos]
y = [u[2] for u in pos]
plot(x,y)
savefig("vdp_ode.png")

vdp_sudials.jl

using Sundials
using Plots 
# 非線形パラメータ
mu = 0.2
# 関数定義
function vdp_sun(t, u, du)
    # dx/dt
    du[1] = u[2]
    # dy/dt
    du[2] = mu * (1.0 - u[1]^2.0) * u[2] - u[1]
    du
end
# パラメータ設定
u0 = [0.2; 0.2]
t = [0.0:0.1:10.0;]
# メインの計算
res = Sundials.cvode(vdp_sun, u0, t); # u0とtの順序に注意（関数定義と逆）
# プロット
x = [a for a in res[:,1]] 
y = [a for a in res[:,2]] 
plot(x,y)
savefig("vdp_sundials.png")