Equações de Laplace e Poisson: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
 
(45 revisões intermediárias por 2 usuários não estão sendo mostradas)
Linha 1: Linha 1:
'''Grupo: Augusto M Giani e Henrique '''
'''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.
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<ref name=garcia>'''GARCIA, Alejandro L. Numerical methods for physics. Englewood Cliffs, NJ: Prentice Hall, 2000.'''</ref>. 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.




== Eletroestática ==
== Eletrostática ==


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


<center><math> \nabla^2\Phi = 0 </math></center>
<center><math> \nabla^2\Phi = 0 </math>.</center>


ou na sua versão em 2 dimensões<ref name=griffiths>GRIFFITHS, David J. Introduction to electrodynamics. New Jersey: Prentice Hall, 1962.</ref>: <math>\frac{\partial^{2}\Phi}{\partial x^2} + \frac{\partial^{2}\Phi}{\partial y^2} = 0</math>
ou na sua versão em 2 dimensões<ref name=griffiths>GRIFFITHS, David J. Introduction to electrodynamics. New Jersey: Prentice Hall, 1962.</ref>:
 
<center><math>\frac{\partial^{2}\Phi}{\partial x^2} + \frac{\partial^{2}\Phi}{\partial y^2} = 0</math>.</center>


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


<center><math>\nabla^2\Phi = \frac{-\rho(x,y)}{\epsilon_0}</math></center>
<center><math>\nabla^2\Phi = - \frac{\rho(x,y)}{\epsilon_0}</math>.</center>
 
ou na sua versão em 2 dimensões<ref name=griffiths></ref>:


ou na sua versão em 2 dimensões<ref name=griffiths></ref>: <math>\frac{\partial ^2\Phi}{\partial x^2} + \frac{\partial^{2}\Phi}{\partial y^2} = \frac{-\rho(x,y)}{\epsilon_0}</math>
<center><math>\frac{\partial ^2\Phi}{\partial x^2} + \frac{\partial^{2}\Phi}{\partial y^2} = - \frac{\rho(x,y)}{\epsilon_0}</math>.</center>


== Método de Relaxação ==
== 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 (<math>t \rightarrow \infty</math>). Esta operação é chamada de 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(<math>t \rightarrow \infty</math>).  


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


<center><math>\frac{\partial f}{dt} = D\cdot \left( \frac{\partial^{2}f}{\partial x^2} + \frac{\partial^{2}f}{\partial y^2} \right)</math></center>
<center><math>\frac{\partial f}{dt} = D \left( \frac{\partial^{2}f}{\partial x^2} + \frac{\partial^{2}f}{\partial y^2} \right)</math>.</center>.


Fazendo <math>t \rightarrow \infty</math>, 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:
Fazendo <math>t \rightarrow \infty</math>, 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:


<center><math>\lim_{t \rightarrow \infty}f(x,y,t) = f(x,y)</math></center>
<center><math>\lim_{t \rightarrow \infty}f(x,y,t) = f(x,y)</math>.</center>
 
Com isso: <math>\frac{\partial f}{\partial t} = 0</math>, 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'''.


Com isso: <math>\frac{\partial f}{\partial t} = 0</math>, 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.
== Métodos Numéricos ==  
Os métodos de Jacobi, Gauss-Seidel e SOR são considerados '''Métodos de Relaxação'''<ref name=garcia>'''GARCIA, Alejandro L. Numerical methods for physics. Englewood Cliffs, NJ: Prentice Hall, 2000.'''</ref>.


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


=== Equação de Laplace ===
=== Método de Jacobi ===
Para equação de Laplace partimos de:
 
<center><math>\frac{\partial \Phi}{\partial t} = \mu \left(\frac{\partial^{2}\Phi}{\partial x^2} + \frac{\partial^{2}\Phi}{\partial y^2} \right)</math></center>
'''Equação de Laplace'''
Para equação de Laplace 2D partimos de:
<center><math>\frac{\partial \Phi}{\partial t} = \mu \left(\frac{\partial^{2}\Phi}{\partial x^2} + \frac{\partial^{2}\Phi}{\partial y^2} \right)</math>.</center>


Discretizando, primeiro chegamos que:  
Discretizando, primeiro chegamos que:  


<center><math> \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)</math></center>
<center><math> \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)</math>.</center>


Seguindo mesmo procedimento do método de FTCS, temos a mesma condição de estabilidade:
onde '''n''' representa o passo no tempo, '''i''' representa o passo em X e '''j''' representa o passo em Y. A constante <math>\mu</math> representa uma constante de difusão.


<center><math>\frac{\mu\Delta t}{(\Delta x)^2} + \frac{\mu\Delta t}{(\Delta y)^2} \leq \frac{1}{2}</math></center>
Seguindo o método de FTCS, temos a condição de estabilidade do método:


No nosso algoritmo ultizamos <math>\Delta x = \Delta y</math> então obtivemos a condição de estabilidade:
<center><math>\frac{\mu\Delta t}{(\Delta x)^2} + \frac{\mu\Delta t}{(\Delta y)^2} \leq \frac{1}{2}</math>.</center>


<center><math>\frac{\mu\Delta t}{(\Delta x)^2} \leq \frac{1}{4}</math></center>
No nosso algoritmo ultizamos <math>\Delta x = \Delta y</math>. Então obtivemos a condição de estabilidade:
 
<center><math>\frac{\mu\Delta t}{(\Delta x)^2} \leq \frac{1}{4}</math>.</center>


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


<center><math>\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)</math></center>
<center><math>\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)</math>.</center>
 
onde '''n''' representa o passo no tempo, '''i''' representa o passo em X e '''j''' representa o passo em Y. A constante <math>\mu</math> 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 ===
'''Equação de Poisson'''


Partindo de:
Partindo de:


<center><math>\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}</math></center>
<center><math>\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}</math>.</center>


chegamos em:  
chegamos em:  


<center><math>\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}</math></center>
<center><math>\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}</math>.</center>


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


<center><math>\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}</math></center>
<center><math>\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}</math>.</center>


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


<center><math>\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)</math></center>
<center><math>\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)</math>.</center>


== Método de Gauss-Seidel ==
=== 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.
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 ===
'''Equação de Laplace'''


Para equação de Laplace, obetemos a seguinte equação de iteração:  
Para equação de Laplace, percorrendo os índices de forma crescente, obtemos a seguinte equação de iteração:  


<center><math>\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)</math></center>
<center><math>\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)</math>.</center>


=== Equação de Poisson ===
'''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 <math>\left(\frac{1}{4}\frac{\Delta x^2 \rho_{i,j}}{\epsilon_0} \right)</math> ao lado direito da equação. Assim temos a seguinte equação de iteração:
Como pode-se notar, o termo que distingue a Equação de Laplace para a Equação de Poisson é apenas o termo que soma <math>\left(\frac{1}{4}\frac{\Delta x^2 \rho_{i,j}}{\epsilon_0} \right)</math> ao lado direito da equação. Assim, temos a seguinte equação de iteração:


<center><math>\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)</math></center>
<center><math>\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)</math>.</center>


== Método SOR (Simultaneous Overrelaxation) ==
=== 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 <math>\Phi^{n+1}_{i,j}</math> é 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 <math>\Phi_{i,j}</math>). 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.  
Como pode-se notar nas equações (é mais intuitivo na forma discretizada da Equação de Laplace), a atualização de um ponto <math>\Phi^{n+1}_{i,j}</math> é 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 <math>\Phi_{i,j}</math>). 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 ===
'''Equação de Laplace'''


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


<center><math>\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)</math></center>
<center><math>\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)</math>.</center>


=== Equação de Poisson ===
'''Equação de Poisson'''


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


<center><math>\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)</math></center>
<center><math>\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)</math>.</center>
 
Onde:


Todas as equações em Azul descritas acima foram implementadas em nossos códigos.
<center><math> w = \frac{2}{1 + \sqrt(1-r^2)} </math>.</center>
 
E <math>r</math> é dado por:
 
<center><math> r = \frac{1}{2} \cos\left( \frac{\pi}{L_x} \right) + \frac{1}{2} \cos\left( \frac{\pi}{L_y} \right)</math>.</center>
 
No nosso caso, temos que <math>L_x=L_y=L</math>, portanto o <math>w_{opt}</math> é:
 
<center><math> w_{opt} = \frac{2}{1 + \sin(\frac{\pi}{L} )} </math>.</center>
 
Como mencionado no início, para registro, as equações em Azul descritas acima foram implementadas em nossos códigos.


== Problemas propostos ==
== Problemas propostos ==


=== Problema da borda carregada ===  
=== Problema da borda carregada ===


Equação de Laplace aplicada a um plano <math>L\times L</math> com uma das quatro "bordas" carregada. Consideramos as dimensões de <math>L_x</math> e <math>L_y</math> como <math>L</math> (um plano quadrado) para facilitar os cálculos.
Equação de Laplace aplicada a um plano <math>L\times L</math> com uma das quatro "bordas" carregada. Consideramos as dimensões de <math>L_x</math> e <math>L_y</math> como <math>L</math> (um plano quadrado) para facilitar os cálculos.
Linha 126: Linha 146:
Para a solução analítica do problema temos que:  
Para a solução analítica do problema temos que:  


<math>\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)} </math>
<math>\Phi(x,y) = \Phi_{0} \sum^{\infty}_{i=1,3,5,...}\frac{4}{\pi n} \left( \frac{n\pi x}{L} \right) \frac{\sinh(n\pi y/L)}{\sinh(n\pi)} </math>


<br />
<br />


<div><ul>  
<center><div><ul>  
<li style="display: inline-block;"> [[Arquivo:Problema1-image.PNG|thumb|left|437px|Problema da borda carregada eletricamente.]] </li>
<li style="display: inline-block;"> [[Arquivo:Problema1-image.PNG|thumb|left|437px|Problema da borda carregada eletricamente.]] </li>
<li style="display: inline-block;"> [[Arquivo:Escorregador analitico.png|thumb|right|400px|Gráfico da solução analítica.]] </li>
<li style="display: inline-block;"> [[Arquivo:Escorregador analitico.png|thumb|right|400px|Gráfico da solução analítica somando até o termo n=199.]] </li>
</ul></div>
</ul></div></center>


=== Problema do Dipolo Elétrico ===
=== Problema do Dipolo Elétrico ===
Linha 145: Linha 165:
\rho(x = L/2,y = L/4) = 1 \\
\rho(x = L/2,y = L/4) = 1 \\
\rho(x = L/2,y = 3L/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\\
\Phi(x = 0,y) = \Phi(x = L,y) = \Phi(x,y = 0) = \Phi(x,y = L) = 0\\
\end{cases}
\end{cases}
</math>
</math>
Linha 151: Linha 171:
Temos a solução do potencial eletroestático de cada partícula no plano dada por:
Temos a solução do potencial eletroestático de cada partícula no plano dada por:


<math>\Phi(x,y)=\frac{1}{4\pi\epsilon_{0}}\frac{q}{\sqrt{x^2+y^2}}</math>
<math>\Phi(x,y)=\frac{1}{4\pi\epsilon_{0}}\frac{q}{r_{(x,y)}}</math>
 
Onde:
 
<math>r = \sqrt{(x-x_0)^2+(y-y_0)^2}</math>


Onde consideramos as constantes com sendo = 1 para facilitar a modelagem do problema.
Onde <math>x_0</math> e <math>y_0</math> são as coordedanas da carga. Também consideramos as constantes com sendo = 1 para facilitar a modelagem do problema.


<div><ul>  
<center><div><ul>  
<li style="display: inline-block;"> [[Arquivo:Problema2.jpeg|thumb|central|500px|Dipolo Elétrico.]] </li>
<li style="display: inline-block;"> [[Arquivo:Problema2.jpeg|thumb|central|500px|Dipolo Elétrico<ref name=ref_dif>https://noic.com.br/materiais-fisica/cursos/teorico/aula-5-3-energia-potencial-eletrica-trabalho-e-potencial-eletrico/</ref>.]] </li>
<li style="display: inline-block;"> [[Arquivo:Analitico 01.png|thumb|central|350px|Solução analítica para o problema do Dipolo.]] </li>
<li style="display: inline-block;"> [[Arquivo:Analitico 01.png|thumb|central|350px|Solução analítica para o problema do Dipolo.]] </li>
<li style="display: inline-block;"> [[Arquivo:Analitico 02.png|thumb|central|350px|Solução analítica para o problema do Dipolo.]] </li>
</ul></div></center>
</ul></div>


<br />
<br />
Linha 165: Linha 188:
== Implementação  ==
== 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
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 do 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:
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 FTCS integra no espaço (2D) a equação de iteração, mas utilizando métodos de Relaxação é preciso também iterar no tempo estas equações. Um exemplo disso é mostrado abaixo:


<source lang="haskell" line='line'>
<source lang="haskell" line='line'>


### Exemplo da evolução temporal no método de relaxação ###
### Exemplo da evolução temporal no método de relaxação ###
### Exemplo para o algoritmo de jacobi, Equação de Laplace ###
# P é a matriz do potencial no tempo n
# P é a matriz do potencial no tempo n
# Q é a matriz do potencial no tempo n+1
# Q é a matriz do potencial no tempo n+1
Linha 184: Linha 208:
   t = t + td
   t = t + td


plt.plot(x,y,P)
plt.plot(x,y,P) # plotagem dos gráficos


</source>
</source>


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<ref name="metcompc1">Numerical solution of partial differential equations, Dr. Louise Olsen-Kettle. The University of Queensland School of Earth Sciences Centre for Geoscience Computing.</ref> <ref name="sayama">SAYAMA, Hiroki. Introduction to the modeling and analysis of complex systems. Open SUNY Textbooks, 2015.</ref>.
Lembrando que estamos resolvendo o problema em 2D, por isso '''P''' e '''Q''' são matrizes, onde cada elemento dessas matrizes 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<ref name="metcompc1">Numerical solution of partial differential equations, Dr. Louise Olsen-Kettle. The University of Queensland School of Earth Sciences Centre for Geoscience Computing.</ref> <ref name="sayama">SAYAMA, Hiroki. Introduction to the modeling and analysis of complex systems. Open SUNY Textbooks, 2015.</ref>.


== Resultados do problema da Borda Carregada ==
Para sabermos se um método converge para o estado estacionário mais rapidamente, analisamos a diferença entre o valor no tempo '''n''' e no tempo '''n+1''', da seguinte forma:
 
<center><math>diff = P^{n+1}_{i,j} - P^{n}_{i,j}</math>.</center>
 
Comparamos esta diferença '''''diff''''' para o mesmo número de iterações no tempo utilizando os três métodos (Jacobi, Gauss-Seidel e SOR), e assim analisamos quais métodos apresentavam valores menores em relação aos outros. No exemplo de implementação acima <math>Q_{i,j}=P^{n+1}_{i,j}</math>.
 
 
== Resultados e Discussões 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.
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 ===
=== Método de Jacobi ===
Linha 199: Linha 229:
Obtivemos a seguinte solução pelo método explícito:
Obtivemos a seguinte solução pelo método explícito:


[[Arquivo:Escorregador jacobi.png|center|thumb|500px|Solução numérica do problema 1.]]
[[Arquivo:Escorregador jacobi.png|center|thumb|500px|Solução numérica do problema da borda carregada.]]


Comparando com a iteração no passo de tempo anterior, obtivemos as seguintes diferenças entre o passo n+1 e o passo n:
Comparando com a iteração no passo de tempo anterior, obtivemos as seguintes diferenças entre o passo n+1 e o passo n:
Linha 205: Linha 235:
[[Arquivo:Jacobi 01.png|center|thumb|500px|Maior diferença percentual ao longo de todas iterações.]]
[[Arquivo:Jacobi 01.png|center|thumb|500px|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.
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.
 
Sobre o método de relaxação observa-se que nas primeiras iterações há uma variação muito grande dos valores da solução. Nas primeiras 60 iterações, os valores do potencial <math>\Phi</math> variaram em até 2500% em relação ao passo anterior. Percebe-se que a partir de 160 iterações começa a ter uma variação menor que 10%, e a partir de 1200 passos já atingimos uma variação menor que 0.1%. Decidimos comparar os métodos com 1500 passos, a partir disso já não notamos mais variação significativa.
 
[[Arquivo:Erro jacobi 01.png|center|thumb|600px|Erro relativo médio para a solução de Jacobi variando o y fixando x]]
 
O gráfico acima representa a comparação entre o método numérico e a solução analítica do problema. A maneira que escolhemos para verificar a acurácia do método foi a seguinte: fixamos o valor de X e comparamos todos os valores da solução numérica e analítica percorrendo em Y. Nessa comparação fizemos a média dos valores absolutos percentuais do erro entre as soluções, com isso temos a média do erro ao longo do eixo X, para diferentes iterações. Esta análise foi repetida para todas as implementações seguintes.
 
Fizemos a tentativa de ao invés de plotar o erro médio, considerar o erro máximo, porém analisamos que alguns valores máximos (perto das bordas) destoavam muito da média dos valores, não representando corretamente o comportamento que gostaríamos de apresentar nesta Wiki.
 
Na região que apresenta mais erro (mais distante da borda carregada), analisamos este comportamento da solução apresentar mais erro em relação à equação analítica, porém este comportamento não identificamos como sendo uma observação geral, acreditamos que seja devido à geometria do problema específico mesmo, pois como será apresentado nos resultados do dipolo elétrico, o maior erro foi onde tinha as cargas, e não longe delas, diferente do verificado neste problema.


<div><ul>
Abaixo temos a comparação entre a solução numérica e analítica:
<li style="display: inline-block; vertical-align: top;"> [[Arquivo:Jacobi 02.png|thumb|310px|Maior diferença percentual entre iteração temporal 0 e 100.]] </li>
<li style="display: inline-block; vertical-align: top;"> [[Arquivo:Jacobi 03.png|thumb|313px|Maior diferença percentual entre iteração temporal 100 e 400.]] </li>
<li style="display: inline-block; vertical-align: top;"> [[Arquivo:Jacobi 04.png||thumb|320px|Maior diferença percentual entre iteração temporal 400 e 1500 (estado final).]] </li>
</ul></div>


[[Arquivo:erro 3d jacobi.png|center|thumb|600px|Erro relativo entre a solução de Jacobi e a solução analítica para 1500 interações.]]


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 ===
=== Método de Gauss-Seidel ===
Linha 220: Linha 256:
Obtivemos a seguinte solução pelo método de Gauss-Seidel:
Obtivemos a seguinte solução pelo método de Gauss-Seidel:


[[Arquivo:Escorregador gauss.png|center|thumb|500px|Solução numérica do problema 1 utilizando método de Gauss-Seidel]]
[[Arquivo:Escorregador gauss.png|center|thumb|500px|Solução numérica do problema da borda carregada 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:
Comparando com a iteração no passo de tempo anterior, obtivemos as seguintes diferenças entre o passo n+1 e o passo n:


[[Arquivo:Gauss 01.png|center|thumb|500px|Maior diferença percentual ao longo de todas iterações]]
[[Arquivo:Gauss 01.png|center|thumb|500px|Maior diferença percentual ao longo de todas iterações.]]
 
A implementação numérica do Método de Gauss-Seidel atingiu o mesmo comportamento que a solução analítica, assim como apresentado pelo Método de Jacobi. Observando os gráficos da diferença dos valores de <math>\Phi</math> por iteração (mostrados na seção de '''Resultados''''), pode-se notar que a variação menor que 10% encontra-se a partir da 90ª iteração (já no método de Jacobi observou-se esta variação somente pelo 160º passo). Já esperávamos que o método de Gauss-Seidel convergisse mais rápido, pelo fato da equação de iteração utilizar dois termos calculados no passo '''n''' e dois termos do passo '''n-1''', enquanto o método de jacobi considera os quatro termos do passo '''n-1'''.


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
<center><div><ul>
<li style="display: inline-block; vertical-align: top;"> [[Arquivo:Erro gauss 01.png|thumb|center|500px|Erro relativo médio para a solução de Gauss-Seidel para várias iterações.]] </li>
</ul></div></center>
 
Os gráficos acima representam o mesmo cálculo de erro médio percentual absoluto utilizado para comparar o Método de Jacobi. Porém para este método ressaltamos o resultado das soluções acima de 1000 passos (gráfico da direita). Verificamos que para o número de iterações menores o erro foi ligeiramente menor, o que não era esperado.
 
Abaixo temos a comparação entre a solução numérica e analítica:
 
[[Arquivo:erro 3d gauss.png|center|thumb|600px|Erro relativo entre a solução de Gauss-Seidel e a solução analítica para 1500 interações.]]


<div><ul>
<li style="display: inline-block; vertical-align: top;"> [[Arquivo:Gauss 02.png|thumb|310px|Maior diferença percentual entre iteração temporal 0 e 100.]] </li>
<li style="display: inline-block; vertical-align: top;"> [[Arquivo:Gauss 03.png|thumb|313px|Maior diferença percentual entre iteração temporal 100 e 400.]] </li>
<li style="display: inline-block; vertical-align: top;"> [[Arquivo:Gauss 04.png|thumb|320px|Maior diferença percentual entre iteração temporal 400 e 1500 (estado final).]] </li>
</ul></div>


=== Método SOR (Simultaneous Overrelaxation) ===
=== Método SOR (Simultaneous Overrelaxation) ===
Linha 238: Linha 279:
Obtivemos a seguinte solução aprimorando o método de Gauss-Seidel para ''Simultaneous OverRelaxation'':
Obtivemos a seguinte solução aprimorando o método de Gauss-Seidel para ''Simultaneous OverRelaxation'':


[[Arquivo:Escorregador sor.png|center|thumb|500px|Solução numérica do problema 1 utilizando método SOR]]
[[Arquivo:Escorregador sor.png|center|thumb|500px|Solução numérica do problema da borda carregada utilizando método SOR]]


Comparando com a iteração no passo anterior, podemos fazer a seguinte análise:
Comparando com a iteração no passo anterior, podemos fazer a seguinte análise:
Linha 244: Linha 285:
[[Arquivo:Sor 01.png|center|thumb|500px|Maior diferença percentual ao longo de todas iterações]]
[[Arquivo:Sor 01.png|center|thumb|500px|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).  
A implementação numérica do Método de Superrelaxação atingiu o mesmo comportamento que a solução analítica. Este método obteve a maior variação, chegando à 25000% nas primeiras 10 iterações. Porém, a partir do 90º passo a variação já era menor que 0.1%, o que demonstra ser o método mais rápido para convergir à uma solução estacionária. Também verificamos que a partir de 400 iterações o método não apresenta mais diferença entre valores do passo '''n+1''' e o passo '''n'''.
 
Um fenômeno não esperado para este método que pudemos observar foram picos de variação nas iterações múltiplas de 50, que é valor do tamanho do nosso plano XY.
 
No gráfico abaixo observamos o fato do erro apresentar valores menores para um número de iterações menores, porém para o método de SOR aconteceu um valor de passos muito menor do que para o método de Gauss-Seidel, abaixo de 50 iterações já notamos que convergimos à solução da Equação de Laplace.
 
[[Arquivo:Erro sor 01.png|center|thumb|500px|Erro relativo médio para a solução de SOR variando o y fixando x.]]
 
Abaixo temos a comparação entre a solução numérica e analítica:


<div><ul>
[[Arquivo:erro 3d sor.png|center|thumb|600px|Erro relativo entre a solução de SOR e a solução analítica para 1500 interações.]]
<li style="display: inline-block; vertical-align: top;"> [[Arquivo:Sor 02.png|thumb|230px|Maior diferença percentual entre iteração temporal 0 e 60.]] </li>
<li style="display: inline-block; vertical-align: top;"> [[Arquivo:Sor 03.png|thumb|230px|Maior diferença percentual entre iteração temporal 30 e 80.]] </li>
<li style="display: inline-block; vertical-align: top;"> [[Arquivo:Sor 04.png|thumb|230px|Maior diferença percentual entre iteração temporal 70 e 100.]] </li>
<li style="display: inline-block; vertical-align: top;"> [[Arquivo:Sor 05.png|thumb|230px|Maior diferença percentual entre iteração temporal 90 e 160.]] </li>
</ul></div>


== Resultados do problema do Dipolo Elétrico ==
== Resultados e Discussões 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.
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.
Uma observação importante é que é esperado que seja infinito o valor do potencial eletríco exatamente nos pontos onde fica as cargas, porém, para os gráficos desse trabalho, colocamos um valor médio nesses pontos somento para representação. O que analisamos para os erros foi a partir dos valores vizinhos.


=== Método de Gauss-Seidel ===
=== Método de Gauss-Seidel ===


[[Arquivo:Sol-gauss-poisson.jpeg|center|thumb|500px|legendas eq poisson - gauss]]
A implementação foi realizada conforme descrito na seção de '''Métodos Numéricos'''.
[[Arquivo:Sol-SOR-poisson.jpeg|center|thumb|500px|legendas eq poisson - -SOR]]
 
[[Arquivo:Erro-poisson1.jpeg|center|thumb|500px|legendas eq poisson - comparação das soluções]]
<center><div><ul>
[[Arquivo:Erro-relativo-poisson.jpeg|center|thumb|500px|legendas eq poisson - erro médio com analitica]]
<li style="display: inline-block; vertical-align: top;"> [[Arquivo:Dipolo gauss 01.png|thumb|500px|Solução para o potencial elétrico do dipolo com o método de Gauss-Seidel.]] </li>
<li style="display: inline-block; vertical-align: top;"> [[Arquivo:Dipolo gauss 02.png|thumb|500px|Solução para o potencial elétrico do dipolo com o método de Gauss-Seidel.]] </li>
</ul></div></center>
 
=== Método de SOR ===
 
Assim como a implementação do método de Gauss-Seidel, obtivemos os seguintes resultados com o método de SOR.
 
<center><div><ul>
<li style="display: inline-block; vertical-align: top;"> [[Arquivo:Dipolo sor 01.png|thumb|500px|Solução para o potencial elétrico do dipolo com o método de superrelaxamento.]] </li>
<li style="display: inline-block; vertical-align: top;"> [[Arquivo:Dipolo sor 02.png|thumb|500px|Solução para o potencial elétrico do dipolo com o método de superrelaxamento.]] </li>
</ul></div></center>
 
Observando os gráficos acima é possível concluir que os métodos implementados atingiram o mesmo comportamento da solução analítica.
 
[[Arquivo:Erro todos corte.png|center|thumb|500px|Comparação entre os métodos numéricos e a solução analítica.]]
 
Porém, também observamos que apesar de alcançar a solução estacionária, a equação de iteração não consegue alcançar com precisão a equação analítica. A divergência da solução analítica é representada corretamente (a plotagem no gráfico na verdade representa a divergência daquela maneira). Porém para os métodos analíticos precisamos selecionar um valor de carga pontual como densidade de carga, por isso não é tratada como uma divergência.
 
== Considerações Finais ==
 
=== Problema da Borda Carregada ===
 
Primeiro começamos a implementação pelo Método de Jacobi, modelando um problema conhecido para compararmos a precisão de nossos algoritmos com a solução analítica. Ao resolver o potencial elétrico nas condições do problema da borda carregada obtivemos o comportamento da curva igual ao esperado analiticamente, porém com certo erro, conforme foi apresentado anteriormente.
 
As discussões de cada implementação foram descritas em suas devidas seções em '''Resultados'''.
 
Por fim, analisamos a precisão de cada método e é possível observar que o Método de Gauss-Seidel apresentou os menores erros, porém necessita mais iterações comparado ao método SOR, que tem um erro ligeiramente maior, porém é cerca de 20 vezes mais rápido.
 
[[Arquivo:Erro todos 01.png|center|thumb|500px|Erro médio por método, ao longo do eixo X. Todos os três métodos foram comparados com 1500 iterações.]]
 
Uma observação importante é que por mais que para interações menores alguns métodos ficaram com um erro menor comparado a analítica, decidimos comparar os métodos com uma interação grande (1500) pois quando estes forem usados para resolver um problema sem solução, não vamos ter essa possibilidade de comparar com a analítica para saber qual é interação ideal para parar a simulação. Assim temos uma noção geral de quão próximo da solução real está cada método.
 
=== Problema do Dipolo Elétrico ===


== Discussão ==
Analisamos a precisão de cada método implementado, não foi observado diferenças significativas entre os métodos. Observa-se que para este problema tivemos mais erros onde se localizam as cargas, diferente dos maiores erros encontrados no problema da borda carregada.


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.
[[Arquivo:Erro todos.png|center|thumb|500px|Erro médio relativo dos dois métodos implementados ao longo do eixo Y.]]


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.
Por fim, um ponto que gostaríamos de ressaltar é a facilidade de modelar um problema complexo de eletroestática com estas implementações. Muitas vezes é difícil de chegar à uma equação analítica para uma densidade de carga complexa, porém em código é possível com algumas restrições modelar problemas com geometrias bem específicas.


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
== Link para Códigos ==


[[Arquivo:Discuss1.jpeg|center|thumb|500px|Erro médio por método, ao longo do eixo x]]
Fizemos no ambiente Colab em ''.ipynb'', segue link do github:[https://github.com/padovanih/equacao-de-laplace]  


== Referências ==  
== Referências ==  


<references/>
<references/>

Edição atual tal como às 13h31min de 6 de abril de 2021

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 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 () de uma determinada região num espaço que não possui nenhuma densidade de carga elétrica (corpo carregado):

.

ou na sua versão em 2 dimensões[2]:

.

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:

.

ou na sua versão em 2 dimensões[2]:

.

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

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

.

.

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

.

Com isso: , 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:

.

Discretizando, primeiro chegamos que:

.

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

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

.

No nosso algoritmo ultizamos . Então obtivemos a condição de estabilidade:

.

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

.

Equação de Poisson

Partindo de:

.

chegamos em:

.

Para nosso problema , então multiplicando os dois lados por , chegamos em:

.

E finalmente, aplicando a condição de estabilidade e cancelando os termos :

.

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, percorrendo os índices de forma crescente, obtemos a seguinte equação de iteração:

.

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 ao lado direito da equação. Assim, temos a seguinte equação de iteração:

.

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

.

Equação de Poisson

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

.

Onde:

.

E é dado por:

.

No nosso caso, temos que , portanto o é:

.

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

Condições iniciais e de contorno:

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


  • Problema da borda carregada eletricamente.
  • Gráfico da solução analítica somando até o termo n=199.

Problema do Dipolo Elétrico

Equação de Poisson considerando uma dipolo elétrico como densidade de carga :

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

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

Onde:

Onde e são as coordedanas da carga. Também consideramos as constantes com sendo = 1 para facilitar a modelagem do problema.

  • Dipolo Elétrico[3].
  • 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 do 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 FTCS integra no espaço (2D) a equação de iteração, mas utilizando métodos de Relaxação é preciso também iterar no tempo estas equações. Um exemplo disso é mostrado abaixo:

### Exemplo da evolução temporal no método de relaxação ###
### Exemplo para o algoritmo de jacobi, Equação de Laplace ###
# 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) # plotagem dos gráficos

Lembrando que estamos resolvendo o problema em 2D, por isso P e Q são matrizes, onde cada elemento dessas matrizes 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[4] [5].

Para sabermos se um método converge para o estado estacionário mais rapidamente, analisamos a diferença entre o valor no tempo n e no tempo n+1, da seguinte forma:

.

Comparamos esta diferença diff para o mesmo número de iterações no tempo utilizando os três métodos (Jacobi, Gauss-Seidel e SOR), e assim analisamos quais métodos apresentavam valores menores em relação aos outros. No exemplo de implementação acima .


Resultados e Discussões 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 da borda carregada.

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.

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.

Sobre o método de relaxação observa-se que nas primeiras iterações há uma variação muito grande dos valores da solução. Nas primeiras 60 iterações, os valores do potencial variaram em até 2500% em relação ao passo anterior. Percebe-se que a partir de 160 iterações começa a ter uma variação menor que 10%, e a partir de 1200 passos já atingimos uma variação menor que 0.1%. Decidimos comparar os métodos com 1500 passos, a partir disso já não notamos mais variação significativa.

Erro relativo médio para a solução de Jacobi variando o y fixando x

O gráfico acima representa a comparação entre o método numérico e a solução analítica do problema. A maneira que escolhemos para verificar a acurácia do método foi a seguinte: fixamos o valor de X e comparamos todos os valores da solução numérica e analítica percorrendo em Y. Nessa comparação fizemos a média dos valores absolutos percentuais do erro entre as soluções, com isso temos a média do erro ao longo do eixo X, para diferentes iterações. Esta análise foi repetida para todas as implementações seguintes.

Fizemos a tentativa de ao invés de plotar o erro médio, considerar o erro máximo, porém analisamos que alguns valores máximos (perto das bordas) destoavam muito da média dos valores, não representando corretamente o comportamento que gostaríamos de apresentar nesta Wiki.

Na região que apresenta mais erro (mais distante da borda carregada), analisamos este comportamento da solução apresentar mais erro em relação à equação analítica, porém este comportamento não identificamos como sendo uma observação geral, acreditamos que seja devido à geometria do problema específico mesmo, pois como será apresentado nos resultados do dipolo elétrico, o maior erro foi onde tinha as cargas, e não longe delas, diferente do verificado neste problema.

Abaixo temos a comparação entre a solução numérica e analítica:

Erro relativo entre a solução de Jacobi e a solução analítica para 1500 interações.


Método de Gauss-Seidel

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

Solução numérica do problema da borda carregada 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.

A implementação numérica do Método de Gauss-Seidel atingiu o mesmo comportamento que a solução analítica, assim como apresentado pelo Método de Jacobi. Observando os gráficos da diferença dos valores de por iteração (mostrados na seção de Resultados'), pode-se notar que a variação menor que 10% encontra-se a partir da 90ª iteração (já no método de Jacobi observou-se esta variação somente pelo 160º passo). Já esperávamos que o método de Gauss-Seidel convergisse mais rápido, pelo fato da equação de iteração utilizar dois termos calculados no passo n e dois termos do passo n-1, enquanto o método de jacobi considera os quatro termos do passo n-1.

  • Erro relativo médio para a solução de Gauss-Seidel para várias iterações.

Os gráficos acima representam o mesmo cálculo de erro médio percentual absoluto utilizado para comparar o Método de Jacobi. Porém para este método ressaltamos o resultado das soluções acima de 1000 passos (gráfico da direita). Verificamos que para o número de iterações menores o erro foi ligeiramente menor, o que não era esperado.

Abaixo temos a comparação entre a solução numérica e analítica:

Erro relativo entre a solução de Gauss-Seidel e a solução analítica para 1500 interações.


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 da borda carregada 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

A implementação numérica do Método de Superrelaxação atingiu o mesmo comportamento que a solução analítica. Este método obteve a maior variação, chegando à 25000% nas primeiras 10 iterações. Porém, a partir do 90º passo a variação já era menor que 0.1%, o que demonstra ser o método mais rápido para convergir à uma solução estacionária. Também verificamos que a partir de 400 iterações o método não apresenta mais diferença entre valores do passo n+1 e o passo n.

Um fenômeno não esperado para este método que pudemos observar foram picos de variação nas iterações múltiplas de 50, que é valor do tamanho do nosso plano XY.

No gráfico abaixo observamos o fato do erro apresentar valores menores para um número de iterações menores, porém para o método de SOR aconteceu um valor de passos muito menor do que para o método de Gauss-Seidel, abaixo de 50 iterações já notamos que convergimos à solução da Equação de Laplace.

Erro relativo médio para a solução de SOR variando o y fixando x.

Abaixo temos a comparação entre a solução numérica e analítica:

Erro relativo entre a solução de SOR e a solução analítica para 1500 interações.

Resultados e Discussões 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.

Uma observação importante é que é esperado que seja infinito o valor do potencial eletríco exatamente nos pontos onde fica as cargas, porém, para os gráficos desse trabalho, colocamos um valor médio nesses pontos somento para representação. O que analisamos para os erros foi a partir dos valores vizinhos.

Método de Gauss-Seidel

A implementação foi realizada conforme descrito na seção de Métodos Numéricos.

  • Solução para o potencial elétrico do dipolo com o método de Gauss-Seidel.
  • Solução para o potencial elétrico do dipolo com o método de Gauss-Seidel.

Método de SOR

Assim como a implementação do método de Gauss-Seidel, obtivemos os seguintes resultados com o método de SOR.

  • Solução para o potencial elétrico do dipolo com o método de superrelaxamento.
  • Solução para o potencial elétrico do dipolo com o método de superrelaxamento.

Observando os gráficos acima é possível concluir que os métodos implementados atingiram o mesmo comportamento da solução analítica.

Comparação entre os métodos numéricos e a solução analítica.

Porém, também observamos que apesar de alcançar a solução estacionária, a equação de iteração não consegue alcançar com precisão a equação analítica. A divergência da solução analítica é representada corretamente (a plotagem no gráfico na verdade representa a divergência daquela maneira). Porém para os métodos analíticos precisamos selecionar um valor de carga pontual como densidade de carga, por isso não é tratada como uma divergência.

Considerações Finais

Problema da Borda Carregada

Primeiro começamos a implementação pelo Método de Jacobi, modelando um problema conhecido para compararmos a precisão de nossos algoritmos com a solução analítica. Ao resolver o potencial elétrico nas condições do problema da borda carregada obtivemos o comportamento da curva igual ao esperado analiticamente, porém com certo erro, conforme foi apresentado anteriormente.

As discussões de cada implementação foram descritas em suas devidas seções em Resultados.

Por fim, analisamos a precisão de cada método e é possível observar que o Método de Gauss-Seidel apresentou os menores erros, porém necessita mais iterações comparado ao método SOR, que tem um erro ligeiramente maior, porém é cerca de 20 vezes mais rápido.

Erro médio por método, ao longo do eixo X. Todos os três métodos foram comparados com 1500 iterações.

Uma observação importante é que por mais que para interações menores alguns métodos ficaram com um erro menor comparado a analítica, decidimos comparar os métodos com uma interação grande (1500) pois quando estes forem usados para resolver um problema sem solução, não vamos ter essa possibilidade de comparar com a analítica para saber qual é interação ideal para parar a simulação. Assim temos uma noção geral de quão próximo da solução real está cada método.

Problema do Dipolo Elétrico

Analisamos a precisão de cada método implementado, não foi observado diferenças significativas entre os métodos. Observa-se que para este problema tivemos mais erros onde se localizam as cargas, diferente dos maiores erros encontrados no problema da borda carregada.

Erro médio relativo dos dois métodos implementados ao longo do eixo Y.

Por fim, um ponto que gostaríamos de ressaltar é a facilidade de modelar um problema complexo de eletroestática com estas implementações. Muitas vezes é difícil de chegar à uma equação analítica para uma densidade de carga complexa, porém em código é possível com algumas restrições modelar problemas com geometrias bem específicas.

Link para Códigos

Fizemos no ambiente Colab em .ipynb, segue link do github:[1]

Referências

  1. GARCIA, Alejandro L. Numerical methods for physics. Englewood Cliffs, NJ: Prentice Hall, 2000.
  2. 2,0 2,1 GRIFFITHS, David J. Introduction to electrodynamics. New Jersey: Prentice Hall, 1962.
  3. https://noic.com.br/materiais-fisica/cursos/teorico/aula-5-3-energia-potencial-eletrica-trabalho-e-potencial-eletrico/
  4. Numerical solution of partial differential equations, Dr. Louise Olsen-Kettle. The University of Queensland School of Earth Sciences Centre for Geoscience Computing.
  5. SAYAMA, Hiroki. Introduction to the modeling and analysis of complex systems. Open SUNY Textbooks, 2015.