解析力学

【解析力学:放物線上に束縛された運動(3)】運動方程式をPythonで解く。

こんにちは(@t_kun_kamakiri)(‘◇’)ゞ

本記事は前回の続きとなります。
前回は以下のお題を出し、放物線上を運動する球の運動方程式を導出しました。

内容をさらっと振り返っておきましょう。

内容の振り返り

2つの球を初期の高さを変えて同時に静かに離すと、どちらが先に最下端に到達するか?

  • 球は質点とみなす
  • 空気抵抗は考えない
  • 摩擦もなし
  • 放物線$y=\frac{1}{2}x^2$上を運動

運動方程式は以下となりました。

\begin{align*}
\ddot{x}=-\frac{x}{1+x^2}\dot{x}^2-\frac{g}{1+x^2}x\tag{1}\end{align*}

これを$[x_{0},x]$まで積分すると、

\begin{align*}
t=\pm\frac{1}{\sqrt{g}}\int^{x}_{x_{0}} \sqrt{\frac{1+x^2}{x_{0}^2-x^2}}\, dx\tag{2}
\end{align*}
となります。

右辺の積分の結果は第2種楕円関数$E(x|k^2)$を使うことで以下

\begin{align*}
t=\pm\frac{1}{\sqrt{g}}E\bigg(\sin^{-1}\big(\frac{x}{x_{0}}\big)|-x_{0}^2\bigg)\tag{3}
\end{align*}

前回のおさらいはここまでです。
本記事では数値計算を使って放物線上に束縛された質点の運動方程式を数値計算で解きます

本記事の内容
  1. (1)の運動方程式を数値計算で解く
  2. (2)の運動方程式の解を積分して数値的に解く

数値計算のプログラミング言語はPythonを用います。

【Pythonの使用環境について】
Anaconda
Google Colaboratory
当ブログはこれをメインに使った計算結果です。

関連記事

当ブログではPythonの基礎文法や数値解析の記事をそろえていますので参考にしてください。

(1)の運動方程式を数値計算で解く

放物線上を運動する質点の運動方程式がわかっているので、数値計算で求めてやれば解の振る舞いを知ることができますね。

運動方程式をおさらいしておきます。

\begin{align*}
\ddot{x}=-\frac{x}{1+x^2}\dot{x}^2-\frac{g}{1+x^2}x\tag{1}\end{align*}

積分は4次のルンゲクッタを使います。

初期位置$x_{0}=1.0[m]$テスト計算

試しに初期位置$x_{0}=1.0[m]$として計算してみます。

【結果】

$x$方向を行ったり来たりしている結果ですね。
放物線上を運動するので正しい結果に見えます。

初期位置から最下端の時刻の関係

では、今度は初期位置と最下端の時刻の関係をグラフ化します。

最下端までにかかる時間なので積分区間は$[x_{0},0]$までとなります。

\begin{align*}
t_{bottom}=\pm\frac{1}{\sqrt{g}}\int^{0}_{x_{0}} \sqrt{\frac{1+x^2}{x_{0}^2-x^2}}\, dx\tag{2}
\end{align*}

30秒くらい時間がかかります。

【結果】

距離を見ると青球と赤球で最下端に来る時間には0.3秒くらいの差しかないことがわかります。

カマキリ

0.3秒=300msは瞬きするくらいの時間間隔ですね

名古屋市科学館で実験すると確かに初期高さが低い球の方が早くしたに到達しました。
実際の球は大きさがある球体のため回転の運動方程式も考える必要がありますが、今回の数値計算は質点としてモデル化して考えました。

(2)の運動方程式の解を積分して数値的に解く

\begin{align*}
\ddot{x}=-\frac{x}{1+x^2}\dot{x}^2-\frac{g}{1+x^2}x\tag{1}\end{align*}

これを初期位置から最下端までの位置$[x_{0},0]$で積分すると、

\begin{align*}
t_{bottom}=\pm\frac{1}{\sqrt{g}}\int^{0}_{x_{0}} \sqrt{\frac{1+x^2}{x_{0}^2-x^2}}\, dx\tag{2}
\end{align*}
になることが前回の記事で示したので、(2)を数値計算で求めて(1)の数値計算結果と比較してみます。
間違いがなければ一致するはずですね。

(2)の積分にはscipyのライブラリを使用しています。

【結果】

  • 赤線:(1)の運動方程式の数値計算
  • 黒線:(2)の運動方程式の解を積分した数値計算

重なっているので(2)の積分は間違っていなさそうです。

まとめ

放物線上を運動する質点の運動方程式の振る舞いを数値計算で確認しました。

今回解いた方程式はこちらです。

\begin{align*}
\ddot{x}=-\frac{x}{1+x^2}\dot{x}^2-\frac{g}{1+x^2}x\tag{1}\end{align*}

これを初期位置から最下端までの位置$[x_{0},0]$で積分すると、

\begin{align*}
t_{bottom}=\pm\frac{1}{\sqrt{g}}\int^{0}_{x_{0}} \sqrt{\frac{1+x^2}{x_{0}^2-x^2}}\, dx\tag{2}
\end{align*}
となります。

(2)をもう少し頑張って積分した結果が以下となることがわかっています。

\begin{align*}
t=\pm\frac{1}{\sqrt{g}}\sqrt{x_{0}^2-x^2}E\bigg(\sin^{-1}\big(\frac{x}{x_{0}}\big)|-x_{0}^2\bigg)\tag{3}
\end{align*}

右辺の積分の結果は第2種楕円関数$E(x|k^2)$を使えばよく、追々頑張って導出したいと思います(‘ω’)ノ

解析力学・量子論 第2版

解析力学・量子論 第2版

須藤 靖
3,080円(11/28 20:54時点)
Amazonの情報を掲載しています

【プロフィール】

カマキリ
(^^)

大学の専攻は物性理論で、Fortranを使って数値計算をしていました。
CAEを用いた流体解析は興味がありOpenFOAMを使って勉強しています。

プロフィール記事はこちら

 

大学学部レベルの物理の解説をします 大学初学者で物理にお困りの方にわかりやすく解説します。

このブログでは主に大学以上の物理を勉強して記事にわかりやすくまとめていきます。

  • ・解析力学
  • ・流体力学
  • ・熱力学
  • ・量子統計
  • ・CAE解析(流体解析)
  • note
    noteで内容は主に「プログラミング言語」の勉強の進捗を日々書いています。また、「現在勉強中の内容」「日々思ったこと」も日記代わりに書き記しています。
  • youtube
    youtubeではオープンソースの流体解析、構造解析、1DCAEの操作方法などを動画にしています。
    (音声はありません_(._.)_)
  • Qiita
    Qiitaではプログラミング言語の基本的な内容をまとめています。
関連記事もどうぞ