Equações de Laplace e Poisson

De Física Computacional
Ir para navegação Ir para pesquisar

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.


Eletroestática

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

2Φ=0

ou na sua versão em 2 dimensões[1]: 2Φx2+2Φy2=0

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

2Φ=ρ(x,y)ϵ0

ou na sua versão em 2 dimensões[1]: 2Φx2+2Φy2=ρ(x,y)ϵ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). 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:

fdt=D(2fx2+2fy2)

Fazendo t, 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:

limtf(x,y,t)=f(x,y)

Com isso: ft=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

Equação de Laplace

Para equação de Laplace partimos de:

Φt=μ(2Φx2+2Φy2)

Discretizando, primeiro chegamos que:

Φi,jn+1=Φi,jn+μΔt(Φi1,jn2Φi,jn+Φi+1,jn(Δx)2+Φi,j1n2Φi,jn+Φi,j+1n(Δy)2)

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

μΔt(Δx)2+μΔt(Δy)212

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

μΔt(Δx)214

Para o algoritmo de Jacobi (Relaxação) escolhemos o valor de 14 e com isso resulta na equação final:

Φi,jn+1=14(Φi1,jn+Φi+1,jn+Φi,j1n+Φi,j+1n)

onde n representa o passo no tempo, i representa o passo em X e j representa o passo em Y. A constante μ 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

Partindo de:

Φt=2Φx2+2Φy2+ρ(x,y)ϵ0

chegamos em:

Φi,jn+1=Φi,jn+Δt(Φi1,jn2Φi,jn+Φi+1,jn(Δx)2+Φi,j1n2Φi,jn+Φi,j+1n(Δy)2)+Δtρi,jϵ0

Para nosso problema Δx=Δy, então multiplicando os dois lados por Δx2Δt, chegamos em:

Φi,jn+1Δx2Δt=Δx2ΔtΦi,jn4Φi,jn+Φi1,jn+Φi+1,jn+Φi,j1n+Φi,j+1n+Δx2ρi,jϵ0

E finalmente, aplicando a condição de estabilidade ΔtΔx2=14 e cancelando os termos Φi,jn:

Φi,jn+1=14(Φi1,jn+Φi+1,jn+Φi,j1n+Φi,j+1n+Δx2ρi,jϵ0)

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:

Φi,jn+1=14(Φi1,jn+1+Φi+1,jn+Φi,j1n+1+Φi,j+1n)

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 (14Δx2ρi,jϵ0) ao lado direito da equação. Assim temos a seguinte equação de iteração:

Φi,jn+1=14(Φi1,jn+1+Φi+1,jn+Φi,j1n+1+Φi,j+1n+Δx2ρi,jϵ0)

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 Φi,jn+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 Φ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:

Φi,jn+1=(1ω)Φi,jn+ω4(Φi1,jn+1+Φi+1,jn+Φi,j1n+1+Φi,j+1n)

Equação de Poisson

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

Φi,jn+1=(1ω)Φi,jnω4(Φi1,jn+1+Φi+1,jn+Φi,j1n+1+Φi,j+1n+Δx2ρi,jϵ0)

Todas 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×L com uma das quatro "bordas" carregada. Consideramos as dimensões de Lx e Ly como L (um plano quadrado) para facilitar os cálculos.

Condições iniciais e de contorno:

{Φ(x=0,y)=Φ0=1Φ(x=L,y)=Φ(x,y=0)=Φ(x,y=L)=0

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

Φ(x,y)=Φ0i=1,2,3,...4πn(nπxL)sinh(nπy/L)sinh(nπ)


  • 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 ρ(x,y):

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

{ρ(x=L/2,y=L/4)=1ρ(x=L/2,y=3L/4)=1Φ(x=0,y)=Φ(x=L,y)=Φ(x,y=0)Φ(x,y=L)=0

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

Φ(x,y)=14πϵ0qx2+y2

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.


Método de Jacobi

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.

Método de Gauss-Seidel

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

Método SOR (Simultaneous Overrelaxation)

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.

Resultados do problema do Dipolo Elétrico

Para o Problema do Dipolo Elétrico (através da Equação de Poisson) compararmos e analisarmos a implementação dos métodos focando no erro proporcionado pelos métodos numéricos.

Método de Gauss-Seidel

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.

UMA COISA QUE EU QUERO FALAR É SOBRE A FACILIDADE DE MODELAR UM PROBLEMA COMPLEXO DE ELETROESTÁTICA SEM TER QUE RESOLVER CONTAS, APENAS APLICAR A EQUACAO DIFERENCIAL E DEU ENTENDEU? si pá vey se pá vei si pa

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

Referências

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