Introdução
Descrição do Modelo
O modelo de Gray-Scott descreve uma reação autocatalítica. Sejam duas substâncias químicas cujas concentrações em um dado ponto do espaço são dadas pelas variáveis
e
, a reação pode ser representada como
Isso significa que uma molécula da substância
é transformada em uma molécula da substância
por meio da ação de outras duas moléculas da substância
, ou seja,
é um catalisador de sua própria produção (daí o termo autocatálise). Além dessa reação, ambas substâncias se difundem pelo meio (por isso esse modelo pertence à classe mais geral de modelos reativos-difusivos) e, portanto, as concentrações
e
mudam com o tempo e diferem em cada ponto. Por simplicidade, assume-se que a reação reversa (i.e.,
) não ocorre.
O comportamento geral do sistema pode ser descrito pelas equações abaixo:
A primeira das equações acima pode ser interpretada da seguinte forma. Em um dado ponto, a variação na concentração
aumenta proporcionalmente ao laplaciano de
naquele ponto, i.e., quando a concentração
na vizinhança desse ponto é alta; e proporcionalmente à taxa
de reposição de
(taxa de alimentação, ou feed rate). A concentração
diminui com o termo reativo
, que representa a reação
.
De outro lado, na segunda das equações acima, a concentração
aumenta com o termo reativo
e também proporcionalmente ao laplaciano de
naquele ponto, mas diminui com a remoção de
a uma taxa
, mais rápida, portanto, do que a reposição de
.
O sistema pode ser imaginado como consistindo em duas substâncias
e
, envoltas por uma membrana semipermeável e imersas em um meio em que essas mesmas duas substâncias estão presentes. A membrana permite a entrada da substância
, mas não da substância
, e permite a saída da substância
, mas não da substância
[1].
Análise de estabilidade
Nota: A análise em toda esta seção pressupõe sempre que os parâmetros e coeficientes de difusão são positivos.
Soluções estacionárias sem difusão
O modelo de Gray-Scott depende dos parâmetros
e dos coeficientes de difusão
das espécies químicas. Ignorando em um primeiro momento os termos de difusão, percebe-se que, por inspeção, o sistema possui uma solução estacionária em
para quaisquer valores dos parâmetros. Esse ponto, no entanto, não é a única solução estacionária do sistema; para encontrar as outras, é necessário impor
nas equações do sistema. Fazendo isso e dispensando os termos de difusão (
), obtém-se o seguinte sistema de equações:
![{\displaystyle {\begin{aligned}-&uv^{2}+F(1-u)=0\\&uv^{2}-(F+k)v=0\quad (2)\\\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7b11de91d0874cfcde3a857b086e88f857e34838)
Somando essas duas equações, relacionamos as variáveis
e
:
![{\displaystyle F(1-u)-(F+k)v=0\Rightarrow u=1-\gamma v\quad (3)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a1c3afcc002e504243e90ee4d55780147e04f780)
onde definiu-se o parâmetro auxiliar
.
Substituindo
na segunda equação do sistema (2) (e reescrevendo
), ficamos com:
![{\displaystyle \left(1-\gamma v\right)v^{2}-\gamma Fv=0\Rightarrow -\gamma v^{3}+v^{2}-\gamma Fv=0\quad (4)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8b30cf254adc9bd565736d0adc73d586f6ab6a94)
Evidentemente,
é solução dessa equação, implicando em
, como já havíamos inspecionado. Alternativamente, considerando
, podemos dividir (4) por
, ficando com
. Resolvendo esta equação quadrática, obtemos duas novas soluções estacionárias para
:
![{\displaystyle v_{\pm }={\frac {1}{2\gamma }}(1\pm {\sqrt {1-4\gamma ^{2}F}})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/06ac57f2e011f42220fa2aff4dcb770028e3f237)
Disso, pela relação (3), temos que os valores correspondentes para
são:
![{\displaystyle u_{\mp }={\frac {1}{2}}(1\mp {\sqrt {1-4\gamma ^{2}F}})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f071866f2e40aed74a63563ce20347da9ec74226)
É necessário apontar que, para que as duas últimas soluções (não-triviais) existam — isto é, sejam números reais — o fator dentro da raiz quadrada tem de ser positivo (
). Por consequência:
, para que existam as soluções não-triviais.
Portanto, há três soluções estacionárias
do sistema:[2]
![{\displaystyle {\begin{aligned}&u_{0}^{*}=1&v_{0}^{*}&=0&\\&u_{1}^{*}={\frac {1}{2}}(1+{\sqrt {1-4\gamma ^{2}F}})&v_{1}^{*}&={\frac {1}{2\gamma }}(1-{\sqrt {1-4\gamma ^{2}F}})&\quad (5)\\&u_{2}^{*}={\frac {1}{2}}(1-{\sqrt {1-4\gamma ^{2}F}})&v_{2}^{*}&={\frac {1}{2\gamma }}(1+{\sqrt {1-4\gamma ^{2}F}})&\\\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/edb5faf0f9cde6fc3fac994762d20562c4df7881)
Estabilidade dos estados estacionários (sem difusão)
Para avaliar a estabilidade das soluções acima, faz-se necessário obter a matriz Jacobiana dos termos de reação,
. Explicitamente, analisando o sistema (1) de equações, temos que
e
. A matriz Jacobiana do sistema é então dada por:
![{\displaystyle J_{R}(u,v)={\begin{bmatrix}\partial _{u}R_{1}&\partial _{v}R_{1}\\\partial _{u}R_{2}&\partial _{v}R_{2}\end{bmatrix}}={\begin{bmatrix}-v^{2}-F&-2uv\\v^{2}&2uv-(F+k)\end{bmatrix}}\quad (6)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ae792c7eeb16dac2e422908aee72a26455c1c399)
Analisemos a estabilidade para os três pares
de soluções estacionárias:
- Para
:
![{\displaystyle J_{R}(u_{0}^{*},v_{0}^{*})={\begin{bmatrix}-F&0\\0&-(F+k)\end{bmatrix}}\quad (7)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0bc3914603b782ac3d8520aa09ba163dab45feda)
- Por essa ser uma matriz diagonal, os autovalores
são justamente as entradas das diagonais; ou seja,
e
. Uma vez que
e
são parâmetros positivos, os dois autovalores são reais e negativos, e portanto o ponto
é sempre estável.
- Para
, podemos utilizar uma estratégia que simplifica as contas. Em particular, nota-se que os dois pontos obedecem à segunda equação do sistema (2) com
. Desse modo, se dividirmos tal equação por
, percebemos que ambos os pontos obedecem a:
![{\displaystyle u_{i}^{*}v_{i}^{*}-(F+k)=0}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6a585eec4fbe7c3d30af86014ff779898af68c9f)
- Dessa equação, podemos calcular as entradas da segunda coluna da matriz jacobiana com facilidade:
![{\displaystyle -2u_{i}^{*}v_{i}^{*}=-2(F+k)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/550e18bb317488d949cfa911a2af0cb249adb60f)
![{\displaystyle 2u_{i}^{*}v_{i}^{*}-(F+k)=(F+k)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d54f976c7beccf7fbc4d4e7dcbef16fe2b3492d4)
- Assim, a matriz jacobiana desses pontos fica:
![{\displaystyle J_{R}(u_{i}^{*},v_{i}^{*})={\begin{bmatrix}-(v_{i}^{*})^{2}-F&-2(F+k)\\(v_{i}^{*})^{2}&(F+k)\end{bmatrix}},\quad i=1,2\quad (8)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d1da2ff0f2bedf72e3d40f1e9f2484af890be038)
- Sabemos que o produto dos autovalores dessa matriz é igual ao seu determinante. Calculando-o, obtém-se:
![{\displaystyle \Delta _{i}:=\operatorname {det} \left(J_{R}(u_{i}^{*},v_{i}^{*})\right)=(F+k)\left[(v_{i}^{*})^{2}-F\right]\quad (9)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b05a6366d96f61f3af1fd4518bd4e80e12fdc964)
- Dividindo por
:
![{\displaystyle {\frac {\Delta _{i}}{F(F+k)}}={\frac {(v_{i}^{*})^{2}}{F}}-1=\left[{\frac {1\mp {\sqrt {1-4(\gamma {\sqrt {F}})^{2}}}}{2(\gamma {\sqrt {F}})}}\right]^{2}-1=\left[{\frac {1\mp {\sqrt {1-4a^{2}}}}{2a}}\right]^{2}-1}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a6dc4e89a50fb14385e69f7ae2244be3caa56cfa)
- onde se definiu
(observação: este é o
definido no Gros). Nota-se que a condição de existência
para os dois pontos não-triviais é equivalente a
. Expandindo os termos, é possível mostrar que a expressão acima pode ser reescrita como:
![{\displaystyle {\frac {\Delta _{i}}{F(F+k)}}={\frac {\sqrt {1-4a^{2}}}{2a^{2}}}\left({\sqrt {1-4a^{2}}}\mp 1\right)\quad (10)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0c9bdc51367898d2353ba8e1397c9427dd709050)
- Para o caso
(sinal negativo em (10)), temos a cota superior
. Portanto,
para todo
que satisfaça a condição de existência. Como o determinante é negativo, sabemos que os autovalores são reais (comentário: como as entradas da matriz são reais, se os autovalores fossem complexos, seriam também conjugados, de modo que o produto deles fosse igual ao módulo ao quadrado de qualquer um, que seria um valor positivo). Ademais, como seu produto é negativo, eles têm sinais opostos; isto é, um deles é positivo, de modo que o ponto
nunca seja estável. Depreendemos desse raciocínio que o determinante da matriz jacobiana de entradas reais ser positivo é uma condição necessária para que haja estabilidade do ponto.
- Já para
(sinal positivo em (10)), temos sempre que
. Para verificar a estabilidade, temos que agora calcular o traço da matriz jacobiana, pois o traço é a soma dos autovalores: se os autovalores são reais, eles têm o mesmo sinal por seu determinante ser positivo, de modo que o traço compartilhe o sinal com os dois autovalores; se os autovalores
são complexos, eles serão conjugados e o traço será
, de modo que a parte real dos autovalores tenha o mesmo sinal do traço. Assim, basta que o traço seja negativo para que o ponto seja estável, e que seja positivo para que seja instável.
- No caso, temos que
. Esse traço é negativo quando
e positivo quando
; ou seja,
é estável quando
e instável quando
(lembrando que
para todo
e
). Desse modo, pode-se caracterizar uma transição de estabilidade quando
.
(...)
Estabilidade dos estados estacionários (com difusão)
Precisamos agora analisar a estabilidade dos pontos estacionários na presença de difusão, como prescreve o sistema de equações (1), que descreve o modelo. Para isso, é necessário levar em consideração, para cada um dos estados de equilíbrio, os autovalores da matriz
, em que
é a matriz diagonal cujas entradas são
e
:[3]
![{\displaystyle D={\begin{bmatrix}D_{u}&0\\0&D_{v}\end{bmatrix}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c4d0b7564740832b583788bb1bcd2ff40f92e96c)
Se escrevermos, genericamente, que
, teremos a seguinte matriz jacobiana de reação-difusão:
![{\displaystyle \left(J_{R}-D\omega ^{2}\right){\Bigg |}_{(u_{i}^{*},v_{i}^{*})}={\begin{bmatrix}a-D_{u}\omega ^{2}&b\\c&d-D_{v}\omega ^{2}\end{bmatrix}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e679937a32845a44e2f7e93931f5e95732f7fcec)
Como já detalhado acima, para que o ponto seja estável, tal matriz tem que ter a parte real de todos os seus autovalores negativa, de modo que seu determinante seja positivo (
) e seu traço negativo (
).[4] Impondo tais condições à matriz acima, obteremos, após manipulações:[3]
![{\displaystyle {\begin{aligned}aD_{v}\omega ^{2}+dD_{u}\omega ^{2}-D_{u}D_{v}\omega ^{4}&<\operatorname {det} (J_{R})\\D_{u}\omega ^{2}+D_{v}\omega ^{2}&>\operatorname {Tr} (J_{R})\\\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4620fefc725d1d2718953874e9c0fc019086e37f)
Se o traço é negativo
Se agora incluímos os termos de difusão , deve-se levar em consideração a matriz
.
Aqui,
é a matriz jacobiana dos termos de reação,
é a matriz diagonal dos termos de difusão e
é o parâmetro que determina a frequência espacial das perturbações. A demonstração da validade desse método pode ser encontrada na referência[5]. Aplicando ao modelo de Gray-Scott em
:
Para que o estado de equilíbrio
seja estável é necessário que o determinante da matriz acima seja positivo e o traço seja negativo. Obtém-se então
Ambas desigualdades são imediatamente satisfeitas para quaisquer valores de
, e
. Portanto, o estado de equilíbrio
permanece estável no modelo de Gray-Scott mesmo após a inclusão dos coeficientes de difusão, sejam quais forem os valores desses coeficientes (lembrando que estamos nos restringindo a valores positivos dos parâmetros e coeficientes).
Esse é um resultado à primeira vista surpreendente. Em geral, o surgimento de padrões complexos e não homogêneos em sistemas reativos-difusivos está relacionado à desestabilização de um ou mais estados de equilíbrio homogêneo causada pela introdução dos coeficientes de difusão (conhecida como instabilidade de Turing)[6].
Entretanto, no caso do modelo de Gray-Scott, o surgimento de padrões complexos e não homogêneos não decorre da instabilidade de Turing, uma vez que o surgimento de padrões não triviais nesse modelo ocorre mesmo quando apenas o estado de equilíbrio trivial
está presente [2].
Implementação
Será usado o método FTCS (Foward Time Central Space) para integrar as equações do modelo. Como existem explicações do método em toda literatura e em outras entradas da Wiki (ver, por exemplo, Modelo de Turing), a explicação aqui será sucinta.
O método consiste em discretizar a derivada parcial em relação ao tempo para frente e discretizar as derivadas parciais de segunda ordem em relação ao espaço centralmente. Para uma função
:
A partir das duas últimas equações acima é fácil mostrar que o laplaciano em duas dimensões, como será usado no presente trabalho, pode ser escrito como
Fazendo
, pode-se simplificar a discretização do laplaciano para
Usando a notação
é possível então escrever as equações do modelo de forma discretizada:
Utilizou-se uma rede quadrada de tamanho
. O estado do inicial do sistema é aquele em que todos os pontos estão no estado de equilíbrio estável trivial
, exceto o ponto central, em que é introduzida uma perturbação com
. Foram usadas condições de fronteira conforme a Figura 1.
Figura 1 - Grid para exemplificar as condições de fronteira usadas na simulação.
Cada elemento na matriz tem quatro vizinhos que são denominados por U (Up), D (Down), L (Left), R (Right). Na Figura 1, o elemento
, possui os vizinhos
,
; o elemento
possui como vizinhos
,
e
; o elemento
tem vizinhos
,
e
; e, finalmente, os vizinhos de
são
,
,
e
.
Essas condições de fronteira e a condição inicial explicada acima buscam reproduzir as mesmas condições usadas na simulação de Sayama[7].
Resultados e discussão
Modelo de Gray-Scott com
|
Concentração de ![{\displaystyle u}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c3e6bb763d22c20916ed4f0bb6bd49d7470cffd8) para ![{\displaystyle (F,k)=(0.015,0.055)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b26db19fe765ecd8f27fe7d2d38f2bf743681eba) , de t=0 até t=2000.
|
Concentração de ![{\displaystyle u}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c3e6bb763d22c20916ed4f0bb6bd49d7470cffd8) para ![{\displaystyle (F,k)=(0.02,0.05)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/72e2f4ce55dc247e885e450b96a9d4be321356a7) , de t=0 até t=2000.
|
Programa
Simulação do Modelo de Gray-Scott
Referências
Bibliografia
- C. Gros, "Complex and Adaptive Dynamical Systems". Springer-Verlag, Berlim, 2015.
- H. Sayama, "Introduction to the Modeling and Analysis of Complex Systems". Open SUNY Textbooks, Geneseo, NY, 2015.