Modelo de Gray-Scott

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 $u$ e $v$, a reação pode ser representada como

$u+2v\to 3v$

Isso significa que uma molécula da substância $u$ é transformada em uma molécula da substância $v$ por meio da ação de outras duas moléculas da substância $v$, ou seja, $v$ é 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 $u$ e $v$ mudam com o tempo e diferem em cada ponto. Por simplicidade, assume-se que a reação reversa (i.e., $3v\to u+2v$) não ocorre. Há reposição de $u$ a uma taxa $F$ (taxa de alimentação, feed rate) e remoção de $v$ a uma taxa ligeiramente mais rápida do que a reposição de $u$.

O comportamento geral do sistema pode ser descrito pelas equações abaixo:

$\begin{array}{rl}\frac{\partial u}{\partial t}& =-u{v}^2+F(1-u)+D_u{\mathrm{\nabla }}^{2}u\\ \frac{\partial v}{\partial t}& =u{v}^2-(F+k)v+D_v{\mathrm{\nabla }}^{2}v(1)\\ \end{array}$

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 $F,k$ e dos coeficientes de difusão $D_u,D_v$ 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 $(u^{*},v^{*})=(1,0)$ 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 $\partial u/\partial t=\partial v/\partial t=0$ nas equações do sistema. Fazendo isso e dispensando os termos de difusão ($D_u=D_v=0$), obtém-se o seguinte sistema de equações:

$\begin{array}{rl}-& u{v}^2+F(1-u)=0\\ & u{v}^2-(F+k)v=0(2)\\ \end{array}$

Somando essas duas equações, relacionamos as variáveis $u$ e $v$:

$F(1-u)-(F+k)v=0\Rightarrow u=1-\gamma v$

onde definiu-se o parâmetro auxiliar $\gamma =\frac{F+k}{F}$.

Substituindo $u$ na segunda equação do sistema (2) (e reescrevendo $F+k=\gamma F$), ficamos com:

$(1-\gamma v)v^2-\gamma Fv=0\Rightarrow -\gamma v^3+v^2-\gamma Fv=0$

Evidentemente, $v=0$ é solução dessa equação, implicando em $u=1$, como já havíamos inspecionado. Alternativamente, considerando $v\ne 0$, podemos dividir a expressão acima por $v$, ficando com $-\gamma v^2+v-\gamma F=0$. Resolvendo esta equação quadrática, obtemos duas novas soluções estacionárias para $v$:

$v_{\pm }=\frac{1}{2\gamma }(1\pm \sqrt{1-4{\gamma }^{2}F})$

Disso, pela relação $u=1-\gamma v$, temos que os valores correspondentes para $u$ são:

$u_{\mp }=\frac{1}{2}(1\mp \sqrt{1-4{\gamma }^{2}F})$

É 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 ( $1-4{\gamma }^{2}F\ge 0$ ). Por consequência:

$4{\gamma }^{2}F\le 1\Rightarrow 4{\left(\frac{F+k}{F}\right)}^{2}F\le 1\Rightarrow F\ge 4(F+k{)}^2$, para que existam as soluções não-triviais.

Portanto, há três soluções estacionárias $(u_i^{*},v_i^{*})$ do sistema:

$\begin{array}{rlr}& 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{array}$

Estabilidade dos estados estacionários

Para avaliar a estabilidade das soluções acima, faz-se necessário obter a matriz Jacobiana dos termos de reação, $R_i(u,v)$. Explicitamente, analisando o sistema (1) de equações, temos que $R_1(u,v)=-u{v}^2+F(1-u)$ e $R_2(u,v)=u{v}^2-(F+k)v$. A matriz Jacobiana do sistema é então dada por:

$J_R(u,v)=\left(\begin{array}{rr}{\partial }_u{R}_1& {\partial }_v{R}_1\\ {\partial }_u{R}_2& {\partial }_v{R}_2\end{array}\right)$

---

Esse estado de equilíbrio é estável porque a matriz jacobiana $J{|}_{(u^{*},v^{*})=(1,0)}=\left(\begin{array}{cc}-F& 0\\ 0& -F-k\end{array}\right)$ possui traço negativo e determinante positivo.

---

Se agora incluímos os termos de difusão $D_u$ e $D_v$, deve-se levar em consideração a matriz $(J-D{\omega }^2){|}_{f=f_{eq}}$. Aqui, $J$ é a matriz jacobiana dos termos de reação, $D$ é a matriz diagonal dos termos de difusão e $\omega$ é 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. Aplicando ao modelo de Gray-Scott em $(u^{*},v^{*})=(1,0)$:

$(\left(\begin{array}{cc}-F-v^2& -2uv\\ v^2& -F-k+2uv\end{array}\right)-\left(\begin{array}{cc}{D}_u& 0\\ 0& D_v\end{array}\right){\omega }^2){|}_{(u^{*},v^{*})=(1,0)}=\left(\begin{array}{cc}-F-D_u{\omega }^2& 0\\ 0& -F-k-D_v{\omega }^2\end{array}\right)$

Para que o estado de equilíbrio $(u^{*},v^{*})=(1,0)$ seja estável é necessário que o determinante da matriz acima seja positivo e o traço seja negativo. Obtém-se então

$(F+D_u{\omega }^2)(F+k+D_v{\omega }^2)>0$

$-2F-(D_u+D_v){\omega }^2-k<0$

Ambas desigualdades são imediatamente satisfeitas para quaisquer valores de $F,k,D_u$, e $D_v$. Portanto, o estado de equilíbrio $(u^{*},v^{*})=(1,0)$ 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).

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 $(u^{*},v^{*})=(1,0)$ está presente .

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 Título do link, 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 $f(x,y,t)$:

$\frac{\partial f}{\partial t}\approx \frac{f(x,y,t+\mathrm{\Delta }t)-f(x,y,t)}{\mathrm{\Delta }t}$

$\frac{{\partial }^{2}f}{\partial x^2}\approx \frac{f(x+\mathrm{\Delta }x,y,t)-2f(x,y,t)+f(x-\mathrm{\Delta }x,y,t)}{\mathrm{\Delta }{x}^2}$

$\frac{{\partial }^{2}f}{\partial y^2}\approx \frac{f(x,y+\mathrm{\Delta }y,t)-2f(x,y,t)+f(x,y-\mathrm{\Delta }y,t)}{\mathrm{\Delta }{y}^2}$

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

${\mathrm{\nabla }}^{2}f(x,y,t)\approx \frac{f(x+\mathrm{\Delta }x,y,t)-2f(x,y,t)+f(x-\mathrm{\Delta }x,y,t)}{\mathrm{\Delta }{x}^2}+\frac{f(x,y+\mathrm{\Delta }y,t)-2f(x,y,t)+f(x,y-\mathrm{\Delta }y,t)}{\mathrm{\Delta }{y}^2}$

Fazendo $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{\Delta }h$, pode-se simplificar a discretização do laplaciano para

${\mathrm{\nabla }}^{2}f(x,y,t)=\frac{f(x+\mathrm{\Delta }h,y,t)+f(x-\mathrm{\Delta }h,y,t)+f(x,y+\mathrm{\Delta }h,t)+f(x,y-\mathrm{\Delta }h,t)-4f(x,y,t)}{\mathrm{\Delta }{h}^2}$

Usando a notação $f(x\pm \mathrm{\Delta }h,y\pm \mathrm{\Delta }h,t\pm \mathrm{\Delta }t)\equiv f_{i\pm 1,j\pm 1}^{n\pm 1}$ é possível então escrever as equações do modelo de forma discretizada:

$u_{i,j}^{n+1}=u_{i,j}^n+[-u_{i,j}^n(v_{i,j}^n{)}^2+F(1-u_{i,j}^n)+D_u\frac{u_{i+1,j}^n+u_{i-1,j}^n+u_{i,j+1}^n+u_{i,j-1}^n-4u_{i,j}^n}{\mathrm{\Delta }{h}^2}]\mathrm{\Delta }t$

$v_{i,j}^{n+1}=v_{i,j}^n+[u_{i,j}^n(v_{i,j}^n{)}^2-(F+k)v_{i,j}^n+D_v\frac{v_{i+1,j}^n+v_{i-1,j}^n+v_{i,j+1}^n+v_{i,j-1}^n-4v_{i,j}^n}{\mathrm{\Delta }{h}^2}]\mathrm{\Delta }t$

Utilizou-se uma rede quadrada de tamanho $69\times 69$ com condições de contorno periódicas. O estado do inicial do sistema é aquele em que todos os pontos estão no estado de equilíbrio estável trivial $(u,v)=(1,0)$, exceto o ponto central, em que é introduzida uma perturbação com $(u,v)=(0,1)$, como simulado por Sayama.

A formação de padrões no modelo depende fortemente não apenas dos parâmetros e coeficientes de difusão, mas também de $\mathrm{\Delta }h$ e $\mathrm{\Delta }t$.