📂Julia

外力を持つ常微分方程式の数値解法 (ジュリア・微分方程式パッケージ - DifferentialEquations)

説明

DifferentialEquations.jlは微分方程式の数値的解法のためのジュリアパッケージだ。この文章では、DifferentialEquations.jlを利用して外力^{forcing term}のある常微分方程式の数値的解法を紹介する。チュートリアルを読んでくると役に立つだろう。

外力がある1階常微分方程式

以下のような簡単な問題を考えよう。

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

これを積分因子法で解くと次のような正確な解を得ることができる。

$$ u(t) = -\dfrac{1}{2} \left( \sin(t) + \cos(t) - e^{-t} \right) $$

さて、パッケージを使用して上記の方程式の数値解を求めてみよう。まず方程式を$\dfrac{du}{dt} = f(u, p, t)$の形に直す。すると、$f(u, p, t) = u(t) + \sin(t)$であり、コードでは次のように定義する。ここで注意する点は、ソリューション$u(t)$はuのように書く必要があり、パラメータ$p(t)$はp(t)のように書かなければならないということだ。

using DifferentialEquations
using Plots

f(u, p, t) = u + p(t)

性能のためにin-place方式で定義するには次のように書く。

function f!(du, u, p, t)
    du[1] = u[1] + p(t)
end

次に初期値、ドメイン、外力関数を定義しよう。fをin-place方式で定義した場合、初期値もベクトルでu0 = [0.0]のように定義する必要があることに注意しよう。

u0 = 0.0 
tspan = (0.0, 1.0)
p(t) = sin(t)

さて、ODEProblemに代入して問題を定義し、solveで問題を解けば終わりだ。

prob = ODEProblem(f, u0, tspan, p)
sol = solve(prob)

グラフを描いて確認すると次のようになる。

plot(sol, lw = 5, label = "numerical solution")
plot!(LinRange(0, 1, 100), t -> -(1/2)*(sin(t) + cos(t) - exp(t)), lw=3, ls=:dash, lc=:red, label="exact solution")

コード全文

using DifferentialEquations
using Plots

f(u, p, t) = u + p(t)

function g!(du, u, p, t)
    du[1] = u[1] + p(t)
end

u0 = 0.0
tspan = (0.0, 1.0)
p(t) = sin(t)

prob = ODEProblem(f, u0, tspan, p)
sol = solve(prob)

plot(sol, lw = 5, label = "numerical solution", dpi=300)
plot!(LinRange(0, 1, 100), t -> -(1/2)*(sin(t) + cos(t) - exp(t)), lw=3, ls=:dash, lc=:red, label="exact solution")

外力がある2階常微分方程式

次のような強制調和振動子問題を考えよう。

$$ \begin{align*} \ddot{x} + 2\dot{x} + 0.5 x &= \cos(2t), \qquad t\in [0, 40] \\ x(0) &= 0 \\ \dot{x}(0) &= 1 \end{align*} $$

これは2階線形非同次方程式なので、補助解$\phi(t)$と特殊解$X(t)$の合計で解を求めることができる。計算すると、補助解は$\phi(t) = 0.656 e^{(-1+\sqrt{0.5})t}-0.532e^{(-1-\sqrt{0.5})t}$であり、特殊解は$X(t) = -\frac{14}{113}\cos(2t) + \frac{16}{113}\sin(2t)$だ。したがって正確な解は次の通り。

$$ \begin{align*} x(t) &= \phi(t) + X(t) \\ &= 0.656 e^{(-1+\sqrt{0.5})t}-0.532e^{(-1-\sqrt{0.5})t} -\frac{14}{113}\cos(2t) + \frac{16}{113}\sin(2t) \end{align*} $$

では、パッケージを使用してこの問題を解いてみよう。ここで$x_{1} = x$、$x_{2} = \dot{x}$とすると、上記の問題を次のような連立常微分方程式で表すことができる。

$$ \begin{align*} \dot{x}_{1} &= \dot{x} = x_{2} \\ \dot{x}_{2} &= \ddot{x} = - 2\dot{x} - 0.5x + \cos(2t) = - 2x_{2} - 0.5x_{1} + \cos(2t) \end{align*} $$

$$ \implies \begin{align*} \dot{x}_{1} &= x_{2} \\ \dot{x}_{2} &= - 2x_{2} - 0.5x_{1} + \cos(2t) \end{align*} $$

上記の数式は次のようにコードで書かれる。

function forced_harmonic_oscillation!(dx, x, p, t)
    dx[1] = x[2]
    dx[2] = -2x[2] -0.5x[1] + p(t)
end

もちろんp(t)を単にcos(2t)と仮定して解いても問題はない。その場合、ODEProblemの引数にpを入れなければよい。前の例のように初期値、ドメイン、外力関数を定義し、以下のコードで数値解を求めることができる。

x0 = [0.0, 1.0]
tspan = (0.0, 40.0)
p(t) = cos(2t)

prob2 = ODEProblem(forced_harmonic_oscillation!, x0, tspan, p)
sol2 = solve(prob2)

実際の解と一致するかを確認するために描いてみると次のようになる。

x(t) = 0.6564172054223187exp((-1+sqrt(0.5))*t) -0.5325234001125843exp((-1-sqrt(0.5))*t)-(14/113)*cos(2t) + (16/113)*sin(2t)

plot(sol2, lw = 5, label = "numerical solution", dpi=300, idxs=1)
plot!(LinRange(0, 40, 100), t -> x(t), lw=3, ls=:dash, lc=:red, label="exact solution")

コード全文

using DifferentialEquations
using Plots

function forced_harmonic_oscillation!(dx, x, p, t)
    dx[1] = x[2]
    dx[2] = -2x[2] -0.5x[1] + cos(2t)
end

x0 = [0.0, 1.0]
tspan = (0.0, 40.0)
p(t) = cos(2t)

prob2 = ODEProblem(forced_harmonic_oscillation!, x0, tspan)
sol2 = solve(prob2)

x(t) = 0.6564172054223187exp((-1+sqrt(0.5))*t) -0.5325234001125843exp((-1-sqrt(0.5))*t)-(14/113)*cos(2t) + (16/113)*sin(2t)

plot(sol2, lw = 5, label = "numerical solution", dpi=300, idxs=1)
plot!(LinRange(0, 40, 100), t -> x(t), lw=3, ls=:dash, lc=:red, label="exact solution")

環境

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