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

De Física Computacional
Ir para navegação Ir para pesquisar
 
(29 revisões intermediárias por 3 usuários não estão sendo mostradas)
Linha 1: Linha 1:
''' Grupo: Carolina Lenzi, Eric Naiber e Vitória Xavier '''
''' 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'').  
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 1968 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==  
==Modelo de Brusselator==  
Linha 7: Linha 7:
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. <ref>J. Tyson, Some further studies of nonlinear oscillations in chemical systems, J. Chem. Phys. 58 (1973) 3919</ref>. Turing 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<ref>Turing A. M., The chemical basis of morphogenesis. 1953. Bull. Math. Biol. 52, 153, discussion 119 (1990).</ref>. 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 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. <ref>J. Tyson, Some further studies of nonlinear oscillations in chemical systems, J. Chem. Phys. 58 (1973) 3919</ref>. Turing 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<ref>Turing A. M., The chemical basis of morphogenesis. 1953. Bull. Math. Biol. 52, 153, discussion 119 (1990).</ref>. 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<ref>I. Prigogine, R. Lefever, Symmetries breaking instabilities in dissipative systems II. J. Phys. Chem. 48, 1695–1700 (1968)</ref> proposto por Prigogine (1970) é dado por:  
O mecanismo de Brusselator proposto por Prigogine<ref>I. Prigogine, R. Lefever, Symmetries breaking instabilities in dissipative systems II. J. Phys. Chem. 48, 1695–1700 (1968)</ref> é dado por:  


:<math>A \to U</math>        (1.a)
:<math>A \to U</math>        (1.a)
Linha 14: Linha 14:
:<math>U \to E</math> (1.d)
:<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 oscilações de ciclo limite<ref>R. Lefever and G. Nicolis, Chemical instabilities and sustained oscillations, J. Theor. Biol. 30 (1971)</ref>.
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 estado de 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 oscilações de ciclo limite<ref>R. Lefever and G. Nicolis, Chemical instabilities and sustained oscillations, J. Theor. Biol. 30 (1971)</ref>.


As equações diferenciais parciais associadas com o sistema Brusselator<ref>G. Adomian, The diffusion-Brusselator equation. Comput. Math. Appl. 29, 1–3 (1995)</ref> são dadas por:
As equações diferenciais parciais associadas com o sistema Brusselator<ref>G. Adomian, The diffusion-Brusselator equation. Comput. Math. Appl. 29, 1–3 (1995)</ref> são dadas por:
Linha 21: Linha 21:
:<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>
:<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.
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> são as constantes de difusão de U e V.


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.
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.
Linha 29: Linha 29:
===Análise de ponto crítico===
===Análise de ponto crítico===


Considerando o sistema livre de difusão, quando <math>\partial x=\partial y=0</math>:
Considerando o sistema livre de difusão, quando <math>D_u = D_v = 0</math>:


:<math> \frac{du}{dt} = f(u, v) = a + u^2v - (b+1)u </math>
:<math> \frac{du}{dt} = f(u, v) = a + u^2v - (b+1)u </math>
Linha 35: Linha 35:
:<math> \frac{dv}{dt} = g(u, v) = bu - u^2v </math>
:<math> \frac{dv}{dt} = g(u, v) = bu - u^2v </math>


 
Encontramos os pontos críticos do sistema igualando as derivadas à zero. Obtemos que o único ponto crítico é <math>(u^*, v^*) = \left ( a, \frac{b}{a} \right)</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
A matriz jacobiana <math>J^*</math> no ponto crítico <math>(u^*, v^*)</math> é dada por


:<math>
:<math>
J^* = \begin{bmatrix}
J^* = \begin{bmatrix}
\frac{\partial f}{\partial u} \Big |_{u^*, v^*} & \frac{\partial f}{\partial v} \Big |_{u^*, v^*} \\
\frac{\partial g}{\partial u} \Big |_{u^*, v^*} & \frac{\partial g}{\partial v} \Big |_{u^*, v^*}
\end{bmatrix}
=
\begin{bmatrix}
b-1 & a^2 \\
b-1 & a^2 \\
-b & -a^2
-b & -a^2
Linha 48: Linha 53:
Os autovalores de <math>J^*</math> são os valores <math>\lambda</math> que satisfazem a equação caracterísitca
Os autovalores de <math>J^*</math> são os valores <math>\lambda</math> que satisfazem a equação caracterísitca


:<math> \lambda ^2 + (1 - a + b^2)\lambda + b^2 = 0 </math>
:<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 - a+ b^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-a+b^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.
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 e determinam a existência de um ciclo limite. Conforme discutido em Twizell<ref name='t'>Twizell, E. H., Gumel, A. B., & Cao, Q. (1999). Journal of Mathematical Chemistry, 26(4), 297–316 </ref>, quando <math>Re(\lambda) < 0 </math> o sistema é estável e não existe ciclo limite. Já para <math>Re(\lambda) > 0 </math> o sistema é instável e existe ciclo limite.
 
Conclui-se que a curva <math>1-a+b^2 = 0</math> governa a estabilidade do sistema.


===Análise de ponto fixo===
===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:
Ainda em Twizell<ref name='t'>Twizell, E. H., Gumel, A. B., & Cao, Q. (1999). Journal of Mathematical Chemistry, 26(4), 297–316 </ref>, foi feita a análise de ponto fixo. O único ponto fixo encontrado é <math>(u^*, v^*) = \left ( a, \frac{b}{a} \right)</math>. Manipulando a matriz Jacobiana nesse ponto e analisando os autovalores, obtém-se que a solução é atraída para o ponto <math>\left ( a, \frac{b}{a} \right)</math> se <math>1-b+a^2 > 0</math>. Caso contrário, a solução não converge.
 
:<math> u = a </math>
:<math> v = \frac{b}{a} </math>
 
Em Twizell et al.<ref>A second-order scheme for the Brusselator reaction-diffusion system
J. Math. Chem, 26 (1999), pp. 297-316
</ref> a manipulação da matriz jacobiana nos pontos fixos <math>U=A</math> e <math>V=B/A</math> resulta em autovalores em que seus denominadores são positivos quando <math> 1 - a + b^2 > 0</math>. Isso deixa claro que os denominadores dos autovalores são sempre positivos quando 1 − A + B2 > 0. O módulo desses autovalores é menor que 1 sempre que <math>1 − A + B2 = 0</math> for verdade. Portanto, uma condição suficiente para o ponto fixo (B, A/B) atrair a sequência gerada pelo sistema é 1 − A + B² > 0.
 
 
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 132: Linha 115:


:<math>
:<math>
\delta u_{i, j}^{n+1} = \delta u_{i, j}^n + K_u \left( \delta u_{i+1, j}^n + \delta u_{i-1, j}^n \delta u_{i, j+1}^n + \delta u_{i, j-1}^n - 4 \delta u_{i, j}^n  \right) + \Delta t f_u \delta u_{i, j}^n + \Delta t f_v \delta v_{i, j}^n
\delta u_{i, j}^{n+1} = \delta u_{i, j}^n + K_u \left( \delta u_{i+1, j}^n + \delta u_{i-1, j}^n + \delta u_{i, j+1}^n + \delta u_{i, j-1}^n - 4 \delta u_{i, j}^n  \right) + \Delta t f_u \delta u_{i, j}^n + \Delta t f_v \delta v_{i, j}^n
</math>
</math>


:<math>
:<math>
\delta v_{i, j}^{n+1} = \delta v_{i, j}^n + K_v \left( \delta v_{i+1, j}^n + \delta v_{i-1, j}^n \delta v_{i, j+1}^n + \delta v_{i, j-1}^n - 4 \delta v_{i, j}^n  \right) + \Delta t g_u \delta u_{i, j}^n + \Delta t g_v \delta v_{i, j}^n
\delta v_{i, j}^{n+1} = \delta v_{i, j}^n + K_v \left( \delta v_{i+1, j}^n + \delta v_{i-1, j}^n + \delta v_{i, j+1}^n + \delta v_{i, j-1}^n - 4 \delta v_{i, j}^n  \right) + \Delta t g_u \delta u_{i, j}^n + \Delta t g_v \delta v_{i, j}^n
</math>
</math>


:<math>
:<math>
{\delta u^n \choose \delta v^n} = A^n e^{i(k_x i \Delta x + k_y j \Delta y)} {\delta u_0 \choose \delta v_0}
{\delta u^n \choose \delta v^n} = \xi ^n e^{i(k_x i \Delta x + k_y j \Delta y)} {\delta u_0 \choose \delta v_0}
</math>
</math>


Para que o método seja estável, é preciso que <math>|A| < 1</math>. Após inserir as soluções nas equações linearizadas e realizar manipulações algébricas, como feito por Scholz (2009)<ref name='scholz'>Scholz, Christian, Morphology of Experimental and Simulated Turing Patterns, 2009, p. 27-29.</ref>, obtemos
Para que o método seja estável, é preciso que <math>|\xi| < 1</math>. Após inserir as soluções nas equações linearizadas e realizar manipulações algébricas, como feito por Scholz<ref name='s'>Scholz, Christian, Morphology of Experimental and Simulated Turing Patterns, 2009, p. 27-29.</ref>, obtemos
:<math>
:<math>
A_\pm = \frac{1}{2} \left [ \beta \pm \sqrt{\beta ^2 - 4 (\Delta t)^2 a^2 b} \right ]
\xi_\pm = \frac{1}{2} \left [ \beta \pm \sqrt{\beta ^2 - 4 (\Delta t)^2 a^2 b} \right ]
</math>
</math>
onde
onde
Linha 152: Linha 135:
:<math>\gamma = sin^2\left( \frac{k_x \Delta s}{2} \right) + sin^2\left( \frac{k_y \Delta s}{2} \right)</math>
:<math>\gamma = sin^2\left( \frac{k_x \Delta s}{2} \right) + sin^2\left( \frac{k_y \Delta s}{2} \right)</math>


Os valores <math>\Delta t = 0.01</math> e <math>\Delta s = 1</math> satisfazem a condição de estabilidade, como mostrado em <ref name='scholz'>Scholz, Christian, Morphology of Experimental and Simulated Turing Patterns, 2009, p. 27-29.</ref>.
Os valores <math>\Delta t = 0.01</math> e <math>\Delta s = 1</math> satisfazem a condição de estabilidade, como mostrado em Scholz<ref name='s'>Scholz, Christian, Morphology of Experimental and Simulated Turing Patterns, 2009, p. 27-29.</ref>.


===Implementação===
===Implementação===


O método foi implementado em Python, utilizando as constantes definidas na tabela abaixo.
A seguir, o trecho do código relativo à implementação do método na linguagem Python. O código completo encontra-se no Github [https://github.com/RicGary/MetCompC/tree/main/Brusselator]
 
Os valores de <math>b</math> e as condições iniciais foram variados, como descrito na seção de resultados.
 
A seguir, o trecho do código relativo à implementação do método. O código completo encontra-se no Github [https://github.com/RicGary/MetCompC/tree/main/Brusselator]


<source lang="python">
<source lang="python">
# u_n, v_n -> concentracao dos reagentes no tempo n (matriz Nx x Ny)
# u_n, v_n -> concentracao dos reagentes no tempo n (matriz Nx x Ny)
# u_n1, v_n1 -> concentracao dos reagentes no tempo n+1
# u_n1, v_n1 -> concentracao dos reagentes no tempo n+1
# Nx, Ny -> Tamanho do recipiente


    while t < t_max:
while t < t_max:


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


            for j in range(Ny):
        for j in range(Ny):
                j_e = (j - 1) % Ny
            j_e = (j - 1) % Ny
                j_d = (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], b) \
            # calcular u e v no tempo n+1
                            + 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])
            u_n1[i, j] = u_n[i, j] + dt * f(u_n[i, j], v_n[i, j], b) \
                v_n1[i, j] = v_n[i, j] + dt * g(u_n[i, j], v_n[i, j], b) \
                        + 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])
                            + 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])
            v_n1[i, j] = v_n[i, j] + dt * g(u_n[i, j], v_n[i, j], b) \
                        + 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])


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


        t += dt
    t += dt
</source>
</source>


==Resultados==
==Resultados==
Esta seção se dedica ao resultados das simulações realizadas, que mostram o caráter oscilatório do Brusselator ao longo do tempo e também como se dá a difusão dos reagentes em um recipiente em duas dimensões.
Esta seção se dedica ao resultados das simulações realizadas, que mostram o caráter oscilatório do Brusselator ao longo do tempo e também como se dá a difusão dos reagentes em um recipiente em duas dimensões.
As constantes utilizadas foram identificadas abaixo para plotar os gráficos (exceto quando indicado o contrário).
As constantes utilizadas foram identificadas abaixo para plotar os gráficos, exceto quando indicado o contrário.


{| class="wikitable" style="margin-left: auto; margin-right: 0px;" style="text-align:center;"  
{| class="wikitable" style="margin-left: auto; margin-right: 0px;" style="text-align:center;"  
Linha 218: Linha 196:
| <math>b</math>
| <math>b</math>
| style="text-align:left;" | Constante relativa à concentração do reagente B
| style="text-align:left;" | Constante relativa à concentração do reagente B
| -
| 1.7
|-
|-
| <math>D_u</math>
| <math>D_u</math>
Linha 230: Linha 208:
| <math>u_0</math>
| <math>u_0</math>
| style="text-align:left;" | Concentração de u no tempo inicial (t=0)
| style="text-align:left;" | Concentração de u no tempo inicial (t=0)
| -
| 1
|-
|-
| <math>v_0</math>
| <math>v_0</math>
| style="text-align:left;" | Concentração de v no tempo inicial (t=0)
| style="text-align:left;" | Concentração de v no tempo inicial (t=0)
| -
| 1
|-
|-
| <math>dt</math>
| <math>dt</math>
Linha 246: Linha 224:
|}
|}


 
As posições das condições iniciais foram escolhidas arbitrariamente, tendo um total de 6 modos que serão explicados na sequência, sendo eles:
As posições das condições iniciais foram escolhidas arbitrariamente, tendo um total de 6 modos que serão explicados logo à frente, sendo eles:
# Condição em formato do sinal "+".
# Condição em formato do sinal "+".
# Condição de borda.
# Condição de borda.
Linha 255: Linha 232:


===Simulação===
===Simulação===
Nas figuras abaixo temos a variação das concentrações dos reagentes ao longo do tempo e os diagramas de fase. Para <math>a = 1</math> e <math>b = 1.7</math> notamos que a solução do sistema converge para <math>u(t) = 1</math> e <math>v(t) = 1.7</math>. Esse resultado está de acordo com o esperado, pois o ponto <math>\left(1, 1.7\right)</math> é um ponto fixo atrator (satisfaz que <math>1-b+a^2 > 0</math>), conforme demonstrado na seção de análise de estabilidade do Brusselator. Notamos também, pelo diagrama de fase, que mesmo alterando as condições iniciais do problema, a solução converge para o estado estacionário.
Nas figuras abaixo temos a variação das concentrações dos reagentes ao longo do tempo e os diagramas de fase. Para <math>a = 1</math> e <math>b = 1.7</math> notamos que a solução do sistema converge para <math>u(t) = 1</math> e <math>v(t) = 1.7</math>. Esse resultado está de acordo com o esperado, pois o ponto <math>\left(1, 1.7\right)</math> é um ponto fixo atrator (satisfaz que <math>1-b+a^2 > 0</math>), conforme discutido na seção de análise de estabilidade do Brusselator. Notamos também, pelo diagrama de fase, que mesmo alterando as condições iniciais do problema, a solução converge para o estado estacionário.




[[Arquivo:Sol_bru_b1.png|thumb|center|790px|'''Figura 1''': Gráfico de concentrações.]]
[[Arquivo:Sol_bru_b1.png|thumb|center|790px|'''Figura 1''': Solução do Brusselator para a = 1 e b = 1.7.]]




Para <math>a = 1</math> e <math>b = 3</math>, temos que <math>1-b+a^2 < 0</math>, portanto a solução do sistema não converge. O que observamos é um comportamento periódico.
Para <math>a = 1</math> e <math>b = 3</math>, temos que <math>1-b+a^2 < 0</math>, portanto a solução do sistema não converge. O que observamos é um comportamento periódico causado pela existência do ciclo limite.




[[Arquivo:Sol_bru_b3.png|thumb|center|790px|'''Figura 2''': Gráfico de concentrações.]]
[[Arquivo:Sol_bru_b3.png|thumb|center|790px|'''Figura 2''': Solução do Brusselator para a = 1 e b = 3.]]


Podemos pensar nas funções de concentração <math>u</math> e <math>v</math> como ondas que sofrem interferência destrutiva, pois estão fora de fase. A animação abaixo mostra este comportamento, observe o reagente <math>U</math>, representado pela cor amarela, no gráfico da direita. Note que quando o ponto do gráfico da esquerda chega ao mínimo local, o valor de <math>U</math> é alto e o de <math>V</math> muito baixo, fazendo com que o gráfico da direita fique muito amarelo.


Note que quando <math>U</math> diminui, <math>V</math> aumenta rapidamente, mantendo-se um ciclo oscilatório em que a amplitude das ondas de <math>U</math> e <math>V</math> vão diminuindo ao longo do tempo (ou não, dependendo das constantes <math>a</math> e <math>b</math> teremos outros comportamentos), como mostrado nas figuras.


Pensando no gráfico da reação-difusão como ondas, é correto afirmar que ambas juntas formam um padrão de interferência destrutiva. A animação abaixo mostra exatamente este comportamento, observe o reagente <math>U</math> no gráfico da direita. Note que quando o ponto do gráfico da esquerda chega ao mínimo local, o valor de <math>U</math> é alto e o de <math>V</math> muito baixo, fazendo com que o gráfico da direita fique muito amarelo (amarelo significa uma alta concentração de <math>U</math> sendo formada).
[[Arquivo:AnaliseMid.gif|thumb|center|1300px|'''Figura 3''': Comportamento do reagente.]]




[[Arquivo:AnaliseMid.gif|thumb|center|1300px|'''Figura 3:''': Comportamento do reagente.]]
É interessante pensar no gráfico da esquerda como a análise em um ponto específico dentro do recipiente, por exemplo, se <math>Nx=Ny=50</math> estamos analisando o ponto <math>(10,10)</math>. Tendo isso em mente conseguimos analisar os outros 3 cantos da figura, já que, por simetria da posição inicial escolhida, os cantos são iguais.




É interessante pensar no gráfico da esquerda como a análise em um ponto específico dentro do recipiente, por exemplo, se <math>Nx=Ny=50</math> estamos analisando o ponto <math>(10,10)</math>. Tendo isso em mente conseguimos analisar os outros 3 cantos da figura, já que por simetria da posição inicial escolhida, os cantos são iguais.
[[Arquivo:ZoomGraph.jpg|thumb|center|800px|'''Figura 4''': Análise de um ponto.]]
 
 
[[Arquivo:ZoomGraph.jpg|thumb|center|800px|'''Figura 4:''': Análise de um ponto.]]


----
----


===Condições Iniciais===
===Condições Iniciais===
Para observar melhor os resultados, foram feitas diversas condições iniciais que foram separadas em tópicos citados em '''Definindo Constantes'''.
Realizamos as simulações variando espacialmente as condições iniciais, para verificar os diferentes padrões formados pelo sistema.
 
Levando em consideração que u_n é um array de tamanho <math>N_x</math> por <math>N_y</math>.


====Condição em Formato do Sinal "+"====  
====Condição em Formato do Sinal "+"====  
Linha 291: Linha 263:
O reagente <math>u_0</math> é distribuído no centro formando um sinal de "+".
O reagente <math>u_0</math> é distribuído no centro formando um sinal de "+".


[[Arquivo:CruzBorda.gif|thumb|right|338px|'''Figura 4:''': Comportamento em formato de "+" para o reagente u. Para o reagente v foi utilizado a condição de borda.]]
[[Arquivo:CruzBorda.gif|thumb|right|338px|'''Figura 4''': Comportamento em formato de "+" para o reagente u. Para o reagente v foi utilizado a condição de borda.]]


<source lang="python">
<source lang="python">
Linha 324: Linha 296:
O reagente <math>v_0</math> foi distribuído nos 4 cantos do recipiente.
O reagente <math>v_0</math> foi distribuído nos 4 cantos do recipiente.


[[Arquivo:CruzBordaB3.gif|thumb|right|338px|'''Figura 5:''': Analisando u com condição de borda "+". Para o reagente v foi utilizado a condição de borda. Foi utilizado b=3 para mostrar a diferença na simulação.]]
[[Arquivo:CruzBordaB3.gif|thumb|right|338px|'''Figura 5''': Analisando u com condição de borda "+". Para o reagente v foi utilizado a condição de borda. Foi utilizado b=3 para mostrar a diferença na simulação.]]


<source lang="python">
<source lang="python">
Linha 353: Linha 325:
Os reagentes são aleatoriamente distribuídos ao longo do recipiente, a quantidade de focos de reagente também é aleatória para <math>u_0</math> e <math>v_0</math>.
Os reagentes são aleatoriamente distribuídos ao longo do recipiente, a quantidade de focos de reagente também é aleatória para <math>u_0</math> e <math>v_0</math>.


[[Arquivo:AleatorioResize.gif|thumb|right|360px|'''Figura 6:''': Posição aleatória.]]
[[Arquivo:AleatorioResize.gif|thumb|right|360px|'''Figura 6''': Posição aleatória.]]


<source lang="python">
<source lang="python">
Linha 379: Linha 351:
====Condição de nove pontos centrais====
====Condição de nove pontos centrais====


Foi utilizado na figura 3.
No centro divide 9 pontos igualmente espaçados de tamanho <math>N_x/4</math>, para <math>u_0</math>. Foi utilizado para <math>u_0</math> na figura 3.
 
No centro divide 9 pontos igualmente espaçados de tamanho <math>N_x/4</math>, para <math>u_0</math>.


[[Arquivo:9PontosBorda.gif|thumb|right|384px|'''Figura 7:''': Utilizando a condição de borda para v.]]
[[Arquivo:9PontosBorda.gif|thumb|right|384px|'''Figura 7''': Utilizando a condição de borda para v.]]


<source lang="python">
<source lang="python">
Linha 409: Linha 379:


====Condição da borda completa====
====Condição da borda completa====
Foi utilizado na figura 3.


Completa toda a borda do recipiente com reagente <math>v_0</math>.
Completa toda a borda do recipiente com reagente <math>v_0</math>. Foi utilizado para <math>v_0</math> na figura 3.


[[Arquivo:9PontosB3.gif|thumb|right|359px|'''Figura 8:''': Utilizando a condição de borda para v e b=3.]]
[[Arquivo:9PontosB3.gif|thumb|right|359px|'''Figura 8''': Utilizando a condição de borda para v e b=3.]]


<source lang="python">
<source lang="python">

Edição atual tal como às 03h33min de 15 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 1968 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. [1]. Turing 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[2]. 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[3] é dado por:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle A \to U} (1.a)
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle B + U \to V + D} (1. b)
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle 2U + V \to 3U} (1.c)
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle U \to E} (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 estado de 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 oscilações de ciclo limite[4].

As equações diferenciais parciais associadas com o sistema Brusselator[5] são dadas por:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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)}
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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)}

onde Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle u(x, y, t)} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle v(x, y, t)} são as concentrações a serem investigadas em função de tempo e espaço, Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle a} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle b} são constantes relativas às concentrações dos reagentes A e B, e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle D_u} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle D_v} são as constantes de difusão de U e V.

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 Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle D_u = D_v = 0} :

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{du}{dt} = f(u, v) = a + u^2v - (b+1)u }
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{dv}{dt} = g(u, v) = bu - u^2v }

Encontramos os pontos críticos do sistema igualando as derivadas à zero. Obtemos que o único ponto crítico é Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle (u^*, v^*) = \left ( a, \frac{b}{a} \right)} . A matriz jacobiana Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle J^*} no ponto crítico Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle (u^*, v^*)} é dada por

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle J^* = \begin{bmatrix} \frac{\partial f}{\partial u} \Big |_{u^*, v^*} & \frac{\partial f}{\partial v} \Big |_{u^*, v^*} \\ \frac{\partial g}{\partial u} \Big |_{u^*, v^*} & \frac{\partial g}{\partial v} \Big |_{u^*, v^*} \end{bmatrix} = \begin{bmatrix} b-1 & a^2 \\ -b & -a^2 \end{bmatrix} }


Os autovalores de Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle J^*} são os valores Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \lambda} que satisfazem a equação caracterísitca

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \lambda ^2 + (1 - b + a^2)\lambda + a^2 = 0 }

Os autovalores claramente mostram dependência em Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle 1-b+a^2} e no determinante Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Delta = (1 - b + a^2)^2 - 4a^2} . Esses autovalores governam a estabilidade do ponto crítico e determinam a existência de um ciclo limite. Conforme discutido em Twizell[6], quando Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle Re(\lambda) < 0 } o sistema é estável e não existe ciclo limite. Já para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle Re(\lambda) > 0 } o sistema é instável e existe ciclo limite.

Análise de ponto fixo

Ainda em Twizell[6], foi feita a análise de ponto fixo. O único ponto fixo encontrado é Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle (u^*, v^*) = \left ( a, \frac{b}{a} \right)} . Manipulando a matriz Jacobiana nesse ponto e analisando os autovalores, obtém-se que a solução é atraída para o ponto Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \left ( a, \frac{b}{a} \right)} se Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle 1-b+a^2 > 0} . Caso contrário, a solução não converge.

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:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{\partial f(x, y, t)}{\partial t} \to \frac{f(x, y, t + \Delta t) - f(x, y, t)}{\Delta t} + \mathcal{O}(\Delta t ^2) }
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{\partial^2 f(x, y, t)}{\partial x^2} \to \frac{f(x-\Delta x, y, t) - 2f(x, y, t) + f(x+\Delta x, y, t)}{(\Delta x)^2} + \mathcal{O}(\Delta x ^3) }
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{\partial^2 f(x, y, t)}{\partial y^2} \to \frac{f(x, y-\Delta y, t) - 2f(x, y, t) + f(x, y+\Delta y, t)}{(\Delta y)^2} + \mathcal{O}(\Delta y ^3) }


Substituindo nas equações do Brusselator

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{u(x, y, t + \Delta t) - u(x, y, t)}{\Delta t} = f(u, v) + D_u \left( \frac{u(x-\Delta x, y, t) - 2u(x, y, t) + u(x+\Delta x, y, t)}{(\Delta x)^2} + \frac{u(x, y-\Delta y, t) - 2u(x, y, t) + u(x, y+\Delta y, t)}{(\Delta y)^2} \right) }
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{v(x, y, t + \Delta t) - v(x, y, t)}{\Delta t} = g(u, v) + D_v \left( \frac{v(x-\Delta x, y, t) - 2v(x, y, t) + v(x+\Delta x, y, t)}{(\Delta x)^2} + \frac{v(x, y-\Delta y, t) - 2v(x, y, t) + v(x, y+\Delta y, t)}{(\Delta y)^2} \right) }

onde Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle f(u, v)} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle g(u, v)} são as funções que representam a reação sem difusão.


Utilizamos discretização do tipo

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle t = 0, 1\Delta t, 2\Delta t, 3\Delta t, ..., t_{max}}
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle x = 0, 1\Delta x, 2\Delta x, 3\Delta x, ..., N_x}
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle y = 0, 1\Delta y, 2\Delta y, 3\Delta y, ..., N_y}


Utilizando a notação Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle f(i\Delta x, j\Delta y, n\Delta t) = f_{i,j}^n} , assumindo Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Delta x = \Delta y = \Delta s} e rearranjando os termos, reescrevemos as equações como

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle u_{i,j}^{n+1} = u_{i, j}^n + f(u_{i,j}^n, v_{i, j}^n)\Delta t + K_u (u_{i-1, j}^n + u_{i+1, j}^n + u_{i, j-1}^n + u_{i, j+1}^n - 4 u_{i, j}^n) }
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle v_{i,j}^{n+1} = v_{i, j}^n + g(u_{i,j}^n, v_{i, j}^n)\Delta t + K_v (v_{i-1, j}^n + v_{i+1, j}^n + v_{i, j-1}^n + v_{i, j+1}^n - 4 v_{i, j}^n) }


onde Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle K_u = \frac{D_u \Delta t}{(\Delta s)^2}} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle K_v = \frac{D_v \Delta t}{(\Delta s)^2}} .

Análise de estabilidade

A análise de estabilidade do método FTCS pode ser feita com a análise local de von Neumann. Para isso precisamos linearizar as equações do modelo Brusselator e escrever as soluções em modos de Fourier, da seguinte forma:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \delta u_{i, j}^{n+1} = \delta u_{i, j}^n + K_u \left( \delta u_{i+1, j}^n + \delta u_{i-1, j}^n + \delta u_{i, j+1}^n + \delta u_{i, j-1}^n - 4 \delta u_{i, j}^n \right) + \Delta t f_u \delta u_{i, j}^n + \Delta t f_v \delta v_{i, j}^n }
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \delta v_{i, j}^{n+1} = \delta v_{i, j}^n + K_v \left( \delta v_{i+1, j}^n + \delta v_{i-1, j}^n + \delta v_{i, j+1}^n + \delta v_{i, j-1}^n - 4 \delta v_{i, j}^n \right) + \Delta t g_u \delta u_{i, j}^n + \Delta t g_v \delta v_{i, j}^n }
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle {\delta u^n \choose \delta v^n} = \xi ^n e^{i(k_x i \Delta x + k_y j \Delta y)} {\delta u_0 \choose \delta v_0} }

Para que o método seja estável, é preciso que Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle |\xi| < 1} . Após inserir as soluções nas equações linearizadas e realizar manipulações algébricas, como feito por Scholz[7], obtemos

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \xi_\pm = \frac{1}{2} \left [ \beta \pm \sqrt{\beta ^2 - 4 (\Delta t)^2 a^2 b} \right ] }

onde

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \beta = 2 - 4 \gamma (K_u + K_v) + \Delta t (b - 1 - a^2)}

e

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \gamma = sin^2\left( \frac{k_x \Delta s}{2} \right) + sin^2\left( \frac{k_y \Delta s}{2} \right)}

Os valores Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Delta t = 0.01} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Delta s = 1} satisfazem a condição de estabilidade, como mostrado em Scholz[7].

Implementação

A seguir, o trecho do código relativo à implementação do método na linguagem Python. O código completo encontra-se no Github [1] 

# u_n, v_n -> concentracao dos reagentes no tempo n (matriz Nx x Ny)
# u_n1, v_n1 -> concentracao dos reagentes no tempo n+1

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

            # calcular u e v no tempo n+1
            u_n1[i, j] = u_n[i, j] + dt * f(u_n[i, j], v_n[i, j], b) \
                         + 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], b) \
                         + 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])

    # 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

Resultados

Esta seção se dedica ao resultados das simulações realizadas, que mostram o caráter oscilatório do Brusselator ao longo do tempo e também como se dá a difusão dos reagentes em um recipiente em duas dimensões. As constantes utilizadas foram identificadas abaixo para plotar os gráficos, exceto quando indicado o contrário.

Tabela de Constantes e Valores
Símbolo Nome Valor
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle N_x} Dimensão analisada ao longo do eixo x 50
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle N_y} Dimensão analisada ao longo do eixo y 50
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle a} Constante relativa à concentração do reagente A 1
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle b} Constante relativa à concentração do reagente B 1.7
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle D_u} Constante de difusão do reagente U 0.5
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle D_v} Constante de difusão do reagente V 1
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle u_0} Concentração de u no tempo inicial (t=0) 1
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle v_0} Concentração de v no tempo inicial (t=0) 1
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle dt} Passo de tempo entre iterações 0.01
Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle ds} Unidade de avanço dos eixos no espaço 1

As posições das condições iniciais foram escolhidas arbitrariamente, tendo um total de 6 modos que serão explicados na sequência, sendo eles:

  1. Condição em formato do sinal "+".
  2. Condição de borda.
  3. Condição aleatória.
  4. Condição de nove pontos centrais.
  5. Condição da borda completa.

Simulação

Nas figuras abaixo temos a variação das concentrações dos reagentes ao longo do tempo e os diagramas de fase. Para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle a = 1} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle b = 1.7} notamos que a solução do sistema converge para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle u(t) = 1} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle v(t) = 1.7} . Esse resultado está de acordo com o esperado, pois o ponto Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \left(1, 1.7\right)} é um ponto fixo atrator (satisfaz que Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle 1-b+a^2 > 0} ), conforme discutido na seção de análise de estabilidade do Brusselator. Notamos também, pelo diagrama de fase, que mesmo alterando as condições iniciais do problema, a solução converge para o estado estacionário.


Figura 1: Solução do Brusselator para a = 1 e b = 1.7.


Para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle a = 1} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle b = 3} , temos que Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle 1-b+a^2 < 0} , portanto a solução do sistema não converge. O que observamos é um comportamento periódico causado pela existência do ciclo limite.


Figura 2: Solução do Brusselator para a = 1 e b = 3.

Podemos pensar nas funções de concentração Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle u} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle v} como ondas que sofrem interferência destrutiva, pois estão fora de fase. A animação abaixo mostra este comportamento, observe o reagente Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle U} , representado pela cor amarela, no gráfico da direita. Note que quando o ponto do gráfico da esquerda chega ao mínimo local, o valor de Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle U} é alto e o de Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle V} muito baixo, fazendo com que o gráfico da direita fique muito amarelo.


Figura 3: Comportamento do reagente.


É interessante pensar no gráfico da esquerda como a análise em um ponto específico dentro do recipiente, por exemplo, se Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle Nx=Ny=50} estamos analisando o ponto Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle (10,10)} . Tendo isso em mente conseguimos analisar os outros 3 cantos da figura, já que, por simetria da posição inicial escolhida, os cantos são iguais.


Figura 4: Análise de um ponto.

Condições Iniciais

Realizamos as simulações variando espacialmente as condições iniciais, para verificar os diferentes padrões formados pelo sistema.

Condição em Formato do Sinal "+"

O reagente Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle u_0} é distribuído no centro formando um sinal de "+".

Figura 4: Comportamento em formato de "+" para o reagente u. Para o reagente v foi utilizado a condição de borda.
# u_n [Nx, Ny]

# Centro
u_n[int(Nx / 2), int(Ny / 2)] = u0


# Lateral do sinal
#       X                   Y
u_n[int(Nx / 2) + 1, int(Ny / 2)] = u0
u_n[int(Nx / 2) - 1, int(Ny / 2)] = u0


# Altura do sinal
#       X                   Y
u_n[int(Nx / 2), int(Ny / 2) + 1] = u0
u_n[int(Nx / 2), int(Ny / 2) - 1] = u0

Condição de borda

O reagente Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle v_0} foi distribuído nos 4 cantos do recipiente.

Figura 5: Analisando u com condição de borda "+". Para o reagente v foi utilizado a condição de borda. Foi utilizado b=3 para mostrar a diferença na simulação.
# u_n [Nx, Ny]

# Sup Esq
v_n[0, 0] = v0


# Sup Dir
v_n[0, Nx - 1] = v0


# Inf Esq
v_n[Nx - 1, 0] = v0


# Inf Dir
v_n[Nx - 1, Nx - 1] = v0

Condição aleatória

Os reagentes são aleatoriamente distribuídos ao longo do recipiente, a quantidade de focos de reagente também é aleatória para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle u_0} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle v_0} .

Figura 6: Posição aleatória.
from random import *

# u_n [Nx, Ny]

aleatorio = randint(int(Nx / 5)


# Aleatório para u
for _ in range(aleatorio, Nx)):
    u_n[randint(0, size), randint(0, size)] = u0


# Aleatório para v
for _ in range(aleatorio, Nx)):
    v_n[randint(0, size), randint(0, size)] = v0

Condição de nove pontos centrais

No centro divide 9 pontos igualmente espaçados de tamanho Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle N_x/4} , para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle u_0} . Foi utilizado para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle u_0} na figura 3.

Figura 7: Utilizando a condição de borda para v.
# Nove pontos centrais (u0)
mid = int(Nx / 2); mov = int(Nx / 4)

# Meio & direita
u_n[mid, mid] = u0
u_n[mid + mov, mid + mov] = u0
# Meio & esquerda
u_n[mid + mov, mid] = u0
u_n[mid, mid + mov] = u0
# Meio, p/baixo & p/lados
u_n[mid + mov, mid - mov] = u0
u_n[mid - mov, mid - mov] = u0
# Meio, p/cima & p/lados
u_n[mid - mov, mid] = u0
u_n[mid, mid - mov] = u0
u_n[mid - mov, mid + mov] = u0

Condição da borda completa

Completa toda a borda do recipiente com reagente Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle v_0} . Foi utilizado para Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle v_0} na figura 3.

Figura 8: Utilizando a condição de borda para v e b=3.
# v_n = [Nx, Ny]


# Toda a borda (v0)


for i in range(Nx):

    # Completa primeira coluna e primeira linha
    v_n[0, i] = v0
    v_n[i, 0] = v0

    # Completa última coluna e última linha
    v_n[Nx - 1, i] = v0
    v_n[i, Nx - 1] = v0

Referências

  1. ↑ J. Tyson, Some further studies of nonlinear oscillations in chemical systems, J. Chem. Phys. 58 (1973) 3919
  2. ↑ Turing A. M., The chemical basis of morphogenesis. 1953. Bull. Math. Biol. 52, 153, discussion 119 (1990).
  3. ↑ I. Prigogine, R. Lefever, Symmetries breaking instabilities in dissipative systems II. J. Phys. Chem. 48, 1695–1700 (1968)
  4. ↑ R. Lefever and G. Nicolis, Chemical instabilities and sustained oscillations, J. Theor. Biol. 30 (1971)
  5. ↑ G. Adomian, The diffusion-Brusselator equation. Comput. Math. Appl. 29, 1–3 (1995)
  6. ↑ 6,0 6,1 Twizell, E. H., Gumel, A. B., & Cao, Q. (1999). Journal of Mathematical Chemistry, 26(4), 297–316
  7. ↑ 7,0 7,1 Scholz, Christian, Morphology of Experimental and Simulated Turing Patterns, 2009, p. 27-29.
Disponível em “http://fiscomp.if.ufrgs.br/index.php?title=Modelo_Brusselator_de_Reação-Difusão&oldid=7299”