AMD系インスタンスでJuliaを始めよう④【オイラー法編】

Contents

1 はじめに
2 前提
3 コードと結果
4 コード解説
5 終わりに

はじめに

　みなさん、こんにJulia！（この挨拶流行らせません？）

　前回の記事では地球トンネルについて常微分方程式（以下ODE）の形式でもモデリングを行い、手計算でがんばって解いて、その上でJuliaで運動の様子を描画したり、トンネル出口への到達時間を計算しましたね。

　上記のような計算ができたのは、解析解が得られていたからです。皆さん、解を仮定したりなんかの演算子法を使ったり思い思いの方法で解析解を手計算で得ていたんですね。でも、正直しんどいかったです。筆者の場合、だいぶ恣意的な仮定を入れて大幅に計算を短縮したのですが、今はもう手計算するとしたら、ここら辺の微分方程式までで許してほしいな、と思いました。偏微・・・？

　もうODEもコンピュータに解かせちゃいましょう！今回は実際の数値演算ではほとんど使用されておりませんが微分の定義から簡単に導出できるオイラー法という手法でODEを解く手法を解説します。

前提

この記事は下記の条件を満たした方を対象としております（特にインスタンスにJuliaをセットアップしていることを前提としております）。

HGAでインスタンスを作成済であること（推奨インスタスタイプはamd系、筆者はamd1dl使用）
インスタンスにログイン可能であること（方法, VSCodeでの接続も推奨します）
動作確認はJupyter Labを使って行いました（参考記事）
インスタンス内にJulia環境がセットアップされていること（参考記事）
前回の記事の続きです

コードと結果

　今回の検証に用いたコード全体および実行結果を示します。

# オイラー法を使って解を求めましょう（詳しい解説は記事にて）
# ここから始める人のための準備
using PhysicalConstants.CODATA2018: G
# r'' = -(4/3)Gπρr = -Cr
ρ = 5514  # 地球の平均密度 [kg/m^3]
R = 6378137  # 地球半径 [m] 
G_float = 6.6743e-11  # 本当はPhysicalConstantsから値のみを取り出せるようにしたい・・・
C = (4/3)*π*G_float*ρ

## オイラー法使用のための準備
t = 0:1:7000  # プロット範囲
T = length(t) # 7001
h = 1

function delta_r1(t, r1, r2)
    r2
end

function delta_r2(t, r1, r2)
    -C*r1
end

## 初期値の設定やループ内で用いる変数の初期化
t_num = [0.0]
r_num = [float(R)]

r1_tmp = r_num[1]
r2_tmp = 0
for j = 1:T-1
    t_tmp = j*h
    r1_new = r1_tmp + h*delta_r1(t_tmp, r1_tmp, r2_tmp)
    r2_new = r2_tmp + h*delta_r2(t_tmp, r1_tmp, r2_tmp)    
#     println(t_tmp + h," " ,r1_new, " " ,r2_new)
    r1_tmp = r1_new
    r2_tmp = r2_new
    push!(t_num, t_tmp+h)
    push!(r_num, r1_new)
end
println(length(t_num))
println(length(r_num))


## グラフ描画
p1 = plot([t t_num], [r_exact_4plt r_num], label = ["analytic" "numerical"], xlabel="Time [sec]", ylabel="Displacement [m]")
relative_err = ((r_num - r_exact_4plt) ./ r_exact_4plt) * 100  # 相対誤差
p2 = plot(t , relative_err, label = "Relative Error", xlabel="Time [sec]", ylabel="relative error [%]")
plot(p1,p2,layout=(1,2),size=(1000,500))

Jupyterで検証する場合、セルに上記を貼り付けCtrl + Enterキーで下記のような出力が得られます。

　まずは、結果のグラフについて確認しましょう。左側は、オイラー法による数値計算結果と解析解をプロットした結果を重ねてプロットしたものです。一見して見てもほとんど違いがないことがわかると思います。単振動レベルなら、オイラー法でも十分な精度が得られると思います。

　一方、右側は解析解に対する相対誤差を求めプロットしたものです。6000 [sec]までは十分小さい誤差と言えるのではないでしょうか。また、ループ文で加算していくという性質から誤差は線形比例で増加しています。気になるのは、変位が０の位置にエラーが発散していると見られる箇所があることですが、これは相対誤差の定義で解析解を分母にとっているので、ここの箇所が０近傍で発散的なふるまいをしてしまうんですよね・・・。そもそも相対誤差の形からして不定形に近いですしね。どうしても気になる人はローパスフィルタをかませてください。

コード解説

オイラー法について

　オイラー法とは微分方程式を数値演算で解く方法です。ところで（常）微分方程式を解くってどういうことかというと

$\frac{dr(t)}{dt} = f(r(t))$

という形の式が与えられたとき、r(t)を求めるものでした。ODEの場合、初期条件が必要でした。

　上記のようなODEの場合、数値的に解くとしたら①（数値）積分する方法と②少なくとも形式的には積分しない方法の２通りがあります。手計算の時と大体同じですよね！定数（ODEで言えば初期条件など）の差を除けば、微分と積分は逆演算の関係です。ですから、積分すれば済むのであれば、積分するのがストレートフォワードです。

　ところが、数値計算（手計算でも）では、積分の計算はなかなか大変です。少なくとも微分に比べるとずっと。よって汎用性を求めると、積分しない方法を採用することになるでしょう。（数値）積分をしないで微分方程式の解を求める方法として、オイラー法を紹介します。

　オイラー法とは微分の定義（ないしテイラー展開を~~恣意的に~~打ち切った結果）を基に、微分方程式を解こうというアプローチです。これから簡単に要領を解説していきますがその際、微分の定義から、頑張って漸化式得ればループ回していけるぜ！とブツブツつぶやきながら式変形を追えば、大まかな理解が得られて良いと思います。ただし医者に連れていかれても責任は取れません。

　さて、十分小さい数をhを用いると、微分の定義から上式は次のように変形されます：

$\frac{r(t+h)-r(t)}{h} = f(r(t))$

よって

$r(t+h)} = r(t) + f(r(t))h$

ここで、時間0 ~ tをn分割することを考えます。均等にn分割するとして

$h = t_{n} - t_{n-1} \ \ \ \ where\ t=t_{n}$

としてあげれば良いことがわかります。実装上、hはユーザ自身が決める定数になります。なかなか直観的に決めづらいと感じる場合は終端時刻Tと分割数nに検討がついているとして

$h = \frac{T}{n-1}$

という関係から求めるのもいいかもしれません（今回は1 [sec]刻みとエイヤで決めています）。

　上記の関係を用いて式を次のように書き換えることができます：

$r(t_{n}) = r(t_{n-1})+f(t_{n-1})h$

これは時間をn分割した際、隣り合う時刻同士の漸化式を表すことがわかります。実際、nの値をどんどん増やしていくと

$r(t_{1}) = r(t_{0})+f(t_{0})h\\r(t_{2}) = r(t_{1})+f(t_{1})h\\\ldots$

という風に順次計算できることがわかるでしょう。なお、右辺第二項は考える時間領域において常に計算可能であるものとしてます。こうやってnを１インクリメントさせて反復計算するのは大体のプログラミング言語の基本的な反復文でできます。

二階ODEに対するオイラー法の適用

　さて、オイラー法の基礎について解説しましたが、一つ不十分な点があります。それは、解説で例に上げたのは一階ODEなのに対し今回解く必要があるのは二階ODEである点です。したがって、前節の内容をそのまま今回のケースに適用して実装するという訳にはいきません。残念ながら。

　もうひと工夫必要です。そこで、以下のように変数r1,r2を定義してみましょう：

$r_{1} = r,\ r_{2}=\frac{dr}{dt}$

すると

$\frac{d^2r(t)}{dt^2} = -\frac{4}{3}\pi\rho G r(t) = -Cr(t)\\C: const.$

を次のように書き換えることができます：

$\left\{ \begin{array}{l} \frac{dr_{1}}{dt}=r_{2} \\ \frac{r_{2}}{dt}=-Cr \end{array} \right.$

これって、結局オイラー法解説のときに使った式が一本増えただけに感じませんか？まったくその通りです。ただし、実装上（も手計算で解く上でも）初期条件がもう一つ必要になります。具体的には、r2に対する初期条件が必要なのですが、これは変位の時間微分（＝速度）であり、今回は初速０でトンネルに入るという前提なので、たやすく求まりました。次節ではもう少し実装を見ていきましょう。

二階ODEに対するオイラー法の実装

　二階ODEをオイラー法で解いてる部分の実装は次の箇所に当たります：

## オイラー法使用のための準備
t = 0:1:7000  # プロット範囲
T = length(t) # 7001
h = 1

## 連立微分方程式を定義
function delta_r1(t, r1, r2)
    r2
end

function delta_r2(t, r1, r2)
    -C*r1
end

## 初期値の設定やループ内で用いる変数の初期化
t_num = [0.0]
r_num = [float(R)]

## オイラー法に従って計算
r1_tmp = r_num[1]
r2_tmp = 0
for j = 1:T-1
    t_tmp = j*h
    r1_new = r1_tmp + h*delta_r1(t_tmp, r1_tmp, r2_tmp)
    r2_new = r2_tmp + h*delta_r2(t_tmp, r1_tmp, r2_tmp)    
#     println(t_tmp + h," " ,r1_new, " " ,r2_new)
    r1_tmp = r1_new
    r2_tmp = r2_new
    push!(t_num, t_tmp+h)
    push!(r_num, r1_new)
end

　上から解説しますと、まず前節の連立微分方程式をそれぞれ関数として実装しております。次に変位r1と速度r2に対する初期条件をそれぞれ定義し、オイラー法で計算を回すという実装の流れです。ループ内を見て々節で解説した内容がそのままコード化している点に着目してください（ただし、二個分になっています）。要は漸化式をループ分で順々に計算しているんですね。

　もし、まだわかりづらいところがあったらTwitterやメールでお気軽にご質問いただければと思います。とりわけ、ループ文や（数学関数ではない）関数など、プログラミング初心者の方の場合、理論は分かっても実装がわからん！といういて当然だと思いますし、筆者も最初はループ文は使い方がよくわからなかった記憶があります。漸化式やサメンションをループ回して計算するのは、数値計算ではよくある手法なので、理解したら自分でも実装してみて慣れることをおすすめします。

他の部分の解説

　オイラー法でぐるぐる計算した後は、計算内容を解析解と合わせてプロットしたり、相対誤差を計算してプロットしたりしています：

## グラフ描画
p1 = plot([t t_num], [r_exact_4plt r_num], label = ["analytic" "numerical"], xlabel="Time [sec]", ylabel="Displacement [m]")
relative_err = ((r_num - r_exact_4plt) ./ r_exact_4plt) * 100  # 相対誤差
p2 = plot(t , relative_err, label = "Relative Error", xlabel="Time [sec]", ylabel="relative error [%]")
plot(p1,p2,layout=(1,2),size=(1000,500))

注意していただきたいのは４行目、相対誤差を定義している箇所で、配列の要素ごとの計算を行うために./のように頭にドットをつけた演算子を用いてることです。要素毎に計算するときはドット（.）をつけるというのを忘れないようにしましょう。

終わりに

　以上で解説終了です！今回は古典力学の枠組みでモデリングした地球トンネルモデルから、オイラー法と呼ばれる手法で運動の様子を得ることができましたね！そして、おおよその結果ばかりではなく精度についても、（数値演算した）解析解の結果と数値解析解の結果がそれほど悪くないことがと感じられたと思います。地球トンネルモデルに、オイラー法を適用するにあたって、少し技巧的な式変形を~~余儀なくされた~~行いましたが、これは慣れの問題です。この技巧について補足しますと、いわゆる現代制御理論ではよく用いるアプローチです。制御の場合、機械システムであれ電気回路であれ制御対象が高々二階の微分方程式で記述できることが多く、前節のような式変換を行うことで比較的簡単にダイナミクスを記述することができます。ダイナミクスとは一階時間微分のことです。上記のリンクを見てほしいのですが、ベクトルや行列を無視して形式にさえ着目すれば、時間微分の微分方程式として記述され、ものの時間変化を取り扱っているのが一目瞭然です。

　最後に、オイラー法および微分方程式の数値解法について少し補足します。今回はオイラー法と手法を用いて、それなりに良い精度だと実感できる結果を得ることができましたが、よくよく結果を見ると時間を減るにつれ（初期状態たる原点から遠ざかるにつれ）精度が悪くなっているのがわかります（発散箇所については一旦忘れてください）。これは、ループを回して漸化式を数値演算していくという性質から、誤差についてもどんどん加算されていくという性質のためです。直線状に誤差が増えているのがまさにその証左で、具体的にはループ内の足し算という演算の結果誤差が蓄積されていると結論づけることができます。この性質が実際の数値計算分野で使用されることはまずありません。例えばMATLABユーザは分野にもよりますがかなりの方がODE４５を選択されていると思います。これはルンゲクッタ法というものをベースにしたソルバ（微分方程式を数値的に解く手法）です。このように、微分方程式を数値計算で解くとは、精度（そして高速性、更に収束性も）との飽くなき戦いの連続であり、取り扱う方程式の形により、用いるアルゴリズムも様々であることがおいおい実感できると思います。そしてループ箇所や配列を並列処理対応にしてGPUリソースで計算したいと思ったあなた！今すぐHGAに登録しましょう！（公式ブログらしく宣伝もばっちりですね）

[HGAoffcial]

[dot]

[/dot]