Диффузионные реакции (реакции с диффузией)

Реакции с диффузией очень широко используются в компьютерной графике и других областях. Однако, найти по ним информацию (особенно на русском языке) достаточно сложно. Все источники либо сосредотачиваются на химических аспектах (а они весьма замысловаты), либо освещают вопрос очень поверхностно.

Здесь я хочу привести работающий, законченный кусок кода, которым можно поиграться. Кроме того, реализация алгоритмов реакций с диффузией позволяет показать всё мощь numpy.

Немного теории

Мы будем рассматривать очень простую модель (самую простую). В ней есть несколько упрощений, делающие её не очень физичной, но сильно упрощающих все выкладки.

Итак.

• Мы будем рассматривать смесь двух веществ. Концентрации мы обозначим A и B.
• Для каждого вещества мы введём коэффициент диффузии.
• Введём правило: если встречаются три молекулы B-B-A, то происходит реакция и A превращается в B (получаем B-B-B).
• Отток A и приток B мы будем компенсировать искусственными слагаемыми в формулах.

Математически это можно записать так:

A' = (Da · ∆A - A·B² + F(1 - A)) · dt
B' = (Db · ∆B + A·B² + (K + F) · B) · dt

Эти выражения описывают, как изменятся A и B через время dt:

• Первое слагаемое: диффузия
• Второе: вероятность встречи B-B-A (не совсем строго, но приемлемо)
• Третье: искусственный отток и приток

Код

#!/usr/bin/python3

# ffmpeg -i foo%04d.png -c:v libx264 -r 30 -pix_fmt yuv420p out.mp4
# или со звуком:
# ffmpeg -ar 48000 -ac 2 -f s16le -i /dev/zero -i foo%04d.png -shortest -c:a aac -c:v libx264 -r 30 -pix_fmt yuv420p -strict experimental out.mp4

import scipy.ndimage as nd
import numpy as np
import matplotlib
import matplotlib.pyplot as plt

def main():
    size = 300
    field_a = np.ones([size, size], dtype=float)
    field_b = np.zeros([size, size], dtype=float)
    t = size / 2 - 10
    field_b[t:size-t,t:size-t] = 1
    D_a = 0.2
    D_b = 0.1
    F = 0.045 # feed
    K = 0.062 # kill
    dt = 1
    for step in range(300):
        print("step=%d" % step)
        for _ in range(50):
            t = field_a * field_b * field_b
            field_a += ( D_a * nd.filters.laplace(field_a) - t + F * (1.0 - field_a) ) * dt
            field_b += ( D_b * nd.filters.laplace(field_b) + t - (K + F) * field_b ) * dt
        plt.figure(figsize=(16, 9))
        plt.axis('equal')
        plt.contourf(field_b)
        plt.colorbar()
        plt.title("step=%d" % step)
        plt.savefig('foo%04d.png' % step, dpi=40)
        plt.clf()

if __name__ == '__main__':
    main()

Этот код создаёт набор картинок из которых можно легко собрать видео (см. комментарий в коде и результат ниже).

Обратите внимание, что чтобы реализовать тот алгоритм нам понадобилось бы четыре вложенных цикла (два — для обхода матриц с плотностями, два — для вычисления лапласиана). Numpy позволило не писать ни одного из них.

И конечно, numpy позволяет в десятки (а то и сотни) раз поднять производительность. Здесь мы создаём 300 кадров с шагом в 50 итераций между кадрами. Это 15000 итераций и они выполняются за считанные секунды!

Оттолкнувшись от этого кода, вы можете самостоятельно поиграться с разными значениями параметров. Коэффициенты диффузии не сильно влияют на картину. Они влияют на размер элементов узора. А вот K и F влияют очень сильно. Попробуйте чуть-чуть их поменять и вы увидите, как изменится картина.

Результат