Skip to content

simulation

Directory actions

More options

Directory actions

More options

Latest commit

 

History

History
 
 

simulation

Folders and files

NameName
Last commit message
Last commit date

parent directory

..
fig
fig
 
 
README.md
README.md
 
 
index.html
index.html
 
 

README.md

数値シミュレーション

[Up] [Repository]

  • 空気抵抗がない場合の弾道計算
  • 空気抵抗がある場合の弾道計算
  • 反応拡散方程式(Gray-Scott Model)
  • 余談：パーソナルスーパーコンピュータ

はじめに

互いに重力で引き合う二つの星の軌道の形が楕円となることは知っているであろう。ではなぜ楕円となるか、答えられるだろうか？様々な答えはあろうが、一つの答えは「距離の逆二乗に比例する力で互いに引き合う二つの物質の運動方程式を解くと、その軌道が楕円となるから」である。さて、運動方程式とはなんだったか覚えているだろうか。最も簡単な運動方程式は$F=ma$である。これは、物質にかかる加速度と力が比例し、その比例係数が物質の質量であることを主張する。さて、加速度とは何か？それは速度の時間変化である。したがって、先ほどの運動方程式は、より正確に書くと

$$ m \frac{d v}{dt} = F $$

となる。さらに速度とは、単位時間当たりの位置の変化分であった。したがって、それもちゃんと書くと、

$$ \begin{aligned} \frac{d v}{dt} &= \frac{F}{m} \\ \frac{d r}{dt} &= v \end{aligned} $$

これが運動方程式である。すなわち、運動方程式とは時間に関する連立常微分方程式である。上記の式は一次元で書いたが、これを三次元の式にして、さらに二つの物質の間に距離の二乗に反比例する力を仮定すると、互いに重力で引き合う二つの星の運動を表す方程式となる。これを解くと、軌道が楕円になることや、面積速度一定則、調和の法則など、いわゆるケプラーの三法則が導かれる。

物理学とは、この宇宙を記述する学問である。そして(なぜかはわからないが)この宇宙は微分方程式で記述されている(らしい)。したがって、物理学とは、微分方程式を解く学問である。幸か不幸か、ほとんどの微分方程式は解析的に解くことができない。しかし、方程式さえわかれば、それを数値的に解くことは可能である。何かの現象に着目し、それを記述する方程式を支配方程式(Governing Equation)と呼ぶ。支配方程式を数値的に解いて、その振る舞いを調べることを数値シミュレーションと呼ぶ。というわけで、今回は数値シミュレーションをやってみよう。

空気抵抗がない場合の弾道計算

ものを斜め上に投げると、どの角度で投げると最も遠くに飛ぶだろうか。ただし空気抵抗は無視するものとする。答えが投げる強さ(初速)に依らず45度であることは知っているであろう。では逆に、初速と的までの距離が決まっているときに、何度で投げれば的に当てることができるだろうか？例えば初速100m/sで、500m先にある的に当てたいときの角度は？空気抵抗を無視すれば、二次方程式を解くだけだが、すぐに暗算するのは難しいであろう。この「重力下で物に初速を与えて飛ばす」という設定は、物理としてはかなり簡単な部類に入るが、戦争においては極めて重要な問題設定であった。例えば敵までの距離がわかっているときに、迫撃砲の角度を何度にすればよいかを「すぐに」決めなくてはならない。当然だが戦闘中にいちいち方程式を解く暇はなく、実際には空気的もあるために距離と角度の関係は難しい。そこで、あらかじめ弾の種類と距離の応じて「射表」と呼ばれる距離と角度の関係表を作っていた。実際に射出して着弾距離を調べることももちろん行われたが、数値計算も行われた。最初期の電子計算機であるENIACは、もともと砲撃の射表の作成のために作られたものだ。ENIACは微分方程式を解くことができ、これが本格的な数値シミュレーションの始まりである。ENIACはその後「マンハッタン計画」にも用いられたことからもわかるように、計算機は軍事利用と深い関係にあり、スーパーコンピュータは半ば「兵器」として扱われた。強力な計算機を保有することが軍事的に優位に立つのに必要であり、実は現在もその名残が様々なところにみられるのだが、それはさておく。

さて、いま皆さんの目の前にあるのは、ちょっと前のスーパーコンピュータなみの計算能力を持つ計算機である。それを使って、簡単なシミュレーションをしてみよう。まずは空気抵抗がない場合である。二次元の場合を考えよう。運動方程式は以下のようになる。

$$ \begin{aligned} \dot{v_x} &= 0 \\ \dot{v_y} &= - g \\ \dot{r_x} &= v_x \\ \dot{r_y} &= v_y \end{aligned} $$

さて、この式は厳密に解けるのだが、その厳密解を知らないとしよう。そこで、シミュレーションで近似的に解を求めたい。ここで「解」とは、重力下で角度$\theta$で物体を投げた時の物体の軌道である。まず、数値計算で扱いやすいように、時間を離散化しよう。離散化とは、本来連続である時間や空間を「小さな単位」に分解して近似することだ。例えばテレビは、画像をピクセルに分解して表現しており(空間の離散化)、また1秒間に30フレームを表示している(時間の離散化)。このように、時間も空間も十分小さい単位に区切って表現すれば、もとの連続空間での微分方程式の解の性質を精度よく表現できると考えられる。

今、時間に依存するある量$f(t)$の時間微分の表式$f'(t)$が与えられているとしよう。時間に関する離散化とは、ある小さな時間刻みを$h$に対して、$f(t)$の値からなんとかして$f(t+h)$の値を推定することである。$f(t+h)$を$t$の周りでテイラー展開して、一次までとると、

$$ \begin{aligned} f(t + h) &= f(t) + f'(t)h + O(h^2) \\ &\sim f(t) + f'(t)h \end{aligned} $$

となる。この式は、右辺、すなわち時刻$t$における$f(t)$の値と、その微係数$f'(t)$がわかっていれば、左辺、すなわち時刻$t+h$における値$f(t+h)$は、$f(t) + f'(t)h$で近似できる、ということを主張する。このように、時間や空間を離散化し、微分を離散単位で近似することを 差分化 と呼ぶ。

先の運動方程式に差分化を適用すると、

$$ \begin{aligned} v_x(t+h) &= v_x(t) \\ v_y(t+h) &= v_y(t) - g h\\ r_x(t+h) &= r_x(t) + v_x(t) h \\ r_y(t+h) &= r_y(t) + v_y(t) h\\ \end{aligned} $$

となる。これを実装してみよう。

課題1:抵抗の無い弾道計算

以下のプログラムを三つのセルに分けて入力し、実行せよ。

import matplotlib.pyplot as plt
import numpy as np
from math import pi, cos, sin
def run(theta):
  rx, ry = 0.0, 0.0
  vx, vy = cos(theta), sin(theta)
  ax, ay = [], []
  g = 1.0
  h = 0.01
  while ry >= 0.0:
    rx += vx * h
    ry += vy * h
    # 残りの速度の更新部分を実装せよ
    ax.append(rx)
    ay.append(ry)
  nx = np.array(ax)
  ny = np.array(ay)
  return nx, ny
nx, ny = run(pi*0.25)
plt.plot(nx,ny, label="A")
plt.legend()
plt.show()

以下のような放物線が表示されたら成功である。

放物線

課題2:厳密解との比較

さて、先程の運動方程式は厳密に解くことができる。

$$ \begin{aligned} v_x(t) &= \cos{\theta} \\ v_y(t) &= \sin{\theta} - g t\\ r_x(t) &= t \cos{\theta} \\ r_y(t) &= t \sin{\theta} - \frac{g t^2}{2}\\ \end{aligned} $$

数値計算において厳密解がある場合は、必ず厳密解と比較すべきである。厳密解を返す関数を作ってみよう。runを入力したセルの次に新たにセルを作り、以下を入力せよ。

def exact(theta):
  rx, ry = 0.0, 0.0
  ax, ay = [], []
  t = 0.0
  g = 1.0
  h = 0.01
  while ry >= 0.0:
    t += h
    rx = t * cos(theta)
    ry = t * sin(theta) - g * t**2/2.0
    ax.append(rx)
    ay.append(ry)
  return np.array(ax), np.array(ay)

上記を実行したら、4つ目のセルを以下のように修正し、実行せよ。

nx, ny = run(pi*0.25)
plt.plot(nx,ny, label="Simulation")
ex, ey = exact(pi*0.25)
plt.plot(ex,ey, label="Exact")
plt.legend()
plt.show()

ほぼ重なった二つのグラフが表示されたら成功である。次に、runの時間刻みhの値を0.01から0.001にして、結果がどう変わるか比較せよ。

空気抵抗がある場合の弾道計算

さて、先程は空気抵抗が無い場合を考えた。次は空気抵抗がある場合を考えよう。一口に空気抵抗と言っても様々な要因があり、その扱いは難しいのだが、ここは単純に「速度に比例して向きと反対向きに『ひっぱられる』力」がかかっていると考えよう。仰角(地面となす角度)$\theta$で打ち上げて、しばらくしてその角度が$\phi$になった状況を考える。物体にかかるのは速度の向きに逆方向の抵抗であり、それを$x$軸、$y$軸に分解するとこんな感じになる。

物体にかかる力

この図をもとに、速度に関係する運動方程式を書き下すとこうなる。

$$ \begin{aligned} \dot{v_x} &= - \gamma v \cos(\phi) \\ \dot{v_y} &= - \gamma v \sin(\phi) - g \\ \end{aligned} $$

さて、三角関数の定義から

$$ \begin{aligned} \cos(\phi) &= \frac{v_x}{v}\\ \sin(\phi) &= \frac{v_y}{v}\\ \end{aligned} $$

である。これを代入すると、最終的に以下の運動方程式を得る。

$$ \begin{aligned} \dot{v_x} &= - \gamma v_x \\ \dot{v_y} &= - \gamma v_y - g \\ \dot{r_x} &= v_x \\ \dot{r_y} &= v_y \end{aligned} $$

これを実装してみよう。

課題3:空気抵抗のある運動方程式の差分化

現在、セルの状況は

  1. import
  2. run
  3. exact
  4. plot

となっているはずである。この2つ目と3つ目の間に新たにセルを作り、抵抗がある場合の時間発展を計算するコードrun_rを実装せよ。

def run_r(theta):
  rx, ry = 0.0, 0.0
  vx, vy = cos(theta), sin(theta)
  ax, ay = [], []
  g = 1.0
  h = 0.01
  gamma = 0.5
  while ry >= 0.0:
    rx += vx * h
    ry += vy * h
    # 残りの速度の更新部分を実装せよ
    ax.append(rx)
    ay.append(ry)
  nx = np.array(ax)
  ny = np.array(ay)
  return nx, ny

入力したら、5つ目のセルを以下のように入力し、実行せよ。

nx, ny = run(pi*0.25)
plt.plot(nx,ny, label="Without resistance")
nx2, ny2 = run_r(pi*0.25)
plt.plot(nx2,ny2, label="With resistance")
plt.legend()
plt.show()

空気抵抗の無い場合と比べて、どのような変化があっただろうか？

また、空気抵抗がある場合に、最も遠くまで飛ぶ角度はどのあたりか確認せよ。5つ目のセルを以下のように修正し、「A」と「B」のどちらが遠くまで飛ぶかを比較しながら、最も遠くまで飛ぶ角度を探してみよ。

nx, ny = run_r(pi*0.25)
plt.plot(nx,ny, label="A")
nx2, ny2 = run_r(pi*0.24)
plt.plot(nx2,ny2, label="B")
plt.legend()
plt.show()

任意課題：空気抵抗がある場合の厳密解

一般に空気抵抗などがあると運動方程式の厳密解を求めるのは困難であるが、上記の場合は厳密解を求めることができる。

まず、速度に関する運動方程式は以下で与えられていた。

$$ \begin{aligned} \dot{v_x} &= - \gamma v_x \\ \dot{v_y} &= - \gamma v_y - g \end{aligned} $$

これはそのまま求積することができる。$v_x(0) = \cos \theta, v_y(0) = \sin \theta$であることから、解は以下のように求まる。

$$ \begin{aligned} v_x(t) &= \mathrm{e}^{-\gamma t} \cos{\theta} \\ v_y(t) &= \mathrm{e}^{-\gamma t} \sin{\theta} - \frac{g}{\gamma}(1 - \mathrm{e}^{-\gamma t} ) \end{aligned} $$

このような式を得たら、必ず$t=0$や$t\rightarrow \infty$の場合にどうなるかを確認すること。特に、十分時間が経過した後は$v_x \rightarrow 0, v_y \rightarrow -g/\gamma$という一定値に収束することがわかる。重力下で抵抗がある場合、最終的に速度が一定に落ち着くが、その落ち着いた速度を 終端速度 (terminal verlocity) と呼ぶ。

さて、速度が求まったので、これを積分すれば任意の時刻による位置も求めることができる。

$$ \begin{aligned} r_x(t) &= \frac{(1 - \mathrm{e}^{-\gamma t})\cos{\theta} }{\gamma} \\ r_y(t) &= \frac{(1 - \mathrm{e}^{-\gamma t})\sin{\theta} }{\gamma} + \frac{g (1 - \mathrm{e}^{-\gamma t} - \gamma t)}{\gamma^2} \end{aligned} $$

これにより、空気抵抗がある場合のシミュレーションと厳密解を比較することができる。空気抵抗がある場合の厳密解exact_rを適当なセルに入力して実行せよ。

def exact_r(theta):
  ax, ay = [], []
  x, y = 0, 0
  g = 1.0
  dt = 0.01
  t = 0.0
  gamma = 0.5
  while y >= 0.0:
    t += dt
    x =  (1.0 - exp(-gamma*t))*cos(theta)/gamma
    y =  (1.0 - exp(-gamma*t))*sin(theta)/gamma
    y += g*(1 - exp(-gamma *t) - gamma *t )/(gamma**2)
    ax.append(x)
    ay.append(y)
  nx = np.array(ax)
  ny = np.array(ay)
  return nx, ny

厳密解と数値計算を比較してみよう。

nx, ny = run_r(pi*0.25)
plt.plot(nx, ny, label="Simulation")
ex, ey = exact_r(pi*0.25)
plt.plot(ex, ey, label="Exact")
plt.legend()
plt.show()

誤差の範囲で一致すれば成功である。ずれが気になる場合は、run_rの時間刻みを細かく(例えばh=0.001)にしてみよ。逆に粗く(h = 0.1)したらどれくらいずれるだろうか。

反応拡散方程式(グレイ・スコット模型)

先ほどまでは弾道計算を行った。飛翔体の計算は奥が深いのだが、結果がやや地味なのは否めない。せっかくシミュレーションをするので、もう少し結果が派手な計算をしてみよう。そのような例として、反応拡散系(Diffusion-Reaction System)を取り上げる。

等モルの塩酸と水酸化ナトリウムを混ぜると、食塩水ができることは知っているであろう。この反応は一方的であり、混ざって食塩ができておしまいである。しかし、ある種類の化合物を混ぜて反応させると、ある物質ができたり消えたりを繰り返すことがある。最も有名な例はBZ反応(ベロウゾフ・ジャボチンスキー反応, Belousov-Zhabotinsky Reaction)であろう。これは、ある溶液を混ぜると、その後しばらく溶液の色が周期的に変動する現象である。非常に雑に説明すると、反応を記述する方程式が時間の二階微分方程式になり、振動解が出てくるのがこの現象の本質である。

さて、BZ反応は時間的に変動する現象であるが、これが「拡散(Diffusion)」と結びつくと、時間的な変動が空間的に伝播していく。これにより複雑な模様ができあがる。化学反応(Reaction)と拡散(Diffusion)が組み合わさった現象であるから反応拡散系(Diffusion-Reaction System)と呼ばれる。反応拡散系は様々な例が知られているが、筆者の趣味でグレイ・スコットモデル(Gray-Scott model)を取り上げる。

グレイ・スコットモデルは、以下のような連立偏微分方程式で記述される。

$$ \begin{aligned} \frac{\partial u}{\partial t} &= D_u \Delta u - uv^2 + F(1-u) \\ \frac{\partial v}{\partial t} &= D_v \Delta v + uv^2 - (F+k)v \end{aligned} $$

右辺の第一項が拡散項、その後ろにあるのが反応を記述する力学系である。三次元を考えることもできるが、ここでは二次元空間を考える。 $\Delta$はラプラシアンと呼ばれる二階微分演算子で、二次元なら以下で定義される。

$$ \Delta \equiv \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} $$

次に空間を離散化しよう。今は二次元空間を考えているが、まずは一次元の世界$f(x)$を考える。微分方程式に二階微分が含まれているので、二階微分を近似したい。そのために、$f(x+h)$と$f(x-h)$をそれぞれ二次までテイラー展開してみよう。

$$ \begin{aligned} f(x + h) &= f(x) + f'(x)h + \frac{h^2}{2} f''(x) + O(h^3) \\ f(x - h) &= f(x) - f'(x)h + \frac{h^2}{2} f''(x) + O(h^3) \end{aligned} $$

両辺足すと$f'(x)$の項が消えるので、整理して

$$ f''(x) = \frac{f(x+h) - 2 f(x) + f(x-h)}{2h} $$

を得る。これを 中央差分 と呼ぶ。

TODO: 差分から見たラプラシアンの物理的意味

全く同様にして、二変数関数の二階微分(ラプラシアン)は、以下のように表現できる。

$$ \Delta f(x,y) \sim \frac{f(x+h,y)+f(x-h,y)+f(x,y+h)+f(x,y-h)-4f(x,y)}{4h} $$

先ほどの微分方程式は、時間に関して一階微分、空間に関して二階微分であったので、これで必要な式がそろった。他の項は「現時刻のその場所の値」で計算可能なので、これで次のステップの値が求まる。

差分計算の流れ

課題4: ラプラシアンの実装

では実際にグレイ・スコットモデルのシミュレーションを実装してみよう。新しいノートブックを開き、最初のセルにimport文を描こう。

import matplotlib.pyplot as plt
import numpy as np

次のセルに、二次元配列sの、座標(ix,iy)におけるラプラシアンの値を返す関数laplacianを実装せよ。ただし、空間差分の単位$h$は1とする。

def laplacian(ix, iy, s):
    ts = 0.0
    # ここを実装せよ
    return ts

実装したらテストしてみよう。三つ目のセルにlaplacianに食わせるテスト用のNumPy二次元配列を作成する。

a = np.arange(9).reshape(3,3)
a

実行すると以下のような表示になるはずである。

array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])

さて、中央に着目しよう。中央の値は4であり、1,3,5,7に囲まれている。周りの4つの数字の平均は4であり、中央の値と一致するため、この地点でのラプラシアンの値は0になるはずである。確認してみよう。三つ目のセルを以下のように書き換えよ。

a = np.arange(9).reshape(3,3)
laplacian(1,1,a) # => 0.0

これは、aという配列の(1,1)地点におけるラプラシアンの中央差分の値を求めよ、という意味である。実行結果として0.0が返ってくれば正しい。

これだけではテストとして不安なので、値を少しずらしてみよう。三つ目のテスト用セルを以下のように修正して実行せよ。

a = np.arange(9).reshape(3,3)
a[0,1] = 0
a

まず、入力するリストが以下のような形になる。

[[0 0 2]
 [3 4 5]
 [6 7 8]]

4の上が1であったのが0になった。したがって、中央は平均よりも高い値になっているから、ラプラシアンは負になるはずである。実際に計算してみよう。周りの数字の合計は15であり、その平均は3.75である。中央の値との差は0.25であるから、その差は-0.25になるはずだ。計算してみよう。

a = np.arange(9).reshape(3,3)
a[0,1] = 0
laplacian(1,1,a) # => 0.25

正しく計算されただろうか？これも正しく計算できていれば、ラプラシアンは正しく実装できているであろう。三つ目のセルは消去してよい。

シミュレーションコードの実装

三つ目、四つ目、のセルに以下を入力せよ。

def calc(u, v, u2, v2):
    (L, _) = u.shape
    dt = 0.2
    F = 0.04
    k = 0.06075
    lu = np.zeros((L, L))
    lv = np.zeros((L, L))
    for ix in range(1, L-1):
        for iy in range(1, L-1):
            lu[ix, iy] = 0.4 * laplacian(ix, iy, u)
            lv[ix, iy] = 0.2 * laplacian(ix, iy, v)
    cu = -v*v*u + F*(1.0 - u)
    cv = v*v*u - (F+k)*v
    u2[:] = u + (lu+cu) * dt
    v2[:] = v + (lv+cv) * dt
def simulation():
    L = 32
    u = np.zeros((L, L))
    u2 = np.zeros((L, L))
    v = np.zeros((L, L))
    v2 = np.zeros((L, L))
    h = L//2
    u[h-6:h+6, h-6:h+6] = 0.9
    v[h-3:h+3, h-3:h+3] = 0.7
    for i in range(1000):
        if i % 2 == 0:
            calc(u, v, u2, v2)
        else:
            calc(u2, v2, u, v)
    return v

ここまで入力したら、セルの状況は以下のようになっているはずである。

  1. import
  2. laplacian
  3. calc
  4. simulation

その状態で、5番目のセルに以下を入力して実行してみよう。

v = simulation()
plt.imshow(v)

シミュレーション結果

こんな白黒の図が現れたら成功である。

JITによる高速化

これは32x32のグリッドを1000ステップ計算するシミュレーションであるが、実行に少し待たされた気がするであろう。まず時間を測定してみよう。5つ目のセルを以下のようにして実行してみよ。

%time v = simulation()
plt.imshow(v)

おそらくこんな結果が表示されたものと思う。

CPU times: user 3.25 s, sys: 1.11 ms, total: 3.25 s
Wall time: 3.26 s

これは、計算に3秒ちょっとかかったよ、という意味だ。32x32のグリッドを2二種類、1000ステップなので、計算量としては、おおざっぱに100万ちょっと、つまり「メガ」のオーダーである。今、目の前にある計算機の処理能力はだいたい「ギガ」、つまり一秒間にメガの千倍の計算を処理する能力があるはずで、メガの量の計算に秒単位でかかるのは遅い気がしてくるであろう。そこでJITで高速化しよう。

TODO: JITの説明

まず、一つ目のセルを以下のように修正せよ。

import matplotlib.pyplot as plt
import numpy as np
from numba import jit # この行を追加

次に、laplacianやcalc、simulationの関数定義の前に@jitと記述せよ。

@jit
def laplacian(ix, iy, s):
@jit
def calc(u, v, u2, v2):
@jit
def simulation():

この状態で、全てのセルを実行してから、5つめのセルを実行しよう。面倒なら「ランタイム」から「再起動してすべてのセルを実行」を選んでもよい。おそらく20～30ms程度の実行時間になったはずである。

模様の確認

さて、それなりの速度になったところで、もう少し大きな系を長時間計算してみよう。まず、システムサイズを64x64にして、時間ステップも10倍の1000にしよう。具体的には、4つ目のセルのsimulationを以下のように書き換えよ。

@jit
def simulation():
    L = 64 # ここを32から64に
    u = np.zeros((L, L))
    u2 = np.zeros((L, L))
    v = np.zeros((L, L))
    v2 = np.zeros((L, L))
    h = L//2
    u[h-6:h+6, h-6:h+6] = 0.9
    v[h-3:h+3, h-3:h+3] = 0.7
    for i in range(10000): # ここを1000から10000に
        if i % 2 == 0:
            calc(u, v, u2, v2)
        else:
            calc(u2, v2, u, v)
    return v

その上で再度5つ目のセルを実行してみよ。不思議な模様が浮かび上がれば成功である。

せっかく模様が出てきたのに、モノクロというのも味気ない。色を付けてみよう。5つ目のプロットのところで、カラーマップを指定してみよう。

%time v = simulation()
plt.imshow(v, cmap="inferno") # ここを修正

カラーマップには他にもviridisや、plasma、magmaなどが用意されているので、興味があれば試してみよ。

任意課題：アニメーション

せっかくの時間発展シミュレーションなので、時間発展の様子をアニメーションで見てみよう。

まず、一つ目のセルに以下を追加せよ。

import matplotlib.pyplot as plt
import numpy as np
from numba import jit
from matplotlib import animation, rc # この行を追加

次に、simulationで、最後の結果だけではなく、途中の経過もまとめてリストで返すようにしよう。

@jit
def simulation():
    L = 64
    u = np.zeros((L, L))
    u2 = np.zeros((L, L))
    v = np.zeros((L, L))
    v2 = np.zeros((L, L))
    h = L//2
    u[h-6:h+6, h-6:h+6] = 0.9
    v[h-3:h+3, h-3:h+3] = 0.7
    r = [] # この行を追加
    for i in range(10000):
        if i % 2 == 0:
            calc(u, v, u2, v2)
        else:
            calc(u2, v2, u, v)
        if i % 100 == 0: # このif文を追加
          r.append(u.copy())
    return r # uでなくrを返すよう修正

その上で、6つ目のセルに以下を入力、実行せよ。

imgs = simulation()
fig = plt.figure()
def update(img):
  plt.cla()
  plt.imshow(img,cmap="inferno")

ani = animation.FuncAnimation(fig, update, frames=imgs,interval=50)
rc('animation', html='jshtml')
ani

アニメーションが表示されたら成功である。

TODO: 左下に余計な図が出る原因を調べること。

余談：パーソナルスーパーコンピュータ

パソコンとは「パーソナルコンピュータ」の略、つまり「個人向け計算機」という意味だ。もともと計算機は貴重品かつ大型であり、組織に一つしかないものだった。それが徐々に小型化し、オフィスに一つ(オフコン)になり、さらに個人で独占して利用できるものになった。パソコンが普及するにつれて、もともと「組織に一つ」しかないような巨大な計算機は「スーパーコンピュータ(スパコン)」と呼ばれ、パソコンと区別されるようになった。スパコンは、安くて一億、高ければ数十、数百億円といったその価格もさることながら、その維持も大変である。計算するのには莫大な電気が必要で、かつ使った電気はすべて熱となるからそれを冷却するシステムも必要である。したがって、本来「パソコン」と「スパコン」は相容れない概念のはずだが、スパコンを個人で所有することでその二つを悪魔合体させ、「パーソナルスーパーコンピュータ」という狂った概念を生み出した人がいる。桑原邦郎氏である。彼は流体力学を専門とする研究者で、親から受け継いだ莫大な財産をすべてスパコンに突っ込んだ。東京大学工学部物理工学科の助手、宇宙科学研究所の助教授を経て、自宅に計算流体力学研究所という研究所を作り、そこにスパコンを購入して運用した。最盛期は七台のスパコンがフル稼働し、電気代だけで月に2000万円かかったという。1980年代後半から1990年代にかけて、計算流体力学を専門とする人はほとんど彼のパーソナルスパコンにお世話になったと思われる。自動車メーカも技術者を派遣していたそうだ。また、米国の諜報機関が「軍事目的に使っているのではないか」と疑ったとのエピソードもある。彼は親から受け継いだ莫大な遺産をすべてスパコンに突っ込み、それを惜しげもなくいろんな人に使わせた。

それから紆余曲折あって、計算流体力学研究所はスパコンを手放し、技術コンサルやパソコンの組み立て、販売をする会社となった。私が大学院に進学した際に与えられたパソコンは、この計算流体力学研究所で購入したAlpha21164のマシンであった。指導教員の「せっかくだから組み立てさせてもらったら？」の言葉に甘え、目黒に行ってアルバイトのお兄さんと一緒に自分の研究に使うマシンを組み立てた。そこに社長である桑原氏が様子を見にやってきた。僕が物工の学生と知ると、興味をもっていろいろ聞いてきた。僕はまさか目の前の社長さんが元物工の助手だったなんて知らなかったので「物工のご出身なんですか？どこの研究室ですか？」と的外れな質問をした。彼はただ笑って何も答えなかったのを思い出す。その時は青二才で何もわからなかった私だが、スパコンを使って研究をするようになった今なら、彼からいろいろ興味深い話が聞けたのではないかと残念に思う。桑原氏は2008年、その豪快な生涯を閉じた。「親の遺産をもっとも有効に学術に活かした」と評されている。