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
![{\displaystyle u+2v\to 3v}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ec1c9bbcb954b57f08fa68235b4ebda25c1a26bd)
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. Há reposição de
a uma taxa
(taxa de alimentação, feed rate) e remoção de
a uma taxa ligeiramente mais rápida do que a reposição de
.
O comportamento geral do sistema pode ser descrito pelas equações abaixo:
![{\displaystyle {\begin{aligned}{\frac {\partial {u}}{\partial {t}}}&=-uv^{2}+F(1-u)+D_{u}\nabla ^{2}u\\{\frac {\partial {v}}{\partial {t}}}&=uv^{2}-(F+k)v+D_{v}\nabla ^{2}v\quad (1)\\\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f1d6d5474c9e3e54de79afe9bf25dad47886cefd)
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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/79e84a3598f6173e5258260920d549acbca0610b)
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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c0b43014a275d7a46ccecde177526bb9d09f4ceb)
Evidentemente,
é solução dessa equação, implicando em
, como já havíamos inspecionado. Alternativamente, considerando
, podemos dividir a expressão acima 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
, 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:[1]
![{\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}})\\&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/c7d9419a34823054dba14aae9b66920f56f08547)
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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b69a5983ff66647bf1d19c85262f3d112dd02f68)
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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/067ba08696c52ddc306666b4b197a7373df14775)
- 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.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/081bbdf733869e9d15630507ec2148aea8180262)
- 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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9a2e3b106e85b5b1c63f4da4d91d7f6bc83d24c2)
- 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)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/76ba804b2a3a7fa00a66d9aa098da69541a0d296)
- Para o caso
(sinal negativo), 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), 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
; ou seja,
é estável quando ![{\displaystyle (v_{2}^{*})^{2}>k}](https://wikimedia.org/api/rest_v1/media/math/render/svg/db0132fcf77b0485a105d68c706164da3b7b2e49)
(...)
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
:[2]
![{\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 (
).[3] Impondo tais condições à matriz acima, obteremos, após manipulações:[2]
![{\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[4]. 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)[5].
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 [1].
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[2].
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.