Sistemas do Tipo Difusão-Reação e Preservação de Pontos Singulares

RESUMO

Motivados por aplicações recentes em computação gráfica, este trabalho apresenta um estudo teórico e computacional de sistemas de difusão-reação baseados no Gradient Vector Flow (GVF), com foco no comportamento do GVF em relação às singularidades do campo inicial. O estudo teórico parte de uma análise local, independente de condições de fronteira. Em seguida, supõe-se condição de fronteira no infinito e usa-se análise de Fourier para estabelecer condições suficientes para preservação do ponto singular. Finalmente, supõe-se um domínio compacto, com geometria retangular, e analisa-se a preservação de um ponto singular em relação à condição de fronteira usando um método de solução de equações diferenciais parciais (EDPs) baseado em wavelets de Haar. Desenvolvemos também uma implementação de um método direto para a equação estacionária do GVF baseado em diferenças finitas (DF) para comparar com a solução tradicional do Euler explícito, no que diz respeito a singularidade. É discutida a influência da vorticidade no problema de interesse usando a função de linhas de corrente e equação de Helmholtz. Nos experimentos computacionais, consideramos duas condições de fronteira, dois tipos de singularidades e os três métodos numéricos (Euler explícito, diferenças finitas para a equação estacionária, e wavelets) para verificar os resultados teóricos obtidos.

Palavras-chave:
difusão-reação; Gradient Vector Flow; singularidades; wavelets de Haar

ABSTRACT

Motivated by recent computer graphics applications, this work presents a theoretical and computational study of diffusion-reaction systems based on Gradient Vector Flow (GVF), focusing on the behavior of GVF concerning the singularities of the initial vector field. The theoretical study starts from a local analysis, regardless of boundary conditions. Then, the boundary condition at infinity is assumed, and Fourier analysis is used to establish sufficient conditions for preserving the singular point. Finally, a compact domain with rectangular geometry is assumed. The preservation of a singular point concerning the boundary condition is analyzed using a method for solving partial differential equations (PDEs) based on Haar wavelets. We have also developed an implementation of a direct method for the GVF stationary equation based on finite differences (DF) to compare with the traditional explicit Euler solution with respect to singularity. The influence of vorticity on the problem of interest is discussed using the streamlines function and the Helmholtz equation. In the computational experiments, we consider two boundary conditions, two types of singularities, and the three numerical methods (explicit Euler, finite differences for the stationary equation, and wavelets) to verify the theoretical results obtained.

Keywords:
diffusion-reaction; Gradient Vector Flow; singularities; Haar wavelets

1 INTRODUÇÃO

Sistemas de difusão-reação baseados em equações diferenciais parciais (EDPs) oferecem modelos teóricos para o estudo de fenômenos em física e engenharia ⁶. Um caso particular destes modelos é o Gradient Vector Flow (GVF) proposto originalmente em ¹⁵ para resolver problemas na área de processamento de imagens. O GVF é uma equação vetorial de difusão-reação, cuja solução estacionária é usada em ¹⁵ para guiar modelos deformáveis em algoritmos de segmentação de imagens digitais (veja também ⁴ e referências citadas).

Mais recentemente, em ⁸, ⁹ e ¹⁰ é discutida a utilização do GVF em outra aplicação na área de computação gráfica: inicialização de simulações de fluidos para animação computacional. Nestes trabalhos, permite-se que um usuário forneça um rascunho do comportamento do fluido, de maneira que ele possa informar graficamente as características desejadas para o escoamento. Este rascunho, denominado sketch, deverá passar por etapas de processamento, compostas por filtragem de curvas e técnicas de difusão-reação para campos vetoriais em ℝ² (ou ℝ³), com o objetivo de gerar um campo inicial, dado pela solução estacionária do GVF, com as propriedades matemáticas necessárias para ser usado no modelo de simulação do fluido.

Assim, uma questão fundamental nesta aplicação é a preservação dos pontos singulares (pontos onde o campo de velocidades se anula) definidos pelo usuário, na solução estacionária da equação de difusão-reação. Este fato é importante visto que tais singularidades são fundamentais para caracterizar as diferentes topologias dos campos vetoriais em ℝ^{n (13}.

Estudos envolvendo singularidades e modelos de difusão-reação são conhecidos no contexto da equação de Fourier ⁵, teoria de perturbação ¹⁴, cinemática de sistemas mecânicos ¹⁴, dentre outros. Porém, não temos conhecimento de trabalhos analisando preservação de singularidades em soluções de modelos de EDPs para sistemas de difusão-reação. Desta forma, o estudo teórico e computacional apresentado aqui é a principal contribuição deste trabalho.

Na linha teórica, nosso trabalho parte de uma análise local, independente de condições de fronteira. Em seguida, supõe-se condição de fronteira no infinito e usa-se análise de Fourier para estabelecer condições suficientes para preservação de singularidades. Finalmente, supõe-se um domínio compacto, com geometria retangular, e analisa-se a preservação de um ponto singular no interior do domínio em relação à condição de fronteira. Neste caso, usa-se o método para a solução de EDPs baseado em wavelets de Haar ¹², para parametrizar o comportamento da solução do GVF em relação à condição de fronteira. Ainda no contexto numérico, apresentamos um método baseado em diferenças finitas (DF) para a equação estacionária do GVF. Até onde pudemos verificar na literatura, não existem trabalhos propondo tal solução para o GVF estacionário, embora se trate apenas da aplicação de uma metodologia numérica tradicional para solução de equações diferenciais parciais elípticas.

Em seguida, analisamos a influência da vorticidade na preservação e/ou criação de singularidades. Primeiramente, observamos que a caracterização das singularidades pode ser feita através do estudo dos pontos críticos da função de corrente, caso essa função exista. Esse fato nos leva a uma equação característica envolvendo a vorticidade, permitindo analisar sua influência na geração de novas singularidades na solução do GVF.

Nos experimentos computacionais, utilizamos condições de fronteira homogênea e não homogênea, dois campos vetoriais com diferentes singularidades e três métodos numéricos (Euler explícito, diferenças finitas para a equação estacionária, e wavelets) para verificar numericamente os resultados teóricos obtidos para o GVF, no que diz respeito a preservação de pontos singulares do campo inicial.

O texto está organizado como segue. Na Seção 2 revisamos o GVF e apresentamos uma análise local referente a preservação de singularidades isoladas. Esse problema é retomado na Seção 3 usando recursos da análise de Fourier. Na Seção 4 o problema de preservação de pontos singulares é estudado em domínio compacto usando o método de wavelets. A solução do GVF via DF é apresentada na Seção 5. A Seção 6 discute o relaxamento da preservação da singularidade por questões numéricas. Na Seção 7 discutimos singularidades e vorticidade via função de corrente. Os experimentos computacionais são apresentados na Seção 8. Finalmente, apresentamos as conclusões e trabalhos futuros na Seção 9.

2 MÉTODO DO GVF

O GVF é um modelo de difusão-reação proposto em ¹⁵, segundo a equação:

$\frac{\partial \mathbf{v}}{\partial t}=g∆\mathbf{v}-h\left(\mathbf{v}\mathbf{-}{\mathbf{v}}_0\right),$ (2.1)

onde, no caso do ℝ², a função desconhecida $\mathbf{v}\mathbf{=}\left(v_1\left(x_1,x_2,t\right),v_2\left(x_1,x_2,t\right)\right)$ é um campo vetorial, supostamente suave em uma região $\mathcal{R}\times {\mathrm{ℝ}}^{+}$, com $\mathcal{R}\subset {\mathrm{ℝ}}^2$, e ${\mathbf{v}}_0=\left(f_1\left(x_1,x_2\right),f_2\left(x_1,x_2\right)\right)$ é um campo vetorial estacionário dado, definido em ℛ. Na equação (2.1) o símbolo ∆ denota o operador Laplaciano: $∆={\partial }^2/\partial x^2+{\partial }^2/\partial y^2$. Neste trabalho, vamos considerar g e h constantes reais estritamente positivos. A equação (2.1) fica sujeita as condições iniciais e de contorno, dadas, respectivamente, por:

${\mathbf{v}}_{\mathcal{R}\times \left\{0\right\}}={\mathbf{v}}_0\left(x_1,x_2\right),{\mathbf{v}}_{\partial \mathcal{R}}=\mathit{\vartheta },$ (2.2)

sendo ∂ℛ a fronteira do domínio ℛ, v _∂ℛ a velocidade na fronteira e ϑ uma função vetorial suave em ∂ℛ. A equação (2.1) é do tipo elíptica, a existência e unicidade de sua solução estacionária, com $g,h\in {\mathrm{ℝ}}_{+}^{\star }$ e sujeita as condições dadas pelas expressões em (2.2), são demonstradas na referência ¹⁶.

No texto que segue, o campo vetorial $\mathbf{v}:\mathcal{R}\subset {\mathrm{ℝ}}^2⟶{\mathrm{ℝ}}^2$ representa um campo de velocidades. Suponha que o campo inicial v ₀ tem um ponto singular isolado na origem, ou seja:

${\mathbf{v}}_0\left(0,0\right)=\left(0,0\right).$ (2.3)

A pergunta fundamental neste trabalho é: a solução estacionária do problema definido pelas equações (2.1)-(2.2) preserva o ponto singular isolado? Ou seja, a solução v do problema de contorno:

$-g∆\mathbf{v}+h\left(\mathbf{v}\mathbf{-}{\mathbf{v}}_0\right)=\mathbf{0},$ (2.4)

sujeito a:

${\mathbf{v}}_{\partial \mathcal{R}}=\mathit{\vartheta },$ (2.5)

satisfaz a condição $\mathbf{v}\left(0,0\right)=\left(0,0\right)$?

Os pontos singulares isolados de um campo vetorial são fundamentais para a caracterização das propriedades topológicas do mesmo. Assim, se o usuário desenha um campo vetorial com um ponto singular isolado em um ponto $\mathbf{p}\in \mathcal{R}$, no contexto da metodologia apresentada em ⁸, a solução estacionária v do problema definido pelas equações (2.4)-(2.5) deve ter também ponto singular em p ou em sua vizinhança, relaxando um pouco a pergunta acima. Em relação a esse fato, uma vez que $\mathbf{v}\mathbf{=}\left(v_1\left(x_1,x_2\right),v_2\left(x_1,x_2\right)\right)$ é solução da equação (2.4) para todo $\left(x_1,x_2\right)\in {\mathrm{ℝ}}^2$, então, particularmente para $\left(x_1,x_2\right)=\left(0,0\right)$ teremos:

$-g∆\mathbf{v}\left(0,0\right)+h\left(0,0\right)·\mathbf{v}\left(0,0\right)=h\left(0,0\right)·{\mathbf{v}}_0\left(0,0\right).$ (2.6)

Uma vez que $h\in {\mathrm{ℝ}}_{+}^{\star }$, podemos escrever:

$\mathbf{v}\left(0,0\right)-{\mathbf{v}}_0\left(0,0\right)=\frac{g∆\mathbf{v}\left(0,0\right)}{h\left(0,0\right)}.$ (2.7)

Assim, se $∆\mathbf{v}\left(0,0\right)=0$ então $\mathbf{v}\left(0,0\right)={\mathbf{v}}_0\left(0,0\right)$ e a singularidade é preservada pela solução do GVF. Apesar dessa análise local ser limitada, uma vez que coloca condições sobre a solução do problema, ela será explorada nos resultados computacionais.

Por outro lado, a possibilidade de criação de novas singularidades isoladas deve ser investigada, pois, dependendo de como isto ocorrer, o GVF pode fornecer uma inicialização totalmente diferente daquela indicada pelo usuário. Este tema está relacionado com o comportamento da vorticidade (ξ), definida em ℝ³ pela expressão:

$\xi \equiv \nabla \times v=\left(\frac{\partial v_3}{\partial x_2}-\frac{\partial v_2}{\partial x_3},\frac{\partial v_1}{\partial x_3}-\frac{\partial v_3}{\partial x_1},\frac{\partial v_2}{\partial x_1}-\frac{\partial v_1}{\partial x_2}\right),$ (2.8)

onde ∇× representa o operador rotacional. Se $\nabla \times {\mathbf{v}}_0=\mathbf{0}$ (campo irrotacional), mas a solução do GVF for um campo com vorticidade não-nula (campo rotacional), então poderemos ter criado singularidades no domínio ℛ. Estas questões são importantes para a fundamentação teórica da metodologia descrita em ⁸. As seções a seguir discutem estes problemas.

3 ANÁLISE NO DOMÍNIO DA FREQUÊNCIA

Nesta seção, é apresentada uma análise utilizando a transformada de Fourier que permite obter condições suficientes sobre o campo v ₀ para que a singularidade seja preservada. Primeiramente, notemos que a equação vetorial (2.4) gera duas equações desacopladas, uma para v ₁ e outra para v ₂. Supondo condição de fronteira no infinito para a equação (2.4) e que existem as transformadas de Fourier de v e v ₀, podemos aplicar a transformada de Fourier nos dois membros da equação (2.4) para encontrar a expressão a seguir:

$g\left({\omega }_1^2{V}_i\left({\omega }_1,{\omega }_2\right)+{\omega }_2^2{V}_i\left({\omega }_1,{\omega }_2\right)\right)+h{V}_i\left({\omega }_1,{\omega }_2\right)-h{F}_i\left({\omega }_1,{\omega }_2\right)=0,i=1,2,$ (3.1)

onde $V_1\left({\omega }_1,{\omega }_2\right),V_2\left({\omega }_1,{\omega }_2\right)$ são as transformadas de Fourier para v ₁, v ₂, e F ₁, F ₂ são as transformadas de Fourier para f ₁, f ₂, respectivamente. Portanto, podemos isolar V _i na equação (3.1) e obter:

$V_i\left({\omega }_1,{\omega }_2\right)=\frac{F_i\left({\omega }_1,{\omega }_2\right)}{\frac{g}{h}\left({\omega }_1^2+{\omega }_2^2\right)+1},i=1,2,$ (3.2)

o que mostra que a solução estacionária do GVF, com coeficientes g e h constantes, é obtida por um filtro passa-baixa aplicado sobre o campo inicial. Esse filtro, no domínio da frequência, é dado pela expressão:

$H\left({\omega }_1,{\omega }_2\right)=\frac{1}{\frac{g}{h}\left({\omega }_1^2+{\omega }_2^2\right)+1}.$ (3.3)

Portanto, a solução estacionária do GVF no domínio do espaço pode ser obtida calculando a transformada inversa de Fourier para a expressão (3.2), ou seja:

$v_i\left(x_1,x_2\right)=\int \int H\left({\omega }_1,{\omega }_2\right)·F_i\left({\omega }_1,{\omega }_2\right)exp\left[2\pi j\left(x{\omega }_1+y{\omega }_2\right)\right]d{\omega }_{1}d{\omega }_2.$ (3.4)

Consequentemente, se desejamos obter $v_1\left(0,0\right)=v_2\left(0,0\right)=0$, então a condição:

$v_i\left(0,0\right)=\int \int H\left({\omega }_1,{\omega }_2\right)·F_i\left({\omega }_1,{\omega }_2\right)d{\omega }_{1}d{\omega }_2=0,i=1,2,$ (3.5)

deve ser satisfeita, onde o domínio de integração é o suporte do integrando.

Por exemplo, se ${\mathbf{v}}_0=\left(f_1\left(x_1,x_2\right),f_2\left(x_1,x_2\right)\right)$ é derivado de um potencial, ou seja, existe uma função f tal que:

$f_1\left(x_1,x_2\right)=\frac{\partial f\left(x_1,x_2\right)}{\partial x},f_2\left(x_1,x_2\right)=\frac{\partial f\left(x_1,x_2\right)}{\partial y},$ (3.6)

então, pelas propriedades da transformada de Fourier, podemos escrever:

$F_i\left({\omega }_1,{\omega }_2\right)=j{\omega }_{i}F\left({\omega }_1,{\omega }_2\right),i=1,2,$ (3.7)

onde $F\left({\omega }_1,{\omega }_2\right)$ é a transformada de Fourier do potencial f . Substituindo as equações dadas pela expressão (3.7) na equação (3.5) obtemos:

$v_i\left(0,0\right)=j\int \int H\left({\omega }_1,{\omega }_2\right)·{\omega }_{i}F\left({\omega }_1,{\omega }_2\right)d{\omega }_{1}d{\omega }_2,i=1,2.$ (3.8)

Por exemplo, seja o potencial Gaussiano:

$g_{\sigma }\left(x_1,x_2\right)=\frac{1}{4\pi \sigma }exp\left(-\frac{x_1^2+x_2^2}{4\sigma }\right).$ (3.9)

onde σ é um parâmetro real positivo. Neste caso, $f=g_{\sigma }$ nas equações em (3.6).

Usando as expressões em (3.7) bem como o fato da transformada de Fourier da Gaussiana (3.9) ser dada por $F\left({\omega }_1,{\omega }_2\right)=exp\left[-\sigma \left({\omega }_1^2+{\omega }_2^2\right)\right]$ (ver ³), podemos escrever as expressões em (3.7), para o potencial Gaussiano (3.9), como:

$F_i\left({\omega }_1,{\omega }_2\right)=j{\omega }_{i}exp\left[-\sigma \left({\omega }_1^2+{\omega }_2^2\right)\right],i=1,2.$ (3.10)

Logo, a partir da equação (3.10) pode concluir-se que a expressão (3.8) é nula, uma vez que $F\left({\omega }_1,{\omega }_2\right)$, bem como $H\left({\omega }_1,{\omega }_2\right)$ dada pela equação (3.3), são funções simétricas em relação à origem. Em outras palavras, a singularidade na solução do GVF é preservada.

4 GVF EM DOMÍNIO COMPACTO

Nessa seção resolveremos o GVF utilizando o método de solução de EDPs via wavelets apresentado em ¹². Seja a equação estacionária do GVF dada pela expressão (2.4), no domínio $\mathcal{R}=[-1,1]\times [-1,1]$. Uma vez que a expressão vetorial (2.4) resulta em duas EDPs desacopladas, vamos tratar sua forma geral, dada por:

$-g\left(\frac{{\partial }^{2}v}{\partial x_1^2}+\frac{{\partial }^{2}v}{\partial x_2^2}\right)+hv=hf,$ (4.1)

onde $g,h\in {\mathrm{ℝ}}_{+}^{\star },f:\mathcal{R}\subset {\mathrm{ℝ}}^2⟶\mathrm{ℝ}$ é uma função dada, e $v:\mathcal{R}\subset {\mathrm{ℝ}}^2⟶\mathrm{ℝ}$ está sujeita as condições de contorno:

$v\left(x_1,-1\right)={\gamma }_1^1\left(x_1\right),v\left(x_1,1\right)={\gamma }_2^1\left(x_1\right),v\left(-1,x_2\right)={\gamma }_1^2\left(x_2\right),v\left(1,x_2\right)={\gamma }_2^2\left(x_2\right),$ (4.2)

onde ${\gamma }_i^j$ são funções escalares suaves em ∂ℝ.

Assim, seguindo a metodologia apresentada em ¹², vamos construir inicialmente as funções de Haar para representar a solução do problema de contorno. Seja um intervalo genérico $\left[A,B\right]\subset \mathrm{ℝ}$, com $A\ne B$, e as sequências $j=0,1,\dots ,J,k=0,1,\dots ,\left(2^j-1\right),i=2^j+k+1$. Além disto, consideremos as expressões:

${\xi }_1\left(j,k\right)=A+\frac{k}{2^j}\left(B-A\right),{\xi }_{s+1}\left(j,k\right)={\xi }_s\left(j,k\right)+\frac{\left(B-A\right)}{2^{j+1}},s=1,2.$

Desta forma, podemos definir a função de escala (h ₁) e as wavelet de Haar $(h_i,i>1)$ em $\left[A,B\right]\subset \mathrm{ℝ}$ como:

$h_1\left(t\right)=\left\{\begin{array}{l}1\text{para}t\in [A,B],\\ 0\text{caso contrário}.\end{array},h_i\left(t\right)=\right.\left\{\begin{array}{c}\begin{array}{l}1\text{para}t\in [{\xi }_1\left(j,k\right),{\xi }_2\left(j,k\right)),\\ -1\text{para}t\in [{\xi }_2\left(j,k\right),{\xi }_3\left(j,k\right)],\\ 0\text{caso contrário}.\end{array}\end{array}\right.$ (4.3)

Para $i=2,3,···,2^{J+1}$.

O próximo passo utilizado em ¹² assume que a solução v do problema de valor de contorno definido pelas expressões (4.1)-(4.2) satisfaz:

$\frac{{\partial }^{4}v}{\partial x_1^2{x}_2^2}=\sum _{i,l=1}^{2^{J+1}}{a}_{il}{h}_i\left(x_1\right)h_l\left(x_2\right),$ (4.4)

onde as funções h _i são funções de Haar definidas em (4.3), particularizadas para o intervalo [−1, 1] uma vez que o domínio de interesse é a região $\mathcal{R}=[-1,1]\times [-1,1]$, e a _il são coeficientes reais a determinar. Assim, integrando a equação (4.4) duas vezes em relação às variáveis x ₁ e x ₂ obtemos:

$\frac{{\partial }^{2}v}{\partial x_j^2}=\sum _{i,l}{a}_{il}{h}_i\left(x_j\right)p_{2l}\left(x_k\right)+\left(x_k+1\right)\frac{d^2{\phi }_1^j\left(x_j\right)}{d{x}_j^2}+\frac{d^2{\phi }_2^j\left(x_j\right)}{d{x}_j^2},$ (4.5)

onde $(j,k)\in \left\{\right(1,2),(2,1\left)\right\},{\phi }_1^j,{\phi }_2^j$ são funções que surgem no processo de integração, e p _2l é dado por ¹¹:

$p_{2l}\left(z\right)={\int }_{-1}^z\left[{\int }_{-1}^x{h}_l\left(t\right)dt\right]dx.$ (4.6)

Integrando duas vezes novamente obtém-se:

$$\begin{gathered} v\left(x_1,x_2\right)= \\ \sum _{i,l=1}^{2^{J+1}}{a}_{il}{p}_{2i}\left(x_1\right)p_{2l}\left(x_2\right)+\left(x_2+1\right){\phi }_1^1\left(x_1\right)+{\phi }_2^1\left(x_1\right)+\left(x_1+1\right){\phi }_1^2\left(x_2\right)+{\phi }_2^2\left(x_2\right). \end{gathered}$$ (4.7)

Vamos agora impor sobre a expressão (4.7) as condições de contorno dadas pela expressão (4.2). A primeira restrição $v\left(x_1,-1\right)={\gamma }_1^1\left(x_1\right)$ gera a equação:

$\sum _{i,l}{a}_{il}{p}_{2i}\left(x_1\right)p_{2l}\left(-1\right)+{\phi }_2^1\left(x_1\right)+\left(x_1+1\right){\phi }_1^2\left(-1\right)+{\phi }_2^2\left(-1\right)={\gamma }_1^1\left(x_1\right).$

Porém, pela expressão (4.6) é imediato que $p_{2l}\left(-1\right)=0$, logo, isolando ${\phi }_2^1\left(x_1\right)$ na equação acima, obtemos:

${\phi }_2^1\left(x_1\right)=-\left(x_1+1\right){\phi }_1^2\left(-1\right)-{\phi }_2^2\left(-1\right)+{\gamma }_1^1\left(x_1\right).$ (4.8)

Analogamente, usando agora a condição de contorno $v\left(-1,x_2\right)={\gamma }_1^2\left(x_2\right)$ na expressão (4.7), podemos isolar ${\phi }_2^2\left(x_2\right)$ obtendo:

${\phi }_2^2\left(x_2\right)=-\left(x_2+1\right){\phi }_1^1\left(-1\right)-{\phi }_2^1\left(-1\right)+{\gamma }_1^2\left(x_2\right).$ (4.9)

Substituindo em (4.9) o valor de ${\phi }_2^1\left(-1\right)$ computado pela expressão (4.8) obtemos:

${\phi }_2^2\left(x_2\right)=-\left(x_2+1\right){\phi }_1^1\left(-1\right)+{\phi }_2^2\left(-1\right)+{\gamma }_1^2\left(x_2\right)-{\gamma }_1^1(-1).$ (4.10)

Seja agora a função (4.7) restrita à condição de fronteira para $x_2=1$:

$v\left(x_1,1\right)=\sum _{i,l}{a}_{il}{p}_{2i}\left(x_1\right)p_{2l}\left(1\right)+2{\phi }_1^1\left(x_1\right)+{\phi }_2^1\left(x_1\right)+\left(x_1+1\right){\phi }_1^2\left(1\right)+{\phi }_2^2\left(1\right)={\gamma }_2^1\left(x_1\right).$

A ideia agora é isolar o termo ${\phi }_1^1\left(x_1\right)$ na expressão acima, mas sem a dependência da função ${\phi }_2^1\left(x_1\right)$. Para isso, substituímos na expressão acima o termo ${\phi }_2^1\left(x_1\right)$ dado pela função (4.8). Desta forma, e considerando que a função em (4.10) satisfaz $-{\phi }_2^2\left(1\right)+{\phi }_2^2(-1)=2{\phi }_1^1(-1)-{\gamma }_1^2\left(1\right)+{\gamma }_1^2(-1)$, obtemos:

${\phi }_1^1\left(x_1\right)=\frac{1}{2}\left(-\sum _{i,l}{a}_{il}{p}_{2i}\left(x_1\right)p_{2l}\left(1\right)+(x_1+1)\left({\phi }_1^2(-1)-{\phi }_1^2\left(1\right)\right)+2{\phi }_1^1(-1)+{\gamma }_2^1\left(x_1\right)-{\gamma }_1^1\left(x_1\right)-{\gamma }_1^2\left(1\right)+{\gamma }_1^2(-1)\right).$ (4.11)

Aplicando a última restrição $v\left(1,x_2\right)={\gamma }_2^2\left(x_2\right)$, da condição de fronteira em (4.2), na função (4.7), obtêm-se:

$\sum _{i,l}{a}_{il}{p}_{2i}\left(1\right)p_{2l}\left(x_2\right)+\left(x_2+1\right){\phi }_1^1\left(1\right)+{\phi }_2^1\left(1\right)+\left(1+1\right){\phi }_1^2\left(x_2\right)+{\phi }_2^2\left(x_2\right)={\gamma }_2^2\left(x_2\right).$

Analogamente ao desenvolvimento anterior, vamos isolar o termo ${\phi }_1^2\left(x_2\right)$. A expressão (4.10) permite evitar a dependência de ${\phi }_1^2\left(x_2\right)$ em relação a ${\phi }_2^2$. A partir do resultado obtido, e usando a função (4.11), obtêm-se:

$$\begin{gathered} {\phi }_1^2\left(x_2\right)=\frac{1}{2}(-\sum _{i,l}{a}_{il}\left[p_{2i}\left(1\right)\left(p_{2l}\left(x_2\right)-\frac{(x_2+1)}{2}{p}_{2l}\left(1\right)\right)\right](x_2+1\left)\right(-{\phi }_1^2(-1)+{\phi }_1^2\left(1\right))+ \\ 2{\phi }_1^2(-1)+{\gamma }_2^2(x_2)-{\gamma }_1^1(1)+{\gamma }_1^1(-1)-{\gamma }_1^2(x_2)+\frac{(x_2+1)(-{\gamma }_2^1(1)+{\gamma }_2^1(-1)+{\gamma }_1^1(1)-{\gamma }_1^1(-1\left)\right)}{2}). \end{gathered}$$ (4.12)

Observando as equações (4.8), (4.10), (4.11), (4.12) é direto mostrar que, para $j\in \{1,2\}$:

$\frac{d^2{\phi }_1^j\left(x_j\right)}{d{x}_j^2}=\frac{1}{2}\left(-\sum _{i,l}{a}_{il}{h}_i\left(x_j\right)p_{2l}\left(1\right)+\frac{d^2{\gamma }_2^j\left(x_j\right)}{d{x}_j^2}-\frac{d^2{\gamma }_1^j\left(x_j\right)}{d{x}_j^2}\right),$ (4.13)

$\frac{d^2{\phi }_2^j\left(x_j\right)}{d{x}_j^2}=\frac{d^2{\gamma }_1^j\left(x_j\right)}{d{x}_j^2}.$ (4.14)

Substituindo os resultados (4.13)-(4.14) nas equações (4.5) obtemos:

$\frac{{\partial }^{2}v}{\partial x_j^2}=\sum _{i,l}{a}_{il}{h}_i\left(x_j\right)\left(p_{2l}\left(x_k\right)-\frac{(x_k+1)}{2}{p}_{2l}\left(1\right)\right)+\frac{(x_k+1)}{2}\left({\gamma }_2^{j\text{'}\text{'}}\left(x_j\right)-{\gamma }_1^{j\text{'}\text{'}}\left(x_j\right)\right)+{\gamma }_1^{j\text{'}\text{'}}\left(x_j\right),$ (4.15)

onde $(j,k)\in \left\{\right(1,2),(2,1\left)\right\}$. O sistema envolvendo os coeficientes a _il como incógnitas é obtido substituindo as expressões em (4.15) para ${\partial }^{2}v/\partial x_1^2,{\partial }^{2}v/\partial x_2^2$ e a expressão (4.7) para v na equação (4.1). Em seguida, calcula-se a expressão resultante para uma malha de pontos $\left(x_{1,r},x_{2,s}\right)$ no domínio ℛ com resolução $2^{J+1}\times 2^{J+1}$, de acordo a metodologia utilizada para compor as expressões (4.3), para obter um sistema linear do tipo:

$$\begin{gathered} \sum _{i,l}{a}_{ij}\{g\left[h_i\left(x_1\right)\left(p_{2l}\left(x_2\right)-\frac{x_2+1}{2}{p}_{2l}\left(1\right)\right)+\left(p_{2i}\left(x_1\right)-\frac{x_1+1}{2}{p}_{2i}\left(1\right)\right)h_l\left(x_2\right)\right] \\ +h\left[\left(p_{2i}\left(x_1\right)-\frac{x_1+1}{2}{p}_{2i}\left(1\right)\right)\left(p_{2l}\left(x_2\right)-\frac{x_2+1}{2}{p}_{2l}\left(1\right)\right)\right]\}=-hf(x_1,x_2) \\ -g[\frac{x_2+1}{2}\left({\gamma }_2^{1\text{'}\text{'}}\left(x_1\right)-{\gamma }_1^{1\text{'}\text{'}}\left(x_1\right)\right)+\frac{x_1+1}{2}\left({\gamma }_2^{2\text{'}\text{'}}\left(x_2\right)-{\gamma }_1^{2\text{'}\text{'}}\left(x_2\right)\right)+{\gamma }_1^{1\text{'}\text{'}}(x_1)+{\gamma }_1^{2\text{'}\text{'}}(x_2\left)\right] \\ -h[\frac{x_2+1}{2}\left({\gamma }_2^1\left(x_1\right)-{\gamma }_1^1\left(x_1\right)-{\gamma }_1^2\left(1\right)+{\gamma }_1^2(-1)\right)+\frac{x_1+1}{2}\left({\gamma }_2^2\left(x_2\right)-{\gamma }_1^2\left(x_2\right)-{\gamma }_1^1\left(1\right)+{\gamma }_1^1(-1)\right) \\ +\frac{(x_1+1)(x_2+1)}{4}\left(-{\gamma }_2^1\left(1\right)+{\gamma }_2^1(-1)+{\gamma }_1^1\left(1\right)-{\gamma }_1^1(-1)\right)+{\gamma }_1^1(x_1)+{\gamma }_1^2(x_2)-{\gamma }_1^1(-1\left)\right]. \end{gathered}$$ (4.16)

Note que o sistema linear (4.16) pode gerar uma matriz cheia e má condicionada, portanto, além do aumento do custo computacional, para grandes sistemas, o método pode gerar erros numéricos significativos.

Retomando a discussão sobre a equação 2.7, vamos utilizar a condição extra:

$∆v\left(0,0\right)=\frac{{\partial }^{2}v\left(0,0\right)}{\partial x_1^2}+\frac{{\partial }^{2}v\left(0,0\right)}{\partial x_2^2}=0,$ (4.17)

suficiente para garantir a preservação da singularidade no ponto (0, 0).

Computando a equação (4.17) usando as expressões em (4.15) obtêm-se:

$$\begin{gathered} \sum _{i,l}{a}_{il}\left[h_i\left(0\right)\left(p_{2l}\left(0\right)-\frac{1}{2}{p}_{2l}\left(1\right)\right)+\left(p_{2i}\left(0\right)-\frac{1}{2}{p}_{2i}\left(1\right)\right)h_l\left(0\right)\right] \\ +\frac{1}{2}\left({\gamma }_2^{1\text{'}\text{'}}\left(0\right)-{\gamma }_1^{1\text{'}\text{'}}\left(0\right)\right)+\frac{1}{2}\left({\gamma }_2^{2\text{'}\text{'}}\left(0\right)-{\gamma }_1^{2\text{'}\text{'}}\left(0\right)\right)+{\gamma }_1^{1\text{'}\text{'}}\left(0\right)+{\gamma }_1^{2\text{'}\text{'}}\left(0\right)=0. \end{gathered}$$ (4.18)

Essa expressão mostra como as condições de fronteira, representadas aqui pelas derivadas das funções ${\gamma }_1^1\left(x_1\right),{\gamma }_2^1\left(x_1\right),{\gamma }_1^2\left(x_2\right),{\gamma }_2^2\left(x_2\right)$, em $x_1=x_2=0$, influenciam na preservação da singularidade. Numericamente, podemos garantir a preservação do ponto singular impondo a condição (4.18) ao sistema (4.16). Porém, nesse caso teremos mais equações do que incógnitas, não havendo como garantir a existência da solução (ver Seção 9). O desenvolvimento desta seção considerou a singularidade em (0, 0). Porém, devemos notar que os cálculos são válidos para qualquer ponto singular no interior do domínio ℛ.

5 APROXIMAÇÃO POR DIFERENÇAS FINITAS

Nos trabalhos ¹⁵ e ² são propostas soluções para o GVF baseadas em Euler explícito, no caso de ¹⁵, e DF semi-implícita, explícita e implícita no caso de ², para a sua equação de evolução (2.1). Como esse modelo do GVF é fundamentado em uma equação vetorial linear desacoplada e o interesse deste trabalho é somente em sua solução $(\partial \mathbf{v}/\partial \mathbf{x}=0)$, é possível obter soluções mais eficientes para o GVF considerando o modelo estacionário (2.4).

Assim, diferentemente das soluções clássicas descritas em ¹⁵ e ², apresentamos um método de diferenças finitas para a solução estacionária. Sendo a equação (2.4) desacoplada é suficiente fornecer uma solução para a expressão (4.1) com as condições de contorno (4.2), utilizadas na Seção 4. Desta forma escolhemos o esquema de discretização:

$-g\left({\partial }_{x_1}^{-}{\partial }_{x_1}^{+}+{\partial }_{x_2}^{-}{\partial }_{x_2}^{+}\right)v^{i,j}+h{v}^{i,j}=h{v}_0^{i,j};i=0,···,N-1,j=0,···,M-1,$ (5.1)

que resulta em uma matriz pentadiagonal esparsa, onde $v^{i,j}\equiv v\left((x_{1,0}+(\Delta x_1)i,x_{2,0}+(\Delta x_2)j\right)$, sendo $\left(x_{1,0},x_{2,0}\right)$ o ponto inicial no domínio ℛ discreto, $\Delta x_1=x_{1,1}-x_{1,0}$ e $\Delta x_2=x_{2,1}-x_{2,0}$ nas direções x ₁ e x ₂, respectivamente. Para impor as condições de contorno (4.2) para v em ∂ℛ fazemos: $v^{0,j}={{\gamma }_1^2}^j;v^{N-1,j}={{\gamma }_2^2}^j;v^{i,0}={{\gamma }_1^1}^i;v^{i,M-1}={{\gamma }_2^1}^i$.

Sendo a expressão (5.1) uma discretização da equação estacionária, não temos o mesmo problema de convergência como do método de Euler apresentado em ¹⁵. Porém, a qualidade da solução depende do condicionamento da matriz resultante da discretização escolhida. Um bom condicionamento é obtido se a matriz resultante de (5.1) for estritamente diagonal dominante, isto é, o módulo de cada valor da diagonal for estritamente maior que a soma em módulo dos valores da linha ou coluna correspondente. Pode-se mostrar que essa condição é satisfeita para o sistema gerado pela expressão (5.1) se:

$\left|\left(2(\Delta x^2+\Delta y^2)+\Delta x^2\Delta y^2\frac{h}{g}\right)\right|>|-\Delta x^2|+|-\Delta y^2|+|-\Delta y^2|+|-\Delta x^2|=2(\Delta x^2+\Delta y^2),$

o que sempre ocorre para a discretização escolhida, uma vez que $g,h\in {\mathrm{ℝ}}_{+}^{\star }$.

6 RELAXANDO A PRESERVAÇÃO DA SINGULARIDADE

Considerando o domínio contínuo, se p (0) é o ponto singular do campo inicial v ₀ então, por continuidade, podemos supor que existe ${\epsilon }_t>0$ tal que a solução v da equação de evolução (2.1) satisfaz:

$0\le t\le {\epsilon }_t\Rightarrow \mathbf{v}\left(x\left(t\right),y\left(t\right),t\right)=\left(0,0\right).$ (6.1)

Porém, em termos da aplicação de interesse, o desejável é que exista uma curva $\mathit{p}\mathbf{}:{\mathrm{ℝ}}_{+}⟶\mathcal{R}$ tal que:

$||\mathit{p}\left(t=0\right)-\mathit{p}\left(t=\infty \right)||<{\epsilon }_{p}e\mathit{v}\left(\mathit{p}\left(t=\infty \right)\right)=\left(0,0\right),$ (6.2)

onde ║ ║ é uma norma conveniente e ${\epsilon }_p>0$ é um valor que pode ser considerado pequeno no contexto da aplicação; ou seja, o ponto singular sofre apenas um pequeno deslocamento na solução do GVF. A condição (6.2) é mais forte do que a condição (6.1) e, em geral, não pode ser garantida apenas por argumentos de continuidade. Indo agora para o contexto numérico, podemos (de fato devemos) relaxar a condição (6.2) para:

$||\mathit{p}\left(t=0\right)-\mathit{p}\left(t=\infty \right)||<{\epsilon }_{p}e||\mathit{v}\left(\mathit{p}\left(t=\infty \right)\right)||<{\epsilon }_v,$ (6.3)

onde ε _v é um valor pequeno. Esse relaxamento é necessário por questões relacionadas à solução numérica do problema e à aritmética de ponto-flutuante, como será verificado nos experimentos da Seção 8. Dentro do contexto deste trabalho, uma forma que poderia acomodar os pontos de vista contínuo e numérico em um modelo consistente, seria aplicar a metodologia usada na seção 4 para tratar numericamente a equação elíptica (2.1) do GVF, seguindo técnicas aplicadas em ¹². Em seguida, poderíamos verificar experimentalmente se a solução satisfaz as condições (6.3) e a influência das condições de fronteira, uma vez que essas aparecem explicitamente na solução (vide Seção 7 da referência ¹ como exemplo). Neste caso, a expressão (2.7) gera a condição:

$-{\mathit{v}}_0\left(\mathit{p}\left(t=\infty \right)\right)=\frac{g}{h}\Delta \mathbf{v}\left(\mathit{p}\left(t=\infty \right)\right),$ (6.4)

que poderia ser explorada também neste contexto.

7 SINGULARIDADES, FUNÇÃO DE CORRENTE E VORTICIDADE

Seja um campo estacionário $\mathbf{v}:\mathcal{R}\subset {\mathrm{ℝ}}^2⟶{\mathrm{ℝ}}^2$ que satisfaz as condições de incompressibilidade $(\mathbf{\nabla }\mathbf{·}\mathbf{v}=0)$ em ℛ e não-escorregamento $(\mathbf{v}{|}_{\partial \mathcal{R}}=\mathbf{0})$ na fronteira ∂ℛ do domínio. Nestas condições ⁵, este campo pode ser representado por uma função escalar $\varphi :U\to \mathrm{ℝ}$, denominada função de corrente, com as seguintes propriedades:

1. As isolinhas de ϕ são linhas de corrente do fluxo de v;

2. O campo vetorial v pode ser computado por:

$\mathbf{v}\left(x_1,x_2\right)=\left(\frac{\partial \varphi }{\partial x_2},-\frac{\partial \varphi }{\partial x_1}\right).$ (7.1)

Usando a equação (7.1) e a definição (2.8) é direto mostrar que a vorticidade pode ser computada usando a função de corrente pela equação:

$\mathit{\xi }=\left(0,0,-\left(\frac{{\partial }^2\varphi }{\partial x_1^2}+\frac{{\partial }^2\varphi }{\partial x_2^2}\right)\right)=\left(0,0,-∆\varphi \right).$ (7.2)

Primeiramente, observemos que, pela expressão (7.1), um ponto singular de v ocorre em um ponto crítico de ϕ, pois se $\mathbf{v}\left(s_1,s_2\right)=\left(0,0\right)$ então $\nabla \varphi \left(s_1,s_2\right)=\left(\frac{\partial \varphi }{\partial x_1}\left(s_1,s_2\right),\frac{\partial \varphi }{\partial x_2}\left(s_1,s_2\right)\right)=\left(0,0\right)$. Uma vez que as linhas de corrente do campo são as isolinhas de ϕ, podemos caracterizar a topologia do campo v, ou seja, seus pontos singulares, através dos pontos críticos de ϕ. A matriz Hessiana da função de corrente é dada pela expressão:

$H_{\mathbf{x}}\varphi =\left[\begin{array}{cc}\frac{{\partial }^2\varphi }{\partial x_1^2}& \frac{{\partial }^2\varphi }{\partial x_2\partial x_1}\\ \frac{{\partial }^2\varphi }{\partial x_1\partial x_2}& \frac{{\partial }^2\varphi }{\partial x_2^2}\end{array}\right],$ (7.3)

portanto, assumindo que $\frac{{\partial }^2\varphi }{\partial x_2\partial x_1}=\frac{{\partial }^2\varphi }{\partial x_1\partial x_2}$, sua equação $(\det \left(H_{\mathbf{x}}\varphi -tI\right)=0)$ de autovalores é dada por:

$t^2-t\left(\frac{{\partial }^2\varphi }{\partial x_1^2}+\frac{{\partial }^2\varphi }{\partial x_2^2}\right)+\left(\frac{{\partial }^2\varphi }{\partial x_1^2}\frac{{\partial }^2\varphi }{\partial x_2^2}-{\left(\frac{{\partial }^2\varphi }{\partial x_1\partial x_2}\right)}^2\right)=0.$ (7.4)

Observemos que, se a vorticidade (expressão (7.2)) é nula, os autovalores da matriz Hessiana são:

${\lambda }_1=+\sqrt{{\left(\frac{{\partial }^2\varphi }{\partial x_1\partial x_2}\right)}^2-\frac{{\partial }^2\varphi }{\partial x_1^2}\frac{{\partial }^2\varphi }{\partial x_2^2}},e{\lambda }_2=-\sqrt{{\left(\frac{{\partial }^2\varphi }{\partial x_1\partial x_2}\right)}^2-\frac{{\partial }^2\varphi }{\partial x_1^2}\frac{{\partial }^2\varphi }{\partial x_2^2}}.$ (7.5)

Uma vez que a matriz Hessiana (7.3) é simétrica e real, seus autovalores são reais também. Assim, as únicas possibilidades para a expressão (7.5) serão dois autovalores nulos $({\lambda }_1={\lambda }_2=0)$ ou ${\lambda }_1·{\lambda }_2<0$. O primeiro caso é degenerado e o segundo caso corresponde a um ponto de sela (Figura 2.(a)).

Figura 1:
(a) Linhas de corrente do campo vetorial (8.1). (b) Linhas de corrente do campo vetorial (8.2).

Figura 2:
(a) Campo vetorial (8.1). (b) Linhas de corrente da solução do GVF via Euler explícito. (c) Solução do GVF usando DF para a equação estacionária. (d) Visualização da solução do GVF usando método de wavelets.

Por outro lado, se a vorticidade for diferente de zero no ponto singular, então a solução da equação (7.4) poderá ter dois autovalores positivos $({\lambda }_1\ge 0,{\lambda }_2\ge 0)$, negativos ${\lambda }_1\le 0,{\lambda }_2\le 0)$ ou com sinais trocados $({\lambda }_1·{\lambda }_2<0)$. Nos dois primeiros casos, se algum autovalor for nulo, teremos um caso degenerado. Caso contrário o ponto crítico será um mínimo local (autovalores estritamente positivos) ou um máximo local (autovalores estritamente negativos). Em qualquer uma destas situações, as isolinhas de ϕ (linhas de corrente) em uma vizinhança de (s ₁ , s ₂) são curvas fechadas contendo (s ₁ , s ₂) em seu interior. Portanto, teremos uma configuração topologicamente equivalente a um centro. No último caso, temos um ponto de sela, como já comentado anteriormente.

Embora a análise acima seja limitada em função das condições necessárias para a existência da função de corrente, ela indica que a vorticidade está diretamente relacionada com outra questão: geração de novas singularidades na solução estacionária do GVF. Especificamente, seja um campo irrotacional v ₀ com uma singularidade do tipo sela. Porém, se algum processo gerar vorticidade, então o novo campo pode passar para a situação descrita pela expressão (7.4) com novas singularidades em seu interior.

No caso do GVF, a geração de vorticidade pode ser formalmente analisada aplicando o operador rotacional nos dois membros da equação (2.4) obtendo:

$-g∆\left(\mathbf{\nabla }\mathbf{\times }\mathbf{v}\right)+h\left(\nabla \times \mathbf{v}\right)=\nabla \times {\mathbf{v}}_0.$ (7.6)

Supondo que o campo v ₀ seja irrotacional, a expressão (7.6) simplifica para:

$g∆\mathit{\xi }+h\mathit{\xi }=0,$ (7.7)

que é uma equação de Helmholtz homogênea, onde $\mathit{\xi }=\left(\nabla \times \mathrm{v}\right)$.

Dependendo da condição de fronteira para a vorticidade, um campo inicial v ₀ irrotacional pode gerar uma solução da equação do GVF com vorticidade não-nula. Como consequência da discussão acima, podemos ter singularidades na solução da equação (2.4) que não estavam originalmente presentes em v ₀.

8 EXPERIMENTOS COMPUTACIONAIS

Nesta seção, apresentaremos alguns testes computacionais para discutir o comportamento da solução numérica do GVF em relação à preservação das singularidades do campo inicial, para o caso do domínio compacto. Os experimentos foram realizados utilizando implementações próprias na linguagem Python3 (v3.6.7), biblioteca SciPy (v1.5.0) para solução de sistemas lineares, em um notebook pessoal com processador Intel Core i7 8550U, 8GB de memória RAM e placa de vídeo GeForce MX150. Em cada experimento, apresentamos os tempos de CPU dos métodos considerados para permitir a comparação dos mesmos em termos do seu custo computacional. Nestes testes, estaremos considerando os seguintes campos vetoriais, no domínio $\mathcal{R}=[-1,1]\times [-1,1]$.

• Campo vetorial com Ponto de Sela:

${\mathbf{v}}_{01}(x_1,x_2)=(-(x_1-s_1),x_2-s_2),$ (8.1)

• Campo vetorial não-linear:

${\mathbf{v}}_{02}(x_1,x_2)=\left(\right(x_1-s_1{)}^2,(x_2-s_2{)}^2-(x_1-s_1)/2),$ (8.2)

As linhas de corrente destes campos podem ser visualizadas nas Figuras 1.(a) e 1.(b), onde utilizamos o software descrito em ⁷ para gerar as visualizações. Primeiramente, vamos utilizar o campo irrotacional tradicional, definido pela equação (8.1), que possui Laplaciano nulo e apenas um ponto singular. Uma verificação direta mostra que o campo (8.1) é solução da equação (2.4) do GVF em ℝ². Pela equação (2.7), desconsiderando efeitos da fronteira, é imediato que neste caso a solução preserva o ponto singular. Esse fato deve ser verificado numericamente em função dos efeitos da borda do domínio. Em seguida, utilizaremos o campo inicial mais geral dado pela expressão (8.2), o qual possui vorticidade não nula, com o valor do Laplaciano sendo diferente de zero na singularidade.

Quanto as condições de fronteira, estaremos interessados em verificar a preservação da singularidade tanto para condições de fronteira homogênea quanto para não homogênea, bem como verificar o surgimento de novos pontos singulares. Assim, vamos utilizar as seguintes condições na fronteira ∂ℛ do domínio $\mathcal{R}=[-1,1]\times [-1,1]$.

• (1) Condição de Fronteira Homogênea:

${\mathbf{v}}_{\partial \mathcal{R}}=(0,0),$ (8.3)

• (2) Condição de Fronteira Não-Homogênea:

$\left\{\begin{array}{l}\mathbf{v}(-1,y)=(0,0),\mathbf{v}(1,y)=(0,0),\\ \mathbf{v}(x,-1)=(0,-(x^2-1)/2),\mathbf{v}(x,1)=(0,-(x^2-1)/2).\end{array}\right.$ (8.4)

Quanto aos métodos numéricos, utilizaremos o método de Euler explícito usado em ¹⁵, um método de DF para a equação estacionária (Seção 5) e o método baseado em wavelets de Haar (Seção 4), cujas respectivas soluções serão denotadas por v _E, v _I e v _W. Os valores dos parâmetros do GVF são $g=0.1$ e $h=0.8$, enquanto as malhas usadas nas discretizações são regulares, com resolução $(N_x,N_y)=(16,16)$, e células com tamanho $\Delta x=\Delta y=2/15$. Foram realizados experimentos adicionais com $(N_x,N_y)=(32,32)$ e $(N_x,N_y)=(64,64)$ para verificar a influência da discretização escolhida. Como os resultados não mudaram qualitativamente com as novas malhas, optamos por apresentar nesta seção apenas os resultados para a malha com resolução $16\times 16$. Assim, os campos vetoriais serão representados numericamente por matrizes no espaço ℝ^16×16×2 . Dados dois campos vetoriais (discretos) $\mathit{F}=\left[F_1,F_2\right],\mathit{G}=\left[G_1,G_2\right]$, as distâncias entre eles serão calculadas pela norma do máximo através da expressão:

${||\mathit{F}-\mathit{G}||}_{\infty }=\underset{\mathit{x}\in \Omega }{\max }\left\{\left|F_1\left(\mathit{x}\right)-G_1\left(\mathit{x}\right)\right|,\left|F_2\left(\mathit{x}\right)-G_2\left(\mathit{x}\right)\right|\right\}.$ (8.5)

Nos testes desta seção, os valores de s ₁ e s ₂ nas expressões (8.1)-(8.2) são escolhidos de modo que a singularidade do campo vetorial inicial coincida com um dos pontos de colocação da solução por wavelets. Como esse método necessita de 2^J+1 pontos em cada uma das direções x e y, como descrito na Seção 4, posicionamos o ponto singular em $(s_1,s_2)=(1/15,1/15)$, que é o ponto central da malha utilizada.

É importante observar que os campos vetoriais das expressões (8.1)-(8.2) não satisfazem as condições de fronteira (8.3)-(8.4). Assim, temos uma descontinuidade do campo vetorial na vizinhança da fronteira do domínio. Embora este fato tenha efeitos teóricos e numéricos que não foram considerados neste trabalho, este se justifica em função da aplicação de interesse em ⁸, onde esta restrição não é satisfeita.

8.1 Campo vetorial 8.1 e condição de fronteira homogênea 8.3

O problema de valor de contorno será formado pela equação estacionária do GVF (2.4), e pela condição de fronteira (8.3). A Figura 1.(a) mostra o campo vetorial v ₀₁ e sua singularidade no domínio considerado. Primeiramente, é direto mostrar que:

1. $∆{\mathbf{v}}_{01}(x_1,x_2)=(0,0)$ em ℝ²,

2. O campo v ₀₁ é solução do GVF em ℝ², ou seja, $\mathbf{v}={\mathbf{v}}_{01}$ ,

3. $\nabla \times {\mathbf{v}}_{01}({\mathrm{x}}_1,{\mathrm{x}}_2)=(0,0)$ em ℝ²,

4. Como estamos usando condição de fronteira homogênea para o GVF, devemos considerar a mesma condição de fronteira para a equação de Helmholtz (expressão (7.7)). Desta forma o problema de contorno correspondente é homogêneo, implicando em solução trivial (nula) para o mesmo.

Os fatos (1)-(2) implicam que a solução v do GVF em ℝ² satisfaz $∆\mathbf{v}({\mathrm{s}}_1,{\mathrm{s}}_2)=(0,0)$. Além disso, a observação (4) implica que não teremos geração de vorticidade. Assim, não esperamos mudanças significativas do comportamento qualitativo entre a solução v do problema de valor de contorno (2.4)-(8.3) e o campo v ₀₁. Porém, em (2) não estamos considerando domínio compacto e, além disso, o campo v ₀₁ não satisfaz a condição de contorno homogênea. Assim, é importante verificar numericamente o efeito da condição de fronteira no resultado do GVF. Os campos mostrados nas Figuras 2.(b)-(d) foram obtidos pela solução do problema de valor de contorno (2.4)-(8.3) com os três métodos numéricos considerados, comprovando a hipótese de pouca mudança qualitativa.

Em termos quantitativos, a expressão (3.2) é um filtro passa-baixa, tendo a característica de diminuir a intensidade do campo. Porém, os resultados da Seção 3 também não consideram condições de fronteira. Assim, o problema de valor de contorno (2.4)-(8.3) foi resolvido usando os três métodos numéricos considerados. As estatísticas mostradas na Tabela 1 permitem a comparação numérica entre os campos mostrados na primeira coluna. A quarta e quinta colunas da Tabela 1 confirmam a redução na intensidade do campo de velocidades nas soluções numéricas.

Tabela 1
Tempo de CPU (T(s)), velocidade na singularidade $(\Vert \mathbf{v}({\mathrm{s}}_1,{\mathrm{s}}_2\left){\Vert }_2\right)$, média (µ _v) e desvio padrão (σ _v ) da intensidade da velocidade, módulo do Laplaciano no ponto singular $\left(\right|∆\mathbf{v}({\mathrm{s}}_1,{\mathrm{s}}_2)\left|\right)$, média (µ _∆) e desvio padrão (σ _∆) do módulo do Laplaciano para o campo inicial v ₀ e as soluções numéricas v _E, v _I e v _W.

×	T(s)	$\Vert \mathbf{v}({\mathrm{s}}_1,{\mathrm{s}}_2){\Vert }_2$	µ v	σ v	$|∆\mathbf{v}({\mathrm{s}}_1,{\mathrm{s}}_2\left)\right|$	µ Δ	σ Δ
v 0	×	0.0000	0.5488	0.3847	0.0000	26.4057	33.1545
v E	0.9740	0.0161	0.1867	0.1200	0.1263	5.1313	3.3770
v I	0.0020	0.0159	0.1867	0.1199	0.1271	5.1304	3.3767
v W	15.2553	0.0147	0.1880	0.1245	0.1271	4.6855	2.9503

Primeiramente, fizemos uma busca em torno do ponto (s ₁ , s ₂) para as soluções v _E, v _I, v _W, e identificamos que o mínimo do módulo da velocidade, representado genericamente por $\Vert \mathbf{v}({\mathrm{s}}_1,{\mathrm{s}}_2){\Vert }_2$ na terceira coluna da Tabela 1, ocorre no próprio (s ₁ , s ₂), como esperado teoricamente. Comparando ainda o módulo da velocidade no ponto (s ₁ , s ₂) com o intervalo estatístico $I_s=\left[{\mathrm{\mu }}_{\mathrm{v}}-{\sigma }_{\mathrm{v}},{\mathrm{\mu }}_{\mathrm{v}}+{\sigma }_{\mathrm{v}}\right]$ de cada método, que satisfaz $I_s\approx \left[0.07,0.31\right]$, pode-se dizer que o valor $\Vert \mathbf{v}(s_1,s_2){\Vert }_2$ é pequeno em relação ao intervalo I _s . Como comentado na Seção 6, e explicitado na expressão (6.3) é preciso relaxar o conceito de singularidade aqui, em função de aproximações inerentes à solução numérica e à aritmética de ponto-flutuante. Consequentemente, os métodos numéricos preservaram a singularidade, como visualmente observado nas Figuras 2.(b)-(d).

Neste contexto, o valor do campo de velocidade gerado neste ponto $(\Vert \mathbf{v}(s_1,s_2\left){\Vert }_2\right)$ é mais próximo de zero na solução via wavelets, embora essa diferença aconteça apenas na terceira casa decimal. Em termos relativos, considerando agora a média e o desvio padrão da intensidade do campo de velocidades, nota-se que estas estatísticas são idênticas até a segunda casa decimal, confirmando assim a pequena vantagem do método de wavelets em relação aos demais. Em relação ao valor do Laplaciano no ponto singular $(∆{\mathbf{v}}_0(s_1,s_2\left)\right)$, teoricamente deveria ser nulo pela expressão (2.7). De fato, os valores mostrados na sexta coluna da Tabela 1 podem ser considerados pequenos, levando em conta que os intervalos estatísticos dos métodos numéricos satisfazem $[{\mathbf{\mu }}_{∆}-{\mathbf{\sigma }}_{∆},{\mathbf{\mu }}_{∆}+{\mathbf{\sigma }}_{∆}]\subset [1.75,6.64]$. Porém, observa-se uma vantagem estatisticamente pouco significativa para o método de Euler explícito: 0.1263 contra 0.1271 para os demais métodos, sendo que os limites dos intervalos estatísticos são da mesma ordem de grandeza.

O tempo de CPU de cada método numérico, mostrando na segunda coluna da Tabela 1, indica que em nossas implementações, o método de wavelets é bem mais caro que os demais. Isso se justifica pela complexidade necessária para montar e resolver o sistema linear (4.16). Visualmente, as diferenças entre as linhas de campo nas Figuras 2.(b)-(d) não são perceptíveis. Esse fato se comprova numericamente, observando que $\Vert {\mathbf{v}}_{\mathbf{E}}-{\mathbf{v}}_{\mathbf{I}}{\Vert }_{\infty }=0.0002,\Vert {\mathbf{v}}_{\mathbf{E}}-{\mathbf{v}}_{\mathbf{W}}{\Vert }_{\infty }=0.02135$ e $\Vert {\mathbf{v}}_{\mathbf{I}}-{\mathbf{v}}_{\mathbf{W}}{\Vert }_{\infty }=0.02129$, onde $\Vert \mathbf{a}\mathbf{-}\mathbf{b}{\Vert }_{\infty }$ é calculado pela equação (8.5). Esses valores podem ser considerados pequenos levando em conta que o intervalo estatístico I _s da cada método para a função $\Vert \mathbf{v}(s_1,s_2){\Vert }_2$ satisfaz $I_s\approx \left[0.07,0.31\right]$, como comentado acima.

8.2 Campo vetorial (8.1) e condição de fronteira heterogênea (8.4)

Agora, vamos utilizar o mesmo campo vetorial do exemplo anterior (expressão (8.1)), mas usando a condição (8.4) na fronteira do domínio $\mathcal{R}=[-1,1]\times [-1,1]$, para resolver a equação estacionária do GVF (2.4) Os fatos (1)-(3) da Seção 8.1 continuam válidos, mas agora a condição de fronteira para a equação de Helmholtz (expressão (7.7)) é dada por:

$\left\{\begin{array}{r}\mathbf{\xi }(-1,y)=(0,0),\\ \mathbf{\xi }(1,y)=(0,0),\\ \mathbf{\xi }(x,-1)=(0,0,-x),\\ \mathbf{\xi }(x,1)=(0,0,-x),\end{array}\right.$ (8.6)

onde estas expressões são obtidas aplicando a definição (2.8) para vorticidade nas equações (8.4), que definem a condição de fronteira para a velocidade. Portanto, o problema de contorno para a vorticidade não terá solução trivial como ocorreu no exemplo da Seção 8.1. Pela discussão apresentada na Seção 7 poderemos ter geração de novos pontos singulares na solução do GVF, o que se confirma na visualização das soluções apresentadas nas Figuras 3.(b)-(d). Nestas figuras, além do ponto singular original, observamos a convergência das linhas de corrente para um novo ponto, localizado na posição (−0.733, 0.067).

Figura 3:
(a) Campo vetorial (8.1). (b) Linhas de corrente da solução do GVF via Euler explícito. (c) Solução do GVF usando DF para a equação estacionária. (d) Visualização da solução do GVF usando método de wavelets.

A Tabela 2 mostra as estatísticas dos módulos do campo de velocidade (Laplaciano) e seus valores no ponto singular original, em fonte usual, e na nova singularidade, em negrito. Apesar da criação de uma nova singularidade em (−0.733, 0.067), o ponto singular original foi preservado pelos três métodos numéricos, uma vez que o intervalo estatístico I _s para a função $\Vert \mathbf{v}(s_1,s_2){\Vert }_2$ nos três métodos satisfaz $I_s\approx \left[0.1,0.37\right]$. Porém, a solução via wavelets fornece um resultado mais próximo de zero que as demais, o que é verificado comparando o valor do módulo da velocidade nas singularidades (terceira coluna) nos três métodos numéricos e levando em conta que os intervalos estatísticos são muito próximos.

Tabela 2
Tempo de CPU (T(s)), velocidade na singularidade $(\Vert \mathbf{v}(s_1,s_2\left){\Vert }_2\right)$, média (µ _v) e desvio padrão (σ _v ) da intensidade da velocidade, módulo do Laplaciano no ponto singular $\left(\right|∆\mathbf{v}(s_1,s_2)\left|\right)$, média (µ _∆) e desvio padrão (σ _∆) do módulo do Laplaciano para o campo inicial v ₀ e as soluções numéricas v _E, v _I e v _W.

×	T(s)	$\Vert \mathbf{v}(s_1,s_2){\Vert }_2$	µ G	σ G	$|∆\mathbf{v}(s_1,s_2\left)\right|$	µ Δ	σ Δ
v0	×	0.0000	0.5877	0.3455	0.0000	26.6846	34.8954
		0.8000			0.0000
vE	1.1492	0.0327	0.2333	0.1339	0.2640	5.3942	4.2628
		0.0770			5.7872
vI	0.0020	0.0329	0.2333	0.1339	0.2631	5.3939	4.2624
		0.0769			5.7866
vW	15.5518	0.0120	0.2339	0.1384	0.6470	4.7602	3.5738
		0.0313			6.4243

Quanto ao valor do módulo do Laplaciano $|∆\mathbf{v}(s_1,s_2\left)\right|$, mostrado na sexta coluna da Tabela 2, observamos que o método de wavelets foi o menos preciso em relação à conservação desta quantidade na singularidade original $(s_1,s_2)=(1/15,1/15)$. Isto se confirma observando o valor de $|∆\mathbf{v}(s_1,s_2\left)\right|$ neste ponto e os intervalos estatísticos dos métodos, pois ${\mu }_{∆}-{\sigma }_{∆}\simeq 1.2$ em todos os casos. Assim como no caso anterior, o tempo de CPU do método de wavelets foi bem maior que para os demais métodos (segunda coluna).

Com relação ao novo ponto singular, na posição $(s_1,s_2)=\left(-0.733,0.067\right)$, o método de wavelets forneceu o melhor resultado para o valor da velocidade neste ponto. Em relação ao Laplaciano neste ponto, podemos ver pelas estatísticas fornecidas na tabela que seu valor está acima da média em todos os casos, mostrando que ele não se anulou no novo ponto. Observando a equação (2.7) esse fato faz sentido, uma vez que o ponto (−0.733, 0.067) não é um ponto singular do campo v ₀.

8.3 Campo vetorial (8.2) e condição de fronteira heterogênea (8.4)

Neste último experimento usamos o campo vetorial da expressão (8.2), utilizando a condição de fronteira não homogênea descrita em 8.4. O campo em questão não atende as condições (1), (2) e (3) presentes na Seção 8.1. Além disso, neste caso o campo inicial é rotacional (vorticidade não-nula), o que pode ser verificado pela aplicação direta da definição (2.8) no campo em questão. Assim, do ponto de vista da vorticidade, estamos na situação descrita pela equação de Helmholtz não-homogênea, dada pela expressão (7.6), com condição de fronteira não-homogênea também.

Como o campo inicial não é solução do GVF em ℝ², a análise local não pode ser usada diretamente como foi feito na Seção 8.1. Nossa discussão será conduzida pela visualização do resultado do GVF e pelas estatísticas da Tabela 3. A Figura 4.(a) mostra as linhas de corrente do campo v ₀₂ em uma vizinhança da singularidade e as Figuras 4.(b)-(d) mostram as soluções do GVF, obtidas pelos métodos numéricos. É importante observar que a singularidade é não-hiperbólica, uma vez que a representação linear do campo ¹³ no ponto singular (s ₁ , s ₂) possui dois autovalores nulos, o que justifica o comportamento complexo das linhas de campo na vizinhança do ponto singular na Figura 4.(a). Porém, nas Figuras 4.(b)-(d) observamos um padrão bem mais simples, visualmente indicando fluxo tubular na região onde se localizava a singularidade inicial. As estatísticas mostradas na Tabela 3 concordam com essa análise. Por exemplo, no caso da solução via DF para a equação estacionária (Figura 4.(c)), observamos que a intensidade média do campo de velocidades e o correspondente desvio padrão fornecem um intervalo estatístico para o módulo da velocidade na solução do GVF dado por (0.1035, 0.3943) enquanto a velocidade na posição do ponto singular tem módulo $\Vert \mathbf{v}(s_1,s_2){\Vert }_2=0.1138$. Considerando o intervalo estatístico, essa velocidade não pode ser considerada próxima de zero. Uma análise do campo v na vizinhança de (s ₁ , s ₂) mostra um comportamento semelhante, indicando a inexistência de singularidades nesta região.

Tabela 3
Tempo de CPU (T(s)), velocidade na singularidade $(\Vert \mathbf{v}(s_1,s_2\left){\Vert }_2\right)$, média (µ _v) e desvio padrão (σ _v ) da intensidade da velocidade, módulo do Laplaciano no ponto singular $\left(\right|∆\mathbf{v}(s_1,s_2)\left|\right)$, média (µ _∆) e desvio padrão (σ _∆) do módulo do Laplaciano para o campo inicial v ₀ e as soluções numéricas v _E, v _I e v _W.

×	T(s)	$\Vert \mathbf{v}(s_1,s_2){\Vert }_2$	µ G	σ G	$|∆\mathbf{v}(s_1,s_2\left)\right|$	µ Δ	σ Δ
v 0	×	0.0000	0.4520	0.3503	|2.8284	22.395	29.292
v E	1.0783	0.1140	0.2489	0.1454	1.4414	3.5898	3.1675
v I	0.0020	0.1138	0.2489	0.1454	1.4415	3.5894	3.1671
v W	13.502	0.1476	0.2724	0.1559	1.5805	3.4479	2.9585

Figura 4:
(a) Campo vetorial (8.2). (b) Linhas de corrente da solução do GVF via Euler explícito. (c) Solução do GVF usando DF para a equação estacionária. (d) Visualização da solução do GVF usando método de wavelets.

Com relação à intensidade do Laplaciano da solução do GVF, via DF, para a equação estacionária, vemos que o intervalo estatístico neste caso é (0.4223, 6.7565) enquanto $\Vert ∆\left(\mathbf{v}\right(s_1,s_2\left)\right){\Vert }_2=1.4415$ não podendo ser considerado nulo também. Pela expressão (2.7) da análise local, vem que o ponto singular não é preservado na solução do GVF, estando coerente com o exposto acima. Análises idênticas podem ser realizadas para a solução via Euler explícito, uma vez que as estatísticas obtidas são muito próximas até a terceira casa decimal.

No caso da solução por wavelets, os números mudam um pouco em relação aos dois primeiros métodos, mas obtêm-se as mesmas conclusões, mantendo também o padrão das linhas de corrente, como pode ser verificado comparando as Figuras 4.(c),(d). Com relação à vorticidade do campo v, sua intensidade média é 0.354 e desvio padrão 0.290, tanto para o DF para a equação estacionária quanto para Euler explícito. No caso do método de wavelets, a média é 0.394 e o desvio padrão 0.322. Assim, não se trata de um campo irrotacional, o que não é surpresa visto que o campo (8.2) é rotacional e o problema de valor de contorno para a vorticidade é não-homogêneo.

9 CONCLUSÕES FINAIS

Nesse trabalho, nosso foco foi a preservação de singularidades do campo inicial na equação de difusão-reação do GVF. Localmente, a preservação da singularidade implica que o Laplaciano da solução seja nulo, e vice-versa. A análise no domínio da frequência mostra que se o campo inicial for derivado de um potencial f cuja transformada de Fourier F seja simétrica em relação à origem do espaço de frequências a singularidade inicial é preservada. Quando supomos domínio compacto, usando o método de solução baseado em wavelets, foi possível explicitar a influência das condições de fronteira no problema de interesse. Numericamente, impor a preservação da singularidade implica em um sistema super-determinado, sem garantias de solução. Do ponto de vista numérico, as soluções do GVF via DF para a equação estacionária e via wavelets foram comparadas com a tradicional, baseada no método de Euler explícito. Os resultados mostram concordância entre as soluções, bem como em relação à preservação (ou não) da singularidade.

No que diz respeito a aplicação do GVF na metodologia apresentada em ⁸, os experimentos realizados mostram que a condição de fronteira homogênea é mais adequada. O objetivo em ⁸ é utilizar a solução do GVF para obter um campo vetorial suave, com topologia semelhante à condição inicial v ₀ gerada com atuação do usuário. Por fim, esta solução é utilizada como condição inicial em animação de fluidos. Os experimentos apresentados em ⁸ mostraram também preservação de singularidades isoladas na solução do GVF com condição de fronteira homogênea, concordando com os resultados da Seção 8.1.

Como trabalhos futuros, pretendemos explorar a restrição (4.18), de uma forma mais fraca que não torne o problema mal-posto, e estudar a preservação do tipo de singularidade na solução do GVF. Serão também analisados aspectos numéricos e teóricos referentes ao relaxamento da condição de singularidade (Seção 6).

Agradecimentos

O presente trabalho foi realizado com apoio da Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Código de Financiamento 001, e da Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ).

REFERÊNCIAS

• ¹
  I. Aziz, A. Al-fhaid & A. Shah. A numerical assessment of parabolic partial differential equations using Haar and Legendre wavelets. Applied Mathematical Modelling, 37 (2013), 9455-9481.
• ²
  D. Boukerroui. Efficient numerical schemes for gradient vector flow. In “ICIP” (2009).
• ³
  C. Chui. “An Introduction to Wavelets”. New York Academic Press (1992).
• ⁴
  H. Cui & Y. Xia. Automatic Coronary Centerline Extraction Using Gradient Vector Flow Field and Fast Marching Method From CT Images. IEEE Access, 6 (2018), 41816-41826.
• ⁵
  V. Girault & P.A. Raviart. “Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms”. Springer-Verlag (1986).
• ⁶
  G. Hariharan. “Wavelet Solutions for Reaction-Diffusion Problems in Science and Engineering”.Springer (2019).
• ⁷
  J.D. Hunter. Matplotlib: A 2D graphics environment. Computing in Science & Engineering, 9(3) (2007), 90-95. doi:10.1109/MCSE.2007.55.
  » https://doi.org/10.1109/MCSE.2007.55
• ⁸
  S. Judice & G. Giraldi. Fluid Animation Using Sketching, Diffusion-Reaction and Lattice Boltzmann Models. IEEE Latin America Transactions, 18(03) (2020), 514-521.
• ⁹
  S.F. Judice & G. Giraldi. Sketching Fluid Flows - Combining Sketch-based Techniques and Gradient Vector Flow for LBM Initialization. In “GRAPP 2012 - International Conference on Computer Graphics Theory and Applications” (2012).
• ¹⁰
  S.F. Judice , J.G. Mayworm, P. Azevedo & G. Giraldi . Perspectives for Sketching Fluids Using Sketch- Based Techniques and Gradient Vector Flow for 3D LBM Initialization. In “Computer Vision, Imaging and Computer Graphics. Theory and Application”. Springer Berlin Heidelberg (2013), pp. 127-141.
• ¹¹
  U. Lepik. Numerical solution of evolution equations by the Haar wavelet method. Applied Mathematics and Computation, 185(1) (2007), 695-704. doi:10.1016/j.amc.2006.07.077.
  » https://doi.org/10.1016/j.amc.2006.07.077
• ¹²
  U. Lepik . Solving PDEs with the aid of two-dimensional Haar wavelets. Computers & Mathematics with Applications, 61(7) (2011), 1873-1879.
• ¹³
  J. Sotomayor. “Lições de Equações Diferenciais Ordinárias”. Projeto Euclides, Gráfica Editora Hamburgo Ltda., São Paulo (1979).
• ¹⁴
  G.S. Weiss & G. Zhang. Existence of a Degenerate Singularity in the High Activation Energy Limit of a Reaction-Diffusion Equation. Communications in Partial Differential Equations, 35 (2009), 185-199.
• ¹⁵
  C. Xu & J.L. Prince. Snakes, shapes, and gradient vector flow. IEEE Transactions on Image Processing, 7(3) (1998), 359-369. doi:10.1109/83.661186.
  » https://doi.org/10.1109/83.661186
• ¹⁶
  C. Xu & J.L. Prince . Global Optimality of Gradient Vector Flow. In “Proc. of 34th Annual Conference on Information Sciences and Systems”. Princeton University (2000).

Datas de Publicação

• Publicação nesta coleção
  02 Jul 2021
• Data do Fascículo
  May-Aug 2021

Histórico

• Recebido
  25 Jun 2020
• Aceito
  11 Jan 2021