27 Julia で ODE

(こちらもしばらくは工事中。

DifferentialEquations.jl: Scientific Machine Learning (SciML) Enabled Simulation and Estimation

Julia だけでなく、Python, R でも使えるらしい。

Julia で使えるようにするには

using Pkg
Pkg.add("DifferentialEquations")
using DifferentialEquations

(結構長い時間がかかる。)

“Solving ODEs in Julia” by Chris Rackauckas

using DifferentialEquations
f(u,p,t) = 1.01*u
u0 = 1/2
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

using Plots
plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
     xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!") # legend=false
plot!(sol.t, t->u0*exp(1.01t),lw=3,ls=:dash,label="True Solution!")
savefig("diff.png")

こちらは , reltol=1e-8, abstol=1e-8 としてある。 それも込めて、プログラムの書き方を見てほっとする。 最初の例からこういうのを見せるべきだ。

ソルバーの指定が出来るそうで、, Tsit5() というのがそれに相当する (Tsitouras 5/4 Runge-Kutta method)。 どういうソルバーが良いかは、 「ODE Solvers」 に詳しく書いてある。 MATLAB だと ode45 というのがあるけれど、 翻訳は DP5() だとか (どういう意味だろう)。 それよりは大抵の場合 Tsit5() の方が効率的だとか。 Tsitouras というのはギリシャの人か。 バリバリ論文を出しているんだ (名前が珍しいので、検索して論文探すの簡単)。 知らなかった。

(2022/1/9追記)      少しツッコミを入れる。 1.01 のようなマジック・ナンバーを入れる (2箇所別々) のは拙い。 パラメーターを渡す仕掛けがあり、それを使うべきだと思う。

using DifferentialEquations
f(u,p,t) = p*u
u0 = 1/2
tspan = (0.0,1.0)
a=1.01
prob = ODEProblem(f,u0,tspan,a)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

using Plots
plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
     xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!") # legend=false
plot!(sol.t, t->u0*exp(a*t),lw=3,ls=:dash,label="Exact Solution!")
savefig("diff.png")

多次元の場合の例として、Loenz 方程式があげられている。

using DifferentialEquations
function lorenz!(du,u,p,t)
    σ,ρ,β = p
    du[1] = σ*(u[2]-u[1])
    du[2] = u[1]*(ρ-u[3]) - u[2]
    du[3] = u[1]*u[2] - β*u[3]
end
u0 = [1.0,0.0,0.0]
p = (10,28,8/3)
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan,p)
sol = solve(prob)

using Plots
plot(sol,vars=(1,2,3))

何と MATLABDiffEq.jl というのもあるらしい (MATLAB をインストールしてあることが必要だとか)。

(某月某日) あ、面白いの見つけた。

「数値計算法基礎 (2020)」