反面教師あり学習

*/

反面教師あり学習

Negative Supervised Learning

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

やったこと

  • 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