Potencial gravitacional de discos achatados usando o teorema de convolução

Gravitational potential of flattened disks using the convolution theorem

Resumo

O potencial gravitacional é muito utilizado nas áreas como astrofísica e física de plasmas, por exemplo, no estudo de formação e evolução de sistemas estelares e de sistemas planetários, em que a dinâmica de suas partículas constituintes é governada principalmente pelo campo gravitacional (aceleração) gerado pelo próprio sistema, que é determinado a partir do seu potencial gravitacional. Alguns desses sistemas têm estruturas relativamente achatadas, em que a espessura da estrutura é bastante pequena em comparação com o diâmetro, por exemplo o disco galáctico de uma galáxia espiral ou anéis planetários, nessas situações é válida a aproximação do sistema como uma estrutura totalmente achatada, é dizer como uma estrutura bidimensional. Neste trabalho, apresentamos uma forma numérica de determinar o potencial gravitacional gerado por um sistema extremadamente achatado em forma de disco fino, em um ponto sobre o plano do disco. O método consiste em dividir o domínio espacial do sistema em pequenas seções retangulares, em que em cada uma das seções a densidade de massa é considerada constante. Determinamos analiticamente a expressão para o potencial gravitacional gerada por uma seção retangular de densidade constante. Logo, devido ao princípio de superposição, o potencial gerado pelo sistema no centro de uma seção em particular será a soma de todos os potenciais geradas por cada uma das seções que formam o sistema. Para realizar essa somatória dos potenciais de uma forma mais eficiente usaremos o teorema de convolução da transformada discreta de Fourier. Finalmente analisaremos a precisão do método comparando nossos resultados com os resultados analíticos encontrados para o disco de Kalnajs e para o disco uniforme.

Palavras-chave:
Potencial gravitacional; astrofísica; métodos numéricos; sistemas contínuos; transformada discreta de Fourier

Abstract

Gravitational potential is widely used in areas such as astrophysics and plasma physics, for example, in the study of the formation and evolution of stellar systems and planetary systems, in which the dynamics of their constituent particles is governed mainly by the gravitational field (acceleration) generated by the system itself, which is determined from its gravitational potential. Some of these systems have relatively flattened structures, where the thickness of the structure is quite small compared to the diameter, for example the galactic disc of a spiral galaxy or planetary rings, in these situations the approach of the system is valid as a totally flattened structure, is to say as a two-dimensional structure. In this work, we present a numerical way to determine the gravitational potential generated by an extremely flattened, thin disc-shaped system, at a point on the disc plan. The method consists of dividing the spatial domain of the system into small rectangular sections, where in each section the mass density is considered constant. We analytically determine the expression for the gravitational potential generated by a rectangular section of constant density. Therefore, due to the principle of superposition, the potential generated by the system in the center of a particular section will be the sum of all the potentials generated by each of the sections that form the system. To perform this summation of potentials in a more efficient way, we will use the convolution theorem of the discrete Fourier transform. Finally, we will analyze the accuracy of the numerical method by comparing our results with the analytical results found for the Kalnajs disk and the uniform disk.

Keywords:
Gravitational potential; astrophysics; numerical methods; continuous systems; discrete Fourier transform

1. Introdução

A determinação do potencial Newtoniano ou potencial gravitacional é um problema recorrente em astronomia, astrofísica, eletrostática, magnetostática, planetologia e física de plasmas. Seja um sistema em que a distribuição de massa (carga) seja contínua, logo o potencial gravitacional (potencial elétrico) Φ(x) em alguma posição $\overrightarrow{x}$ é definido como [¹,²,³],

(1)$$\mathrm{\Phi }(x)=-G\int d^3{\overrightarrow{x}}^{\prime }\frac{\rho ({\overrightarrow{x}}^{\prime })}{|{\overrightarrow{x}}^{\prime }-\overrightarrow{x}|}$$

em que G é a constante gravitacional (para o caso do potencial elétrico a constante G deve ser substituida por −K, em que K é a constate de Coulomb), $\rho ({\overrightarrow{x}}^{\prime })$ é a distribuição de massa (carga) na posição ${\overrightarrow{x}}^{\prime }$, e a integração inclui todos os pontos do sistema. Aplicando o operador laplaciano ∇² ≡ ∇⋅∇ à equação (1), é dizer aplicando a gradiente e em seguida a divergência, e logo usando o teorema de divergência para transformar a integral de volume em uma integral de superfície [¹,²,³] é possível obter a equação de Poisson

(2)$${\nabla }^2\mathrm{\Phi }=4\pi G\rho .$$

Esta é uma equação diferencial parcial que pode ser resolvida para Φ(x) dada ρ(x) e uma condição de contorno apropriada. Para um sistema isolado, a condição de contorno é Φ → 0 quando |x| → ∞. O potencial dado pela equação (1) satisfaz automaticamente essa condição de contorno. No caso especial em que ρ = 0 (em um ponto do espaço onde exista ausência de massa ou carga) a equação de Poisson se torna a equação de Laplace,

(3)$${\nabla }^2\mathrm{\Phi }=0.$$

Existem diversas técnicas para resolver a equação (2), por exemplo, a técnica baseada nos harmônicos esféricos [¹], que funciona melhor para sistemas que não são nem muito achatados e nem muito alongados. Para sistemas que têm superfícies de equidensidade em forma elipsoidal ou esferoidal concêntricas, a equação de Poisson é resolvida usando coordenadas esferoidais, logo fazendo tender a zero um dos eixos de um elipsoide, mantendo constante os outros eixos, o elipsoide se torna um disco fino (um disco totalmente achatado) e, assim, pode-se obter os potenciais de discos finos como um caso limite de elipsoide [¹]. Outra forma de resolver a equação de Poisson para sistemas não muito alogados e nem muito achatados é usando métodos numéricos, por exemplo, o método de diferença finita [¹].

Na literatura existem muitos modelos de discos finos, em que são apresentados os pares potencial-densidade, por exemplo, os discos de Kalnajs [⁴], os discos e anéis de Kuzmin-Toomre [⁵], os discos e anéis de Morgan & Morgan [⁶], e disco de densidade uniforme [⁷]. A importância de determinar o potencial de um disco achatado é devido a que, por exemplo, em uma galáxia espiral típica a maior parte da luz emitida vem de um disco achatado [¹], [⁸] é dizer, que uma fração substancial da massa visível da galáxia está concentrada no disco e, portanto, uma determinação eficiente do potencial gravitacional de um disco achatado pode ser muito importante para estudar a influência desse componente galáctico sobre as outras componentes da galáxia. Por exemplo, na Via Láctea a maioria das estrelas estão em uma estrutura achatada e aproximadamente axissimétrica conhecida como disco galáctico. Em noites claras e escuras, a luz cumulativa da miríade de estrelas do disco é visível como uma banda luminosa que se estende pelo céu, que é a fonte do nome “Via Láctea” para a nossa galáxia. O plano médio desse disco é chamado de plano galáctico e serve como equador das coordenadas galácticas. Na direção perpendicular ao plano galáctico a densidade das estrelas cai exponencialmente com a distância ao plano [¹],

(4)$$\rho (R,\phi ,z)=\rho (R,\phi ,0){\text{e}}^{-|z|/z_d(R)},$$

em que R, φ e z são as coordenadas cilíndricas sobre o disco galáctico, com o eixo z perpendicular ao disco galáctico, e z_d(R) é uma certa medida da espessura do disco galáctico na posição R. A espessura z_d(R) do disco galáctico depende da idade das estrelas que estão sendo consideradas, sendo maoir para as populações estelares mais antigas. Na vizinhança do sol, a espessura varia de ≲ 100 pc para as estrelas jovens a ≃ 300 pc para as estrelas mais velhas com idades de ordem de 10 Gyr que constituem a maior parte da massa de disco [¹]. Por outro lado o diâmetro do disco galáctico da Via Láctea é em torno de 25 000 pc [¹, ⁸] logo a razão entre a espessura do disco e o diâmetro do disco é da ordem de 10⁻², é dizer, o disco é uma estrutura bastante achatada. Portanto é viável considerar, pelo menos aproximadamente, a esse disco galactico como um disco totalmente achatado, é dizer reescrevendo a densidade (4) como

(5)$$\rho (R,\phi ,z)\approx \mathrm{\Sigma }(R,\phi )\delta (z),$$

em que δ(z) é a função delta de Dirac [⁹] e Σ(R, φ) deve ser a densidade superficial de massa do disco galactico. Com essa aproximação a espessura do disco é desconsiderada e considerando-se simplesmente como uma estrutura bidimensional.

Neste artigo, determinamos numericamente o potencial gravitacional de um sistema em forma de disco totalmente achatado, a partir da sua forma integral (1), primeiramente transformando a integral em uma somatória e logo usando o teorema de convolução da transformada discreta de Fourier para calcular a somatória de uma forma mais rápida. O método é primeiramente testado em um modelo de brinquedo, uma estrutura extremadamente alongada (sistema unidimensional), em que analisamos a eficiência desse método. Finalmente determinamos numericamente o potencial do disco de Kalnajs [⁴] e do disco uniforme [⁷] e analisamos o erro percentual entre o potencial numérico e o potencial analítico para esses discos.

2. Transformada Discreta de Fourier

A transformada discreta de Fourier (TDF) [⁹] transforma uma sequência de N números complexos, {x_n} ≡ {x₁, x₂, ⋯, x_N}, em outra sequência de números complexos, $\{{\widehat{x}}_k\}\equiv \{{\widehat{x}}_1,{\widehat{x}}_2,\mathrm{\cdots },{\widehat{x}}_N\}$, definida como

(6)$${\widehat{x}}_k=\frac{1}{\sqrt{N}}\sum _{n=1}^N{x}_n{\gamma }_{nk},k=1,2,\mathrm{\cdots }N$$

em que

(7)$${\gamma }_{nk}={\text{e}}^{-\frac{2\pi ink}{N}},$$

com as seguintes propriedades

(8)$${\gamma }_{0k}={\gamma }_{Nk}=1$$$${\gamma }_{nk}={\gamma }_{kn}$$$${\gamma }_{n+N,k}={\gamma }_{nk}.$$

A TDF inversa é definida como

(9)$$x_n=\frac{1}{\sqrt{N}}\sum _{k=1}^N{\widehat{x}}_k{\gamma }_{kn}^{*},n=1,2,\mathrm{\cdots }N$$

em que ${\gamma }_{kn}^{*}$ representa o complexo conjugado de γ_kn [¹⁰]. A TDF compartilha muitas das propriedades da transformada contínua de Fourier [⁹], mas para isso é necessário definir as quantidades x_n para n fora do intervalo [1, N] pela regra

(10)$$x_n=x_{n+mN},\text{para todo o número inteiro}\mathrm{}m$$

ou seja, assumir que x_n é periódico com o período N, isso implica a partir da equação (6) que o ${\widehat{x}}_k$ também é periódico com período N. Com as propriedades anteriores pode-se provar o teorema de convolução para a TDF, que consiste em

(11)$$\sum _{n=1}^N{x}_{k-n}{y}_n=\sqrt{N}{z}_k,k=1,2,\mathrm{\cdots }N$$

onde a sequência {z_k} é a TDF inversa da sequência $\{{\widehat{x}}_n\cdot {\widehat{y}}_n\}\equiv \{{\widehat{x}}_1\cdot {\widehat{y}}_1,{\widehat{x}}_2\cdot {\widehat{y}}_2,\mathrm{\cdots },{\widehat{x}}_N\cdot {\widehat{y}}_N\}$, é dizer,

(12)$$z_k=\frac{1}{\sqrt{N}}\sum _{n=1}^N{\widehat{x}}_n\cdot {\widehat{y}}_n{\gamma }_{nk}^{*}.\mathit{}k=1,2,\mathrm{\cdots }N$$

Para que o teorema de convolução seja válido é necessário que as sequências envolvidas, x_k, y_n e z_k sejam períodicas com o mesmo período, é dizer, devem satisfazer a equação (10). A utilidade do teorema de convolução está em que é possível calcular a cada uma das somatórias que aparece na equação (11) de uma forma alternativa, primeiro determinando a TDF da sequência x_n e y_n, é dizer determinando ${\widehat{x}}_k$ e ${\widehat{y}}_k$, respectivamente, em seguida multiplicando um a um os elementos dessas sequências transformadas para obter a sequência $\{{\widehat{x}}_k\cdot {\widehat{y}}_k\}$, e finalmente determinar a TDF inversa dessa última sequência.

A extensão para duas dimensões da TDF pode ser feita da seguinte maneira. Seja a sequência de N⋅M números complexos, {x_nm} ≡ {x₁₁, x₁₂, ⋯, x_1N, x₂₁, x₂₂, ⋯, x_NM}, a TDF dessa sequência será uma outra sequência de números complexos, $\{{\widehat{x}}_{kl}\}\equiv \{{\widehat{x}}_{11},{\widehat{x}}_{12},\mathrm{\cdots }$, ${\widehat{x}}_{1N},{\widehat{x}}_{21},{\widehat{x}}_{22},\mathrm{\cdots },{\widehat{x}}_{NM}\}$, em que

(13)$${\widehat{x}}_{kl}={\displaystyle \frac{1}{\sqrt{N\cdot M}}}{\displaystyle \sum _{n=1}^N}{\displaystyle \sum _{m=1}^M}{x}_{nm}{\gamma }_{nk}{\gamma }_{ml},\mathrm{}k=1,2,\mathrm{\cdots }N\mathrm{}l=1,2,\mathrm{\cdots }M$$

E a TDF inversa será

(14)$$x_{nm}={\displaystyle \frac{1}{\sqrt{N\cdot M}}}{\displaystyle \sum _{k=1}^N}{\displaystyle \sum _{l=1}^M}{\widehat{x}}_{kl}{\gamma }_{kn}^{*}{\gamma }_{lm}^{*},\mathrm{}n=1,2,\mathrm{\cdots }N\mathrm{}m=1,2,\mathrm{\cdots }M$$

Logo o teorema de convolução será

(15)$${\displaystyle \sum _{n=1}^N}{\displaystyle \sum _{m=1}^M}{x}_{k-n,l-m}{y}_{nm}=\sqrt{N\cdot M}{z}_{kl},\mathrm{}k=1,2,\mathrm{\cdots }N\mathrm{}l=1,2,\mathrm{\cdots }M$$

onde a sequência {z_kl} é a TDF inversa da sequência $\{{\widehat{x}}_{nm}\cdot {\widehat{y}}_{nm}\}$, é dizer,

(16)$$z_{kl}=\frac{1}{\sqrt{N\cdot M}}\sum _{n=1}^N\sum _{m=1}^M{\widehat{x}}_{nm}\cdot {\widehat{y}}_{nm}{\gamma }_{nk}^{*}{\gamma }_{ml}^{*}.$$

Para que essa propriedade seja válida é preciso que todas as sequências envolvidas sejam períodicas, é dizer, por exemplo

(17)$$x_{nm}=x_{n+pN,m+qM},\text{para todo número inteiro}\mathrm{}p\mathrm{}\text{e}\mathrm{}q$$

de período N em relação ao primeiro índice e de período M em relação ao segundo subíndice.

3. Determinação Numérica do Potencial

3.1 Sistema unidimensional

Consideremos um sistema extremadamente alongado, algo raro de se encontrar na natureza, mas devido a sua simplicidade não deixa de ser importante como um modelo de brinquedo. Seja assim um sistema finito retilineamente alongado, em que sua densidade de massa ou de carga se anula fora do seu domínio espacial.

Dado que alguns métodos numéricos discretiza o domínio espacial em seções e considera a densidade de massa aproximadamente constante sobre cada uma dessas seções [¹¹,¹²,¹³], também aqui vamos dividir o domínio espacial do sistema em N seções, todas do mesmo comprimento Δx. Em seguida vamos definir a densidade linear de massa λ como uma função constante por seções, é dizer

(18)$$\lambda (x)={\lambda }_i,\text{se}\mathit{}x\in ]x{}_{i-\frac{1}{2}},x{}_{i+\frac{1}{2}}[$$

em que x é a coordenada cartesiana de um ponto sobre o sistema. $x_{i-\frac{1}{2}}=x_i-\mathrm{\Delta }x/2$ e $x_{i+\frac{1}{2}}=x_i+\mathrm{\Delta }x/2$ são respectivamente as fronteiras da esquerda e da direita da i-ésima seção, cujo centro é x_i. Cada λ_i adota um valor constante, assim λ adota um valor constante sobre cada seção.

Em um ponto x sobre um sistema retilineamente alongado se estabelece o seguinte potencial gravitacional [¹³],

(19)$$\mathrm{\Phi }(x)=-G\int \frac{\lambda (x^{\prime })}{|x-x^{\prime }|+\epsilon }\mathit{d}{x}^{\prime },$$

em que a integração é calculada sobre todo o sistema, e ε é o parâmetro que suaviza o potencial gravitacional, e evita divergência do potencial quando é calculado sobre o próprio sistema. Considerando G = 1, por simplicidade, e levando em consideração a definição (18) para λ, a equação (19) torna-se

(20)$$\mathrm{\Phi }(x)=-\sum _{j=1}^N{\lambda }_j{\int }_{x_{j-1/2}}^{x_{j+1/2}}\frac{d{x}^{\prime }}{|x-x^{\prime }|+\epsilon }$$

logo, o potencial no centro da i-ésima seção será

(21)$${\mathrm{\Phi }}_i\equiv \mathrm{\Phi }(x_i)=\sum _{j=1}^N{\lambda }_j{g}_{i,j}$$

em que

(22)$$g_{i,j}=-{\int }_{x_{j-1/2}}^{x_{j+1/2}}\frac{d{x}^{\prime }}{|x_i-x^{\prime }|+\epsilon }$$

Logo de resolver a integral (22) e considerando que x_i − x_j = (i − j)Δx, podemos reescrever

(23)$$g_{i,j}=\{\begin{array}{cc}\hfill 2\ln \left[\frac{2\epsilon }{\mathrm{\Delta }x+2\epsilon }\right],\hfill & \hfill \text{se}i=j,\\ \hfill \ln \left[\frac{\left(|i-j|-1/2\right)\mathrm{\Delta }x+\epsilon }{\left(|i-j|+1/2\right)\mathrm{\Delta }x+\epsilon }\right],\hfill & \hfill \text{se}i\ne j\end{array}$$

Pode-se notar que o valor de g_i,j depende somente da diferença i - j, logo fazendo i − j = k, podemos reescrever a equação (23) para k > 0, resultando

(24)$$g_{i,j}\equiv g_k=\{\begin{array}{cc}\hfill 2\ln \left[\frac{2\epsilon }{\mathrm{\Delta }x+2\epsilon }\right],\hfill & \hfill \text{se}\mathrm{}k=0,\\ \hfill \ln \left[\frac{\left(k-1/2\right)\mathrm{\Delta }x+\epsilon }{\left(k+1/2\right)\mathrm{\Delta }x+\epsilon }\right],\hfill & \hfill \text{se}\mathrm{}k>0\end{array}$$

A partir da equação (23) também notamos que g_{i, j} = g_{j, i}, de onde deduzimos a propriedade g_k = g_−k, essa propriedade vai permitir encontrar os valores de g_k para k < 0, utilizando simplesmente a equação (24). Finalmente podemos reescrever a equação (21) como

(25)$${\mathrm{\Phi }}_i=\sum _{j=1}^N{g}_{i-j}{\lambda }_j$$

A qual tem a aparência do teorema de convolução da transformada discreta de Fourier [¹,⁹,¹⁴], mas para isso é necessário que os três conjuntos numéricos {Φ_i}, {g_i} e {λ_i} tenham o mesmo número de elementos e sejam periódicas, é dizer, precisam satisfazer uma relação do tipo g_i = g_{i + mT}, para todo número inteiro m, em que T deve ser o período. Inicialmente nenhum desses conjuntos têm essa propriedade de periodicidade, e não têm o mesmo número de elementos, pois, a partir da equação (25) pode-se notar que o subíndice em g_{i – j} varia de −N + 1 (quando i = 1 e j = N) a N - 1 (quando i = N e j = 1), assim adotando 2N - 1 valores diferentes, enquanto os subíndices em Φ_i e em λ_j adotam valores de 1 a N, cada um totalizando N valores diferentes. Em resumo {g_i} tem 2N - 1 elementos, mas {Φ_i} e {λ_j} tem somente N elementos cada um. Para que os conjuntos {Φ_i}, {g_i} e {λ_i} adotem as propriedades acima mencionadas, vamos começar por estender o domínio espacial do sistema, incluindo as seções i = −N + 1 a i = 0, e logo redefinindo λ_i como

(26)$${\lambda }_i=\{\begin{array}{ccc}\hfill 0\hfill & \hfill \text{se}\hfill & \hfill -N+1\le i\le 0\hfill \\ \hfill {\lambda }_i\hfill & \hfill \text{se}\hfill & \hfill 1\le i\le N\hfill \end{array}$$

em que λ_i se anula na faixa estendida −N + 1 ≤ i ≤ 0, pois essa faixa inclui seções que estão fora do domínio espacial original do sistema, é dizer fora do domínio de nosso interesse. Alem disso, agora, vamos redefinir o conjunto {λ_i} com a propriedade de periodicidade, é dizer, λ_i = λ_{i + mT}, para todo número inteiro m, com período T = 2N. Agora o subíndice de λ pode assumir qualquer valor inteiro, isso significa que precisamos estender ainda mais o domínio espacial do sistema para um sistema espacial infinito e periódico. Como resultado teremos regiões com λ = 0 e λ ≠ 0 fora do domínio espacial original, e as regiões com λ ≠ 0 terão influência na determinação do potencial nesse domínio original, porém devido a que essas regiões se encontram longe do domínio original, então, suas influências devem ser suficientemente pequenas e, portanto, podem se considerar desprezíveis, pois o potencial newtoniano atenua-se para grandes distâncias. Do mesmo jeito, a faixa de valores para o subíndice em Φ_i também deve ser estendida, logo podemos reescrever a equação (25) como,

(27)$${\mathrm{\Phi }}_i=\sum _{j=-N+1}^N{g}_{i-j}{\lambda }_j,i=-N+1,\mathrm{\dots },N$$

que também deve ter a propriedade de periodicidade com o mesmo período que o conjunto {λ_i}. Agora os conjuntos {λ_i} e {Φ_i} estão conformados por 2N elementos diferentes, mas a partir da equação (27) notamos que, agora, o conjunto {g_i} adota 4N - 1 elementos diferentes, com subíndice que varia de − 2N + 1 a 2N - 1, logo é preciso desconsiderar os elementos com i < − N + 1 e i > N para que os três conjuntos numéricos tenham o mesmo número de elementos, isso não deve modificar a determinação do potencial, pois g_i se atenua quando i assume valores grandes, tal como se pode verificar na equação (24). E para finalizar precisamos redefinir {g_i} como um conjunto períodico de período 2N.

Devido a que é mais confortável trabalhar com subíndices positivos, e levando em consideração a periodicidade de λ, podemos reescrever o conjunto {λ_i} como

(28)$${\lambda }_i=\{\begin{array}{ccc}\hfill {\lambda }_i\hfill & \hfill \text{se}\hfill & \hfill 1\le i\le N\hfill \\ \hfill 0\hfill & \hfill \text{se}\hfill & \hfill N+1\le i\le 2N\hfill \end{array}$$

em que λ_i = λ_{i − 2N} = 0 para N + 1 ≥ i ≥ 2N, tal como se pode verificar usando a equação (26). E para o conjunto {g_i} temos

(29)$$g_i=\{\begin{array}{ccc}{g}_i\hfill & \hfill \text{se}\hfill & \hfill 1\le i\le N\hfill \\ g_{i-2N}\hfill & \hfill \text{se}\hfill & \hfill N+1\le i\le 2N\hfill \end{array}$$

e da mesma forma para o conjunto {Φ_i}. Após o ajuste dos subíndices, a equação (27) pode ser reescrita como

(30)$${\mathrm{\Phi }}_i=\sum _{j=1}^{2N}{g}_{i-j}{\lambda }_j,i=1,2,\mathrm{\cdots },2N$$

agora com a ajuda do teorema de convolução discreta de Fourier ( [¹]), se obtem

(31)$${\widehat{\mathrm{\Phi }}}_k=\sqrt{2N}{\widehat{g}}_k{\widehat{\lambda }}_k,k=1,2,\mathrm{\cdots },2N$$

em que os conjuntos {${\widehat{\mathrm{\Phi }}}_k$}, {${\widehat{g}}_k$} e {${\widehat{\lambda }}_k$} representam as transformadas discretas de Fourier dos conjuntos {Φ_i}, {g_i} e {λ_i}, respectivamente. Logo {Φ_i}, o potencial gravitacional nas seções, será simplesmente a transformada inversa de Fourier de {${\widehat{\mathrm{\Phi }}}_k$}.

Desta forma temos duas formas numericas de determinar o potencial gravitacional, a primeira consiste em usar a soma tradicional, calculando diretamente a somatória que aparece na equação (25), é dizer, para cada valor fixo do subíndice i realizamos a soma dos termos g_{i − j}λ_j, com j variando de 1 a N, a segunda consiste em usar o método de convolução, é dizer, utilizar a equação (31) para determinar ${\widehat{\mathrm{\Phi }}}_k$, logo aplicando a transformada inversa de Fourier podemos obter o potencial gravitacional. A transformada discreta de Fourier é calculada usando o esquema de transformada rápida de Fourier [¹⁵,¹⁶].

Para testar o método numérico consideremos a seguinte densidade de massa de tipo cossenoidal,

(32)$${\lambda }_i=\{\begin{array}{cc}\hfill \frac{\pi }{2}\cos \left(\pi x_i\right),\hfill & \hfill \text{se}|x_i|\le 0.5\hfill \\ \hfill 0,\hfill & \hfill \text{em outro caso}\hfill \end{array}$$

em que x_i representa o centro da i-ésima seção e consideramos ϵ = 0.01 para o parâmetro de suavidade. Consideramos o intervalo [-1,1] como o domínio espacial para a determinação do potencial, logo esse domínio foi dividido em N = m² seções, em que m será um número inteiro. O centro da i-ésima seção é determinado como x_i = − 1 + (i − 0.5)Δx.

A Figura 1 mostra a razão entre o tempo de processamento usando a soma tradicional, t_T, e o tempo de processamento usando o método de convolução, t_C, em função de m. A simulação foi feita usando um computador doméstico, processador Intel(R) Core(TM) i5-9400 CPU @ 2.90GHz, 2904 Mhz, 6 Núcleos, 6 processadores lógicos, e a linguagem de programação utilizada foi o FORTRAN. Os resultados mostram t_T/t_C < 2 para m < 8, isso significa que o método de convolução não te vantagem sobre a soma tradicional, porém para m > 9, a vantagem do método de convolução é significante, por exemplo para m = 10, o tempo gasto pelo método de convolução é aproximadamente a quarta parte do tempo gasto pela soma tradicional. Na mesma Figura, no quadro interior, apresentamos a diferença percentual, e% entre os resultados para o potencial usando a soma tradicional, Φ_T, e o método de convolução, Φ_C, é dizer

Figura 1
Comparação dos tempos de processamento, t_T, gasto pela soma tradicional e, t_C, gasto pelo método de convolução, para diferentes valores de m.

(33)$$e\%=\frac{|{\mathrm{\Phi }}_T-{\mathrm{\Phi }}_C|}{{\mathrm{\Phi }}_T}\times 100\%$$

de onde notamos que existe uma diferença quase nula entre os resultados de ambos os métodos.

Na Figura 2 a curva contínua representa o resultado para o potencial gravitacional usando o método de convolução, e a curva tracejada representa a densidade de massa cossenoidal (32). O centro do poço potencial coincide com a região com maior concentração de massa.

Figura 2
Resultado para o potencial usando o método de convolução, considerando uma densidade cossenoidal.

Como segundo exemplo consideremos uma densidade de massa de tipo modular, é dizer

(34)$${\lambda }_i=\{\begin{array}{cc}\hfill 4|x_i|,\hfill & \hfill \text{se}|x_i|\le 0.5\hfill \\ \hfill 0,\hfill & \hfill \text{em outro caso}\hfill \end{array}$$

o resultado é apresentado na Figura 3, novamente a curva contínua representa o resultado para o potencial gravitacional usando o método de convolução, e a curva tracejada representa a densidade de massa modular. Neste caso, observamos dois poços potenciais centrados na região com maior densidade de massa. O comportamento do tempo de processamento usando a soma tradicional e o método de convolução é similar que para o caso cossenoidal (Figura 1).

Figura 3
Resultado para o potencial usando o método de convolução, considerando uma densidade modular.

3.2 Sistema bidimensional

Consideremos um sistema finito e extremadamente achatado em forma de um plano, cuja densidade superficial de massa σ seja uma função constante por seções, é dizer

(35)$$\sigma (x,y)={\sigma }_{i,j},\text{se}\mathit{}(x,y)\in {𝒞}_{i,j}$$

em que (x, y) são as coordenadas cartesianas de um ponto sobre o sistema, ${𝒞}_{i,j}=\mathit{}]x{}_i-\frac{\mathrm{\Delta }x}{2},x{}_i+\frac{\mathrm{\Delta }x}{2}[\mathit{}\times \mathit{}]y{}_j-\frac{\mathrm{\Delta }y}{2},y{}_j+\frac{\mathrm{\Delta }y}{2}[\mathit{}$ é a i-ésima seção sobre o sistema na direção do eixo x e j-ésima seção na direção do eixo y, cujo centro é (x_i, y_j) e os σ_i,j são constantes em cada seção. Consideramos que as seções têm as mesmas dimensões espaciais, é dizer, Δx é o comprimento na direção x e Δy é o comprimento na direção y.

Para sistemas totalmente achatados (planos) pode-se estabelecer o seguinte potencial gravitacional no ponto interior ao sistema (x, y),

(36)$$\mathrm{\Phi }(x,y)=-G\int \frac{\sigma (x^{\prime },y^{\prime })}{\sqrt{{(x-x^{\prime })}^2+{(y-y^{\prime })}^2}}\mathit{d}{x}^{\prime }\mathit{d}{y}^{\prime },$$

em que a integração é calculada sobre todo o domínio espacial do sistema. Note que, neste caso não é necessário inserir o parâmetro de suavidade ε. Novamente, por simplicidade consideramos G = 1. Considerando que σ é uma função constante por seções, a equação (36) torna-se

(37)$$\mathrm{\Phi }(x,y)=\mathrm{}-{\displaystyle \sum _{k=1}^N}{\displaystyle \sum _{l=1}^M}{\sigma }_{k,l}{\displaystyle {\int }_{y_l-\frac{\mathrm{\Delta }y}{2}}^{y_l+\frac{\mathrm{\Delta }y}{2}}}{\displaystyle {\int }_{x_k-\frac{\mathrm{\Delta }x}{2}}^{x_k+\frac{\mathrm{\Delta }x}{2}}}$$$$\times {\displaystyle \frac{d{x}^{\prime }d{y}^{\prime }}{\sqrt{{(x-x^{\prime })}^2+{(y-y^{\prime })}^2}}}$$

em que N e M representam, respectivamente, os números de seções na direção de x e na direção y. O número total de seções será N × M. Logo, o potencial no centro da seção i, j_𝒞, será

(38)$${\mathrm{\Phi }}_{i,j}\equiv \mathrm{\Phi }(x_i,y_j)=\sum _{k=1}^N\sum _{l=1}^M{\sigma }_{k,l}{g}_{i,j,k,l}$$

em que

(39)$$g_{i,j,k,l}=-{\displaystyle {\int }_{y_l-\frac{\mathrm{\Delta }y}{2}}^{y_l+\frac{\mathrm{\Delta }y}{2}}}{\displaystyle {\int }_{x_k-\frac{\mathrm{\Delta }x}{2}}^{x_k+\frac{\mathrm{\Delta }x}{2}}}{\displaystyle \frac{d{x}^{\prime }d{y}^{\prime }}{\sqrt{{(x_i-x^{\prime })}^2+{(y_j-y^{\prime })}^2}}}$$

Fazendo as seguintes mudanças de variáveis, x′ = x_k + u e y′ = y_l + v, na integral anterior, e considerando que x_i − x_k = [i − k]Δx e y_j − y_l = [j − l]Δy, obtemos

(40)$$g_{i,j,k,l}$$$$\mathrm{}={\displaystyle {\int }_{-\frac{\mathrm{\Delta }y}{2}}^{\frac{\mathrm{\Delta }y}{2}}}{\displaystyle {\int }_{-\frac{\mathrm{\Delta }x}{2}}^{\frac{\mathrm{\Delta }x}{2}}}{\displaystyle \frac{-dudv}{\sqrt{{([i-k]\mathrm{\Delta }x-u)}^2+{([j-l]\mathrm{\Delta }y-v)}^2}}}$$

em que tal como no caso unidimensional, a expressão g_i,j,k,l depende somente das diferenças, i - k e j - l, então podemos redefinir simplesmente como g_{i, j, k, l} = g_{i − k, j − l}, e resolvendo a integral (40), obtemos

(41)$$g_{r,s}=-f(\frac{\mathrm{\Delta }x}{2}-r\mathrm{\Delta }x,\frac{\mathrm{\Delta }y}{2}-s\mathrm{\Delta }y)$$$$+f(\frac{\mathrm{\Delta }x}{2}-r\mathrm{\Delta }x,-\frac{\mathrm{\Delta }y}{2}-s\mathrm{\Delta }y)$$$$+f(-\frac{\mathrm{\Delta }x}{2}-r\mathrm{\Delta }x,\frac{\mathrm{\Delta }y}{2}-s\mathrm{\Delta }y)$$$$-f(\frac{-\mathrm{\Delta }x}{2}-r\mathrm{\Delta }x,-\frac{\mathrm{\Delta }y}{2}-s\mathrm{\Delta }y)$$

em que r = i − k, s = j − l e

(42)$$f(x,y)=y\ln \left[x+\sqrt{x^2+y^2}\right]+x\ln \left[y+\sqrt{x^2+y^2}\right]-y$$

Agora podemos reescrever a equação (38) como

(43)$${\mathrm{\Phi }}_{i,j}=\sum _{k=1}^N\sum _{l=1}^M{g}_{i-k,j-l}{\sigma }_{k,l}$$

podemos notar que também esta equação tem a aparência do teorema de convolução, mas para isso é preciso que os conjuntos {Φ_i,j}, {g_j,j} e {σ_i,j} sejam periódicos e tenham o mesmo número de elementos. Para que isso seja possível, podemos seguir um procedimento semelhante que no caso unidimensional. Assim, obtemos para σ

(44)$${\sigma }_{i,j}=\{\begin{array}{ccc}\hfill {\sigma }_{i,j}\hfill & \hfill \text{se}\hfill & \hfill 1\le i\le N\text{e}1\le j\le M\hfill \\ \hfill 0\hfill & \hfill \text{se}\hfill & \hfill N+1\le i\le 2N\text{ou}M+1\le j\le 2M\hfill \end{array}$$

e para g temos

(45)$$g_{i,j}=g_{i^{\prime },j^{\prime }}$$

em que

(46)$$i^{\prime }=\{\begin{array}{ccc}i\hfill & \hfill \text{se}\hfill & \hfill 0\le i\le N\hfill \\ i-2N\hfill & \hfill \text{se}\hfill & \hfill N+1\le i\le 2N\hfill \end{array}$$

e

(47)$$j^{\prime }=\{\begin{array}{ccc}j\hfill & \hfill \text{se}\hfill & \hfill 1\le j\le M\hfill \\ j-2M\hfill & \hfill \text{se}\hfill & \hfill M+1\le j\le 2M\hfill \end{array}$$

e Φ_i,j segue a mesma regra que g_i,j. Lembrando que agora os conjuntos {Φ_i,j}, {g_j,j} e {σ_i,j} são periódicos, com períodos 2N e 2M em relação aos subíndices i e j, respectivamente. Logo finalmente, a equação (43) pode ser reescrita como

(48)$${\mathrm{\Phi }}_{i,j}=\sum _{j=1}^{2N}\mathit{}\sum _{j=1}^{2M}{g}_{i-k,j-l}{\sigma }_{k,l}$$

e com a ajuda do teorema de convolução discreta de Fourier, obtemos

(49)$${\widehat{\mathrm{\Phi }}}_{\alpha ,\beta }=\sqrt{2N\cdot 2M}{\widehat{g}}_{\alpha ,\beta }{\widehat{\lambda }}_{\alpha ,\beta },\alpha =1,2,\mathrm{\cdots },2N\mathrm{}\beta =1,2,\mathrm{\cdots },2M$$

em que os novos conjuntos {${\widehat{\mathrm{\Phi }}}_{\alpha ,\beta }$}, {${\widehat{g}}_{\alpha ,\beta }$} e {${\widehat{\lambda }}_{\alpha ,\beta }$} representam as transformadas discretas de Fourier dos conjuntos {Φ_i,j}, {g_i,j} e {λ_i,j}, respectivamente. Logo o potencial gravitacional será simplesmente a transformada inversa de Fourier de {${\widehat{\mathrm{\Phi }}}_{\alpha ,\beta }$}.

A exatidão do método numérico aqui apresentado é verificada comparando o resultado numérico Φ_C e o resultado analítico Φ se houver. Assim definimos o erro percentual, e%, de Φ_C como

(50)$$e\%=\frac{|{\mathrm{\Phi }}_C-\mathrm{\Phi }|}{{\mathrm{\Phi }}_C}\times 100\%$$

Para testar o método apresentado devemos considerar os modelos de discos achatados com extensão finita, é dizer que tenham um raio finito, além disso que tenham um potencial analítico conhecido. Para começar consideremos o primeiro membro da família de discos de Kalnajs [⁴], cuja densidade de massa é dada como

(51)$$\sigma (x,y)=\frac{3}{2\pi }\sqrt{1-x^2-y^2}$$

em que o raio do disco é a unidade, e seu potencial é dado como

(52)$$\mathrm{\Phi }(x,y)=\frac{3\pi }{8}\left(x^2+y^2-2\right)$$

essa expressão para o potencial é somente válida sobre o próprio disco, é dizer, para x² + y² ≤ 1. Fora dessa região não temos uma expressão analítica conhecida para o potencial.

Para a simulação consideramos N = M = m², em que m é um número natural. O σ_i,j é determinado avaliando o valor da função σ do modelo no centro da célula 𝒞_{i, j}. Para o disco de Kalnajs a função σ é dada pela equação (51). Na seguinte tabela apresentamos a comparação do tempo de processamento do método em função de m para o disco de Kalnajs,

m	3	4	5	6	7
tT/tC	0.9	3.2	9.1	27.4	90.5

em que t_T e t_C são, respectivamente, os tempos de processamento usando a soma tradicional e o método de convolução. De onde notamos que t_T/t_C cresce quase exponencialmente com m, e para m ≥ 4 a vantagem do método de convolução é indiscutível.

Na Figura 4 apresentamos a evolução do erro percentual e% do potencial numérica via método de convolução em função de m, em que consideramos a expressão (52) como o potencial analítico. As comparações de Φ_C e Φ são feitos somente no interior do disco, dado que fora do disco não existe um valor analítico para comparar. As regiões de cor mais escura são as regiões com maior erro percentual e as regiões de cor mais clara são as regiões com baixo erro percentual. A escala no lado direito de cada gráfico indica o valor do erro percentual em porcentagem. Podemos notar claramente que com o aumento de m o valor de e% se atenua rapidamente e que as regiões com maior valor para e% sempre ficam na fronteira do disco, por exemplo para m = 6, o erro percentual é sempre menor que 0.15%, o erro percentual é ainda menor na região central do disco, e somente na fronteira adota valores próximos de 0.15%.

Figura 4
Evolução do erro percentual do potencial via método de convolução no interior do disco de Kalnajs em função de m.

A Figura 5 mostra o potencial numérico usando o método de convolução com m = 6 para o disco de Kalnajs. Nesse caso o método também permite encontrar o potencial fora do disco, é dizer para x² + y² ≥ 1, região na qual não é conhecida uma expressão analítica para o potencial.

Figura 5
Resultado para o potencial do disco de Kalnajs usando o método de convolução com m = 6.

Outro modelo de disco considerado para testar nosso método numérico é o disco uniforme [⁷], cuja densidade é dada como

(53)$$\sigma (x,y)=\{\begin{array}{cc}\hfill \frac{1}{\pi },\mathrm{}\text{se}\mathrm{}{x}^2+y^2\le 1,\hfill & \\ & \\ \hfill 0,\mathrm{}\text{caso contrário}\hfill & \end{array}$$

em que o raio do disco é a unidade, e seu potencial analítico é dado como

(54)$$\mathrm{\Phi }(x,y)=-\frac{2}{\pi }\left[(1+r)E(\eta )+(1-r)K(\eta )\right]$$

em que $r=\sqrt{x^2+y^2}$, η = 4r/(1 + r)², K(η) é a integral elíptica completa de primeira classe, e E(η) é a integral elíptica completa de segunda classe.

Na Figura 6 apresentamos a evolução do erro percentual e% do potencial numérica em função de m para o disco uniforme, em que consideramos a expressão (54) como o potencial analítico. Neste caso as comparações de Φ_C e Φ são feitos tanto no interior e fora do disco. Novamente, as regiões de cor mais escura são as regiões com maior erro percentual e as regiões de cor mais clara são as regiões com baixo erro percentual. Também notamos que o aumento de m diminui o valor de e% e que as regiões com maior erro percentual também ficam na fronteira do disco, por exemplo para m = 6, o erro percentual na fronteira do disco é próximo de 2.5%, porém no interior e fora do disco os valores são menores.

Figura 6
Evolução do erro percentual do potencial via método de convolução para o disco uniforme em função de m.

Para cada m, no disco uniforme se apresenta maior erro percentual que no disco de kalnajs, é possível que isso se deve a que a densidade de massa do disco uniforme é uma função descontínua na fronteira do disco e do disco de Kalnajs é uma função contínua.

A Figura 7 mostra o potencial numérico usando o método de convolução com m = 6 para o disco uniforme.

Figura 7
Resultado para o potencial usando o método de convolução com m = 6, para o disco uniforme.

4. Conclusão

Foi determinado numericamente o potencial gravitacional para diferentes estruturas, para o caso unidimensional consideramos dois casos, a densidade cossenoidal e a densidade modular, nesse casos não existem resultados analíticos para fazer as comparações. Foi comparado o tempo de processamento, usando a soma tradicional (soma direta) e o método da convolução (usando a transformada discreta de Fourier), e notamos uma vantagem no uso do método de convolução para m > 7, é dizer para número de seções maiores que 2⁷ = 128.

Para o caso bidimensional encontramos analiticamente o potencial gerado por uma plataforma retangular de densidade constante (41), e logo determinamos o potencial newtoniano numérico para o disco de Kalnajs [⁴] e para o disco uniforme [⁷], e comparamos com seus resultados analíticos. Para o disco de Kalnajs com m = 6 (2⁶ × 2⁶ seções) e m = 7 (2⁷ × 2⁷ seções) o erro percentual é menor que 0.15% e 0.09%, respectivamente. Apresentando boa exatidão do método numérico. Para o disco uniforme, a concordância entre o resultado numérico e analítico é um pouco menor, para m = 6 e m = 7 o erro percentual máximo é próximo de 2.5% e 1.2%, respectivamente.

Em ambos os discos a concordância entre o valor numérico e valor analítico é mais pobre na vizinhança da fronteira do disco, possivelmente deve-se ao fato de que o método numérico utiliza seções retangulares que tem fronteiras retilíneas e paralelas aos eixos cartesianos e que difere da fronteira dos discos que são arcos de circunferência. Para o caso do disco uniforme o erro percentual na região próxima da fronteira do disco é ainda maior, possivelmente deve-se ao fato de que a função densidade é descontínua na fronteira do disco.

Também para o caso bidimensional foi comparado o tempo de processamento, usando a soma tradicional e o método da convolução, e notamos uma vantagem no uso do método de convolução para m > 3, é dizer para número de seções maiores que 2³ × 2³ = 64.

O método apresentado é um pouco parecido ao método apresentado por Fukushima [¹²], porém o método dele não usa o teorema de convolução que economiza o tempo de processamento computacional na determinação do potencial. Pórem uma limitação de nosso método é que permite determinar o potencial somente sobre um ponto interior ao disco e se for um ponto exterior precisa estar sobre o mesmo plano que contem o disco, é dizer, não permite determinar o potencial em um ponto fora do plano que contem o disco. O método de Fukushima permite determinar o potencial em qualquer ponto do espaço. Da mesma forma que o método de Fukushima, o nosso método também permite determinar o potencial gravitacional de discos que não tenham simétria axial.

Referências

• ^[1] J. Binney e S. Tremaine, Galactic Dynamics (Princeton University Press, New Jersey, 2008), 2 ed.
• ^[2] J.R. Reitz, F.J. Milford e R.W. Christy, Fundamentos da Teoria Eletromagnética (Editora Campus, Rio de Janeiro, 1982).
• ^[3] K.D. Machado, Eletromagnetismo (Editora TODAPALAVRA, Ponta Grossa, 2012), v. 1.
• ^[4] G.A. Gonzáles e J.I. Reina, MNRAS 371, 1873 (2006).
• ^[5] D. Vogt e P.S. Letelier, MNRAS 396, 1487 (2009).
• ^[6] P.S. Letelier, MNRAS 381, 1031 (2007).
• ^[7] H. Lass e L. Blitzer, Celest. Mech. Dyn. Astron. 30, 225 (1983).
• ^[8] K. de Oliveira e M. de Fátima, Astronomia & Astrofísica (Editora Livraria da Física, São Paulo, 2004), 2 ed.
• ^[9] G.B. Arfken e H.J. Weber, Física Matemática: Métodos matemáticos para engenharia e física (Elsevier, Rio de Janeiro, 2007).
• ^[10] J.W. Brown e R.V. Churchill, Variáveis complexas e aplicações (AMGH, Porto Alegre, 2015).
• ^[11] T. Fukushima, MNRAS 456, 3702 (2016).
• ^[12] T. Fukushima, MNRAS 459, 3825 (2016).
• ^[13] C.J. Alarcón LLacctarímay, Rev. Bras. Ensino Fís. 43, e20200351 (2021).
• ^[14] R.W. Hockney e J.W. Eastwood, Computer Simulation Using Particles (IOP Publishing Ltd, New York, 1988).
• ^[15] J.W. Cooley e J.W. Turkey, An Algorithm for the Machine Computation of the Complex Fourier Series, Mathematics of Computation, 19, 297 (1965).
• ^[16] G.D. Bergland, IEEE Spectrum 6, 41 (1969).

Editado por

• Editor-in-Chief:
  Marcello Ferreira https://orcid.org/0000-0003-4945-3169

Datas de Publicação

• Publicação nesta coleção
  07 Abr 2025
• Data do Fascículo
  2025

Histórico

• Recebido
  27 Nov 2024
• Revisado
  25 Fev 2025
• Aceito
  09 Mar 2025