Equações de Laplace e Poisson

Grupo: Augusto M Giani e Henrique Padovani

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

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: ${\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 (Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{\partial f}{dt} = D\cdot \left( \frac{\partial^{2}f}{\partial x^2} + \frac{\partial^{2}f}{\partial y^2} \right)}

Fazendo Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \lim_{t \rightarrow \infty}f(x,y,t) = f(x,y)}

Com isso: Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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.

Discretizações

Método de Jacobi "FTCS"

Equação de Laplace Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \nabla^2\Phi = 0}

Para equação de Laplace partimos de:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Phi^{n+1}_{i,j} = \Phi^{n}_{i,j} + \mu\Delta t \left(\frac{\Phi^{n}_{i-1,j} - 2\Phi^{n}_{i,j} + \Phi^{n}_{i+1,j}}{(\Delta x)^2} + \frac{\Phi^{n}_{i,j-1} - 2\Phi^{n}_{i,j} + \Phi^{n}_{i,j+1}}{(\Delta y)^2} \right)}

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

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{\mu\Delta t}{(\Delta x)^2} + \frac{\mu\Delta t}{(\Delta y)^2} \leq \frac{1}{2}}

No nosso algoritmo ultizamos Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Delta x = \Delta y} então obtivemos a condição de estabilidade:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{\mu\Delta t}{(\Delta x)^2} \leq \frac{1}{4}}

Para o algoritmo de Jacobi (Relaxação) escolhemos o valor de Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{1}{4}} e com isso resulta na equação final:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \color{MidnightBlue}\Phi^{n+1}_{i,j} = \frac{1}{4} \left(\Phi^{n}_{i-1,j} + \Phi^{n}_{i+1,j} + \Phi^{n}_{i,j-1} + \Phi^{n}_{i,j+1} \right)}

onde n representa o passo no tempo, i representa o passo em X e j representa o passo em Y. A constante Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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 Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \nabla^2\Phi = \frac{-\rho(x,y)}{\epsilon_0}}

Partindo de:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Phi^{n+1}_{i,j} = \Phi^{n}_{i,j} + \Delta t \left(\frac{\Phi^{n}_{i-1,j} - 2\Phi^{n}_{i,j} + \Phi^{n}_{i+1,j}}{(\Delta x)^2} + \frac{\Phi^{n}_{i,j-1} - 2\Phi^{n}_{i,j} + \Phi^{n}_{i,j+1}}{(\Delta y)^2} \right) + \frac{\Delta t \rho_{i,j}}{\epsilon_0}}

Para nosso problema Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Delta x = \Delta y} , então multiplicando os dois lados por Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{\Delta x^2}{\Delta t}} , chegamos em:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Phi^{n+1}_{i,j}\frac{\Delta x^2}{\Delta t} = \frac{\Delta x^2}{\Delta t}\Phi^{n}_{i,j} - 4\Phi^{n}_{i,j} + \Phi^{n}_{i-1,j} + \Phi^{n}_{i+1,j} + \Phi^{n}_{i,j-1} + \Phi^{n}_{i,j+1} + \frac{\Delta x^2 \rho_{i,j}}{\epsilon_0}}

E finalmente, aplicando a condição de estabilidade Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \frac{\Delta t} {\Delta x^2} = \frac{1}{4}} e cancelando os termos Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Phi^{n}_{i,j}} :

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \color{MidnightBlue} \Phi^{n+1}_{i,j} = \frac{1}{4}\left(\Phi^{n}_{i-1,j} + \Phi^{n}_{i+1,j} + \Phi^{n}_{i,j-1} + \Phi^{n}_{i,j+1} + \frac{\Delta x^2 \rho_{i,j}}{\epsilon_0}\right)}

Método de Gauss-Seidel

Como pode-se notar, o termo que distingue a Equação de Laplace para a Equação de Poisson é apenas o termo que soma Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \left(\frac{1}{4}\frac{\Delta x^2 \rho_{i,j}}{\epsilon_0} \right)} ao lado direito da equação. Para demonstrar as próximas discretizações, as deduções foram deixadas de lado pelo fato de que são irrelevantes, tendo entendido de onde vem as equações.

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, respectivamente para equação de Laplace e Poisson, utilizamos na nossa implementação:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \color{MidnightBlue}\Phi^{n+1}_{i,j} = \frac{1}{4}\left(\Phi^{n+1}_{i-1,j} + \Phi^{n}_{i+1,j} + \Phi^{n+1}_{i,j-1} + \Phi^{n}_{i,j+1}\right)}

e

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \color{MidnightBlue}\Phi^{n+1}_{i,j} = \frac{1}{4}\left(\Phi^{n+1}_{i-1,j} + \Phi^{n}_{i+1,j} + \Phi^{n+1}_{i,j-1} + \Phi^{n}_{i,j+1} + \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 Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Phi^{n+1}_{i,j}} é 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 Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \color{MidnightBlue}\Phi^{n+1}_{i,j} = \left(1-\omega\right)\Phi^{n}_{i,j} + \frac{\omega}{4}\left(\Phi^{n+1}_{i-1,j} + \Phi^{n}_{i+1,j} + \Phi^{n+1}_{i,j-1} + \Phi^{n}_{i,j+1}\right)}

e

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \color{MidnightBlue}\Phi^{n+1}_{i,j} = \left(1-\omega\right)\Phi^{n}_{i,j} \frac{\omega}{4}\left(\Phi^{n+1}_{i-1,j} + \Phi^{n}_{i+1,j} + \Phi^{n+1}_{i,j-1} + \Phi^{n}_{i,j+1} + \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

Solução analítica do problema 1

Equação de Laplace aplicada a um plano Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle L\times L} com uma das quatro "bordas" carregada. Consideramos as dimensões de Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle L_x} e Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle L_y} como Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle L} (um plano quadrado) para facilitar os cálculos.

Condições iniciais e de contorno:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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:

Falhou ao verificar gramática (MathML com retorno SVG ou PNG (recomendado para navegadores modernos e ferramentas de acessibilidade): Resposta inválida ("Math extension cannot connect to Restbase.") do servidor "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Phi(x,y) = \Phi_{0} \sum^{\infty}_{i=1,2,3,...}\frac{4}{\pi n} \left( \frac{n\pi x}{L} \right) \frac{\sinh(n\pi y/L)}{\sinh(n\pi)} }

Problema 2

Dipolo elétrico localizado ...

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.

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

Não sei o que escrever aqui kkk

[]

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.

[Colocar gráficos único com as 3 curvas dos 3 métodos.]