AMD系インスタンスでJuliaを始めよう③【ODE解析解編】

Contents

1 はじめに
2 前提
3 モデリング
4 解析解の導出
5 Juliaによる運動の様子の描画・到達時間の計算
6 パッケージの導入
- 6.1 運動のプロットと到達時間の算出
7 次回予告：ODEを数値計算で解く

はじめに

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

　突然ですが、小さい頃「地球（日本）の裏側はブラジルだ」という言葉を聞いたことがありませんか？

　筆者はこの言葉を聞いた時、「地面を延々と掘ればブラジルに行けるやん！」なんてことを考えて公園の砂場だかで友達と延々とスコップを掘った記憶があります。ところがその時は~~残念ながら~~車で通りがかった近所の人（多分）に見つかり散々怒られた挙句、元に戻した後車で家まで連行され今度は親に怒られるはめになった記憶があります。こうしてブラジルへの不法入国の企ては未然に処理されました。

　ところで、筆者は小さいころ仮にブラジルまで貫通するトンネルができたとして、それに入ったらどうなるんだろう？みたいな疑問をしばらく抱いていた記憶があります。この疑問、結構皆さん抱いたことあるのではないでしょうか。

　その疑問、2021年の必修科目Juliaが簡単に答えてくれます！

前提

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

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

モデリング

　地球（半径R）を完全な球体とみなし、かつ密度（ρ）は一様分布であるものとしましょう。そして小さい頃の我々は無事地面を重心を通るよう一直線に掘り進め、無事日本－ブラジル間のトンネルを掘り上げるのに成功したとします（~~ならもうシミュレーション必要ないじゃん~~）。もちろん死者は出ていませんし外交問題や法規上の問題は何ら発生していません（~~シミュレーション必要ないじゃん~~）。早速トンネルに入ってみましょう。

　体重m [kg]の人は、大きな地球に対してはただの点にすぎません。したがって人間は質量mの質点として扱うことができるでしょう。質点に働く力は地球（の重心）との万有引力のみであるとしてます。地球の質量は我々の体重に比べ十分すぎるほど大きいため、万有引力の影響でトンネル経路に影響が出ることもありません。また、地球の重心からトンネル内を通って日本に向かう方向に座標軸rを設けます。トンネル内に入る時（時刻t=0とします）、r成分の速度は０であるとします。

　以上を整理して質点についての運動方程式を立てると次のようになります：

$m\frac{d^2r(t)}{dt^2} = -\frac{4}{3}\pi\rho G m r(t)$

（数式の作成にはこちらのページを利用させていただきました。とても便利です！）

整理して以下の常微分方程式が得られます：

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

Cは既知の値（万有引力定数Gや円周率π）から成る定数です。

　ここで、注意したいのは万有引力を考慮する際、地球の方の質量も位置に応じた（つまり、時変の）関数であることです：

$F(t) = -G\frac{M(t)m}{r^2(t)} = -\frac{4}{3}\pi r^3(t) \rho G\frac{m}{r^2(t)} = -\frac{4}{3}\pi r(t) \rho G$

いろいろ話が飛びましたが、地球トンネル内（および境界）の運動についての数学モデルが得られましたね！

解析解の導出

　（常）微分方程式（以降ODE）が得られたので解析解を導出します。ODEの解き方にはいろいろな手法がございますが、ここでは運動の様子を仮定して、それが式を満たしているか確認するという方法をとってみましょう（筆者としてはこの場合ラプラス演算子法をおすすめします）！

　数学モデル（ODE）導出時の仮定をもう一度振り返ってみましょう。初速0でトンネル入り口（r=R）に入ることを考えます。あなたは（筆者は嫌です）万有引力により地球重心に向けどんどん加速されながら向かっていくでしょう。重心につくと、万有引力は（極限として）0となり、あなたには外部からの力が加わっていないことになります。ところが、車は急に止まれないよろしく、慣性のため重心を過ぎてそのままブラジル（r = -R）に向かうことになりますが、再び重心から離れるにつれ万有引力が働きどんどん減速しながらブラジルに向かうことになります。減速によりブラジルまで到達できないか心配される方もいるかもしれませんが、エネルギ保存則から到達できるのは間違いないです。どちらかというと裁判とかそこら辺のことを心配すべきです。

　実は、ブラジルに到達できても、何もしないでそのまま日本に戻ることになります（ちょうどトンネルに入る時と同じ状態になっていることに注目です！）。どうも周期運動になりそうなイメージですね！

　以上の運動の様子を横軸に時刻tをとり、縦軸に変位r(t)をとって表すと下図のようになります。

　上記のような運動をイメージできれば微分方程式の解を仮定することができます：

$r(t) = Rcos\omega t$

実際に微分方程式に代入してチェックすると

$\omega = \sqrt{\frac{4\pi \rho G}{3}}$

であれば、仮定の式は微分方程式を満たすことがわかります。したがって解析解は

$r(t) = Rcos\sqrt{\frac{4\pi \rho G}{3}}t$

となります

Juliaによる運動の様子の描画・到達時間の計算

　ここからがJuliaの出番です！

今回は前節で得られた解析解をプロットし、更にブラジルまでの到達時間を計算しましょう。

なお、次回の記事ではオイラー法を使った微分方程式の数値計算を行ってみようと思います（実はJuliaには微分方程式をスパッと解くパッケージがあります）。

今回は苦労して得た解析解をベースに運動の様子をみてみましょう！

パッケージの導入

　特別なパッケージは必要がありませんが、物理定数の値が知りたいときはPhysicalConstantsというパッケージがおすすめです。導入には

julia

と入力しREPLに切り替えて]と入力します：

すると下図のようなpkg>という、パッケージマネージメント用のモードに切り替わるのでそこで

add PhysicalConstants

と入力し、JupyterのカーネルでJuliaを実行するためのパッケージが導入できます。

パッケージを導入したらJupyterLab上やREPLで下記のように呼び出しを行い、基本的な物理量を得ることができます。

using PhysicalConstants.CODATA2018
G = NewtonianConstantOfGravitation

運動のプロットと到達時間の算出

　ここでは以前の記事でグラフ描画用のパッケージ（plots）を導入していることを前提として、実際にグラフの描画とブラジルまでの到達時間を求めます。

　先に具体的なコードと、その実行結果を載せちゃいます！

# 解析解に基づきブラジルまでの移動時間の算出とマス変位をプロット
using PhysicalConstants.CODATA2018: G
# r'' = -(4/3)Gπρr = -Cr
ρ = 5514  # 地球の平均密度 [kg/m^3]

println(G)
G_float = 6.6743e-11  # 本当はPhysicalConstantsから値のみを取り出せるようにしたい・・・

C = (4/3)*π*G_float*ρ


## 解析解（記事にて解説）
R = 6378137  # 地球半径 [m] 
ω = sqrt(C)  # 角周波数 [Hz]

function r_exact(t)  # 解析解
    R*cos(ω*t)
end

## プロットの準備
using Plots
gr()

## プロット範囲の策定と演算
t = 0:1:7000
r_exact_4plt = r_exact.(t)


## 地球の反対側まで行く時間を求めます（解析解から）
min_dis_tuple = findmin(r_exact_4plt)
println("ブラジルまで", round(t[min_dis_tuple[2]] / 60, digits=1), "[分]で行けます！")  # 地球の反対まで行く時間 [sec] → [min]


## グラフ描画
plot(t, r_exact_4plt, label = "motion of mass(human...)" ,xlabel="Time [sec]", ylabel="Displacement [m]", size=(800, 400))  # 解析解に基づく変位の描画

上記コードを実行すれば、ブラジルまで42.2 [分]で到達でき、（解析解通り）cos運動をすることがわかりました！

コードについて少し解説すると、まず到達時間を求めるのに

function r_exact(t)  # 解析解
    R*cos(ω*t)
end

## プロット範囲の策定と演算
t = 0:1:7000
r_exact_4plt = r_exact.(t)
## 地球の反対側まで行く時間を求めます（解析解から）
min_dis_tuple = findmin(r_exact_4plt)
println("ブラジルまで", round(t[min_dis_tuple[2]] / 60, digits=1), "[分]で行けます！")  # 地球の反対まで行く時間 [sec] → [min]

のように解析解に具体的に時間配列をそれぞれ代入し、r_exact_4plt というcos(t)の値からなる配列を定義します。一周期以内だと、cos(t)の最小値（r=-R）の箇所における時刻がブラジルまでの到達時間となります。Juliaのfindminという関数を使うと(min値, min時のindex値)というタプルが得られるのでブラジル、つまり運動の最小となる時刻はt[min_dis_tuple[2]]で得ることができます。

　また、グラフ描画については時間配列tとtのそれぞれの値を解析解に代入して得られる解析解配列r_exact_4pltから得ることができます。プロットは以下のコードで行うことができます：

plot(t, r_exact_4plt, label = "motion of mass(human...)" ,xlabel="Time [sec]", ylabel="Displacement [m]", size=(800, 400))  # 解析解に基づく変位の描画

基本的にはplot(x軸配列, y軸配列)でプロットができます。

次回予告：ODEを数値計算で解く

　ここまでお読みになられた方、本当にありがとうございました！

　今回は、解析解をベースにグラフを描画し、トンネル内での運動の様子とブラジルへ行く時間を算出しました。この速さに匹敵するのは弾道ミサイルくらいのものでしょう。

　ところで、ODEを手計算で解くのって何かと手間がかかりますよね！解におおよその当たりがついてればそれほど難しくないのですが、ついてない場合は悲惨の一言につきます。解の仮定時に微積演算上の優位性や三角関数も出せる、ネイピア指数関数を用いたり、演算子法を用いるのも良いでしょう。

　でも、一番楽な方法はそもそも微分方程式を手計算で解こうとなんて思わず、コンピュータに解いてもらうことですね！次回はその方法をご紹介します。オイラー法と呼ばれる手法です（実は微分方程式の解を得るパッケージはあるのですが、記事では扱わない予定です）。

　次回もお楽しみに！

[HGAoffcial]

[dot]

[/dot]