第4回　物理のシミュレーションをしてみよう

Last updated at 2024-12-23Posted at 2024-12-23

今回は、第3回までの内容を応用して、ボールの軌道をシミュレーションしてみます。
いよいよシリーズのメインディッシュです！張り切っていきましょう！！

問題設定

$ x_0 = 0$, $y_0 = 333 $, $t=0$, 初速$v_0 = 50 $m/s, 地面からの角度$\theta = 45$度でボールを投げた。地面に付くまでの軌道をプロットせよ。
重力加速度は9.8 m/s^2 とする。

実装の手順としては、

シミュレーションする運動方程式を書き出してみる
時間のパラメタ動かして、座標の時間発展を追いかける
グラフにプロット

という感じでしょうか。

皆さんにはもうこれを実装するだけの基礎が身についているはずです。ぜひ挑戦してから解説を読んでください。

解説

もちろんこのくらいだったら手で解いて$y = f(x)$の式を求めることができますが、あえてそれをせずに、時事刻々$x$, $y$の座標を求めることで軌道をプロットしましょう。

まずは方程式を整理します。今回は$t$のみがパラメタで
$$x=v_0\cos(\theta)t$$ $$y = h + v_0 \sin(\theta) t - \frac{1}{2}gt^2$$
と書けますね。この2式に$t$を代入してどんどん増やしていき、$y\le0$になったらやめれば良さそうです。
まずは初期条件を用意しておきます

t = 0
h = 333
v0 = 10
theta = 45

for文を用いて、十分大きな回数繰り返すのも手ですが、今回は条件を満たす間繰り返すwhile文を使ってみます。while y>0:であれば y>0 である限り繰り返します。
やっていることはfor文と同じです。

import numpy as np

t = 0
h = 333
v0 = 10
theta = 45

while y>0:
  x = v0*np.cos(theta)*t
  y = h + v0*np.sin(theta)*t - 0.5*9.8*t**2

  t = t+0.1

1行目では三角関数を計算するためにライブラリをimportしておきました。
y>0の限り、x, yを求めてからtに0.1を足します。
ただ、三角関数を計算するのに用いたnp.cosなんかは、引き数の単位が度だったりラジアンだったりするので、調べるかテストして確かめておきましょう。

print(np.cos(60))
print(np.cos(60/180*np.pi))

実行結果

-0.9524129804151563
0.5000000000000001

どうやらnp.cos()の中身にはラジアンを入れるようですね。
では、初期条件を書き直しておきます。

import numpy as np

t = 0
h = 333
v0 = 50
theta = 45/180*np.pi

x = 0
y = h
while y>0:
  x = v0*np.cos(theta)*t
  y = h + v0*np.sin(theta)*t - 0.5*9.8*t**2

  t = t+0.1

whileを実行する前に、x, yの初期値を入れておきました。こういう操作を初期化と呼びます。プログラミングでは重要な手順です。

それでは、プロットするためにデータをlistに保存しておきましょう。
「listの最後に要素を追加」は、[list名].append([追加する値])です。

import numpy as np

t = 0 # [s]
h = 333 # [m]
v0 = 50 # [m/s]
theta = 45/180*np.pi  # [rad]

xData = []
yData = []

x = 0
y = h
while y>0:
  x = v0*np.cos(theta)*t
  y = h + v0*np.sin(theta)*t - 0.5*9.8*t**2

  xData.append(x)
  yData.append(y)

  t = t+0.1

print(xData)
print(yData)

なんだかそれっぽいlistが二つ出力されました。これで良さげです。

最後にプロットする部分を追加しましょう。

import numpy as np
import matplotlib.pyplot as plt

t = 0 # [s]
h = 333 # [m]
v0 = 50 # [m/s]
theta = 45/180*np.pi  # [rad]

xData = []
yData = []

x = 0
y = h
while y>0:
  x = v0*np.cos(theta)*t
  y = h + v0*np.sin(theta)*t - 0.5*9.8*t**2

  xData.append(x)
  yData.append(y)

  t = t+0.1

# print(xData)
# print(yData)

plt.plot(xData, yData)
plt.show()

うまくいきました!
実際にコードを書くときは、ちょいちょいprint()で問題ないか確認しつつ進めています。うまくいっていたらコメントアウトすればOKです。

おまけ（アニメーション）

せっかくなので、シミュレーション結果をアニメーション表示してみましょう。若干発展的な内容なので、ここからは理解できなくても問題ありません。

これまた便利なライブラリがあるので、ありがたく使わせていただきます。

import matplotlib.animation as animation

ani = animation.ArtistAnimation(fig, ims, interval=100)

こんなのがあるらしいです。第1引数にプロットする図を、第2引数にはアニメーションにする図を、第3引数に図を送る間隔[ms]を指定するようです。
では、while文で画像をlistに保存しておきましょう。

ims = []

x = 0
y = h
while y>0:
  x = v0*np.cos(theta)*t
  y = h + v0*np.sin(theta)*t - 0.5*9.8*t**2

  im = plt.plot(x, y, 'bo')
  ims.append(im)    # グラフを配列 ims に追加

  t += 0.1

while文中でプロットして、それをimsという配列に追加しています。plot()の第3引数でプロットする点の見た目を指定しています。

あとは先ほどのArtistAnimation()を使うだけですが、これはグラフの下地を用意しておいて、どんどんパラパラ漫画を重ねていくイメージなので、下地は作っておかなければいけません。

fig = plt.subplots()  # Figureの生成

これで準備が整いました。以下のコードを実行してみてください。
（ちょっと重かったので初期値を小さくしていますが、気にしないでください）

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

t = 0 # [s]
h = 50 # [m]
v0 = 10 # [m/s]
theta = 45/180*np.pi  # [rad]

fig, ax = plt.subplots()  # Figureの生成

ims = []

x = 0
y = h
while y>0:
  x = v0*np.cos(theta)*t
  y = h + v0*np.sin(theta)*t - 0.5*9.8*t**2

  im = plt.plot(x, y, 'bo', markersize=4)
  ims.append(im)    # グラフを配列 ims に追加

  t += 0.1


# 軸の設定
ax.set_xlim(0, v0 * np.cos(theta) * t * 1.1)
ax.set_ylim(0, h*2)
ax.set_xlabel('x [m]')
ax.set_ylabel('y [m]')

# 10枚のプロットを 100ms ごとに表示
ani = animation.ArtistAnimation(fig, ims, interval=100)   # intervalは[ms]

ani.save("trajectory.gif")

google colabの画面上ではgifが表示できないようなので、フォルダに保存しています。

うまくいきました!
せっかくなので、Chat GPTに聞いてアニメーションの速度を現実の時間と合わせて、グラフに時刻も表示してみましょう。

適当なことも言ってくるので、何度かやり取りして、僕はこんな感じになりました。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

t = 0 # [s]
h = 50 # [m]
v0 = 10 # [m/s]
theta = 45 / 180 * np.pi  # [rad]

fig, ax = plt.subplots()  # FigureとAxesの生成

ims = []
texts = []

x = 0
y = h
while y > 0:
    x = v0 * np.cos(theta) * t
    y = h + v0 * np.sin(theta) * t - 0.5 * 9.8 * t ** 2

    im = ax.plot(x, y, 'bo', markersize=4)
    ims.append(im[0])  # グラフを配列 ims に追加

    # 時刻を表示
    text = ax.text(0.9, 0.9, f't = {t:.1f} s', transform=ax.transAxes, ha='right', va='top', fontsize=12)
    texts.append(text)

    t += 0.1

# 軸の設定
ax.set_xlim(0, v0 * np.cos(theta) * t * 1.1)
ax.set_ylim(0, h * 2)
ax.set_xlabel('x [m]')
ax.set_ylabel('y [m]')

# フレームレートと1フレームあたりの表示時間の設定
fps = 10  # フレームレート (frames per second)
frame_duration = 1 / fps  # 1フレームあたりの表示時間

def update(frame):
    for im, text in zip(ims, texts):
        im.set_data([v0 * np.cos(theta) * frame], [h + v0 * np.sin(theta) * frame - 0.5 * 9.8 * frame ** 2])
        text.set_text(f't={frame:.1f}s')
    return ims + texts

# アニメーション作成
ani = animation.FuncAnimation(fig, update, frames=np.arange(0, t, 0.1), interval=frame_duration * 1000, blit=True)

ani.save("trajectory.gif")

結構いい感じになりました。
パラメタをいじったり、グラフを整えたりして遊んでみてください。

例題

円運動や単振動など、好きな軌道をプロットせよ

いかがでしたか？
世の中のシミュレーションは、結局こんな感じのことをしています。
なんとなくイメージは掴めてきたでしょうか。

これで本シリーズのメインが終了です。
プログラミングを使ったシミュレーションのイメージがなんとなく掴めたでしょうか。
プログラミングに初めて触れる物理学生にとって楽しい記事になっていたら幸いです。

ここから応用編、番外編と続きますので、是非もう少しだけお付き合いください！

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up

第4回 物理のシミュレーションをしてみよう

問題設定

解説

おまけ（アニメーション）

例題

第4回　物理のシミュレーションをしてみよう