Skip to content

Latest commit

 

History

History
525 lines (367 loc) · 29.6 KB

README.md

File metadata and controls

525 lines (367 loc) · 29.6 KB
Raw

数値シミュレーション

本稿で学ぶこと

  • 差分化
  • 弾道計算
  • 反応拡散方程式

数値シミュレーションとは

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

さて、この世界の空間や時間は(おそらく)連続的であるのに対し、コンピュータは原則として離散的な値しか扱うことができない。数値シミュレーションとはこの世界の出来事をコンピュータの中に再現することであるから、計算するにあたって連続的な値を離散的な値に変換する必要がある。例えば、パソコンやスマホで動画を見ることがあるだろう。動画は、我々が目にする世界のように、空間的にも時間的にも連続的に見えるが、実際には時間方向は静止画像を高速にコマ送りすることで連続的に見せており、静止画像も拡大するとピクセル単位で離散化されている。

時間と空間の離散化

このように、連続的な値をコンピュータで扱うために離散的な値にすることを 離散化(discretization) と呼ぶ。以下では、時間や空間を離散化することで微分方程式を計算機が扱いやすい形にして、その振る舞いを数値シミュレーションで解析してみよう。

差分化

今、ある量$f(t)$の時間微分$f'(t)$が与えられているとしよう。時間に関する離散化とは、ある小さな時間刻み$h$に対して、$f(t)$の値からなんとかして$f(t+h)$の値を推定することである。時間微分$f'(t)$がわかっているので、厳密な表式は積分で与えられる。

$$ f(t+h) = f(t) + \int_t^{t+h} f'(t) dt $$

この式の意味は時刻$t$において$f(t)$の値である時、それに$t$から$t+h$までの時間変化$f'(t)$をすべて積算したものを加えると$f(t+h)$の値になります、ということであり、なんら難しいことは言っていない。

さて、一般にはこの積分を求積することはできないので、なんらかの近似をする。最も単純な近似は、時間刻み$h$が小さいので、その間は$f'(t)$が時間変化しないものとみなすことだろう。すると、$f'(t)$を積分の外に出すことができるので、

$$ f(t+h) \sim f(t) + h f'(t) $$

と近似することができる。$f(t+h)$を$t$の周りでテイラー展開して、1次までとると、

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

となることから、先程の近似は$h$の1次まで正しい。これを1次近似と呼ぶ。この式は、右辺、すなわち時刻$t$における$f(t)$の値と、その微係数$f'(t)$がわかっていれば、左辺、すなわち時刻$t+h$における値$f(t+h)$は、$f(t) + f'(t)h$で近似できる。

時間の差分化

以上から、ある時刻$t_0$における値$f(t_0)$がわかっているとき、

$$ \begin{aligned} f(t_0+h) &= f(t_0) + f'(t_0)h \\ f(t_0+2h) &= f(t_0+h) + f'(t_0+h)h \\ f(t_0+3h) &= f(t_0+2h) + f'(t_0+2h)h \\ & \cdots \end{aligned} $$

と、代入を繰り返すことで任意の時刻の値を求めることができる。このように、微分を離散単位で近似することを差分化と呼び、差分化により微分方程式を扱う手法を 差分法(finite difference method) と呼ぶ。特に、先程のように1次近似の差分法を オイラー法 (Euler's method) と呼ぶ。

ニュートンの運動方程式

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

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

となる。さらに速度とは、単位時間当たりの位置の変化分であった。それもあわせて書くと、以下の式が得られる。

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

すなわち、運動方程式とは時間に関する連立常微分方程式である。ここで、2つの物体の間に距離の2乗に反比例する力を仮定すると、宇宙空間内で互いに重力で引き合う2つの星の運動を表す方程式となる。これを解くと、軌道が楕円になることや、面積速度一定則、調和の法則など、いわゆるケプラーの三法則が導かれる。 自然科学において法則とは、実験や観測などで再現可能な自然の振る舞いである。物理学とは、より基本的な原理から、観測事実である「法則」を導き、説明しようとする営みである。以下では、運動方程式の簡単な場合として、重力下での自由運動、すなわち弾道計算を行ってみよう(課題1)。

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

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

さて、いま皆さんの目の前にあるのは、ちょっと前のスーパーコンピュータなみの計算能力を持つ計算機である。それを使って、簡単なシミュレーションをしてみよう。2次元の場合を考えよう。速度ベクトルを$\vec{v} = (v_x, v_y)$、位置ベクトルを$\vec{r} = (r_x, r_y)$とし、運動方程式をそれぞれの要素について書き下すと以下のようになる。

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

さて、この式は厳密に解けるのだが、その厳密解を知らないものとし、シミュレーションで近似的に解を求めることにする。求めたいものは重力下で角度$\theta$で物体を投げた時の物体の軌道である。まず、数値計算で扱いやすいように、時間を離散化しよう。

先の運動方程式に1次の差分化(オイラー法)を適用すると、

$$ \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} $$

となる。この計算を一度行うことで、時刻$t$の物理量から時刻$t+h$の物理量が得られる。あとはこれを繰り返せば(時間刻み$h$の精度で)任意の時刻の位置と速度がわかることになる。

厳密解も求めておこう。まず、運動方程式から速度に関する部分を抜き出すと以下のようになる。

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

これは簡単に求積できる。初期条件として、仰角$\theta$、速度$v_0$で投げたとすると、$v_x(0) = v_0 \cos \theta, v_y(0) = v_0 \sin \theta$であるから、

$$ \begin{aligned} v_x(t) &= v_0 \cos \theta \\ v_y(t) &= -gt + v_0 \sin \theta \end{aligned} $$

となる。これをさらに時間積分したものが座標であるから、初期位置が原点、すなわち$(r_x, r_y) = (0,0)$ であったとすると、

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

これが求めたい軌道であった。着弾するまでの時間は$r_y(t_f) = 0$となる$t_f$であるから、

$$ t_f = 0, \frac{2 v_0 \sin \theta}{g} $$

である。$x$方向の速度はずっと$v_0 \cos \theta$であるから、着弾までに飛んだ水平距離$l(\theta)$は

$$ l(\theta) = v_0 t_f \cos \theta = \frac{2 v_0^2 \sin \theta \cos \theta}{g} $$

これを最大にする角度は、$l'(\theta)=0$を満たす$\theta$であり、$v_0$や$g$の値によらず$\theta = \pi/2$、すなわち45度であることがわかる。

高低差がある場合の弾道計算

高低差がある場合の弾道

さて、先程は、水平な地面から物を投げて、そのまま地面にぶつかるまでの距離を計算した。この場合、最も遠くに飛ぶ角度は45度である。では、台の上に乗って物を投げた場合はどうだろうか？なんとなく一番遠くに飛ぶ角度が45度からずれることは想像がつくであろう。では、上向きにずれるだろうか？それとも下向きにずれるだろうか？

いま、初速$v_0$、仰角$\theta$で投げるのは同じとして、最初に高さ$y_0$から物体を投げたとしよう。物体が従う運動方程式は地面から投げた場合と同じだが、初期条件が異なる。まず、速度については解は変わらない。

$$ \begin{aligned} v_x(t) &= v_0 \cos \theta \\ v_y(t) &= -gt + v_0 \sin \theta \end{aligned} $$

しかし、位置については、$r_y$に$y_0$が加わる。

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

違いはたったこれだけだが、急に式が面倒になる。まず、着弾までの時刻$t_f$は$r_y(t_f)=0$の解であるから、

$$ t_f = \frac{v_0 \sin \theta + \sqrt{v_0^2 \sin^2 \theta + 2 g y_0}}{g} $$

ただし、複号は正の解をとった。着弾までに進む水平距離は、

$$ l(\theta) = v_0 t_f \cos\theta $$

これを$\theta$で微分して$l'(\theta)=0$となる$\theta$を探せばよいが、実際に計算すると、解が初等的には求まらないことがわかるであろう。これを数値計算で求めてみよう。

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

先程はニュートンの運動方程式を数値的に解くことで弾道計算を行った。この時、時間のみ離散化し、空間は連続のままだった。次は時間と同時に空間も離散化してみよう。そんな系の例として 反応拡散系(diffusion-reaction system) を取り上げる(課題2)。

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

BZ反応

さて、BZ反応は時間的に変動する現象であるが、これが拡散と結びつくと、時間的な変動が空間的に伝播していく。これにより複雑な模様ができあがる。化学反応(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} $$

右辺の第1項が拡散項、その後ろにあるのが反応を記述する力学系である。3次元を考えることもできるが、ここでは2次元空間を考える。

Gray-Scott模型

まず、空間を差分化して扱うことにしよう。本来連続的な空間を、グリッドに分割することで離散的な表現に落とす。

$\Delta$はラプラシアンと呼ばれる2階微分演算子で、2次元なら以下で定義される。

$$ \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)}{h^2} $$

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

$$ \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)}{h^2} $$

簡単のため、$h=1$としよう。空間の離散化により$f(m, n)$の値を2次元配列s[m][n]で表現する。すると、位置$(m,n)$におけるラプラシアンは、この配列sを用いて

s[m+1][n] + s[m-1][n] + s[m][n+1] + s[m][n-1] - 4 * s[m][n]

と表現できる。これは、前後左右のセルの値の合計から、自分の値の4倍を引いたものであり、自分が「まわりの平均」よりも大きい時に負、小さい時に正となる。このように、ラプラシアンは「なるべくまわりに合わせよう、全体を平均化しよう」という働きをする。

差分計算の流れ

もともとの式のうち、時間微分は1次の差分を取ることで計算する。その時間発展にはラプラシアンという空間微分が含まれているが、それについては中央差分で近似した。これで数値計算に必要な式がすべて揃った。

数値シミュレーション:課題

課題1:運動方程式

まずは物を投げた場合の運動を、運動方程式を解くことで数値的に追ってみよう。新しいノートブックを開き、ballistic.ipynbという名前で保存せよ。

1. ライブラリのインポート

最初のセルで、必要なライブラリのインポートをしよう。

import matplotlib.pyplot as plt
from math import pi, cos, sin
from pylab import rcParams

2. プロットサイズの変更

デフォルトサイズでは、ややグラフが小さいので、少しサイズを大きくしておこう。

rcParams['figure.figsize'] = 10, 5

3. 運動方程式の数値解法

時間方向について差分化した運動方程式を数値的に解くルーチン、throwを実装しよう。

def throw(theta, y0):
    rx, ry = 0.0, y0
    vx, vy = cos(theta), sin(theta)
    ax, ay = [], []
    g = 1.0
    h = 0.001
    while ry >= 0.0:
        rx += vx * h
        ry += vy * h
        vy -= g * h
        ax.append(rx)
        ay.append(ry)
    return ax, ay

ここで、thetaは投げる角度(仰角)、y0は最初に立っている「台」の高さである。

4. プロット

投げる角度と台の高さを受け取って、軌道をグラフにする関数plotを実装しよう。比較のため、複数の角度をリストで受け取り、それぞれの角度で投げた場合の軌道を重ねて表示している。

def plot(angles, y0=0.0):
    for theta in angles:
        nx, ny = throw(theta / 180.0 * pi, y0)
        plt.plot(nx, ny, label=theta)
    plt.legend()
    plt.show()

5. 地面から投げた場合

まず、地面から投げた場合(y0=0の場合)を見てみよう。40度、45度、50度で投げてみた軌道を重ねて表示する。

angles = [40, 45, 50]
plot(angles)

予想通り、45度がもっとも遠くに届いているだろうか？

6. 台の上から投げた場合

次に、台の上から投げてみよう。最初に高さ$1$の台に乗ってものを投げた場合にどうなるか調べてみる。遠くに届く角度が45度から増えるか減るか、予想してから実行すること。

angles = [40, 45, 50]
plot(angles, 1.0)

実際に一番遠くに届く角度は45度から増えただろうか、減っただろうか。入力する角度を5度刻みで変化させ、最も遠くに届く角度を調べよ。

発展課題：初速依存性

先程の運動方程式では、初速を$1$に固定していた。関数throwの初速を与える箇所を、以下のように変更すると、初速を$2$に変えることができる。

    vx, vy = 2.0*cos(theta), 2.0*sin(theta)

初速を大きくした場合、「台から投げた時」の「最も遠くに届く角度」は増えるだろうか、減るだろうか？初速を$10$にしたらどうなるだろう？初速が大きい極限でどうなるか考察せよ。

課題2:反応拡散方程式

反応拡散方程式を数値的に解いてみよう。新しいノートブックを開き、gs.ipynbという名前で保存せよ。

1. ライブラリのインポート

最初のセルで、必要なライブラリのインポートをしよう。

import matplotlib.pyplot as plt
import numpy as np
from numba import jit
from matplotlib import animation, rc

2. ラプラシアンの実装

2次元配列s[m][n]の、位置$(m,n)$におけるラプラシアンは、

s[m+1][n] + s[m-1][n] + s[m][n+1] + s[m][n-1] - 4 * s[m][n]

と表現できる。これを用いて、2次元配列sと位置m,nを受け取って、その位置でのラプラシアンの値を返す関数laplacianを実装したい。以下はそれを途中まで実装したものだ。

@jit
def laplacian(m, n, s):
    ts = 0.0
    ts += s[m+1][n]
    ts += s[m-1][n]
    ts += s[m][n+1]
    ts += s[m][n-1]
    # ここを埋めよ
    return ts

上記がラプラシアンの値となるように、未完成の部分を完成させよ。

なお、関数定義の前に@jitとあるのは、「この関数をJITコンパイルせよ」という指示(デコレータ)である。JITはJust in Timeの略で、実行時にコードをコンパイルすることでコードの実行を加速する。ここではJIT、デコレータともに詳細には触れない。

3. ラプラシアンのテスト

laplacianを実装したらテストしてみよう。以下はすべて三つ目のセルで実行せよ。

まず、laplacianに食わせるテスト用に、NumPyの2次元配列を作成する。

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)

これは、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になった。したがって、中央は平均よりも高い値になっているから、ラプラシアンは負になるはずである。実際に計算してみよう。

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

上記の実行結果が-1.0になれば正しく計算されている。

4. 時間発展

2つの配列を受け取り、1ステップだけ時間を進める関数calcを実装しよう。以下を4つ目のセルに入力せよ。

@jit
def calc(u, v, u2, v2):
    L, _ = u.shape
    dt = 0.2
    F = 0.04
    k = 0.06075
    Du = 0.1
    Dv = 0.05
    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] = Du * laplacian(ix, iy, u)
            lv[ix, iy] = Dv * 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

最初の@jitデコレータを忘れないこと。なお、関数の冒頭でL, _ = u.shapeとあるのは、u.shapeが2要素のタプルを返す際、そのうち最初の要素をLで受け取るが、2つ目は使わないことをアンダースコア「_」で表現している。宣言だけしてあとで使わない変数はバグの元なので、あとで使わないことがわかっている変数に対しては積極的に_を使って使わないことを明示するようにしたい。

5. シミュレーションループ

1ステップ時間を進める関数が書けたら、それを何度も呼び出すことで時間発展をさせよう。また、初期条件として模様の「種」を作る。5つ目のセルに以下を入力せよ。

@jit
def simulation(L, loop):
    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(loop):
        calc(u, v, u2, v2)
        u, u2, v, v2 = u2, u, v2, v
        if i % 100 == 0:
          r.append(v.copy())
    return r

これも、最初の行の@jitを忘れないこと。

6. シミュレーションの実行

imgs = simulation(64, 10000)
n = len(imgs)
for i in range(4):
    im = plt.imshow(imgs[n // 4 * i])
    plt.show()

ここまで正しく入力されていれば、不思議な模様が4枚現れたはずである。

7. アニメーション

せっかくシミュレーションしたので、アニメーションを表示させよう。そのための準備をする。

fig = plt.figure()
im = plt.imshow(imgs[-1])
def update(i):
    im.set_array(imgs[i])

実行後、不思議な模様が出力されるはずである。

8. アニメーションの表示

ではアニメーションを表示してみよう。以下を実行せよ(少し時間がかかる)。画像の下に操作パネルが出たら、再生ボタン(右向きの三角)を押してみよう。アニメーションが表示されたら成功である。

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

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

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

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