Evolução temporal da função de distribuição de sistemas unidimensionais

Temporal evolution of the distribution function of one-dimensional systems

Resumos

A dinâmica no espaço de fase cumpre papel importante em sistemas de muitas partículas que caminham para o equilíbrio. Neste trabalho, estudamos tal dinâmica em alguns sistemas unidimensionais: (i) pêndulos não lineares e não interagentes; e (ii) muitas partículas interagindo gravitacionalmente. A evolução temporal de sistemas com muitas partículas, em regime não colisional, é governada pela equação de Vlasov, a qual resolvemos numericamente. As equações originais são separadas em equações de transporte via método splitting de Lie. Usa-se o método numérico fluxo positivo e conservativo, com um limitador de inclinação de tipo monotonized central-difference (MC). Apesar de apresentarmos soluções de problemas comuns da área, as discussões detalhadas não são frequentes. Logo este trabalho pode ser explorado no ensino de física.

Palavras-chave
física estatística; mistura no espaço de fase; métodos numéricos; sistemas contínuos; equação de Vlasov

Phase space dynamics plays an important role in systems of many particles that move towards equilibrium. In this work, we study such dynamics in some one-dimensional systems: (i) non-linear and non-interacting pendulums; and (ii) many particles interacting gravitationally. The temporal evolution of systems with many particles, in a non-collisional regime, is governed by the Vlasov equation, which we solve numerically. The original equations are separated into transport equations via Lie’s splitting method. The positive and flux conservative numerical method is used, with a monotonized central-difference (MC) slope limiter. Although we present solutions to common problems in the area, detailed discussions are not frequent. Soon this work can be explored in the teaching of physics.

Keywords
statistical physics; phase space mixing; numerical methods; continuous systems; Vlasov equation

1. Introdução

Para sistemas formados por um número elevado de partículas, em que o potencial de interação seja de longo alcance (como sistemas estelares), a dinâmica de qualquer partícula é governada principalmente pelo potencial médio do sistema. Em outras palavras, a influência da interação de dois corpos entre as partículas é muito baixa. Este fato permite fazer uma importante simplificação: tanto a distribuição das partículas no espaço, quanto o potencial do sistema podem ser aproximados por funções suaves da posição (aproximação contínua). Isto permite definir a função de distribuição (FD) f, tal que $f(\overrightarrow{x},\overrightarrow{v},t)$d${}^3\overrightarrow{x}$d${}^3\overrightarrow{v}$, seja a probabilidade que no tempo t, uma partícula aleatoriamente escolhida se encontre no volume hexa-dimensional d${}^3\overrightarrow{x}$d${}^3\overrightarrow{v}$ do espaço de fase, em torno da posição $\overrightarrow{x}$ e velocidade $\overrightarrow{v}$ [¹]. A evolução temporal desta FD é governada pela equação de continuidade, conhecida como a equação de Vlasov ou a equação de Boltzmann sem colisões,

(1) $$\frac{\partial f}{\partial t}+\overrightarrow{v}\cdot \frac{\partial f}{\partial \overrightarrow{x}}-\frac{\partial \mathrm{\Phi }}{\partial \overrightarrow{x}}\cdot \frac{\partial f}{\partial \overrightarrow{v}}=0,$$

em que Φ denota o potencial médio do sistema. Se existir um potencial externo, então este também deve estar incluido em Φ.

Segundo o teorema de Liouville, a FD f permanece constante ao longo do movimento no espaço de fase de uma partícula. Isso significa que, mesmo se a FD é inicialmente uma função suave, esta ficará irregular quando as partículas forem embaralhadas no espaço de fase [²]. A FD fica tão irregular que na prática qualquer medição de f torna-se impossível, ou seja, a resolução finita ou discreta do sistema quebra a validade do limite contínuo. Nesse caso, o sistema é melhor descrito por uma média local de f, conhecida como função de distribuição de grão grosso$F(\overrightarrow{x},\overrightarrow{v},t)$ [³, ⁴], definida como

(2) $$F(\overrightarrow{x},\overrightarrow{v},t)=\frac{1}{\mathrm{\Delta }\mu }{\int }_{\mathrm{\Delta }\mu }f({\overrightarrow{x}}^{\prime },{\overrightarrow{v}}^{\prime },t){\text{d}}^3{\overrightarrow{x}}^{\prime }{\text{d}}^3{\overrightarrow{v}}^{\prime },$$

em que $\mathrm{\Delta }\mu ={\mathrm{\Delta }}^3\overrightarrow{x}{\mathrm{\Delta }}^3\overrightarrow{v}$ é o volume de uma região pequena, porém finita, do espaço de fase, centrada no ponto $(\overrightarrow{x},\overrightarrow{v})$. Esta região é usualmente chamada de macro-célula. Ou seja, a FD F de grão grosso é o resultado da mistura (média aritmética) dos diferentes valores da FD de grão fino f, sobre uma determinada macro-célula. Como consequência, a FD de grão grosso cumpre o papel de suavizar as irregularidades da FD de grão fino. A este processo denomina-se mistura no espaço de fase.

Também, em sentido estrito, a FD de grão fino não atinge uma distribuição de equilíbrio, mas desenvolve filamentos cada vez mais estreitos. Entretanto, espera-se que a FD de grão grosso atinja um estado de equilíbrio, como produto da mistura. Geralmente é sugerido que esse estado de equilíbrio seja uma solução estacionária específica da equação de Vlasov [⁵].

Como uma forma de quantificar a mistura no espaço de fase, define-se os funcionais-H de Boltzmann [⁴] (Intuitivamente, uma funcional pode-se considerar como uma função de uma função)

(3) $$H=-\int C(F){\text{d}}^3\overrightarrow{x}{\text{d}}^3\overrightarrow{v},$$

em que C deve ser uma função convexa com C(0) = 0. Por exemplo, C(F) = F² e C(F) = Fln⁡F. Assim, dizemos que uma função F₂ é mais misturada que F₁, se para qualquer funcional-H, cumpri-se H[F₂] > H[F₁]. Na evolução temporal de sistemas não colisionais todos os funcionais-H sempre aumentam com o tempo [⁶, ⁴], ou seja a FD de grão grosso evolui de modo que ela fique cada vez mais misturada.

A mistura oferece uma vantagem ao suavizar a FD, porém, este processo leva inevitavelmente à perda de informação sobre os detalhes mais finos da FD. Ou seja, um processo de mistura leva ao aumento de grau de ignorância (entropia) sobre a FD.

O funcional-H para C(F) = Fln⁡F é igual a k_B (constante de Boltzmann) vezes a entropia, e para sistemas colisionais é a única funcional-H que cresce durante uma evolução temporal [³]. Assim, para este tipo de sistemas, quando sua FD for mais misturada será também maior sua entropia (maior grau de ignorância).

Denomina-se relaxação violenta à evolução temporal de um sistema governado pelo seu próprio potencial, em que este potencial sofre mudanças rápidas e violentas durante a evolução do sistema. Em uma relaxação violenta, o sistema tende para o estado de equilíbrio, mas não consegue atingi-lo. As flutuações do potencial desaparecem antes que o processo de relaxamento seja concluído [⁶]. Um estado de equilíbrio, se existir, deve maximizar os funcionais-H de Boltzmann.

Em dinâmica galáctica o principal objetivo é estudar a evolução de um sistema estelar auto-gravitante, resolvendo a equação de Vlasov (1) e a seguinte integral de Poisson

(4) $$\mathrm{\Phi }(\overrightarrow{x},t)=-GM\int \frac{f({\overrightarrow{x}}^{\prime },\overrightarrow{v},t)}{|\overrightarrow{x}-{\overrightarrow{x}}^{\prime }|}{\text{d}}^3{\overrightarrow{x}}^{\prime }{\text{d}}^3\overrightarrow{v},$$

em que G é a constante gravitacional e M é a massa total do sistema. Contudo, destaca-se que soluções analíticas são restritas a soluções estacionárias ou quase estacionárias. Para isso, vários conceitos teóricos foram desenvolvidos, como o teorema de Jeans e a teoria das perturbações [¹]. A dinâmica galáctica fora do equilíbrio, por outro lado, como em uma fusão ou colapso de galáxias, é quase inteiramente tratada com simulações de N-corpo [¹, ³], que são programas computacionais que seguem o movimento de um grande número de elementos de massas (partículas), sujeitos a sua atração gravitacional mútua. Entretanto, o método de N-corpo funciona somente para um número limitado de partículas, e os sistemas realistas geralmente contêm muito mais partículas do que é possível seguir em um computador. As galáxias tipicamente contêm aproximadamente 10¹¹ partículas (estrelas) [¹]. Atualmente, mesmo os maiores computadores não podem funcionar eficientemente com mais de 10¹⁰ partículas [¹]. O potencial gravitacional em uma simulação de N-corpo não é geralmente uma função suave, e como consequência pode aparecer o relaxamento de dois corpos, o que não é desejável para sistemas sem colisões [⁷].

Outro problema inerente às simulações de N-corpo é o aparecimento das dispersões com grandes ângulos de dispersão, causado pelos encontros com baixo parâmetro de impacto, o que pode ser evitado introduzindo o softening gravitacional. Para superar essas deficiências das simulações de N-corpo, várias abordagens alternativas foram exploradas [⁷].

A abordagem comum para simulações numéricas dos sistemas autogravitantes sem colisão seria a integração direta da equação de Vlasov, combinada com a equação de Poisson [⁷, ⁸]. A integração direta da equação de Vlasov para problemas em uma ou duas dimensões espaciais é aplicada sem inconvenientes, porém a resolução em mais dimensões requer memória e tempo computacional elevados. Entretanto, avanços computacionais recentes permitiram desenvolvimentos nessa área. Os supercomputadores tornaram possível simular sistemas autogravitantes sem colisão no espaço de fase hexa-dimensional, integrando numericamente as equações de Vlasov-Poisson com resolução cientificamente significativa [⁷, ⁹].

Neste trabalho, estudaremos a evolução de alguns sistemas unidimensionais, resolvendo a equação de Vlasov e a integral de Poisson. O estudo desses sistemas pode contribuir com a compreensão da física estatística para sistemas com interação de longo alcance, e facilitar seu ensino na graduação. Esse é um assunto muito pouco tratado na graduação. Para resolver a equação de Vlasov usaremos o método de fluxo positivo e conservativo (PFC, do inglês: positive and flux conservative ) [⁸], que será apresentado na Seção ².

Na primeira parte da Seção ³, estudaremos a evolução de um sistema de pêndulos. Escolhemos este tipo de sistema por ser simples e apresentar na evolução da sua FD a mistura no espaço de fase. Aqui, entenderemos em detalhes como acontece a mistura e como influencía na evolução do sistema.

Na segunda parte da Seção ³, estudaremos a evolução de um sistema formado por um número muito grande de partículas que interagem entre si por um potencial quase gravitacional, neste caso deve acontecer a relaxação violenta. O estudo deste tipo de sistema nos permitirá entender como um sistema evolui para um estado de equilíbrio que seja autoconsistente com seu potencial médio. Também analisaremos se ao impor diferentes condições iniciais, o sistema tende a alcançar o mesmo estado de equilíbrio. É importante mencionar que, existem muito poucos trabalhos sobre o assunto na língua portuguesa.

A vantagem de estudar sistemas unidimensionais consiste na possibilidade de analisar sua FD sem complicações e observar todos os detalhes da sua evolução temporal. Para sistemas com dimensões superiores isso não é possível. Todas as simulações serão feita usando um computador ordinário.

2. Método numérico para a solução da equação de Vlasov

Para estudar a evolução de sistemas em uma dimensão, é necessário resolver a versão unidimensional da equação de Vlasov,

(5) $$\frac{\partial f}{\partial t}+v\frac{\partial f}{\partial x}-a\frac{\partial f}{\partial v}=0,$$

em que a = ∂⁡Φ/∂⁡x e nesse caso f = f(x,v,t). A solução numérica consiste em duas etapas: separar a Equação (5) em duas equações de transporte, e resolvê-las.

2.1. Método separador temporal

A Equação (5) é uma equação diferencial de primeira ordem com três variáveis: x, y, e t. O método separador temporal (ou the time-split method [¹⁰]) consiste em transformar (5) em duas equações de transporte unidimensionais. A fundamentação matemática deste método, começa por reescrever a Equação (5) como

(6) $$\frac{\partial f}{\partial t}=(A+B)f,f(0)=g,$$

em que g é a FD inicial, A e B são os operadores diferenciais

(7) $$A=-v\frac{\partial }{\partial x}\mathit{}\text{e}\mathit{}B=a\frac{\partial }{\partial v}.$$

Formalmente, a solução exata da Equação (6) será

(8) $$f(t)={\text{e}}^{t(A+B)}g={\text{e}}^{tA+tB}g.$$

Nessa parte, se realiza a seguinte aproximação

(9) $${\text{e}}^{tA+tB}\approx {\text{e}}^{tA}{\text{e}}^{tB},$$

que denomina-se Separação de Lie [¹⁰]. Esta é uma aproximação de primeira ordem em t, cujo erro local é

(10) $${\text{e}}^{tA}{\text{e}}^{tB}-{\text{e}}^{tA+tB}=\frac{t^2}{2}(AB-BA)+0(t^3).$$

Logo, o erro pode ser desprezível quando se considera um valor pequeno para t. Com essa aproximação, a solução da Equação (6) é

(11) $$f(t)={\text{e}}^{tA}{\text{e}}^{tB}g.$$

Agora, é possível reconhecer a expressão e^tBg como a solução da seguinte equação

(12) $$\frac{\partial f}{\partial t}=Bf\mathit{}\text{ou}\mathit{}\frac{\partial f}{\partial t}-a\frac{\partial f}{\partial v}=0,$$

com a condição inicial g. Finalmente, f(t) = e^tA(e^tBg) deve ser a solução da equação

(13) $$\frac{\partial f}{\partial t}=Af\mathit{}\text{ou}\mathit{}\frac{\partial f}{\partial t}+v\frac{\partial f}{\partial x}=0,$$

com a “condição inicial” e^tBg. As Equações (12) e (13) são equações de transporte unidimensionais. Em resumo, resolver a equação de Vlasov (5) implica resolver as Equações mais simples (12) e (13), uma depois da outra. Essa é a essência do método separador temporal.

2.2. Solucionador da equação de transporte

As Equações (12) e (13) têm a forma

(14) $$\frac{\partial f}{\partial t}+c\frac{\partial f}{\partial z}=0,$$

em que c é uma constante na maioria dos casos, z representa as variáveis x ou v, e f = f(t,z).

Para resolver numericamente a Equação (14), é necessário definir o domínio da variável z, que será −L/2≤z≤L/2. Logo, vamos dividir esse domínio em N células, denotadas como 𝒞_i (i = 1, 2, …, N). Todas as células terão larguras iguais, Δz = L/N, e o centro da i-ésima célula é

(15) $$z_i=-\frac{L}{2}+\frac{2i-1}{2}\mathrm{\Delta }z.$$

O tempo t é discretizado considerando os valores t^k = kΔt, em que k = 0, 1, 2, … e Δt é um pequeno intervalo de tempo. O valor médio da função f sobre a célula 𝒞_i no tempo t^k será

(16) $$f_i^k=\frac{1}{\mathrm{\Delta }z}{\int }_{z_{i-\frac{1}{2}}}^{z_{i+\frac{1}{2}}}f(t^k,z)\mathit{d}z,$$

em que $z_{i-\frac{1}{2}}$ e $z_{i+\frac{1}{2}}$ representam as fronteiras da célula. Denominamos $f_i^k$ como o valor da função na célula 𝒞ⁱ no tempo t^k. Desta forma, $f_i^k$ desempenha o papel da FD de grão grosso.

Vamos descrever brevemente o método PFC[⁸, ⁷]. Supomos que os valores de $f_i^k$ no tempo t^k são conhecidos para todas as células. Dado que a equação de transporte é uma equação conservativa, então

(17) $${\int }_{z_{i-1/2}}^{z_{i+1/2}}f(t^{k+1},z)\mathit{d}z={\int }_{Z(t^k,t^{k+1},z_{i-1/2})}^{Z(t^k,t^{k+1},z_{i+1/2})}f(t^k,z)\mathit{d}z,$$

em que t^{k + 1} = t^k + Δt, Z(t₁,t₂,z) é o valor da coordenada-z da curva característiva no tempo t₁ e que no tempo t₂ origina z como a coordenada-z. Agora, usando a propriedade da integral definida

$${\displaystyle {\int }_{Z(t^k,t^{k+1},z_{i-1/2})}^{Z(t^k,t^{k+1},z_{i+1/2})}}={\displaystyle {\int }_{Z(t^k,t^{k+1},z_{i-1/2})}^{z_{i-1/2}}}+{\displaystyle {\int }_{z_{i-1/2}}^{z_{i+1/2}}}-{\displaystyle {\int }_{Z(t^k,t^{k+1},z_{i+1/2})}^{z_{i+1/2}}},$$

podemos reescrever a Equação (17) como

$$f_i^{k+1}={\displaystyle {\int }_{Z(t^k,t^{k+1},z_{i-1/2})}^{z_{i-1/2}}}f(t^k,z)\mathit{d}z+f_i^k-{\displaystyle {\int }_{Z(t^k,t^{k+1},z_{i+1/2})}^{z_{i+1/2}}}f(t^k,z)\mathit{d}z.f_i^{k+1}$$

A seguir, denotamos

(18) $${\mathrm{\Phi }}^{-}=\frac{1}{\mathrm{\Delta }z}{\int }_{Z(t^k,t^{k+1},z_{i-1/2})}^{z_{i-1/2}}f(t^k,z)\mathit{d}z.$$

e

(19) $${\mathrm{\Phi }}^{+}=\frac{1}{\mathrm{\Delta }z}{\int }_{Z(t^k,t^{k+1},z_{i+1/2})}^{z_{i+1/2}}f(t^k,z)\mathit{d}z,$$

Assim, pode-se reescrever a expressão para $f_i^{k+1}$ como

(20) $$f_i^{k+1}=f_i^k+{\mathrm{\Phi }}_i^{-}-{\mathrm{\Phi }}_i^{+},$$

${\mathrm{\Phi }}_i^{+}$ e ${\mathrm{\Phi }}_i^{-}$ se calcula por interpolação dos valores da função nas células. Este é a base do método PFC. Pode-se interpretar Φ⁺ e Φ⁻ como o fluxo de f à direita e à esquerda da célula 𝒞_i, respectivamente.

O próximo passo considera a seguinte aproximação da função f para z ∈ 𝒞_i

(21) $$f(t,z)\approx f_i(t)+{\sigma }_i(z-z_i),$$

em que f_i(t) é o valor médio de f em C_i, e σ_i será o limitador de inclinação MC limiter [¹⁰, ¹¹, ¹², ¹³, ¹⁴], dado por

(22) $${\sigma }_i=\frac{1}{\mathrm{\Delta }z}\times \text{minmod}(\frac{f_{i+1}-f_{i-1}}{2},\text{minmod}(2[f_i-f_{i-1}],2[f_{i+1}-f_i])),$$

em que a função minmod definem-se, para dois quaisquer número real, α e β, como

(23) $$\text{minmod}[\alpha ,\beta ]=\{\begin{array}{cc}\text{min}(|\alpha |,|\beta |)\cdot \text{sgn}(\alpha ),\hfill & \hfill \text{se}\alpha \beta >0\hfill \\ & \\ 0,\hfill & \hfill \text{se}\alpha \beta \le 0\hfill \end{array}$$

O papel do limitador de inclinação é evitar as oscilações não desejadas na solução numérica. Logo, pode-se determinar o valor dos fluxos ${\mathrm{\Phi }}_i^{+}$ e ${\mathrm{\Phi }}_i^{-}$. Para isso, devemos primeiro encontrar a célula 𝒞_j que contenha o ponto Z(t^k, t^{k + 1}, z_{i + 1/2}), ou seja, tal que Z(t^k, t^{k + 1}, z_{i + 1/2}) ∈ 𝒞_j. Logo, define-se α_i como

(24) $${\alpha }_i=\{\begin{array}{cc}\hfill z_{j+1/2}-Z(t^k,t^{k+1},z_{i+1/2}),\hfill & \hfill \text{se}c>0\hfill \\ & \\ \hfill z_{j-1/2}-Z(t^k,t^{k+1},z_{i+1/2}),\hfill & \hfill \text{se}c<0\hfill \end{array}$$

Reforçando que j é o índice da célula que contém Z(t^k, t^{k + 1}, z_{i + 1/2}). Assim, α_i deve satisfazer a seguinte desigualdade, 0≤|α_i|≤Δz, e o fluxo à direita será

(25) $${\mathrm{\Phi }}_i^{+}=\{\begin{array}{c}\hfill \frac{{\alpha }_i}{\mathrm{\Delta }z}\left[f_j+\frac{{\sigma }_j\mathrm{\Delta }z}{2}\left(1-\frac{{\alpha }_i}{\mathrm{\Delta }z}\right)\right]\hfill \\ \hfill +\sum _{k=j+1}^i{f}_k,\text{se}c>0\hfill \\ \\ \hfill \frac{{\alpha }_i}{\mathrm{\Delta }z}\left[f_j-\frac{{\sigma }_j\mathrm{\Delta }z}{2}\left(1+\frac{{\alpha }_i}{\mathrm{\Delta }z}\right)\right]\hfill \\ \hfill -\sum _{k=i+1}^{j-1}{f}_k,\text{se}c<0\hfill \end{array}.$$

O fluxo a esquerda pode ser determinado a partir do fluxo a direita, fazendo ${\mathrm{\Phi }}_i^{-}={\mathrm{\Phi }}_{i-1}^{+}$, isso quer dizer, que não precisamos calcular os ${\mathrm{\Phi }}_i^{-}$ em separado, e sim somente ${\mathrm{\Phi }}_i^{+}$ para cada célula. Reduzindo o custo computacional.

2.2.1. Determinando j e α_i

As curvas características para a equação de transporte são uma família de retas de inclinação c sobre o plano tz, cuja equação pode ser escrita como z = z_o + ct, em que z_o é uma constante que representa o intercepto-z da reta. Por definição, os pontos (Z(t^k, t^{k + 1}, z_{i + 1/2}), t^k) e (z_i+1/2, t^{k + 1}) devem pertencer a uma mesma curva característica. Logo, deve-se cumprir

$$\begin{array}{l}Z(t^k,t^{k+1},z_{i+1/2})=z_o+c{t}^k,\\ z_{i+1/2}=z_o+c{t}^{k+1}.\end{array}$$

Combinando as equações anteriores obtemos

(26) $$Z(t^k,t^{k+1},z_{i+1/2})=z_{i+1/2}-c\mathrm{\Delta }t.$$

Agora, seja 𝒞_j a célula que contenha o ponto Z(t^k, t^{k + 1}, z_{i + 1/2}), então

$$z_{j-1/2}\le Z(t^k,t^{k+1},z_{i+1/2})<z_{j+1/2},$$ $$z_{j-1/2}\le z_{i+1/2}-c\mathrm{\Delta }t<z_{j+1/2}.$$

Dado que z_i = z₁ + (i−1)Δz, obtemos a desigualdade

(27) $$j\le i+1-s<j+1,$$

em que s é que o número de Courant, definido por

(28) $$s=c\mathrm{\Delta }t/\mathrm{\Delta }z$$

Finalmente, da desigualdade anterior podemos deduzir que

(29) $$j=\text{floor}(i+1-s),$$

em que floor(★) retorna o maior valor inteiro menor do que ★, o que é compatível com a desigualdade (27).

Para determinar α_i, substituimos a Equação (26) na Equação (24), e obtemos

(30) $${\alpha }_i=\{\begin{array}{cc}(j-i+s)\mathrm{\Delta }z,\hfill & \hfill \text{se}c>0\hfill \\ & \\ (j-i-1+s)\mathrm{\Delta }z,\hfill & \hfill \text{se}c<0\hfill \end{array}$$

em que s é novamente o número de Courant. Umas das vantagens do método PFC é que o número de Courant não precisa satisfazer a condição de CFL (Courant-Friedrichs-Lewy) [¹⁵]. O detalhe do código em FORTRAN é apresentado no material suplementar.

2.2.2. Determinando o Δt apropriado

A escolha de N grande sempre permitirá que a solução numérica seja precisa e tenha boa resolução. Entretanto, existem duvidas sobre o valor de Δt. Um Δt pequeno precisará de mais iterações para completar a evolução temporal, e o erro acumulado em todo o processo pode ser considerável. Um Δt grande aumentará o erro de truncamento do método separador temporal. Como método de checagem, pode-se comparar soluções analíticas e numéricas para a determinação de Δt.

A solução analítica da equação de Vlasov para um sistema contínuo de partículas que realizam movimentos harmônicos simples, todas com frequência angular igual a unidade, é dada por

(31) $$f(x,v,t)=g(x\cos t-v\text{sen}t,x\text{sen}t+v\cos t)$$

em que g(x,v) é a FD inicial. Neste caso, vamos considerar a seguinte condição inicial

(32) $$g(x,v)=\{\begin{array}{cc}\frac{{\pi }^2}{4{\delta }_x{\delta }_v}\cos \left(\frac{\pi x}{{\delta }_x}\right)\cos \left(\frac{\pi v}{{\delta }_v}\right),\hfill & \text{se}\mathit{}\{\begin{array}{c}\hfill |x|<\frac{{\delta }_x}{2},\hfill \\ \\ \hfill |v|<\frac{{\delta }_v}{2}\hfill \end{array}\hfill \\ & \\ 0\hfill & \text{em outro caso},\hfill \end{array}$$

em que δ_x = 4 é a largura da distribuição de probabilidade na direção x e δ_v = 2 na direção v.

Para todas as simulações apresentadas neste artigo, serão adotados L = 8 e N = 800, tanto para x e v. Por simplicidade, trabalharemos com unidades adimensionais.

A Figura 1 apresenta a evolução temporal do erro global e_g, para diferentes valores de Δt. Inicialmente, a simulação apresenta menor erro global a medida que aumentamos o valor de Δt, mas apartir de Δt = 0.005π, a simulação voltar apresentar erros globais altos. O erro global calcula-se usando a relação [⁷]

Figura 1
Evolução temporal do erro global, e_g, em uma simulação de um sistema contínuo de partículas que realizam movimentos harmônicos simples. A simulação apresenta menor e_g para Δt = 0,005π.

(33) $$\begin{array}{c}{e}_g(t)=\frac{{\displaystyle \sum _{i,j}|f_{i,j}(t)-f_e(x_i,v_j,t)|\text{Δ}x\text{Δ}v}}{{\displaystyle \int f_e(x,v,t)dxdv}}\\ =\text{Δ}x\text{Δ}v{\displaystyle \sum _{i,j}|f_{i,j}(t)-f_e(x_i,v_j,t)|}\end{array}$$

em que f_i,j(t) é a solução númerica para a célula centrada em (x_i,v_j) e f_e(x_i,v_j,t) representa sua respectiva solução analítica. A simplificação na ultima linha acontece devido a FD sempre ser uma função normalizada. A somatoria deve ser realizada sobre todas as células.

Fazendo outros testes, verificamos que o Δt que minimiza e_g depende da FD inicial, porém geralmente fica próximo de Δt = 0.01π. Por esta razão, consideramos este valor de Δt como referência para as simulações. Adicionalmente, vamos considerar 800×800 células para discretizar o espaço de fase bidimensional.

3. Aplicações em sistemas unidimensionais

3.1. Pêndulo não linear

Um pêndulo consiste de um pequeno objeto de massa m, que se encontra pendurado a um suporte por um fio inextensível de comprimento l, e que oscila em torno de um ponto fixo. A equação de movimento do pêndulo é dada por

(34) $$\ddot{\phi }+{\omega }^2\text{sen}\phi =0,$$

em que φ é o ângulo formado pelo fio e a vertical, e $\omega =\sqrt{g/l}$, em que g é a aceleração numa queda livre. Os dois pontos em cima de φ indicam a segunda derivada em relação ao tempo. Por simplicidade vamos considerar g = m = l = 1, logo ω = 1.

Quando φ não é restrito a valores pequenos, este sistema denomina-se pêndulo não linear. Para o caso em que a velocidade ângular inicial do pêndulo seja nula, este oscila com um período T, dado por [¹⁶]

(35) $$T=4\cdot 𝒦\left[{\text{sen}}^2\left({\phi }_o/2\right)\right],$$

em que φ_o é o deslocamento angular inicial e a função 𝒦 é a integral elíptica de primeira classe. Nota-se que T é uma função crescente de φ_o, tal como se vê na Figura 2. Isto significa que, um pêndulo que oscila com grandes amplitudes vai demorar mais para completar um ciclo, que um pêndulo com menor amplitude ângular. Esta propriedade faz com que um sistema de pêndulos seja interessante para o estudo de misturas no espaço de fase. Por sua vez, sabe-se que, para oscilações com pequena amplitude, o período de qualquer pêndulo tende a T_o = 2π, independentemente da condição inicial.

Figura 2
O período T = 4 ⋅ 𝒦 [sen² (φ_o/2)], para um pêndulo não linear, que parte do repouso. Em que φ_o é o deslocamento angular inicial.

Na ausência de forças dissipativas, a energia individual E de um pêndulo é uma constante de movimento, e sua dinâmica é governada pelo campo gravitacional. A equação das trajetórias que segue no espaço de fase será

(36) $$\frac{1}{2}{\stackrel{.}{\phi }}^2-\cos \phi =E.$$

Para E < 1 essas trajetórias são curvas fechadas. Os pêndulos com menor energia completam suas trajetórias mais rapidamente que aqueles mais energéticos, como consequência da Equação (35).

Consideremos agora um sistema contínuo de pêndulos não lineares, não interagentes, todos idênticos, porém com diferentes condições iniciais. Podemos definir uma FD para este sistema de pêndulos, como $f(\phi ,\stackrel{.}{\phi },t)$, e que deve satisfazer a seguinte equação de Vlasov unidimensional

(37) $$\frac{\partial f}{\partial t}+\stackrel{.}{\phi }\frac{\partial f}{\partial \phi }-\text{sen}\phi \frac{\partial f}{\partial \stackrel{.}{\phi }}=0,$$

em que −senφ é o campo externo que age sobre os pêndulos, especificamente uma componente do campo gravitacional.

Aqui também vamos considerar uma condição inicial da forma (32), porém, neste caso φ deve desempenhar o papel de x e $\stackrel{.}{\phi }$ de v, também δ_x = 3π/2 e δ_v = 1.

Daqui em diante, vamos nos referir à FD de grão grosso, simplesmente como FD, e o denotaremos simplesmente como f, em vez de F. A Figura 3 apresenta a evolução temporal da FD de um sistema de pêndulos nos instantes iniciais, onde a superfície representa o gráfico da FD no espaço de fase. A superfície gira ao redor da origem, no sentido horário, o que é uma propriedade dos sistemas que seguem a equação de Hamilton. Mas nem todas as partes da superfície giram com a mesma “velocidade angular”. Notamos um certo atraso no giro das partes mais afastadas da origem (pêndulos mais energéticos), o que é compatível com a Equação (35). Essa gradiente na “velocidade angular” causa a “deformação” da distribuição de probabilidade no espaço de fase: inicialmente tinha uma forma alongada e logo se transforma em uma “s” invertida, formando dois “braços”. Esses braços vão se enrolando cada vez mais, um sobre o outro (ver curvas de nível da FD na Figura 3).

Figura 3
Evolução temporal da FD de um sistema contínuo de pêndulos. As curvas de nível são apresentadas no plano $\phi \stackrel{.}{\phi }$.

A deformação da distribuição de probabilidade acontece de modo que esta fica mais filamentada e alongada. A Figura 4 apresenta as FDs para tempos maiores. Com o passar do tempo os braços da distribuição de probabilidade ficam mais enrolados, formando uma superfície com muitas “cristas”, ou seja, uma superfície irregular. Depois as cristas se aproximam, de modo que para tempos grandes estas cristas ficam tão próximas que se misturaram. A mistura das cristas começa nas regiões mais afastadas da origem (cristas antigas), e logo a mistura avança na direção à origem. Nas regiões em que as cristas se misturam a FD se suaviza [ver na Figura 4(b)].

Figura 4
A FD de um sistema contínuo de pêndulos para t = 10π e t = 80π. As curvas de nível são apresentadas no plano $\phi \stackrel{.}{\phi }$.

A distribuição de probabilidade na região central gira em torno da origem sem apresentar deformação apreciável, isso deve-se a que nesta região todos os pêndulos oscilam aproximadamente com o mesmo período (aproximação linear). Nesta região a mistura é quase inexistente.

Para quantificar a mistura no espaço de fase, usaremos o funcional-H de Boltzmann com C(f) = fln⁡f, e determinamos o incremento percentual de H usando a seguinte relação,

(38) $$\frac{\mathrm{\Delta }H}{H}=\frac{H(t)-H(0)}{H(0)}\times 100\%$$

Até o tempo 5π o incremento percentual de H (entropia) é no máximo 1%, tal como mostra a Figura 5(a). O pequeno incremento inicial de H indica que a mistura ainda está acontecendo muito de vagar. Para t = 10π e t = 20π o incremento de H é ∼5% e ∼20%, respectivamente. A mistura vai se tornando gradualmente acentuada.

Figura 5
(a) Evolução da funcional H de Boltzmann, para um sistema de pêndulos. (b) Evolução das energias mecânicas do sistema. (c) Relação entre f e E para o estado inicial e para t = 80π.

Quando t ≳ 70π a FD se mistura cada vez menos, tal como pode-se ver na Figura 5(a). O funcional H não consegui alcançar seu valor máximo, devido a que a mistura ainda continua acontecendo, porém mais lentamente.

As energias mecânicas do sistema são determinadas da seguinte maneira,

(39) $$\begin{array}{l}{K}_s=\frac{1}{2}{\displaystyle \int {\dot{\phi }}^{2}f(\phi ,\dot{\phi },t)d\phi d\dot{\phi }}\\ U_s={\displaystyle \int \left[-\cos \phi \right]f(\phi ,\dot{\phi },t)d\phi d\dot{\phi }}\\ E_s=K_s+U_s\end{array}$$

em que K_s, U_s, e E_s são as energias cinética, potencial e total do sistema, respectivamente. Em que −cos⁡φ é o potencial gravitacional que age sobre os pêndulos. A Figura 5(b) mostra que K_s e U_s adotam valores oscilantes com amplitudes que se atenuam com o tempo, porém a energia total do sistema se mantem aproximadamente constante, tal como é esperado, dado que o campo gravitacional que age sobre os pêndulos é conservativo. Entre os tempos 0 e ∼2π, E_s apresenta uma leve oscilação devido a difusão numérica inicial.

A Figura 5(c) mostra a FD em função da energia individual E dos pêndulos. A FD inicial proposta não apresenta nenhuma relação com a energia E, porém, a medida que passa o tempo vai se formando uma relação entre f e E, como consequência da mistura. Por exemplo, para t = 80π, para cada valor de $E\gtrsim -\frac{1}{2}$, corresponde aproximadamente um único valor da FD, assim f→f(E). Para $E\lesssim -\frac{1}{2}$ não existe relação entre f e E. Esta faixa de energia corresponde aos pendulos que se encontram próximos da origem do espaço de fase, assim nessa região a FD se mistura muito pouco.

Segundo o teorema de Jeans [¹], qualquer solução estacionária ou de equilíbrio da equação de Vlasov deve depender apenas das constantes de movimento do sistema. Em nosso caso, dado que f∼f(E) para $E\gtrsim -\frac{1}{2}$, e considerando que E é a única constante de movimento, então a FD determinada, deve ser aproximadamente uma FD de equilíbrio.

3.2. Relaxamento violento

Consideremos um sistema em forma de uma linha reta e finita, com densidade linear de massa λ. Em que a interação interpartículas seja do tipo gravitacional suavizada (para sistemas unidimensionais, uma interação gravitacional “real” gera um potencial médio divergente). Escolhemos o eixo x de modo que coincida com o sistema, logo o potencial médio no tempo t em um ponto x sobre a linha será

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

em que a integração é calculada sobre todo o sistema, e ε é o parâmetro que suaviza a interação gravitacional. Por simplicidade vamos considerar G = 1. A densidade de massa será cálculada em qualquer tempo a partir da FD, usando a relação

(41) $$\lambda (x,t)=M\int f(x,v,t)\mathit{d}v,$$

em que M é a massa total do sistema. Também, por simplicidade vamos considerar M = 1.

Numericamente, a FD será conhecida somente sobre as células discretas do espaço de fase (FD de grão grosso), assim podemos reescrever a Equação (41) como

(42) $$\lambda (x_i,t)=\sum _j{f}_{i,j}(t)\mathrm{\Delta }v\equiv {\lambda }_i(t),$$

em que λ será uma função constante por seções, adotando valor λ_i sobre a célula i-ésima. A Equação (40) torna-se

(43) $$\mathrm{\Phi }(x_i,t)=-\sum _j{\lambda }_j(t){\int }_{x_{j-1/2}}^{x_{j+1/2}}\frac{d{x}^{\prime }}{|x_i-x^{\prime }|+\epsilon }\equiv {\mathrm{\Phi }}_i(t).$$

Calculando separadamente a integral, obtemos

(44) $${\int }_{x_{j-1/2}}^{x_{j+1/2}}\frac{d{x}^{\prime }}{|x_i-x^{\prime }|+\epsilon }=\{\begin{array}{cc}2\ln \left[\frac{\mathrm{\Delta }x+2\epsilon }{2\epsilon }\right],\hfill & \hfill \text{se}i=j,\\ & \\ \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}$$

Assim, fica determinado o potencial numérico, que também adota um valor constante sobre cada célula. O campo gravitacional g_i sobre a célula i será determinado aproximando a derivada pela diferenciação finita,

(45) $$g_i=-\frac{{\mathrm{\Phi }}_{i+1}-{\mathrm{\Phi }}_{i-1}}{2\mathrm{\Delta }x},$$

que também adota valor constante sobre cada célula. Para ε pequeno, nosso potencial se aproxima melhor ao potencial gravitacional real, porém um valor muito pequeno de εpode levar a um potencial divergente. Assim, para nossas simulações escolheremos ε = 0.2.

Quando Φ depende explicitamente do tempo, tal como acontece em uma relaxação violenta, cumpri-se [¹]

(46) $$\frac{dE}{dt}=\frac{\partial \mathrm{\Phi }}{\partial t}.$$

Logo, a energia individual por unidade de massa $E=\frac{1}{2}{v}^2+\mathrm{\Phi }$ das partículas que formam o sistema não serão constantes de movimento. Contudo, quando o sistema alcançar o equilíbrio, Φ será estacionário. Neste caso, E será uma constante de movimento. Em uma relaxação violenta, as energias mecânicas do sistema serão determinadas da mesma forma que no caso dos pêndulos (Equação (39)), porém, a determinação de U_s deve ser modificada, dado que neste caso existe uma interação entre as partículas constituintes do sistema. Neste caso, U_s será

(47) $$U_s=\frac{1}{2}\int \lambda (x)\mathrm{\Phi }(x,t)\mathit{d}x.$$

Interpretamos U_s como a energia potencial armazenada pela estrutura do sistema.

Novamente, vamos considerar uma condição inicial da forma (32), com δ_x = 2 e δ_v = 1. Dado que em um processo de relaxamento violento, o potencial varia muito rapidamente, usaremos neste caso Δt = 0.002 (um valor menor que no caso dos pêndulos).

As Figuras 6 e 7 mostram a evolução temporal da FD. A distribuição de probabilidade gira em torno da origem, em que a porção mais próxima da origem (porção interna) gira mais rapidamente que a porção mais afastada da origem (porção externa).

Figura 6
Evolução da FD para um sistema contínuo de partículas que interagem com um potencial quase-gravitacional. As curvas de nível são apresentadas no plano xv.

Figura 7
Evolução da FD para um sistemas contínuo de partículas que interagem com um potencial quase-gravitacional. As curvas de nível são apresentadas no plano xv.

A porção interna da distribuição de probabilidade apresenta somente uma deformação inicial, para t ≳ 2 esta porção muda pouco. Em t ≳ 10, a porção interna alcança aproximadamente um estado de equilíbrio. Na porção externa aparecem dois braços, esses braços se filamentam e se enrolam ao redor da porção interna, para t ≳ 4 algumas partes desses braços começam a se misturar, e para t=10 as diferentes partes dos braços se misturam quase em sua totalidade. A porção interna forma uma núcleo “denso”, e a porção externa forma um halo extenso e pouco denso.

A Figura 8(a) apresenta a evolução do funcional H de Boltzmann, que quantifica a mistura. Para t < 0.5 o incremento de H é menor que 1%, a distribuição se deforma, mas apresenta pouca mistura. Em t > 1 os braços começam se filamentar consideralvemente, alongando-se cada vez mais, e provocando a mistura nas pontas dos braços (ver a Figura 6(c)). Para t > 2, H começa a crescer rapidamente, devido a que os braços començam a se misturar. Em t > 8 a velocidade de crescimento de H diminui consideralvemente, devido a que a mistura na porção interna e na sua vizinhança já vai acontecendo muito de vagar, e cada vez vai ficando mais lenta. No final, as ultimas misturas acontecem nas partes mais externas da distribuição de probabilidade. Em t > 20 a funcional H permanece quase inalterável indicando que o sistema está perto de alcançar o estado de equilíbrio.

Figura 8
(a) Evolução da funcional H de Boltzmann, para um sistemas de partículas que interagem com um potencial quase-gravitacional. (b) Evolução das energias mecânicas do sistema. (c) Relação entre f e E para o estado inicial e para t = 50.

A diferença do sistema de pêndulos, em uma relaxação violenta, o sistema tende alcançar o equilíbrio mais rapidamente, o que pode ser confirmado pelo comportamento de H.

A Figura 8(b) apresenta a evolução temporal das energias do sistema. K_s e U_s adotam valores oscilantes, e simétricos, que são típicos de sistemas conservativos. Assim, E_s é aproximadamente constante. A oscilação de K_s e U_s se atenua a medida que aumenta o tempo, alcançando um valor aproximadamente estacionário, K_s≈0.22 e U_s≈−1.15.

A Figura 8(c) mostra a relação da FD com a energia individual por unidade de massa $E=\frac{1}{2}{v}^2+\mathrm{\Phi }$ das partículas, para t = 50. Como esperado, a FD tende a ser uma função suave de E, embora que a FD inicial não tenha nenhuma relação com E. Inicialmente a energia individual das partículas se encontra dentro da faixa −2≤E≤−1, aproximadamente. A medida que passa o tempo essa faixa se estende, muitas partículas tendem a adotar energias menores que –2 e algumas poucas partículas adotam energias maiores que –1.

A Figura 9(a) mostra a evolução temporal da densidade de massa do sistema até t = 5. O sistema começa a evolução com um rápido encolhimento, formando um núcleo denso, depois disso a densidade fica oscilante, expulsando para fora do núcleo, a cada certo tempo, algumas partículas. Essas partículas expulsadas acabam formando parte dos braços na Figura 6 e 7. Quando o sistema alcança o estado de equilíbrio, as partículas expulsadas do núcleo formam parte do halo da estrutura. A oscilação da densidade do núcleo se atenua com o tempo. A Figura 9(b) apresenta a densidade de massa do sistema e potencial do sistema para t = 50. Para esse tempo a FD é praticamente estacionária, o que significa que o sistema está muito perto do estado de equilíbrio. Logo esses valores da densidade e do potencial podem ser considerados valores correspondente ao estado de equilíbrio.

Figura 9
(a) Evolução temporal da densidade de massa, λ, para um sistemas de partículas que interagem com um potencial quase-gravitacional. As curvas de nível são apresentadas no plano tx. (b) Densidade de massa e o potencial do sistema, Φ, para t = 50 (próximo do estado de equilíbrio).

Comparando a Figura 8(b) e a Figura 9(a), notamos que o sistema perde energia potencial quando o sistema se encolhe, devido as partículas ficarem mais próximas umas das outras e o potencial de interação interpartícula diminuir. Por outro lado, quando o sistema se estende ganha energia potencial, em contrapartida a energia cinética diminui.

3.2.1. Outras condições iniciais

Nesta seção, verificaremos se um sistema com diferentes condições iniciais se aproxima ao mesmo estado de equilíbrio. Para isso, serão consideradas as seguintes FD iniciais:

1. a FD de tipo I (g_I) será a função (32) com δ_x = 1 e δ_v = 2;

2. a FD de tipo II será

   (48)$$g_{II}(x,v)=\{\begin{array}{cc}\hfill \frac{1}{{\delta }_x{\delta }_v},\hfill & \text{se}\mathit{}|x|<\frac{{\delta }_x}{2}\mathit{}\text{e}\mathit{}|v|<\frac{{\delta }_v}{2}\hfill \\ & \\ \hfill 0,\hfill & \text{em outro caso}\hfill \end{array}$$
   com δ_x = 2 e δ_v = 1;

3. a FD de tipo III será

   (49)$$g_{III}(x,v)=\{\begin{array}{cc}\frac{{\pi }^2}{8}\cos \left[\pi (x-1)\right]\cos \left[\pi v\right],\hfill & \hfill |x-1|,|v|<\frac{1}{2}\\ & \\ \frac{{\pi }^2}{8}\cos \left[\pi (x+1)\right]\cos \left[\pi v\right],\hfill & \hfill |x+1|,|v|<\frac{1}{2}\\ & \\ 0,\hfill & \hfill \text{em outro caso}\end{array}$$

4. e a FD de tipo IV será

   (50)$$g_{IV}(x,v)=\{\begin{array}{cc}\frac{\pi }{4(\pi -2)R^2}\cos \left(\frac{\pi \sqrt{x^2+v^2}}{2R}\right),\hfill & \hfill x^2+v^2<R^2\\ 0,\hfill & \hfill \text{em outro caso}\end{array}$$
   com R = 1.

A Figura 10 apresenta o gráfico das FDs iniciais dos diferentes tipos. A Figura 11 apresenta as FDs no tempo t = 50. Todas as FDs apresentam um núcleo denso e um halo pouco denso. Existe uma certa semelhança entre as FDs associadas as condições iniciais de tipo I e IV. Em geral, diferentes condições iniciais levam também a diferentes FDs de “equilíbrio”, embora que suas curvas de nível sejam um pouco parecidas.

Figura 10
Quatro diferentes FDs inicial, para um sistemas de partículas que interagem com um potencial quase-gravitacional. As curvas de nível são apresentadas no plano xv.

Figura 11
FDs em t = 50, para um sistemas de partículas que interagem com um potencial quase-gravitacional, correspondente a quatro diferentes condições iniciais. As curvas de nível são apresentadas no plano xv.

A Figura 12(a) apresenta a evolução temporal da funcional H, para as diferentes condições iniciais. De onde, notamos que o sistema tende alcançar o equilíbrio mais rapidamente se considerarmos uma condição inicial do tipo I, e com as condições iniciais do tipo II e III o sistema tarda mais para se aproximar ao estado de equilíbrio. A funcional H associada a condição inicial de tipo III é a que sofre maior variação percentual, logo podemos interpretar que a FD associada a esta condição inicial sofre mais mistura, em comparação com as outras.

Figura 12
(a) Evolução do funcional H, correspondente as quatro diferentes condições iniciais. (b) Densidade de massa do sistema em t = 50 (próximo de equilíbrio), correspondente as quatro diferentes condições iniciais.

Na Figura 12(b) apresentamos a densidade de massa aproximadamente estacionária do sistema, associada as diferentes condições iniciais, de onde notamos que, diferentes condições iniciais geram estruturas com diferentes distribuições de massa. A condição inicial de tipo I leva a uma estrutura com um núcleo mais denso e com pouco halo, e uma condição inicial de tipo III leva a uma estrutura com maior halo.

A Figura 13 mostra a relação da FD com a energia individual das partículas, para as diferentes condições iniciais. Para todos os casos, a FD de “equilíbrio” aparenta ser uma função suave da energia individual.

Figura 13
Relação entre f e E em t = 50, para um sistemas de partículas que interagem com um potencial quase-gravitacional, correspondente a quatro diferentes condições iniciais.

4. Conclusão

O método PFC combinado com o método separador temporal é um método númerico robusto para resolver a equação de Vlasov, mostrando boa estabilidade.

Foi possível comprovar mediante simulações que a mistura no espaço de fase ajuda a que um sistema fora de equilíbrio se aproxime a um estado de equilíbrio. E que a medida que passa o tempo, a FD de um sistema se transforma de uma função da energia individual das partículas constituintes do sistema, cumprindo-se o teorema de Jeans. Foi possível também comprovar que a funcional H de Boltzmann quantifica corretamente a mistura no espaço de fase, mostrando com simplicidade as diferentes etapas da evolução da FD.

Conseguimos verificar que as condições iniciais impostas a um sistema influênciam na sua evolução temporal. Existem condições iniciais que permitem que o sistema se aproxime mais rapidamente ao estado de equilíbrio, e outras fazem que esse processo seja mais demorado.

Para abordar o estudo da física estatística de sistemas com interação de longo alcance, é importante a implementação de métodos que resolvam a equação de Vlasov, porém esses métodos são geralmente abordados somente nos cursos de pós-graduação. Para o caso de sistemas unidimensionais a implementação desses métodos é bastante acessível, tal como mostrado neste trabalho, e o custo computacional é baixo, e as simulações podem ser feitas em computadores domésticos. Isso permite pensar que, seria possível fazer um estudo introdutório da física estatística de sistemas com interação de longo alcance também na graduação, mas para isso seria necessário como pré-requisito, implementar uma disciplina que aborde métodos numéricos e computacionais que: (i) resolvam a equação de transporte unidimensional, a qual é a base para a solução da equação de Vlasov, e (ii) resolvam equações diferenciais não-lineares. A equação de Vlasov e suas aplicações poderiam ser estudadas no terceiro ano de um curso de graduação em física, tendo como pré-requisito o curso de equações diferenciais parciais.

Agradecimentos

Ao professor, Dr. Luis Paulo Silveira Machado, pela detalhada leitura deste manuscrito e pelas discussões científicas.

Referências

• ^1.
  J. Binney e S. Tremaine, Galactic Dynamics (Princeton University Press, New Jersey, 2008), 2nd ed.
• ^2.
  T. Fujiwara, Publ. Astron. Soc. Jpn 33, 531 (1981).
• ^3.
  W. Dehnen, Monthly Notices of the Royal Astronomical Society 360, 892 (2005).
• ^4.
  S. Tremaine, M. Hénon e D. Lynden-Bell, Monthly Notices of the Royal Astronomical Society 219, 285 (1986).
• ^5.
  P.H. Chavanis e F. Bouchet, Astron. Astrophys. 430, 771 (2005).
• ^6.
  D. Lynden-Bell, Monthly Notices of the Royal Astronomical Society 136, 101 (1967).
• ^7.
  K. Yoshikawa, N. Yoshida e M. Umemura, Astrophys. J. 762, 116 (2013).
• ^8.
  F. Filbet, E. Sonnendrucker e P. Bertrand, J. Comput. Phys. 172, 166 (2001).
• ^9.
  S. Tanaka, K. Yoshikawa, T. Minoshima e N. Yoshida, Astrophys. J. 849, 76 (2017).
• ^10.
  R.J. Leveque, Time-Split Methods for Partial Differential Equations PhD Dissertation, StanfordUniversity, Palo Alto (1982).
• ^11.
  J.W. Thomas, Numerical Partial Differential Equations: Conservation Laws and Elliptic Equations (Springer, New York, 1999).
• ^12.
  R.J. Leveque, D. Mihalas, E.A. Dorfi e E. Muller, Computational Methods for Astrophysical Fluid Flow (Springer, Berlin, 1998).
• ^13.
  R.J. Leveque, Numerical Methods for Conservation Laws (Springer Basel AG, 1992), 2nd ed.
• ^14.
  P.L. Roe, Lectures in applied mathematics 22, 163 (1985).
• ^15.
  A.R. Mitchell e D.F. Griffiths, The Finite Difference Method in Partial Differential Equations (Wiley, Norwich, 1980).
• ^16.
  A. Beléndez, C. Pascual, D.I. Méndez, T. Beléndez e C. Neipp, Rev. Bras. Ensino Fís. 29, 645 (2007).

Datas de Publicação

• Publicação nesta coleção
  08 Jan 2021
• Data do Fascículo
  2021

Histórico

• Recebido
  20 Ago 2020
• Revisado
  16 Nov 2020
• Aceito
  17 Nov 2020