Welcome, guest | Sign In | My Account | Store | Cart

Reaction-Diffusion Simulation using Gray-Scott Model.

Python, 93 lines
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93``` ```# Reaction-Diffusion Simulation Using Gray-Scott Model # https://en.wikipedia.org/wiki/Reaction-diffusion_system # FB - 20151016 import math import random from PIL import Image, ImageDraw imgx = 256; imgy = 256 # image size image = Image.new("RGB", (imgx, imgy)) draw = ImageDraw.Draw(image) pixels = image.load() steps = 1000 dA = 1.0 # Diffusion Rate A dB = 0.5 # Diffusion Rate B f = 0.0545 k = 0.062 h = 0.01 # time step size # 3x3 Convolution Matrix weights = [[0.05, 0.2, 0.05], [0.2, -1.0, 0.2], [0.05, 0.2, 0.05]] ##arA = [[[1.0 for x in range(imgx)] for y in range(imgy)] for z in range(2)] ##arB = [[[0.0 for x in range(imgx)] for y in range(imgy)] for z in range(2)] ##arB[0][(imgy - 1) / 2][(imgx - 1) / 2] = 1.0 # seed arA = [[[random.random() for x in range(imgx)] for y in range(imgy)] for z in range(2)] arB = [[[random.random() for x in range(imgx)] for y in range(imgy)] for z in range(2)] # simulation p = 0; print "%" + str(p).zfill(2) z = 1 for i in range(steps): z = 1 - z for iy in range(imgy): for ix in range(imgx): lapA = 0.0 lapB = 0.0 for jy in range(3): for jx in range(3): kx = ix + (jx - 1) ky = iy + (jy - 1) if kx < 0: kx = imgx - 1 if kx > imgx - 1: kx = 0 if ky < 0: ky = imgy - 1 if ky > imgy - 1: ky = 0 lapA += arA[z][ky][kx] * weights[jy][jx] lapB += arB[z][ky][kx] * weights[jy][jx] a = arA[z][iy][ix] b = arB[z][iy][ix] abb = a * b * b an = a + (dA * lapA - abb + f * (1.0 - a)) * h bn = b + (dB * lapB + abb - (k + f) * b) * h # prevent negative chemical amounts if an < 0: an = 0.0 if bn < 0: bn = 0.0 arA[1 - z][iy][ix] = an arB[1 - z][iy][ix] = bn pn = 100 * (i + 1) / steps # percent completed if pn != p: p = pn print "%" + str(p).zfill(2) # paint the final state z = 1 - z aMin = arA[z][0][0]; aMax = aMin bMin = arB[z][0][0]; bMax = bMin for iy in range(imgy): for ix in range(imgx): a = arA[z][iy][ix] b = arB[z][iy][ix] if a < aMin: aMin = a if a > aMax: aMax = a if b < bMin: bMin = b if b > bMax: bMax = b if aMin != aMax: for iy in range(imgy): for ix in range(imgx): a = arA[z][iy][ix] cA = int(255 * (a - aMin) / (aMax - aMin)) pixels[ix, iy] = (cA, cA, cA) # label = "f = " + str(f) + " k = " + str(k) # draw.text((0, 0), label, (0, 255, 0)) image.save("ReactionDiffusionSim_A.png", "PNG") if bMin != bMax: for iy in range(imgy): for ix in range(imgx): b = arB[z][iy][ix] cB = int(255 * (b - bMin) / (bMax - bMin)) pixels[ix, iy] = (cB, cB, cB) # label = "f = " + str(f) + " k = " + str(k) # draw.text((0, 0), label, (0, 255, 0)) image.save("ReactionDiffusionSim_B.png", "PNG") ```

#### 1 comment

FB36 (author) 7 years, 10 months ago

Here is a much faster NumPy version:

``````# Reaction-Diffusion Simulation Using Gray-Scott Model
# https://en.wikipedia.org/wiki/Reaction-diffusion_system
# http://www.labri.fr/perso/nrougier/teaching/numpy/numpy.html#
# FB - 20160130
import random
import numpy as np
from PIL import Image, ImageDraw
n = 256
imgx = n; imgy = n # image size
image = Image.new("RGB", (imgx, imgy))
draw = ImageDraw.Draw(image)
steps = 10000

params = []
params.append((0.16, 0.08, 0.035, 0.065)) # Bacteria 1
params.append((0.14, 0.06, 0.035, 0.065)) # Bacteria 2
params.append((0.16, 0.08, 0.060, 0.062)) # Coral
params.append((0.19, 0.05, 0.060, 0.062)) # Fingerprint
params.append((0.10, 0.10, 0.018, 0.050)) # Spirals
params.append((0.12, 0.08, 0.020, 0.050)) # Spirals Dense
params.append((0.10, 0.16, 0.020, 0.050)) # Spirals Fast
params.append((0.16, 0.08, 0.020, 0.055)) # Unstable
params.append((0.16, 0.08, 0.050, 0.065)) # Worms 1
params.append((0.16, 0.08, 0.054, 0.063)) # Worms 2
params.append((0.16, 0.08, 0.035, 0.060)) # Zebrafish
(Du, Dv, F, k) = random.choice(params)

Z = np.zeros((n+2,n+2), [('U', np.double), ('V', np.double)])
U,V = Z['U'], Z['V']
u,v = U[1:-1,1:-1], V[1:-1,1:-1]

r = 20
u[...] = 1.0
U[n/2-r:n/2+r,n/2-r:n/2+r] = 0.50
V[n/2-r:n/2+r,n/2-r:n/2+r] = 0.25
u += 0.05*np.random.random((n,n))
v += 0.05*np.random.random((n,n))

p = 0
for i in xrange(steps):
Lu = (                 U[0:-2,1:-1] +
U[1:-1,0:-2] - 4*U[1:-1,1:-1] + U[1:-1,2:] +
U[2:  ,1:-1] )
Lv = (                 V[0:-2,1:-1] +
V[1:-1,0:-2] - 4*V[1:-1,1:-1] + V[1:-1,2:] +
V[2:  ,1:-1] )
uvv = u*v*v
u += (Du*Lu - uvv +  F   *(1-u))
v += (Dv*Lv + uvv - (F+k)*v    )

pn = 100 * (i + 1) / steps # percent completed
if pn != p:
p = pn
print "%" + str(p).zfill(2)

# paint the final state
vMin=V.min(); vMax=V.max()
for iy in range(imgy):
for ix in range(imgx):
w = V[iy, ix]
c = int(255 * (w - vMin) / (vMax - vMin))
pixels[ix, iy] = (c, c, c)
label = "Du=" + str(Du) + " Dv=" + str(Dv) + " F=" + str(F) + " k=" + str(k)
draw.text((0, 0), label, (0, 255, 0))
image.save("ReactionDiffusionSim.png", "PNG")
``````
 Created by FB36 on Fri, 16 Oct 2015 (MIT)