Simulaçao Belousov-Zhabotinsky

<source lang = "py"> import matplotlib matplotlib.use('TkAgg') from pylab import * import numpy as np

n = 100 Dh = 1./n Dt = 0.001

Du = 1e-5 Dv = 1e-5

eps = 0.2

q = 1e-3

f = 1

def inicia():

   global u, v, nextu, nextv
   u = np.loadtxt('u_file.txt')
   v = np.loadtxt('v_file.txt')

   nextu = zeros([n, n])
   nextv = zeros([n, n])

def update():

   global u, v, nextu, nextv
   for i in range(1, n - 1):
       for j in range(1, n - 1):
           uC = u[i, j]
           uR = u[(i + 1) % n, j]
           uL = u[(i - 1) % n, j]
           uU = u[i, (j + 1) % n]
           uD = u[i, (j - 1) % n]
           vC = v[i, j]
           vR = v[(i + 1) % n, j]
           vL = v[(i - 1) % n, j]
           vU = v[i, (j + 1) % n]
           vD = v[i, (j - 1) % n]

           uLap = (uR + uL + uU + uD - 4 * uC) / (Dh ** 2)
           vLap = (vR + vL + vU + vD - 4 * vC) / (Dh ** 2)
           nextu[i, j] = uC + (uC * (1 - uC) - f * vC * (uC - q) / (uC + q) + Du * uLap) * Dt / eps
           nextv[i, j] = vC + (uC - vC + Dv * vLap) * Dt

   u, nextu = nextu, u
   v, nextv = nextv, v

def observeu():

   global u, v, nextu, nextv
   fig, ax = plt.subplots(figsize=(16, 9))
   ax.matshow(u, cmap='plasma')
   ax.axis(False);

def observev():

   global u, v, nextu, nextv
   fig, ax = plt.subplots(figsize=(16, 9))
   ax.matshow(v, cmap='winter')
   ax.axis(False);

inicia() t = 0 observeu() for t in range(1, 100001):

   update()
   if t % 500 == 0:
       print(t)
   if t % 1000 == 0:
       observeu()
       savefig('u-BZ-{}.png'.format(t))
       observev()
       savefig('v-BZ-{}.png'.format(t))