Modelo Brusselator de Reação-Difusão: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
(Criou página com '''' Grupo: Carolina Lenzi, Eric Naiber e Vitória Xavier ''' O objetivo deste trabalho é implementar o modelo de reação-difusão Brusselator em duas dimensões, frequentem...')
 
Sem resumo de edição
Linha 5: Linha 5:
==Modelo de Brusselator==  
==Modelo de Brusselator==  


O estudo de sistemas químicos e biológicos frequentemente requer o uso de modelos que caracterizem reações de reação-difusão.
O estudo de sistemas químicos e biológicos frequentemente requer o uso de modelos que caracterizam reações de reação-difusão. Um dos modelos mais utilizados é o modelo de Brusselator, que é utilizado para descrever o mecanismo químico de reação-difusão com oscilações não lineares. [J. Tyson, Some further studies of nonlinear oscillations in chemical systems, J. Chem. Phys. 58 (1973) 3919.] Turing [ref?] observou que quando determinadas reações são associadas a difusão, é possível obter um padrão espacial estável, e isso leva a teoria de morfogênese. Além de processos de reação-difusão, o modelo Brusselator é observado em reações enzimáticas e na física de plasma e de lasers.
 
O mecanismo de Brusselator proposto por Prigogine (1970) é dado por [I. Prigogine, R. Lefever, Symmetries breaking instabilities in dissipative systems II. J. Phys. Chem. 48, 1695–1700 (1968)]:
 
:<math>A \to U</math> (1.a)
:<math>B + U \to V + D</math> (1. b)
:<math>2U + V \to 3U</math> (1.c)
:<math>U \to E</math> (1.d)
 
Onde U e V são as espécies químicas de interesse. Assumimos A e B em excesso para que o sistema não atinja o equilíbrio.  Esse sistema químico foi importante para o avanço na área de sistemas complexos porque possibilita o uso de modelos matemáticos de duas dimensões, já que U e V são variáveis dependentes, e admite “limit-cycle oscilations”. [ R. Lefever and G. Nicolis, Chemical instabilities and sustained oscillations, J. Theor. Biol. 30 (1971) 267.].
 
As equações diferenciais parciais associadas com o sistema Brusselator são dadas por(G. Adomian, The diffusion-Brusselator equation. Comput. Math. Appl. 29, 1–3 (1995)):
 
:<math>\frac{\partial u}{\partial t} = a + u^2v - (b+1)u + D_u \left(\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{\partial y^2} \right)</math>
:<math>\frac{\partial v}{\partial t} = bu - u^2v + D_v \left(\frac{\partial ^2 v}{\partial x^2} + \frac{\partial ^2 v}{\partial y^2}\right)</math>
 
onde <math>u(x, y, t)</math> e <math>v(x, y, t)</math> são as concentrações a serem investigadas em função de tempo e espaço, <math>a</math> e <math>b</math> são constantes relativas às concentrações dos reagentes A e B, e <math>D_u</math> e <math>D_v</math> constantes de difusão.
 
A solução analítica do sistema reação-difusão Brusselator ainda não é conhecida e por isso há o interesse de explorá-la numericamente.
 
 
==Análise da estabilidade do sistema==
 
===Análise de ponto crítico===
 
Considerando o sistema livre de difusão, quando <math>\partial x=\partial y=0</math>:
 
:<math> \frac{du}{dt} = f(u, v) = a + u^2v - (b+1)u \qquad t < 0 \quad u(0) = U^0 </math>
 
:<math> \frac{dv}{dt} = g(u, v) = bu - u^2v \qquad \qquad \qquad \, \, t < 0 \quad v(0) = V^0 </math>
 
 
Onde <math>u = u(t)</math> e <math>v = v(t)</math> e <math>a</math> e <math>b</math> são constantes positivas e reais. A matriz jacobiana <math>J^*</math> no ponto crítico <math>(u^*, v^*)</math> é dada por
 
:<math>
J^* = \begin{bmatrix}
b-1 & a^2 \\
-b & -a^2
\end{bmatrix}
</math>
 
 
Os autovalores de <math>J^*</math> são os valores <math>\lambda</math> que satisfazem a equação caracterísitca
 
:<math> \lambda ^2 + (1 - b + a^2)\lambda + a^2 = 0 </math>
 
Os autovalores claramente mostram dependência em <math>1-b+a^2</math> e no determinante <math>\Delta = (1-b+ a^2)^2 - 4a^2</math>. Esses autovalores governam a estabilidade do ponto crítico ou determinam a existência de um ciclo limite. As propriedades de estabilidade ou a existência de um ciclo limite estão sumarizadas na tabela abaixo, em relação a figura 1.
 
[[Arquivo:]]
 
Utilizando a teoria Hopf, é mostrado que o ponto crítico perde sua estabilidade quando A e B movem da região 2 para região 3, na figura 1, atravessando a curva <math>1-b+a^2 = 0</math>. Uma bifurcação Hopf ocorre quando ao passo que essa curva é atravessada e um ciclo limite estável existe para A e B nas regiões 1 e 2, mas não para A e B nas regiões 3 e 4.
 
Conclui-se que a curva <math>1-b+a^2 = 0</math> governa a estabilidade do sistema.
 
 
===Análise de ponto fixo===
 
O estado estacionário do sistema pode ser encontrado igualando o coeficiente de difusão, e portanto as derivadas parciais, a zero. Percebe-se que esse sistema converge para os pontos fixos:
 
:<math> u = a </math>
:<math> v = \frac{b}{a} </math>
 
Em [twizell] a manipulação da matriz jacobiana nos pontos fixos resulta nos seguintes autovalores
 
:<math>\mu_1 = \frac{1+ \frac{1}{2} \ell (1-b+a^2) - \frac{1}{4}\ell^2a^2 }{1+ \frac{1}{2} \ell (1-b+a^2) + \frac{1}{4}\ell^2a^2} </math>
 
 
:<math>\mu_2 = \frac{1 - \frac{1}{2} \ell (1-b+a^2) - \frac{1}{4}\ell^2a^2}{1+ \frac{1}{2} \ell (1-b+a^2) + \frac{1}{4}\ell^2a^2} </math>
 
Onde fica claro que os denominadores dos autovalores são sempre positivos quando <math> 1 - b + a^2 > 0</math> and <math> \ell > 0 </math>. As inequações
 
:<math> |\mu_1| < 1 \quad </math> e <math> \quad |\mu_2| < 1 </math>
 
São verdadeiras sempre que <math> 1 - b + a^2 = 0</math>. Portanto, uma condição suficiente para o ponto fixo <math>\left(a, \frac{b}{a}\right)</math> atrair a sequência gerada pelo sistema é <math> 1 - b + a^2 > 0</math>.
 
Ainda em [Twizell] o modelo reação-difusão Brusselator foi discretizado e a análise dos pontos fixos concluiu que o sistema converge para <math>u(x, t) = a</math> e <math> v(x, t)= \frac{b}{a} </math>, sendo esse o único estado estacionário do sistema.
 
Concluiu que também que o sistema apresenta estado oscilatório quando
 
:<math> 1 - b + a < 0 </math>
 
Estado em que o sistema não converge para nenhum ponto.


(...)


==Método FTCS==
==Método FTCS==
Linha 57: Linha 137:


onde <math>K_u = \frac{D_u \Delta t}{(\Delta s)^2}</math> e <math>K_v = \frac{D_v \Delta t}{(\Delta s)^2}</math>.
onde <math>K_u = \frac{D_u \Delta t}{(\Delta s)^2}</math> e <math>K_v = \frac{D_v \Delta t}{(\Delta s)^2}</math>.
===Análise de estabilidade do método===
==Implementação==
O método foi implementado em Python, considerando <math>\Delta t = 0.1, \Delta s = 1, N_x = N_y = 25</math>, variando as constantes <math>a</math> e <math>b</math> e as condições iniciais do problema.
<source lang="python">
import numpy as np
import matplotlib.pyplot as plt
Nx = Ny = 25
t_max = 70
dt = 0.1
ds = 1
t = 0
a = 1
b = 3
Du = 0.5
Dv = 1
ku = Du * dt / (ds ** 2)
kv = Dv * dt / (ds ** 2)
def f(u, v):
    return a - (b + 1) * u + u * u * v
def g(u, v):
    return b * u - u * u * v
# vetores no tempo n
u_n = np.zeros((Nx, Ny))
v_n = np.zeros((Nx, Ny))
# vetores no tempo n+1
u_n1 = np.zeros((Nx, Ny))
v_n1 = np.zeros((Nx, Ny))
# condicoes iniciais
"""for i in range(Nx):
    for j in range(Ny):
        u_n[i, j] = 2  # + 0.25 * j
        v_n[i, j] = 1  # + 0.8 * i"""
# condicoes novas
for i in range(5):
    u_n[int(Nx / 2) - 1 + i, int(Ny / 2) - 1 + i] = 2
for i in range(5):
    v_n[i - 5, i - 5] = 1
lista_t = []
lista_u = []
lista_v = []
lista_mesh_u = np.zeros((Nx, Ny))
# FTCS
while t < t_max:
    for i in range(Nx):
        i_e = (i - 1) % Nx  # vizinho a esquerda de 0 é o da ultima posicao
        i_d = (i + 1) % Nx  # vizinho a direita da ultima posicao é o zero
        for j in range(Ny):
            j_e = (j - 1) % Ny
            j_d = (j + 1) % Ny
            u_n1[i, j] = u_n[i, j] + dt * f(u_n[i, j], v_n[i, j]) \
                        + ku * (u_n[i_e, j] + u_n[i_d, j] + u_n[i, j_e] + u_n[i, j_d] - 4 * u_n[i, j])
            v_n1[i, j] = v_n[i, j] + dt * g(u_n[i, j], v_n[i, j]) \
                        + kv * (v_n[i_e, j] + v_n[i_d, j] + v_n[i, j_e] + v_n[i, j_d] - 4 * v_n[i, j])
            lista_mesh_u[i, j] = u_n1[i, j]
    lista_u.append(u_n1[0][0])
    lista_v.append(v_n1[0][0])
    # atualizar u_n e v_n
    for i in range(Nx):
        for j in range(Ny):
            u_n[i, j] = u_n1[i, j]
            v_n[i, j] = v_n1[i, j]
    t += dt
    lista_t.append(t)
# concentracoes x tempo
# plt.plot(lista_t, lista_u, color='red')
# plt.plot(lista_t, lista_v, color='blue')
# plt.show()
# concentracao no espaco (xy) em um instante de tempo
plt.imshow(lista_mesh_u, interpolation='none')
plt.show()
</source>

Edição das 23h49min de 9 de março de 2022

Grupo: Carolina Lenzi, Eric Naiber e Vitória Xavier

O objetivo deste trabalho é implementar o modelo de reação-difusão Brusselator em duas dimensões, frequentemente utilizado para estudar sistemas complexos químicos e biológicos. O modelo é um sistema não linear de equações diferenciais parciais e foi proposto em 1970 por Ilya Prigogine e seus colaboradores da Universidade Livre de Bruxelas. Desde então tem sido aplicado para analisar reações oscilatórias e autocatalíticas. O método computacional utilizado para implementar o modelo foi o método FTCS (Forward Time Centered Space).

Modelo de Brusselator

O estudo de sistemas químicos e biológicos frequentemente requer o uso de modelos que caracterizam reações de reação-difusão. Um dos modelos mais utilizados é o modelo de Brusselator, que é utilizado para descrever o mecanismo químico de reação-difusão com oscilações não lineares. [J. Tyson, Some further studies of nonlinear oscillations in chemical systems, J. Chem. Phys. 58 (1973) 3919.] Turing [ref?] observou que quando determinadas reações são associadas a difusão, é possível obter um padrão espacial estável, e isso leva a teoria de morfogênese. Além de processos de reação-difusão, o modelo Brusselator é observado em reações enzimáticas e na física de plasma e de lasers.

O mecanismo de Brusselator proposto por Prigogine (1970) é dado por [I. Prigogine, R. Lefever, Symmetries breaking instabilities in dissipative systems II. J. Phys. Chem. 48, 1695–1700 (1968)]:

(1.a)
(1. b)
(1.c)
(1.d)

Onde U e V são as espécies químicas de interesse. Assumimos A e B em excesso para que o sistema não atinja o equilíbrio. Esse sistema químico foi importante para o avanço na área de sistemas complexos porque possibilita o uso de modelos matemáticos de duas dimensões, já que U e V são variáveis dependentes, e admite “limit-cycle oscilations”. [ R. Lefever and G. Nicolis, Chemical instabilities and sustained oscillations, J. Theor. Biol. 30 (1971) 267.].

As equações diferenciais parciais associadas com o sistema Brusselator são dadas por(G. Adomian, The diffusion-Brusselator equation. Comput. Math. Appl. 29, 1–3 (1995)):

onde e são as concentrações a serem investigadas em função de tempo e espaço, e são constantes relativas às concentrações dos reagentes A e B, e e constantes de difusão.

A solução analítica do sistema reação-difusão Brusselator ainda não é conhecida e por isso há o interesse de explorá-la numericamente.


Análise da estabilidade do sistema

Análise de ponto crítico

Considerando o sistema livre de difusão, quando :


Onde e e e são constantes positivas e reais. A matriz jacobiana no ponto crítico é dada por


Os autovalores de são os valores que satisfazem a equação caracterísitca

Os autovalores claramente mostram dependência em e no determinante . Esses autovalores governam a estabilidade do ponto crítico ou determinam a existência de um ciclo limite. As propriedades de estabilidade ou a existência de um ciclo limite estão sumarizadas na tabela abaixo, em relação a figura 1.

[[Arquivo:]]

Utilizando a teoria Hopf, é mostrado que o ponto crítico perde sua estabilidade quando A e B movem da região 2 para região 3, na figura 1, atravessando a curva . Uma bifurcação Hopf ocorre quando ao passo que essa curva é atravessada e um ciclo limite estável existe para A e B nas regiões 1 e 2, mas não para A e B nas regiões 3 e 4.

Conclui-se que a curva governa a estabilidade do sistema.


Análise de ponto fixo

O estado estacionário do sistema pode ser encontrado igualando o coeficiente de difusão, e portanto as derivadas parciais, a zero. Percebe-se que esse sistema converge para os pontos fixos:

Em [twizell] a manipulação da matriz jacobiana nos pontos fixos resulta nos seguintes autovalores


Onde fica claro que os denominadores dos autovalores são sempre positivos quando and . As inequações

e

São verdadeiras sempre que . Portanto, uma condição suficiente para o ponto fixo atrair a sequência gerada pelo sistema é .

Ainda em [Twizell] o modelo reação-difusão Brusselator foi discretizado e a análise dos pontos fixos concluiu que o sistema converge para e , sendo esse o único estado estacionário do sistema.

Concluiu que também que o sistema apresenta estado oscilatório quando

Estado em que o sistema não converge para nenhum ponto.


Método FTCS

O FTCS (Forward Time Centered Space) é um método de diferença finita que utiliza a derivada à direita ("para frente") no tempo e a derivada segunda centralizada no espaço para discretizar as variáveis. As derivadas no tempo e no espaço bidimensional ficam:


Substituindo nas equações do Brusselator

onde e são as funções que representam a reação sem difusão.


Utilizamos discretização do tipo


Utilizando a notação , assumindo e rearranjando os termos, reescrevemos as equações como


onde e .

Análise de estabilidade do método

Implementação

O método foi implementado em Python, considerando , variando as constantes e e as condições iniciais do problema.

import numpy as np
import matplotlib.pyplot as plt

Nx = Ny = 25
t_max = 70
dt = 0.1
ds = 1

t = 0

a = 1
b = 3

Du = 0.5
Dv = 1
ku = Du * dt / (ds ** 2)
kv = Dv * dt / (ds ** 2)


def f(u, v):
    return a - (b + 1) * u + u * u * v


def g(u, v):
    return b * u - u * u * v


# vetores no tempo n
u_n = np.zeros((Nx, Ny))
v_n = np.zeros((Nx, Ny))

# vetores no tempo n+1
u_n1 = np.zeros((Nx, Ny))
v_n1 = np.zeros((Nx, Ny))

# condicoes iniciais
"""for i in range(Nx):
    for j in range(Ny):
        u_n[i, j] = 2  # + 0.25 * j
        v_n[i, j] = 1  # + 0.8 * i"""

# condicoes novas
for i in range(5):
    u_n[int(Nx / 2) - 1 + i, int(Ny / 2) - 1 + i] = 2

for i in range(5):
    v_n[i - 5, i - 5] = 1

lista_t = []
lista_u = []
lista_v = []
lista_mesh_u = np.zeros((Nx, Ny))

# FTCS
while t < t_max:

    for i in range(Nx):
        i_e = (i - 1) % Nx  # vizinho a esquerda de 0 é o da ultima posicao
        i_d = (i + 1) % Nx  # vizinho a direita da ultima posicao é o zero

        for j in range(Ny):
            j_e = (j - 1) % Ny
            j_d = (j + 1) % Ny

            u_n1[i, j] = u_n[i, j] + dt * f(u_n[i, j], v_n[i, j]) \
                         + ku * (u_n[i_e, j] + u_n[i_d, j] + u_n[i, j_e] + u_n[i, j_d] - 4 * u_n[i, j])
            v_n1[i, j] = v_n[i, j] + dt * g(u_n[i, j], v_n[i, j]) \
                         + kv * (v_n[i_e, j] + v_n[i_d, j] + v_n[i, j_e] + v_n[i, j_d] - 4 * v_n[i, j])

            lista_mesh_u[i, j] = u_n1[i, j]

    lista_u.append(u_n1[0][0])
    lista_v.append(v_n1[0][0])

    # atualizar u_n e v_n
    for i in range(Nx):
        for j in range(Ny):
            u_n[i, j] = u_n1[i, j]
            v_n[i, j] = v_n1[i, j]

    t += dt
    lista_t.append(t)


# concentracoes x tempo
# plt.plot(lista_t, lista_u, color='red')
# plt.plot(lista_t, lista_v, color='blue')
# plt.show()

# concentracao no espaco (xy) em um instante de tempo
plt.imshow(lista_mesh_u, interpolation='none')
plt.show()