テーマ2：反応拡散方程式

目的

化学反応と拡散過程の組み合わせにより、不思議な模様が出現することを見る。

プログラム

先程組んだ拡散方程式とプログラムが重複するので、コピーしてしまおう。もしまだ名前をつけていない場合は、左上の「Untitled」となっているところをクリックして、出現した「Rename Notebook」の「Enter a new notebook name」に「diffusion」と入力して「Rename」ボタンを押す。

次に、現在のプログラムにdiffusionという名前がついている状態で、左上の「File」メニューから「Make a Copy...」を選ぶ。すると、「diffusion-Copy1」という名前のノートブックが作成されるので、名前を「gray-scott」に変えること。そして「Kernel」メニューの「Restart & Clear Output」を実行する。確認ダイアログが出るが、「Restart and Clear All Outputs」をクリックする。すると、全ての結果が消え、セルも未実行の状態となる。

最初のセル

最初のセルは変更しなくて良い。実行だけしておくこと。

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

2番目のセル

二番目のラプラシアンを定義する関数もそのままで良い。実行だけしておくこと。

@jit
def laplacian(ix, iy, s):
    ts = 0.0
    ts += s[ix-1, iy]
    ts += s[ix+1, iy]
    ts += s[ix, iy-1]
    ts += s[ix, iy+1]
    ts -= 4.0*s[ix, iy]
    return ts

3番目のセル

3番目のセルは、一回呼び出すたびに、２つの化学物質uとvの濃度分布を1ステップ変化させる関数calcである。これを以下のように書き換えよ。これは大きく変えるので、全部消して書き直しても良いし、変更点を自分で探して修正しても良い。

@jit
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.1 * laplacian(ix, iy, u)
            lv[ix, iy] = 0.05 * 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

4番目のセル

4番目のセルは、1ステップ時間をすすめる関数calcを何度も呼んで、時間発展を記述する関数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:
            r.append(u.copy())
    return r

先程と同様に写真を撮っているが、100ステップごとに100枚、合計10000ステップの計算を行う。rは写真の配列であり、100枚写真が撮られている。

5番目のセル

5番目のセルで、ただしく計算ができているかの確認をする。三行目で、カラースキーム(値を色に変換する方法)を変更する。vminとvmaxの指定も消しておこう。

imgs = simulation()
fig = plt.figure()
im = plt.imshow(imgs[-1], cmap="GnBu")

上記を実行して、なにか模様が出てくれば成功である。

6番目のセル

得られた100の写真を使ってアニメーションをする6番目のセルは変更しなくて良い。そのまま実行せよ。

def update(i):
  im.set_array(imgs[i])
  return im,

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

しばらくするとコントロールパネルが表示されるので「▶」をクリックせよ。ここまで正しく入力できていれば、アニメーションが表示されるはずである。

より大きい系の計算

せっかくなので、もう少し計算領域を大きくしてみよう。

まず、4番目のセルで、以下の変更をしよう。

システムサイズLを64から128に
計算ステップ数を10000から30000に
写真撮影のタイミングを100ステップから300ステップに

変更した結果、以下のようなコードになるはずである。

@jit
def simulation():
    L = 128  # 64から128に変更
    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(30000): # 10000から30000に
        if i % 2 == 0:
            calc(u, v, u2, v2)
        else:
            calc(u2, v2, u, v)
        if i % 300 == 0:   # 100から300に
            r.append(u.copy())
    return r

4番目、5番目、6番目と順番にセルを実行することで、大きなサイズの計算ができるはずである。

解説

赤い水と青い水を混ぜると紫色の水になる。等モルの塩酸と水酸化ナトリウムを混ぜると、食塩水になる。このように「なにかを混ぜると、一様にまざった何かになる」ことがほとんどである。しかし、うまく化学反応を作ると、濃度がいったりきたりするような系を作ることができる。その例がベロウソフ・ジャボチンスキー反応、通称BZ反応である。BZ反応では、二つの溶液をビーカーに入れてよく混ぜると色が周期的に変化する。さらに面白いのは、二つの溶液をシャーレのような薄い容器に入れて反応させると、濃度が振動しながら拡散することにより螺旋状の不思議な模様が成長していくことだ。

こうした、化学反応と拡散の組み合わせで不思議な模様ができることを示したのはイギリスの天才科学者、アラン・チューリングである。チューリングは生物の模様がこの方程式で説明できるのではないかと考えていたようだ。この説は長らく顧みられることはなかったが、後に熱帯魚の一種の模様がこの機構により作られていることが実験的に確認された。

チューリングによるチューリング・パターンの論文発表は1952年のことだが、それよりかなり前の1933年の日本で、物理学者によるある仮説が議論を呼んだ。その物理学者はキリンの模様に着目した。このまだら模様は泥や粘土が乾く時に見られる割れ目に似ている。実は両者に共通する物理現象が関与しているのではないか？例えばキリンの成長速度に皮膚の成長がついていけず、それでひび割れが起きてあの模様になるのではないか、という説である。これに生物学者が反発する。後に「キリンの斑論争」と呼ばれる論争は、両者が感情的になってしまい、あまり実りある議論にならなかったようだ。実際にはキリンの模様は「ひび割れ」ではなかったようだが、果物のメロンの模様は「ひび割れ」機構で説明できるようである。

物理では、まったく異なる現象でも、同じ微分方程式で記述される系は同じとみなすことができる。たとえばバネにおもりをつけて上下に揺れている系と、コイルとコンデンサで作った回路に流れる電流の振動現象は、片方は力学的、もう一方は電磁気的と現象を引き起こすエネルギーはまったく異なるが、実は同じ微分方程式に従うため、同じ振る舞いをする。こうして、「全く異なる現象に共通する性質を見出していく」のも物理の重要な役割である。

余談だが、チューリング・パターンの提唱者、アラン・チューリングは、数学、論理学、コンピュータ科学など多方面において大きな業績を残した天才である。特に戦争中に、それまで解読不能と言われていたドイツの「エニグマ」という暗号を解読したことで有名である(サイモン・シンの「暗号解読」という本に詳しいので参照されたい)。他にも、「人間であるかAIであるかのテスト」である「チューリング・テスト」の提唱でも有名である。今でこそAIが人間を超える「シンギュラリティ」が話題となっているが、チューリング・テストの発表は1950年である。また、チェスの思考ルーチンも考えるなど、現在では「人工知能の父」とも呼ばれている。彼は同性愛者であったが、当時のイギリスでは同性愛は違法であった。それが直接の原因かはわからないが、その数年後に彼は自宅で青酸カリによって亡くなった。自殺と言われている。その後、彼の業績は再評価され、2019年には彼の業績をたたえて新50ポンド紙幣に採用されることとなった。