プログラミング入門II
2025.07.16
コンピュータシミュレーション

Back to index page


  1. 本日の作業内容

  2. 前回の宿題について

    今回は解答用紙に番号と名前の無いものが出てきてしまいました．注意しましょう．

    num and name: 24_57

    以下は例によって問題のあるものなのですが，今回はグラフのプロットなので，結果の図でお示しします．

    今回一番多かったものが上の例で，tight_layout() を指定しなかったために右Y軸のキャプションが切れているものでした．

    次に多かったのは凡例の不適切な表示で，以下のようないろんなパターンが見られました．

    一番残念なのは以下のもので，今回の主題であるダブルY軸になっていないもので，さすがにこれはダメですね．

  3. 前回の復習

  4. コンピュータシミュレーション

    Wikipedia によりますと，シミュレーションとはラテン語の similis （似ている），simulare （模倣する），simulat （真似た）などの言葉から生まれた概念とされています．対象となるシステムで働いている法則を推定・抽出し，それを真似るようにして組み込んだモデル，模型，コンピュータプログラムなどを用いて行われます．現実のシステムを動かしてその挙動や結果を確かめることが困難，不可能または危険である場合に用いられるものとされます．

    昨年NHK朝ドラ「虎に翼」で話題になった戦前の日本の「総力戦研究所」が行った模擬戦では日本がアメリカに必ず負けるという結果になったことは有名ですね．

    「シュミレーション」と言われることも多いようですがそれは音位転換であり，「シミュレ」よりも「シュミレ」の方が日本人には発音しやすいことに起因します．

    シミュレーションをコンピュータを用いて行う場合にコンピュータシミュレーションと呼ばれます．

    • 解析的に解けない場合の計算を数値計算を用いて行う
    • 物理モデルに基づいて複雑で膨大な計算を行う
    • 乱数を用いてランダムな挙動を取り込んだ計算を行う

    今回は乱数を用いて実際に実験可能な作業を数値化して画面に動画で実行結果を表示することにします．

    • 再描画
      •

      シミュレーションの状況を動的に提示するために，グラフを逐次再描画する方法が使えます．以下では，それを試してみましょう．

      • 進行する平面波

        • ご存じのように1次元の進行波は以下の式で表現されます．

        f (x, t) = sin(kx - ωt)

        ここで，k は波数，x は座標，ω は各振動数，そして t は時間です．乱数を用いて，k を1から3の定数，ω も同じく1-3の整数として，進行波を表示する動画を作製します．以下のソースをご覧ください．

        from random import *
import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(0, 12, 1000)
k = randint(1, 3)
omega = randint(1, 3)
t = 0
txt = f'sin({k}x-{omega}t)'

fig, ax = plt.subplots()

ax.set_xlim(0, 12)
ax.set_ylabel('sin(kx-ωt)')
ax.set_xlabel('x')
ax.text(0, 0.9, txt, transform=ax.transAxes)

line, = plt.plot(x, np.sin((k*x-omega*t)))

while t < 100:
    line.set_data(x, np.sin(k*x-omega*t))

    plt.tight_layout()
    plt.pause(0.01)

    t += 1

        注目すべきはこれまで学習してきた代入とちょっと違う部分があることでしょうか．17行目を見てください．

        line, = plt.plot(x, np.sin((k*x-omega*t)))

        オブジェクト名の後ろにカンマだけついてきています．これまでタプルを利用した多重代入，例えば以下のようなものを経験してきました．

        a, b = c, d

        上の例では変数 a に変数 c の値が代入され，変数 b には変数 d の値が代入されるのですが，

        a, = c, d

        このようになると，変数 c の値が変数 a に代入されるだけとなります．意味の分からない面倒なものに見えますが，代入するオブジェクトが複数の値のセット（リストやタプルなど）になっている状態では，先頭の値だけを取り出すためにこのような書き方になっています．

      • コイン投げにおける表の頻度

        • 前回の演習問題 iii. でコインを10枚投げた時の表の枚数の頻度について作業しました．今回は頻度がリアルタイムで更新されるよう動画にしたものを試してみましょう．このソースを使用します．なお，デフォルトの試行回数は1万回となっていますが，それが終わるまでには数分かかります．

        from random import *
import matplotlib.pyplot as plt

freq = [0 for _ in range(11)]
x = [i for i in range(11)]
ylim = 100

fig, ax = plt.subplots()

ax.set_ylim(0, ylim)
ax.set_ylabel('Frequency')
ax.set_xlabel('Number of head')

text = ax.text(0.05, 0.95, 'Number of cast:', transform=ax.transAxes)

line, = plt.plot(x, freq, 'o-')

for i in range(10000):
    head = 0
    if max(freq) > ylim * 0.75:
        ylim *= 2
        ax.set_ylim(0, ylim)        

    for _ in range(10):
        head += randint(0, 1)

    freq[head] += 1

    line.set_data(x, freq)

    text.set_text(f'Number of cast: {i + 1}')
    fig.canvas.draw
    fig.canvas.flush_events()
    
    plt.tight_layout()
    plt.pause(0.01)

freq_text = ''
for i in range(11):
    freq_text += str(f'{freq[i]:>8}')

ax.text(0, 0.9, freq_text, transform=ax.transAxes)
plt.show()

    • 3次元のグラフ

      • matplotlib を使用すると３次元のグラフも簡単に描画できます．今回は３次元表現の中でも一番単純な散布図について試してみましょう．

      まず３次元グラフとするためには

      ax = fig.add_subplot(projection='3d')

      を記述し，散布図の場合にはプロットに scatter を使用します．リンク先のソースは座標 (5, 5, 5) を中心に標準偏差 1 で散らばる100個の点の分布です． 

      from random import *
import matplotlib.pyplot as plt
import numpy as np

av = 5
sigma = 1

xyz = [[gauss(av, sigma) for _ in range(100)] for _ in range(3)]

fig = plt.figure()
ax = fig.add_subplot(projection='3d')

ax.set_xlim(0, 10)
ax.set_ylim(0, 10)
ax.set_zlim(0, 10)

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

ax.text(5, 1, 1, f'mean = {av}, stdev = {sigma}')

ax.scatter(xyz[0], xyz[1], xyz[2])

plt.show()

    • アニメーション機能使用

      • matplotlib に備わっているアニメーション機能を利用する処理も体験しておきましょう．

      次の作業は while 文の中で乱数を発生させて，その都度グラフにプロット点を追加するものです．それには matplotlib モジュールの animation を利用します．

      from time import *
from random import randint as randi
from matplotlib import pyplot as plt
from matplotlib import animation

def graph_upgrade(data):
    t, y1, y2 = data
    xdata.append(t)
    y1data.append(y1)
    y2data.append(y2)

    line1.set_data(xdata, y1data)
    line2.set_data(xdata, y2data)

    return line1, line2

def graph():
    ax.set_ylim(0, 100)
    ax.set_xlim(0, 20)
    ax.set_xlabel('time [s]')
    ax.set_ylabel('Random numbers')

    line1.set_data(xdata, y1data)
    line2.set_data(xdata, y2data)

    return line1, line2

def make_data():
    t = 0
    while t < 19:
        t = time() - start
        y1 = randi(1, 50)
        y2 = randi(51, 100)

        print(f'{t:.2f},{y1},{y2}')

        yield t, y1, y2

        sleep(1)

start = time()

fig, ax = plt.subplots()
line1, = ax.plot([], [])
line2, = ax.plot([], [])
xdata, y1data, y2data = [], [], []

ani = animation.FuncAnimation(fig, graph_upgrade,
                              make_data, blit = True,
                              interval = 100,
                              cache_frame_data = False,
                              repeat = False, init_func = graph)

plt.show()

      定義する関数は以下の３つになります．

      • リアルタイムでグラフを描く関数 graph_upgrade
      • グラフの軸を決める関数 graph
      • 乱数を作る関数 make_data

      0.17,31,80
1.21,9,99
2.24,35,89
3.27,17,61
4.30,50,52
5.32,11,85
6.35,7,52
7.37,35,56
8.41,4,71
9.44,39,85
10.47,50,81
11.50,16,73
12.52,29,65
13.55,33,69
14.57,15,95
15.60,1,65
16.64,7,76
17.67,35,67
18.69,11,52
19.73,12,61
21.18,16,74


  5. 演習

    今回の演習問題です．

  6. 宿題

    いつものようにMoodleを使用して行います．公開は本日の18:00です．締め切りは7月23日（水）の10:00です．

  7. 宿題の講評

    宿題の講評は特別ページを用意します．締め切り後に閲覧可能にしますので，気になる人は見てください．

目次ページに戻る