logo

줄리아 미분방정식 패키지 DiffetentialEquations 튜토리얼 📂Julia

줄리아 미분방정식 패키지 DiffetentialEquations 튜토리얼

説明

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は、その名の通り常微分方程式問題を定義する関数だ。基本적으로4つの引数を受け取ることができるので、一つずつ見ていこう。

  • 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 ↩︎ ↩︎