やったこと

JuliaでRunge-Kutta法実装した
2階定数係数線形微分方程式をRunge-Kutta法で数値的に解いて解析解と比べた

参考にした文献

『中古』数値計算法 (新コンピュータサイエンス講座)

価格:4,989円
(2019/10/4 02:36時点)
感想(0件)

微分方程式を解析的に解く

以下のような $y = y(t)$ についての2階微分方程式を解く．

$y'' + 2y' + 3y = e^{-2t}, \; y(0) = 0, y'(0) = 0$

まずは解析解をサクッと導出すると，次のようになる．

$y(t) = \frac{1}{3}e^{-2t} + e^{-t} \left( -\frac{1}{3}\cos(\sqrt 2\cdot t) + \frac{1}{3\sqrt 2}\sin(\sqrt 2\cdot t) \right)$

微分方程式を数値的に解く

微分方程式のRunge-Kutta法による解法

微分方程式を数値的に解くには独立変数 $t$ を離散化し， $t_0, t_1, t_2\cdots$ のように書く．離散化された $t$ に対する $y(t)$ の値を $y_n = y(t_n)$ と決める．

次のような1階の微分方程式

$\frac{\mathrm dy}{\mathrm dt} = f(t, y), \; y(t_0) = y_0$

の近似解$y_n$はRunge-Kutta法により，次式のように計算できる．

$y_{n+1} = y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)$

$k_1 = hf(t_n, y_n)$

$k_2 = hf(t_n+\frac{h}{2}, y_n+\frac{k_1}{2})$

$k_3 = hf(t_n+\frac{h}{2}, y_n+\frac{k_2}{2})$

$k_4 = hf(t_n+h, y_n+k_3)$

ここで， $h =t_{n+1} - t_n$ とする．

m元連立微分方程式におけるRunge-Kutta法

上の節ではひとつの微分方程式についてのRunge-Kutta法について述べたが， m元の連立微分方程式に対してもRunge-Kutta法を使うことができる．

変数を $\mathbf y = (y^{(1)}(t), \cdots, y^{(m)}(t))^\top$ ,

$\mathbf y_0 = (y_0^{(1)}, \cdots, y_0^{(m)})^\top$ ,

$\mathbf f = (f_1(t, \mathbf y), \cdots, f_m(t, \mathbf y))^\top$

とおくと，次のようなm元の連立微分方程式

$\frac{\mathrm d\mathbf y}{\mathrm dt} = \mathbf(t, \mathbf y), \; \mathbf y(t_0) = \mathbf y_0$

の近似解$\mathbf y_n$はRunge-Kutta法により，次式のように計算できる．

$\mathbf y_{n+1} = \mathbf y_n + \frac{1}{6}(\mathbf k_1 + 2\mathbf k_2 + 2\mathbf k_3 + \mathbf k_4)$

$\mathbf k_1 = h\mathbf f(t_n, \mathbf y_n)$

$\mathbf k_2 = h\mathbf f(t_n+\frac{h}{2}, \mathbf y_n+\frac{\mathbf k_1}{2})$

$\mathbf k_3 = h\mathbf f(t_n+\frac{h}{2}, \mathbf y_n+\frac{\mathbf k_2}{2})$

$\mathbf k_4 = h\mathbf f(t_n+h, \mathbf y_n+\mathbf k_3)$

ここですごいのは，未知関数が増えてもRunge-Kutta法は同じものがそのまま使えること．すごい．

2階の微分方程式を2元の連立微分方程式に書き直す

2階の微分方程式のままだとRunge-Kutta法が使えないので式変形をする．ちょっとトリッキーな感じがするが，制御工学で状態方程式使うときはよくやる変形だと思う．

$y'' + 2y' + 3y = e^{-2t}$

上式において， $x_1(t) = y(t)$ , $x_2(t) = y'(t)$ と置くと，

$x_2' + 2x_2 + 3x_1 = e^{-2t}$

と書ける．これを $x_2'(t)$ について解くと次式のようになる．

$x_2' = -2x_2 - 3x_1 + e^{-2t}$

あとは， $x_1'(t) = x_2(t)$ と連立させて，

$\begin{cases}x_1'(t) = x_2(t) \\ x_2'(t) = -2x_2 - 3x_1 + e^{-2t} \end{cases}$

とすれば連立微分方程式のできあがり．ここで，

$x(t) = (x_1(t), x_2(t))^\top$

[tex: A = \begin{pmatrix}0 & 1 \\ -3 & -2\end{pmatrix}]

[tex: B = \begin{pmatrix} 0 \\ 1 \end{pmatrix}]

[tex: C = \begin{pmatrix} 1 & 0 \end{pmatrix} ]

（ここの数式の表示のバグ全然直らないのでLaTeXコマンドから読み取ってほしい）

とおくと，

$\begin{cases}\mathbf x'(t) = A\mathbf x(t) + Be^{-2t} \\ y(t) = C\mathbf x(t) \end{cases}$

のように，状態方程式の形で書ける．状態方程式の一つ目の式は「 $\frac{\mathrm d\mathbf y}{\mathrm dt} = f(t, \mathbf y)$ 」の形になっているので，そのままRunge-Kutta法に組み込める．

Juliaプログラム

Runge-Kutta法のJulia 1.0実装は以下の通り．数式がそのままコードになる感じがとても気持ち良い．

using Plots

function main()
    
    # 時間の定義
    t = collect(0:0.005:5)
    h = t[2] - t[1]
    
    # 解析解
    y_analytical = exp.(-2 * t) / 3.0 + exp.(-t) .* ((-1.0/3.0) * cos.(√2 * t) .+ (1.0 / (3.0*√2.0)) * sin.(√2 * t))
    
    # 数値解を突っ込む配列
    x = zeros(2, size(t, 1));
    
    # 初期条件を入れる
    x[:, 1] .= 0.0
    
    # 微分方程式の定義
    A = [0 1; -3 -2];
    B = [0; 1];
    C = [1 0];
    f(t::Float64, x::Array{Float64, 1}) = A * x + B * exp(-2*t)
    
    # Runge-Kutta法を実行する
    for n in 2:size(t, 1)
        k1 = h * f(t[n-1], x[:, n-1])
        k2 = h * f(t[n-1] + h / 2.0, x[:, n-1] + k1 / 2.0)
        k3 = h * f(t[n-1] + h / 2.0, x[:, n-1] + k2 / 2.0)
        k4 = h * f(t[n-1] + h, x[:, n-1] + k3)
        x[:, n] = x[:, n-1] + (k1 + 2*k2 + 2*k3 + k4) / 6.0
    end
    
    # 数値解を取り出す
    y_numerical = C * x
    y_numerical = reshape(y_numerical, size(y_numerical, 2));
    
    # 2つの解をプロットする
    plot(t, y_numerical, line=(:steppre, :dash, 0.5, 4, :blue), label="numerical"), 
    plot!(t, y_analytical, line=(:steppre, :dot, 0.5, 4, :red), label="analytical")
    savefig("runge-kutta.png");
end

main()

実行結果は以下の通り．解析解とよく一致している．

f:id:eqseqs:20181101220458p:plain

（旧）反面教師あり学習

Negative Supervised Learning

2階定数係数線形微分方程式をRunge-Kutta法で数値的に解いて解析解と比べた