Simulaçao Belousov-Zhabotinsky

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);


# In[ ]:


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))