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 eletrostática, descritos pelas equações de Laplace e Poisson^[1]. Os erros serão avaliados em comparação à solução analítica. Como sabemos o Método de Relaxação consiste em fazer uma equação evoluir no tempo até seu estado estacionário, para assim analisarmos o problema que realmente queremos, sabendo disso, também serão avaliados os erros em relação às iterações no tempo para analisarmos o quanto estamos chegando perto da solução que queremos.

Eletrostática

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^[2]:

{\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^[2]:

{\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 propostas não dependem do tempo. O Método consiste em resolver estas equações aplicando algum método numérico (aqui utilizaremos o método FTCS - Forward Time Central Space) em uma equação similar, e fazer a evolução temporal durar tempo suficiente para a solução convergir ao estado estacionário( $t\rightarrow \infty$ ).

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\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$ , assim temos a Equação de Laplace e também possibilitando chegar na discretização da Equação de Poisson da mesma maneira. Então basicamente, utiliza-se da mesma discretização de uma equação de difusão, a evolução temporal vai servir para convergirmos à solução dos problemas envolvendo a Equação de Laplace e de Poisson que propomos. Os métodos de Jacobi, Gauss-Seidel e SOR são considerados Métodos de Relaxação, .

Métodos Numéricos

Foram implementadas as seguintes discretizações para os problemas de eletrostática. As equações em azul foram implementadas em código.

Método de Jacobi

Equação de Laplace Para equação de Laplace 2D 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)

.

onde n representa o passo no tempo, i representa o passo em X e j representa o passo em Y. A constante $\mu$ representa uma constante de difusão.

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

{\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)

.

Equação de Poisson

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 acelera a convergência à solução estacionária, utilizando termos que já foram calculados num passo anterior de tempo para atualizar 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 é observado, há uma variação mais rápida dos valores do plano, quando comparado com o Método de Jacobi.

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.

Equação de Laplace

Aplicando os pesos na forma discretizada da equação de Laplace pelo método de Gauss-Seidel:

\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)

.

Equação de Poisson

Da mesma forma aplicamos os pesos na forma discretizada da equação de Poisson pelo método de Gauss-Seidel:

\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)

.

Onde:

w={\frac {2}{1+{\sqrt {(}}1-r^{2})}}

.

E $r$ é dado por:

r={\frac {1}{2}}cos({\frac {\pi }{L_{x}}})+{\frac {1}{2}}cos({\frac {\pi }{L_{y}}})

.

No nosso caso, temos que $L_{x}=L_{y}=L$ , portanto o $w_{opt}$ é:

w_{opt}={\frac {2}{1+sin({\frac {\pi }{L}})}}

.

Como mencionado no início, para registro, as equações em Azul descritas acima foram implementadas em nossos códigos.

Problemas propostos

Problema da borda carregada

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 )}}$

Problema da borda carregada eletricamente.
Gráfico da solução analítica.

Problema do Dipolo Elétrico

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

Utilizamos como condições iniciais e de contorno do problema:

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

Temos a solução do potencial eletroestático de cada partícula no plano dada por:

$\Phi (x,y)={\frac {1}{4\pi \epsilon _{0}}}{\frac {q}{\sqrt {x^{2}+y^{2}}}}$

Onde consideramos as constantes com sendo = 1 para facilitar a modelagem do problema.

Dipolo Elétrico.
Solução analítica para o problema do Dipolo.

Implementação

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^[3] ^[4].

Resultados do problema da Borda Carregada

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.