Equações de Laplace e Poisson

Grupo: Augusto M Giani e Henrique

O objetivo deste trabalho é implementar os métodos de Relaxação, Gauss-Seidel e SOR (Simultaneous OverRelaxation) em problemas de eletroestática, resolvidos pelas equações de Laplace e Poisson. Também temos como objetivo comparar seus resultados: erro entre os métodos e a solução analítica, tempo para estabilização das soluções.

Problema físico envolvendo as Equações de Laplace e Poisson ^[1]

A Equação de Laplace descreve o Potencial Elétrico ( $\Phi$ ) de uma determinada região num espaço que não possui nenhuma densidade de carga elétrica (corpo carregado):

\nabla ^{2}\Phi =0

ou na sua versão em 2 dimensões: ${\frac {\partial ^{2}\Phi }{\partial x^{2}}}+{\frac {\partial ^{2}\Phi }{\partial y^{2}}}=0$

Quando neste determinado espaço, delimitado pelas condições de contorno, existe uma densidade de carga, o campo $\Phi$ já não se iguala mais à zero, mas sim à densidade de cargas dentro daquela região, sendo descrito agora pela Equação de Poisson:

\nabla ^{2}\Phi ={\frac {-\rho (x,y)}{\epsilon _{0}}}

ou na sua versão em 2 dimensões^[1]: ${\frac {\partial ^{2}\Phi }{\partial x^{2}}}+{\frac {\partial ^{2}\Phi }{\partial y^{2}}}={\frac {-\rho (x,y)}{\epsilon _{0}}}$

Método de Relaxação

Como podemos ver ambas as equações não dependem do tempo, porém podemos usar um truque para resolver estas equações aplicando o método **FTCS** (Forward Time Central Space) em uma equação parecida, e fazer a evoluçao temporal durar tempo sufiente para a solução convergir ( $t\rightarrow \infty$ ). Esta operação é chamada de Método de Relaxação.

O que usamos para convergir à solução da Equação de Laplace foi uma equação de difusão genérica:

{\frac {\partial f}{dt}}=D\cdot \left({\frac {\partial ^{2}f}{\partial x^{2}}}+{\frac {\partial ^{2}f}{\partial y^{2}}}\right)

Fazendo $t\rightarrow \infty$ , para a equação de difusão temos a intuição que dada condição inicial estacionária, a solução não diverge e "relaxa" para uma função que não depende mais do tempo:

\lim _{t\rightarrow \infty }f(x,y,t)=f(x,y)

Com isso: ${\frac {\partial f}{\partial t}}=0$ , e chegando assim à Equação de Laplace e possibilitando chegar na discretização da Equação de Poisson. Então basicamente utiliza-se da mesma discretização de uma equação de difusão, porém a evolução temporal só serve para convergirmos à solução da Equação de Laplace com as condições iniciais que propomos. Os métodos de Jacobi, Gauss-Seidel e SOR são considerados Métodos de Relaxação^[2].

Método de Jacobi FTCS

Equação de Laplace $\nabla ^{2}\Phi =0$

Para equação de Laplace partimos de:

{\frac {\partial \Phi }{\partial t}}=\mu \left({\frac {\partial ^{2}\Phi }{\partial x^{2}}}+{\frac {\partial ^{2}\Phi }{\partial y^{2}}}\right)

Discretizando, primeiro chegamos que:

\Phi _{i,j}^{n+1}=\Phi _{i,j}^{n}+\mu \Delta t\left({\frac {\Phi _{i-1,j}^{n}-2\Phi _{i,j}^{n}+\Phi _{i+1,j}^{n}}{(\Delta x)^{2}}}+{\frac {\Phi _{i,j-1}^{n}-2\Phi _{i,j}^{n}+\Phi _{i,j+1}^{n}}{(\Delta y)^{2}}}\right)

Seguindo mesmo procedimento do método de FTCS, temos a mesma condição de estabilidade:

{\frac {\mu \Delta t}{(\Delta x)^{2}}}+{\frac {\mu \Delta t}{(\Delta y)^{2}}}\leq {\frac {1}{2}}

No nosso algoritmo ultizamos $\Delta x=\Delta y$ então obtivemos a condição de estabilidade:

{\frac {\mu \Delta t}{(\Delta x)^{2}}}\leq {\frac {1}{4}}

Para o algoritmo de Jacobi (Relaxação) escolhemos o valor de ${\frac {1}{4}}$ e com isso resulta na equação final:

\color {MidnightBlue}\Phi _{i,j}^{n+1}={\frac {1}{4}}\left(\Phi _{i-1,j}^{n}+\Phi _{i+1,j}^{n}+\Phi _{i,j-1}^{n}+\Phi _{i,j+1}^{n}\right)

onde n representa o passo no tempo, i representa o passo em X e j representa o passo em Y. A constante $\mu$ somente representava uma similaridade com a equação de difusão para demonstrar que este valor não interfere na equação final, ele sequer aparece (portanto podemos desconsiderá-lo, como faremos na equação de Poisson).

Equação de Poisson $\nabla ^{2}\Phi ={\frac {-\rho (x,y)}{\epsilon _{0}}}$

Partindo de:

{\frac {\partial \Phi }{\partial t}}={\frac {\partial ^{2}\Phi }{\partial x^{2}}}+{\frac {\partial ^{2}\Phi }{\partial y^{2}}}+{\frac {\rho (x,y)}{\epsilon _{0}}}

chegamos em:

\Phi _{i,j}^{n+1}=\Phi _{i,j}^{n}+\Delta t\left({\frac {\Phi _{i-1,j}^{n}-2\Phi _{i,j}^{n}+\Phi _{i+1,j}^{n}}{(\Delta x)^{2}}}+{\frac {\Phi _{i,j-1}^{n}-2\Phi _{i,j}^{n}+\Phi _{i,j+1}^{n}}{(\Delta y)^{2}}}\right)+{\frac {\Delta t\rho _{i,j}}{\epsilon _{0}}}

Para nosso problema $\Delta x=\Delta y$ , então multiplicando os dois lados por ${\frac {\Delta x^{2}}{\Delta t}}$ , chegamos em:

\Phi _{i,j}^{n+1}{\frac {\Delta x^{2}}{\Delta t}}={\frac {\Delta x^{2}}{\Delta t}}\Phi _{i,j}^{n}-4\Phi _{i,j}^{n}+\Phi _{i-1,j}^{n}+\Phi _{i+1,j}^{n}+\Phi _{i,j-1}^{n}+\Phi _{i,j+1}^{n}+{\frac {\Delta x^{2}\rho _{i,j}}{\epsilon _{0}}}

E finalmente, aplicando a condição de estabilidade ${\frac {\Delta t}{\Delta x^{2}}}={\frac {1}{4}}$ e cancelando os termos $\Phi _{i,j}^{n}$ :

\color {MidnightBlue}\Phi _{i,j}^{n+1}={\frac {1}{4}}\left(\Phi _{i-1,j}^{n}+\Phi _{i+1,j}^{n}+\Phi _{i,j-1}^{n}+\Phi _{i,j+1}^{n}+{\frac {\Delta x^{2}\rho _{i,j}}{\epsilon _{0}}}\right)

Método de Gauss-Seidel

O Método de Gauss-Seidel adianta (no tempo) a chegada da solução estacionária, utilizando termos que já foram calculados num passo anterior de tempo para calcular o ponto atual. Isto é, ao fazer a iteração, para calcular o valor de algum ponto são utilizados dois valores que já foram atualizados na mesma iteração e dois que foram calculados na iteração passada. Com isso é de se esperar que há uma variação mais rápida dos valroes do plano, quando comparado com o método anterior.

Equação de Laplace

Para equação de Laplace, obetemos a seguinte equação de iteração:

\color {MidnightBlue}\Phi _{i,j}^{n+1}={\frac {1}{4}}\left(\Phi _{i-1,j}^{n+1}+\Phi _{i+1,j}^{n}+\Phi _{i,j-1}^{n+1}+\Phi _{i,j+1}^{n}\right)

Equação de Poisson

Como pode-se notar, o termo que distingue a Equação de Laplace para a Equação de Poisson é apenas o termo que soma $\left({\frac {1}{4}}{\frac {\Delta x^{2}\rho _{i,j}}{\epsilon _{0}}}\right)$ ao lado direito da equação. Assim temos a seguinte equação de iteração:

\color {MidnightBlue}\Phi _{i,j}^{n+1}={\frac {1}{4}}\left(\Phi _{i-1,j}^{n+1}+\Phi _{i+1,j}^{n}+\Phi _{i,j-1}^{n+1}+\Phi _{i,j+1}^{n}+{\frac {\Delta x^{2}\rho _{i,j}}{\epsilon _{0}}}\right)

Método SOR, Simultaneous Overrelaxation

Como pode-se notar nas equações (é mais intuitivo na forma discretizada da Equação de Laplace), a atualização de um ponto $\Phi _{i,j}^{n+1}$ é feita através de uma espécie de "média" dos pontos, no tempo anterior, ao seu arredor (o ponto acima, à direita, à esquerda e abaixo na matriz dos $\Phi _{i,j}$ ). O método introduz nesse cálculo de "média" (ainda no método de Gauss-Seidel), pesos para a contribuição dos pontos da vizinhança e também um peso para o próprio ponto no tempo anterior. Respectivamente para a Equação de Laplace e para Poisson:

\color {MidnightBlue}\Phi _{i,j}^{n+1}=\left(1-\omega \right)\Phi _{i,j}^{n}+{\frac {\omega }{4}}\left(\Phi _{i-1,j}^{n+1}+\Phi _{i+1,j}^{n}+\Phi _{i,j-1}^{n+1}+\Phi _{i,j+1}^{n}\right)

e

\color {MidnightBlue}\Phi _{i,j}^{n+1}=\left(1-\omega \right)\Phi _{i,j}^{n}{\frac {\omega }{4}}\left(\Phi _{i-1,j}^{n+1}+\Phi _{i+1,j}^{n}+\Phi _{i,j-1}^{n+1}+\Phi _{i,j+1}^{n}+{\frac {\Delta x^{2}\rho _{i,j}}{\epsilon _{0}}}\right)

As equações em Azul descritas acima foram as que implementamos em nossos códigos.

Problemas propostos

Problema 1

Equação de Laplace aplicada a um plano $L\times L$ com uma das quatro "bordas" carregada. Consideramos as dimensões de $L_{x}$ e $L_{y}$ como $L$ (um plano quadrado) para facilitar os cálculos.

Condições iniciais e de contorno:

${\begin{cases}\Phi (x=0,y)=\Phi _{0}=1\\\Phi (x=L,y)=\Phi (x,y=0)=\Phi (x,y=L)=0\\\end{cases}}$

Para a solução analítica do problema temos que:

$\Phi (x,y)=\Phi _{0}\sum _{i=1,2,3,...}^{\infty }{\frac {4}{\pi n}}\left({\frac {n\pi x}{L}}\right){\frac {\sinh(n\pi y/L)}{\sinh(n\pi )}}$

Gráfico da solução analítica do problema 1
Solução analítica do problema 1

Problema 2

Equação de Poisson considerando uma dipolo elétrico como densidade de carga $\rho (x,y)$ :

${\begin{cases}\rho (x=L/2,y=L/4)=1\\\rho (x=L/2,y=3L/4)=-1\\\end{cases}}$

Implementação ^[3] ^[4]

Implementamos as simulações em Python3, no ambiente Colab da Google. Junto com as soluções numéricas também implementamos a solução analítica de um problema para compararmos com a solução numérica

Os códigos se encontram no final desta Wiki, mas uma observação geral é que além de utilizarmos as equações destacadas em Azul para implementar as soluções, é importante lembrar que o resultado final só é atingido quando iteramos as soluções no "tempo", então é preciso iterar os elementos da matriz no espaço, mas também fazer ela evoluir com o tempo, exemplo para o algoritmo de Jacobi, Equação de Laplace:

### Exemplo da evolução temporal no método de relaxação ###
# P é a matriz do potencial no tempo n
# Q é a matriz do potencial no tempo n+1

while t < tmax: # Loop temporal
  
  for i in range(1,L+1):  # Loop em x
    for j in range(1,L+1): # Loop em y
      Q[i][j] = (P[i+1][j] + P[i-1][j] + P[i][j+1] + P[i][j-1])/4 
  
  P = Q.copy()
  t = t + td

plt.plot(x,y,P)

Lembrando que estamos resolvendo o problema em 2D, por isso P e Q são matrizes, onde cada elemento representa um ponto no plano. Como pode-se ver, somente é plotado um gráfico, ou somente se é considerado como resultado final o estado final do vetor P, depois que ele sai do loop while. Esta lógica foi usada para todos os métodos que aplicamos.

Resultados

Apresentamos as seguintes soluções do problema com a Equação de Laplace pois já é um problema bem documentado com resultados mais palpáveis para compararmos e analisarmos a implementação dos métodos.

FTCS, Equação de Laplace

Obtivemos a seguinte solução pelo método explícito:

Solução numérica do problema 1.

Comparando com a iteração no passo de tempo anterior, obtivemos as seguintes diferenças entre o passo n+1 e o passo n:

Maior diferença percentual ao longo de todas iterações.

Dando mais zoom nas iterações futuras podemos ver que as variações (%) se tornam menores, indicando que o método converge para a solução estacionária.

Maior diferença percentual entre iteração temporal 0 e 100.
Maior diferença percentual entre iteração temporal 100 e 400.
Maior diferença percentual entre iteração temporal 400 e 1500 (estado final).

O maior erro aqui considerado no gráfico é do ponto do plano de 50x50 do nosso problema que apresentou maior diferença de valor em relação à iteração temporal anterior. Utilizamos o mesmo padrão de resultado de erro nos outros métodos.

Gauss-Seidel, Equação de Laplace

Obtivemos a seguinte solução pelo método de Gauss-Seidel:

Solução numérica do problema 1 utilizando método de Gauss-Seidel

Comparando com a iteração no passo de tempo anterior, obtivemos as seguintes diferenças entre o passo n+1 e o passo n:

Maior diferença percentual ao longo de todas iterações

Aplicamos zoom nas iterações futuras, onde podemos ver que as variações (%) se tornam tão pequenas que praticamente o método converge para a solução estacionária. Iremos discutir posteriormente que esta solução converge, com erro bem menor, antes do que o método de Jacobi

Maior diferença percentual entre iteração temporal 0 e 100.
Maior diferença percentual entre iteração temporal 100 e 400.
Maior diferença percentual entre iteração temporal 400 e 1500 (estado final).

SOR, Equação de Laplace

Obtivemos a seguinte solução aprimorando o método de Gauss-Seidel para Simultaneous OverRelaxation:

Solução numérica do problema 1 utilizando método SOR

Comparando com a iteração no passo anterior, podemos fazer a seguinte análise:

Maior diferença percentual ao longo de todas iterações

Aplicando zoom nas iterações futuras, onde podemos notar um comportamento anômalo no "erro" em iterações múltiplas de 50 (tamanho da grade).

Maior diferença percentual entre iteração temporal 0 e 60.
Maior diferença percentual entre iteração temporal 30 e 80.
Maior diferença percentual entre iteração temporal 70 e 100.
Maior diferença percentual entre iteração temporal 90 e 160.

Métodos com a Equação de Poisson

mostrar os gráficos de erro para os métodos de FTCS e SOR (já feitos) para o problema de dipolo

legendas eq poisson - analitica

legendas eq poisson - gauss

legendas eq poisson - -SOR

legendas eq poisson - comparação das soluções

legendas eq poisson - erro médio com analitica

Discussão

Importante ressaltar que este método além de convergir mais rápido para uma solução, ele apresenta menos erro com menos iterações do que o método de Gauss-Seidel, isso se deve ao fato que ele leva em consideração o seu mesmo ponto no passo anterior para atualizá-lo no próximo passo, deixando o método mais preciso.

Importante ressaltar que para o método de Gauss uma quantidade menor de iterações no tempo, obtivemos erros menores, comparando com o método explícito, o que mostra uma eficácia maior do método.

Erro médio por método, ao longo do eixo x

Referências

↑ ^1,0 ^1,1 GRIFFITHS, David J. Introduction to electrodynamics. New Jersey: Prentice Hall, 1962.
↑ GARCIA, Alejandro L. Numerical methods for physics. Englewood Cliffs, NJ: Prentice Hall, 2000.
↑ Numerical solution of partial differential equations, Dr. Louise Olsen-Kettle. The University of Queensland School of Earth Sciences Centre for Geoscience Computing.
↑ SAYAMA, Hiroki. Introduction to the modeling and analysis of complex systems. Open SUNY Textbooks, 2015.

[griffiths-1] 1,0 ^1,1 GRIFFITHS, David J. Introduction to electrodynamics. New Jersey: Prentice Hall, 1962.

[garcia-2] GARCIA, Alejandro L. Numerical methods for physics. Englewood Cliffs, NJ: Prentice Hall, 2000.

[metcompc1-3] Numerical solution of partial differential equations, Dr. Louise Olsen-Kettle. The University of Queensland School of Earth Sciences Centre for Geoscience Computing.

[sayama-4] SAYAMA, Hiroki. Introduction to the modeling and analysis of complex systems. Open SUNY Textbooks, 2015.

[1]

[2]

[3]

[4]

Equações de Laplace e Poisson

Índice

Problema físico envolvendo as Equações de Laplace e Poisson ^[1]

Método de Relaxação