logo

Juliaの微分方程式パッケージDifferentialEquationsチュートリアル 📂Julia

Juliaの微分方程式パッケージDifferentialEquationsチュートリアル

説明

DifferentialEquations.jlはSciMLグループに属するパッケージの一つで、微分方程式の数値的解法のために開発された。このパッケージで解ける方程式は次の通りである。

  • Discrete equations (function maps, discrete stochastic (Gillespie/Markov) simulations)
  • 常微分方程式(ODE)
  • Split and Partitioned ODEs (Symplectic integrators, IMEX Methods)
  • 確率微分方程式(SDE)
  • Stochastic differential-algebraic equations (SDAEs)
  • Random differential equations (RODEs or RDEs)
  • Differential algebraic equations (DAEs)
  • Delay differential equations (DDEs)
  • Neutral, retarded, and algebraic delay differential equations (NDDEs, RDDEs, and DDAEs)
  • Stochastic delay differential equations (SDDEs)
  • Experimental support for stochastic neutral, retarded, and algebraic delay differential equations (SNDDEs, SRDDEs, and SDDAEs)
  • Mixed discrete and continuous equations (Hybrid Equations, Jump Diffusions)
  • (Stochastic) partial differential equations ((S)PDEs) (with both finite difference and finite element methods)

チュートリアルを超えた使い方は以下を参考にしよう。

インストール

JuliaのREPLで次のように入力してインストールする。

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

あるいは]を押してパッケージモードに入り、次のように入力してインストールする。

(@v1.10) pkg> add DifferentialEquations

ODEProblem

DifferentialEquations.jlにおける常微分方程式の解法は次のように簡単にまとめられる。

problem = ODEProblem(f, u0, tspan, p)
solution = solve(prob)

この過程は以下で具体的に扱うことにして、まずはODEProblemについて詳しく調べてみよう。関数ODEProblemは名前の通り常微分方程式問題を定義する関数である。基本的に四つの引数を入力できるが、一つずつ見ていこう。

  • f: SciMLBase.ODEFunction型の関数である。常微分方程式を以下のように整理したとき、$f(u, p, t)$に該当する部分を定義すればよい。 $$ \dfrac{du}{dt} = f(u, p, t) $$ このときメモリ効率のためにin-place方式で定義できるが、この場合はf!(du, u, p, t)のように定義すればよい。
  • u0: 初期条件である。$u$がスカラー関数ならu0 = 1.0のように、ベクトル関数ならu0 = [1.0, 1.0, 1.0]のように置けばよい。
  • tspan: 問題で定義された$u$のドメインをタプルで入力する。tspan = (t_start, t_end)
  • p: 解$u$を除いたパラメータ(外力など)を入力する。

では以下のような常微分方程式が与えられたとしよう。

$$ \begin{align*} \dfrac{du(t)}{dt} &= u(t) + \sin(t) \qquad t \in [0, 1] \\ u(0) &= 0.0 \end{align*} $$

すると次のようなコードで問題を定義する。

using DifferentialEquations

f(u, p, t) = u + p(t)
u0 = 0.0
tspan = (0.0, 1.0)
p(t) = sin(t)
prob = ODEProblem(f, u0, tspan, p)

probをsolve関数の引数に入れると方程式の解を得られる。では以下の例を通じてより具体的に調べてみよう。

常微分方程式の解法1

最も簡単な常微分方程式の例として次のような指数成長方程式を考えよう。

$$ \dfrac{du}{dt} = \alpha u $$

この方程式の解はある定数$A$に対して$A e^{\alpha t}$と等しいことをすでに知っているが、例題のために最も簡単な形の微分方程式を持ってきた。具体的には次のような初期値問題を考えよう。

$$ \begin{align*} \dfrac{du(t)}{dt} &= 4 u(t) \qquad t \in [0, 1] \\ u(0) &= 0.5 \end{align*} $$

まず右辺を関数$f(u, p, t)$として定義する。ここで$p$はパラメータparameterであり、普通は各項の係数、この方程式では$\alpha = 2$を意味する。そして初期値と時間の範囲を設定する。

より正確に説明すると$u$を更新する規則を定めることである。$du$がどのように表現されるかを関数として定義し、これをODEProblemの引数に入れるのである。$dt$は入力しなければ自動的に設定される。すなわち微分方程式を$du = f(u,p,t)dt$の形に直したのち$f(u,p,t)$部分を定義してやればよいということである。このようにパラメータが定数の場合には$p$を明示的に関数内部で使う必要はない。

using DifferentialEquations

f(u, p, t) =  4u
u0 = 0.5
tspan = (0.0, 1.0)

これらをODEProblemに引数として入れて、解くべき問題problemを作ることができる。

julia> problem = ODEProblem(f, u0, tspan)
ODEProblem with uType Float64 and tType Float64. In-place: false
timespan: (0.0, 1.0)
u0: 0.5

ではproblemをsolveに代入すると微分方程式の数値解を得られる。ソルバが解を7次元ベクトルとして求めてくれた。

julia> @time solution = solve(problem)
  0.000074 seconds (192 allocations: 28.750 KiB)
retcode: Success
Interpolation: 3rd order Hermite
t: 7-element Vector{Float64}:
 0.0
 0.07581611843444613
 0.21771610504885602
 0.39336627353058856
 0.6101928843972066
 0.8638706427431251
 1.0
u: 7-element Vector{Float64}:
 0.5
 0.581866090044407
 0.7728154714821961
 1.0981042913641448
 1.6942469046842847
 2.8139612733277692
 3.6945246786578037

proertynamesを使うとsolutionの属性を確認できる。

julia> propertynames(solution)
(:u, :u_analytic, :errors, :t, :k, :discretes, :prob, :alg, :interp, :dense, :tslocation, :stats, :alg_choice, :retcode, :resid, :original, :saved_subsystem)

解$u$とドメイン$t$をそれぞれ得たい場合は次のようにできる。

julia> solution.u
7-element Vector{Float64}:
 0.5
 0.581866090044407
 0.7728154714821961
 1.0981042913641448
 1.6942469046842847
 2.8139612733277692
 3.6945246786578037

julia> solution.t
7-element Vector{Float64}:
 0.0
 0.07581611843444613
 0.21771610504885602
 0.39336627353058856
 0.6101928843972066
 0.8638706427431251
 1.0

可視化

solutionに対するレシピがすでに定義されているため、単にplotに代入するとグラフを描ける。よく見るとsolution.uとsolution.tを直接代入して描いたときと違って自動的に値を補間interpolationしてくれることがわかる。またx軸が$t$であることも自動的に表示される。

using Plots
p1 = plot(solution)
p2 = plot(solution.t, solution.u)
plot(p1, p2, layout = (2, 1))

実際の解である$u(t) = 0.5 e^{4t}$と比較してみよう。

plot(solution, label="solver", lw=4)
plot!(LinRange(0,1,100), t->0.5*exp(4t), label="exact", lw=3, ls=:dash, lc=:red, dpi=300)

詳細設定

キーワード引数を直接設定してソルバを選択したり、許容誤差を設定したりできる。たとえば$0.05$タイムステップごとに解を計算したいならsaveat = 0.05save atに設定すればよい。詳しくはリンクで確認できる。

# 솔버를 Tsit5()로 설정
# 매 타임스텝 0.05마다 솔루션을 저장.
# 허용 절대오차는 1e-8로 설정.

julia> solve(problem, Tsit5(), saveat=0.05, abstol=1e-8)
retcode: Success
Interpolation: 1st order linear
t: 21-element Vector{Float64}:
 0.0
 0.05
 0.1
 0.15
 0.2
 0.25
 0.3
 0.35
 0.4
 0.45
 0.5
 0.55
 0.6
 0.65
 0.7
 0.75
 0.8
 0.85
 0.9
 0.95
 1.0
u: 21-element Vector{Float64}:
  0.5
  0.6107013642059927
  0.7459122224392803
  0.9110596077813341
  1.1127694183818455
  1.3591415739400332
  1.6600546629913426
  2.0275986016325755
  2.4765152918979516
  3.0248023813153906
  3.6945244502370813
  4.512485296082109
  5.511506995093366
  6.731841509374868
  8.22227376044907
 10.042494301852722
 12.266093755920075
 14.981974632546809
 18.298554052093934
 22.34964366466662
 27.298573538274052

コード全文

using DifferentialEquations

f(u, p, t) =  4u
u0 = 0.5
tspan = (0.0, 1.0)

problem = ODEProblem(f, u0, tspan)
@time solution = solve(problem)

using Plots
p1 = plot(solution, title="plot(solution)")
p2 = plot(solution.t, solution.u, title="plot(solution.t, solution.u)")
plot(p1, p2, layout = (2, 1), dpi=300)

plot(solution, label="solver", lw=4)
plot!(LinRange(0,1,100), t->0.5*exp(4t), label="exact", lw=3, ls=:dash, lc=:red, dpi=300)

solve(problem, Tsit5(), saveat=0.05, abstol=1e-8)

連立常微分方程式の解法2

連立1階微分方程式の例として次のようなローレンツ方程式を考えよう。

$$ \begin{align*} \dfrac{dx}{dt} &= -\sigma x + \sigma y \\ \dfrac{dy}{dt} &= -xz + \rho x - y \\ \dfrac{dz}{dt} &= xy - \beta z \end{align*} $$

連立微分方程式を解くことも微分方程式を解くことと似ている。$du = (dx, dy, dz)$を更新する規則を定義し、これをODEProblemの引数に入れれば終わりである。パラメータを$\sigma = 10$、$\beta = 8/3$、$\rho = 28$に設定し、$du = f(u,p,t)dt$の形に直そう。

$$ \begin{align*} dx &= (-10 x + 10 y)dt \\ dy &= (-xz + 28 x - y)dt \\ dz &= (xy - \frac{8}{3} z)dt \end{align*} $$

では右辺の式を$dt$を除いて関数として定義する。このときlorenzではなくin-place関数であるlorenz!として定義するのは性能のためである2。この場合には関数の一番目の引数に出力であるduを追加してやればよい。$u = (x,y,z)$、$du = (dx,dy,dz)$と置けば、

$$ \begin{align*} dx &= 10(y-x)dt \\ dy &= \left[x(28-z)- y\right]dt \\ dz &= (xy - \frac{8}{3} z)dt \end{align*} \implies \begin{align*} du_{1} &= 10(u_{2}-u_{1}) \\ du_{2} &= u_{1}(28-u_{3})- u_{2} \\ du_{3} &= u_{1}u_{2} - \frac{8}{3} u_{3} \end{align*} $$

function lorenz!(du, u, p, t)
du[1] = 10(u[2] - u[1])
du[2] = u[1]*(28-u[3]) - u[2]
du[3] = u[1]*u[2] - (8/3)u[3]
end

初期値を$u(0) = (1,1,1)$、ドメインを$[0, 100]$に設定したのち解を計算する。解をプロットするときは引数idx=(1,2,3)を必ず追加してやらなければならない。

u0 = [1.0, 1.0, 1.0]
tspan = (0.0, 100)
prob = ODEProblem(lorenz!, u0, tspan)
sol = solve(prob)

plot(sol, idxs=(1,2,3))

環境

  • OS: Windows11
  • Version: Julia 1.10.0, DifferentialEquations v7.15.0, Plots v1.40.3

  1. https://docs.sciml.ai/DiffEqDocs/stable/getting_started/#Example-1-:-Solving-Scalar-Equations ↩︎

  2. https://docs.sciml.ai/DiffEqDocs/stable/getting_started/#Example-2:-Solving-Systems-of-Equations ↩︎ ↩︎