Grupo1 - Dif em 2D: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
Sem resumo de edição
 
(31 revisões intermediárias pelo mesmo usuário não estão sendo mostradas)
Linha 3: Linha 3:
<math>\Delta u = \frac{\partial ^{2} u }{\partial x^2} + \frac{\partial ^{2} u }{\partial y^2} = -g(x,y)</math>
<math>\Delta u = \frac{\partial ^{2} u }{\partial x^2} + \frac{\partial ^{2} u }{\partial y^2} = -g(x,y)</math>


é uma equação do tipo Elíptica que representa fenômenos físicos estácionarios relacionados a Eletrostatica, Dinâmica de Fluídos e Transferência de Calor, se <math> g(x,y) \equiv 0 </math> a equação passa a ser chamada de Equação de Laplace. Os problemas relacionados a equação de Laplace são estudados pela "Teoria do Potencial".  
é uma equação do tipo Elíptica que representa fenômenos físicos estácionarios relacionados a Eletrostatica, Dinâmica de Fluídos e Transferência de Calor. Se <math> g(x,y) \equiv 0 </math> a equação passa a ser chamada de Equação de Laplace. Os problemas relacionados a equação de Laplace são estudados pela "Teoria do Potencial".  


As soluções <math> u = u(x,y) </math> da Equação de Laplace são denominadas funções Harmônicas. Como já era de se esperar, os problemas mais habituais na vida de um físico, engenheiro ou matemático ao se depararem com uma EDP,  são os problemas com Condições de Contorno, essencialmente será trabalhada a Condição de Dirichlet, que possui fronteiras conhecidas, tendo o seguinte formato:
As soluções <math> u = u(x,y) </math> da Equação de Laplace são denominadas funções Harmônicas. Os problemas mais habituais na vida de um físico, engenheiro ou matemático ao se depararem com uma EDP,  são os problemas com Condições de Contorno em um dominío <math> \Omega \in \mathbb{R}^2</math>, essencialmente será trabalhada a Condição de Dirichlet, que possui fronteiras (<math> \partial \Omega </math>) conhecidas, tendo o seguinte formato:


<math> \begin{cases}
<math> \begin{cases}
Linha 22: Linha 22:




Para tais problemas, estudaremos os métodos de Relaxação e Super-Relaxação para encontrar as soluções da Equação de Laplace no Quadrado de Lado <math>L </math>  <math> ( \{ \Omega = (0,L)\times (0,L) \} ) </math> e faremos (se possível) a solução númerica via Formula de Poisson para a Bola Unitária centrada na origem <math> \{ \Omega = \mathcal{B}_{1}(0) \} </math>.
Para tais problemas, estudaremos os métodos de Relaxação e Super-Relaxação para encontrar as soluções da Equação de Laplace na região de Quadrado de Lado <math>L </math>  <math> ( \{ \Omega = (0,L)\times (0,L) \} ) </math>.
 
== Dominio Quadrado de Lado <math> L </math> ==
== Quadrado de Lado <math> L </math> ==


=== Solução Analítica da Equação de Laplace  ===
=== Solução Analítica da Equação de Laplace  ===
Linha 49: Linha 48:
</math>
</math>


Separamos o problema geral de Dirchlet em 4 problemas "menores" tal que obtemos os problemas desde:
Separamos o problema geral de Dirchlet em 4 problemas "menores", com condições de contorno diferentes de zero em apenas um trecho da fronteira, de modo que obtemos desde:


<math>
<math>
Linha 72: Linha 71:




Podemos então utilizar o Método da Separação de Variáveis para resolver os 4 "probleminhas" e, como a Equação de Laplace é linear, sua soma será a solução completa do Problema de Dirichlet. O método consiste em supor <math> u(x, y) = \phi(x)\theta(y)</math>, para então, ao substituirmos na equação obtermos a seguinte expressão:
Podemos então utilizar o Método da Separação de Variáveis para resolver os 4 problemas e, como a Equação de Laplace é linear, sua soma será a solução completa do Problema de Dirichlet. O método consiste em supor <math> u(x, y) = \phi(x)\theta(y)</math>, para então, ao substituirmos na equação obtermos a seguinte expressão:


<math> \Delta u_{i} = \ddot{\phi_{i}}(x)\theta_{i}(y) + \phi_{i}(x) \ddot{\theta_{i}}(y) = 0 </math>
<math> \Delta u_{i} = \ddot{\phi_{i}}(x)\theta_{i}(y) + \phi_{i}(x) \ddot{\theta_{i}}(y) = 0 </math>


Podemos isolar as funções <math>\phi_{i}</math> e <math>\theta_{i}</math>, de fato ficamos com com duas relações que dependem de suas variações, portanto para elas serem sempre iguais, é necessário que sejam constantes (<math> = \lambda</math>):
Podemos isolar as funções <math>\phi_{i}</math> e <math>\theta_{i}</math>, de fato ficamos com com duas relações que dependem ou apenas de <math> x </math> e <math> y </math> portanto para elas serem sempre iguais, é necessário que sejam constantes (<math> = \lambda</math>):


<math> \frac{\ddot{\phi_{i}}}{\phi_{i}} = -\frac{\ddot{\theta_{i}}}{\theta_{i}} = \lambda </math>
<math> \frac{\ddot{\phi_{i}}}{\phi_{i}} = -\frac{\ddot{\theta_{i}}}{\theta_{i}} = \lambda </math>
Linha 89: Linha 88:
<math>
<math>
\begin{cases}
\begin{cases}
\ddot{\phi_{1}}(x) - \lambda \phi_{1} =0; \\
\ddot{\phi_{1}}(y) - \lambda \phi_{1} =0; \\
\phi_{1} (L) = 0; \\
\phi_{1} (L) = 0; \\
\end{cases}
\end{cases}
Linha 103: Linha 102:


<math> \phi_{1}(x) = A_{1}e^{\sqrt{\lambda} x} + A_{2}e^{- \sqrt{\lambda} x} </math>
<math> \phi_{1}(x) = A_{1}e^{\sqrt{\lambda} x} + A_{2}e^{- \sqrt{\lambda} x} </math>
Utilizando a C.C. vemos que <math> -A_{1}/A_{2} = tanh(\sqrt{\lambda}L)</math>.


Partindo para a segunda equação <math> \theta_{1}(y) </math>,
Partindo para a segunda equação <math> \theta_{1}(y) </math>,
Linha 125: Linha 122:
Ou seja, temos solução <math> \theta_{1}(y) = B_{1}e^{i\sqrt{\lambda} y} + B_{2}e^{-i\sqrt{\lambda} y} </math>  
Ou seja, temos solução <math> \theta_{1}(y) = B_{1}e^{i\sqrt{\lambda} y} + B_{2}e^{-i\sqrt{\lambda} y} </math>  


Utilizando a primeira C.C. obtemos <math>0 = B_{1} + B_{2}</math>
Utilizando a primeira C.C. obtemos <math>B_{1} = - B_{2} = B</math>
 
ou seja, temos que  
ou seja, temos que  


<math>\theta_{1}(y) = B sen(\sqrt{\lambda} y ). </math> Utilizando a segunda C.C. temos <math> 0 = sin(\sqrt{\lambda} y)  \Rightarrow \lambda = \frac{n^2 \pi^2}{L^2} </math>, ou seja, existem infinitos <math>n</math> tal que <math>\theta_{1}</math> é solução.  
<math>\theta_{1}(y) = B sen(\sqrt{\lambda} y ). </math> Utilizando a segunda C.C. temos <math> 0 = sen(\sqrt{\lambda} y)  \Rightarrow \lambda = \frac{n^2 \pi^2}{L^2}, </math>
 
ou seja, existem infinitos <math>n</math> tal que <math>\theta_{1}</math> é solução.  
 
Voltando a <math>\phi_{1}</math>, temos <math> \phi_{1}(x) = senh(\frac{n\pi (L-x)}{L}). </math>


Voltando a <math>\phi_{1}</math>, temos <math> \phi_{1}(x) = senh(\frac{n\pi (L-x)}{L}) </math>
Finalmente unindo as respostas, temos  
Finalmente unindo as respostas, temos  


Linha 147: Linha 148:
A solução completa do problema de Dirichlet no quadrado de Lado <math>L</math> é a soma das quatro soluções parciais: <math> u(x, y) = u_{1}(x, y) + u_{2}(x, y) + u_{3}(x, y) + u_{4}(x, y) </math>.
A solução completa do problema de Dirichlet no quadrado de Lado <math>L</math> é a soma das quatro soluções parciais: <math> u(x, y) = u_{1}(x, y) + u_{2}(x, y) + u_{3}(x, y) + u_{4}(x, y) </math>.


=== Algoritmo de Relaxação ===
=== Método da Relaxação ===
 
O Método da Relaxação é um método iterativo utilizado para obter a solução numérica para a equação de Laplace e Poisson. A ideia do método é de, utilizando a vizinhança iterar os pontos da malha até que convirjam para uma solução.


Discretizando a equação temos <math>x\mapsto x_{i}</math> e <math>y\mapsto y_{j}</math> para <math>i</math>, <math>j = 1, ..., N</math> e <math> h = \Delta x = \Delta y = N/L</math> a função <math> g(x, y) = g_{ij}</math>, nos deparamos com uma matriz <math>\mathcal{M}_{i j} = u(x_{i},y_{j}) = u_{i j}</math> quadrada sendo as bordas <math>\mathcal{M}_{1 j} = f(0, y) = u_{1 j}</math>, <math>\mathcal{M}_{N j} = f(L, y) = u_{N j}</math>, <math>\mathcal{M}_{i 1} = f(x, 0) = u_{i 1}</math> e <math>\mathcal{M}_{i N} = f(x, L) = u_{i N}</math>.
Discretizando a equação temos <math>x\mapsto x_{i}</math> e <math>y\mapsto y_{j}</math> para <math>i</math>, <math>j = 1, ..., N</math> e <math> h = \Delta x = \Delta y = N/L</math> a função <math> g(x, y) = g_{ij}</math>, nos deparamos com uma matriz <math>\mathcal{M}_{i j} = u(x_{i},y_{j}) = u_{i j}</math> quadrada sendo as bordas <math>\mathcal{M}_{1 j} = f(0, y) = u_{1 j}</math>, <math>\mathcal{M}_{N j} = f(L, y) = u_{N j}</math>, <math>\mathcal{M}_{i 1} = f(x, 0) = u_{i 1}</math> e <math>\mathcal{M}_{i N} = f(x, L) = u_{i N}</math>.
Linha 159: Linha 162:
Substituindo na Equação, temos
Substituindo na Equação, temos


<math> \frac{u_{(i+1) j} + u_{(i-1) j} -2u_{i j}}{h^2}  = - \frac{u_{i (j+1)} + u_{i (j-1)} -2u_{i j}}{h^2} + g_{ij} </math>,  
<math> \frac{u_{(i+1) j} + u_{(i-1) j} -2u_{i j}}{h^2}  + \frac{u_{i (j+1)} + u_{i (j-1)} -2u_{i j}}{h^2} = - g_{ij} </math>,  
ou seja:
ou seja:


Linha 168: Linha 171:
<math> u_{i j} = \frac{(\Delta y)^{2}(u_{(i-1) j} + u_{(i+1) j}) + (\Delta x)^2(u_{i (j-1)} + u_{i (j+1)}) + (\Delta y \Delta x)^2 g_{ij}}{2((\Delta x)^2 + (\Delta y)^2)}, </math>
<math> u_{i j} = \frac{(\Delta y)^{2}(u_{(i-1) j} + u_{(i+1) j}) + (\Delta x)^2(u_{i (j-1)} + u_{i (j+1)}) + (\Delta y \Delta x)^2 g_{ij}}{2((\Delta x)^2 + (\Delta y)^2)}, </math>


para <math> i, j = 2, ..., N-1 </math>
para <math> i, j = 2, ..., N+1 </math>


Como condição de parada, foi convencionado tomar o Erro Relativo entre as iterações <math> k </math> e <math> k+1 </math>, para estimar o erro se faz:
Para condição de parada, foi convencionado tomar o erro relativo entre as iterações <math> k </math> e <math> k+1 </math>, para estimar o erro, se optou por tomar como valores, o ponto médio da malha (já que é o ultimo ponto a ser alcançado nas iterações, portanto, quando sua variação diminuir é sinal de que a solução já está convergindo), para observarmos a evolução em relação a outros pontos que variam desde o inicio, foram utilizados as diagonais interiores, tal que o erro relativo é:


<math> \epsilon = \Big|\frac{v^{k} - v^{k+1}}{v^{k}}\Big|</math>
<math> \epsilon = \Big|\frac{v^{k} - v^{k+1}}{v^{k}}\Big|</math>


<math> v^{k} = \frac{1}{8} (u_{2 2} + u_{2 (M+1)} + u_{(M+1) 2} + 4u_{\frac{(M+2)}{2}  \frac{(M+2)}{2}}+u_{(M+1) (M+1)}) </math>
fazendo a média ponderada com peso 4 para o ponto médio:


<math> v^{k} = \frac{1}{8} (u_{2 2} + u_{2 (N+1)} + u_{(N+1) 2} + 4u_{\frac{(N+2)}{2}  \frac{(N+2)}{2}}+u_{(N+1) (N+1)}) </math>
=== Método da Super Relaxação (SOR) ===
Podemos também realizar uma média entre os valores já calculados e os ainda não calculados na iteração, o método da Super Relaxação ou Sobrerrelaxação (SOR) é da seguinte forma:
<math> u_{ij}^{k+1} = u_{i j}^{k}(1 - \omega) + \omega u_{ij}^{R} </math>
tal que <math> u_{ij}^{R} </math> é o valor calculado através do método da Relaxação e <math> \omega \in [0:2]</math> é o fator de relaxamento, se <math> \omega = 1 </math> temos a Relaxação normal.


=== Estabilidade ===
=== Estabilidade ===


[[Arquivo:erro2.png|thumb|400px| Erro relativo ao método de Relaxação em função da quantidade de iterações]]
==== Relaxação ====
[[Arquivo:erro2.png|thumb|400px| Erro relativo (método da Relaxação) em função da quantidade de iterações]]


A relaxação é um método Iterativo sobre os pontos vizinhos que pode ser feita de 2 modos, pelo Algoritmo de Jacobi, e pelo de Gauss-Seidel.  
A relaxação é um método Iterativo sobre os pontos vizinhos que pode ser feita de 2 modos, pelo Algoritmo de Jacobi, e pelo de Gauss-Seidel.  


O algoritmo de Jacobi pega valores "antigos" para a iteração e possui convergencia muito lenta, por isso não é muito utilizado. Já o algoritmo de Gauss-Seidel pega os valores "novos" (que ja foram calculados) e os "antigos" (que não foram calculados), possui convergencia mais rapida, porém ainda é lenta.  
O algoritmo de Jacobi pega valores "antigos" para a iteração e possui convergencia muito lenta, por isso não é muito utilizado. Já o algoritmo de Gauss-Seidel pega os valores "novos" (que ja foram calculados na iteração) e os "antigos" (que não foram calculados na iteração <math> k </math>), possui convergência mais rapida, porém ainda é lenta.  
 
Como a média definida anteriormente foi feita utilizando o ponto médio do domínio, o erro cresce após decair, pois é quando efetivamente ocorrem variações maiores no ponto médio.
 
Para a relaxação, o algoritmo de Jacobi faz o seguinte cálculo:
 
<math> u_{i j}^{k+1} = \frac{u_{(i-1) j}^{k} + u_{(i+1) j}^{k} + u_{i (j-1)}^{k} + u_{i (j+1)}^{k} + h^{2}g_{ij}}{4}</math>,
 
já o algoritmo de Gauss-Seidel faz:
 
<math> u_{i j}^{k+1} = \frac{u_{(i-1) j}^{k+1} + u_{(i+1) j}^{k} + u_{i (j-1)}^{k+1} + u_{i (j+1)}^{k} + h^{2}g_{ij}}{4}</math>,


Algoritmos iterativos tendem a convergir para solução unica, se a matriz que as representa for Diagonal Dominante, ou seja:
Algoritmos iterativos tendem a convergir para solução unica, se a matriz que as representa for Diagonal Dominante, ou seja:


<math> |a_{ii}| \ge \sum_{j=1 (i\neq j)}^{N} |a_{ij}| </math>
<math> |a_{ii}| \ge \sum_{j=1 (i\neq j)}^{N} |a_{ij}| </math>
Sendo <math> a_{ii} = u_{ii} </math> e  <math> a_{ij} = u_{ij} </math>, então  <math> |u_{ii}| = 4 </math> e <math> |u_{ij}| = 1 </math> e
<math>  \sum_{j=1 (i\neq j)}^{N+2} |a_{ij}| = 4 </math>
Como <math> 4 = 4 </math>, a desigualdade vale e o método converge.


De fato, podemos ver que a equação de Laplace respeita tal desigualdade.
De fato, podemos ver que a equação de Laplace respeita tal desigualdade.




Caso façamos um retangulo <math> L_{x} \neq L_{y} </math> com <math> \Delta x \neq \Delta y </math>, obtemos o erro da imagem a seguir:
Caso façamos um retangulo <math> L_{x} \neq L_{y} </math> com <math> \Delta x \neq \Delta y </math>, obtemos os erros da imagem a seguir, feito utilizando <math> L_{x} =1 </math> e <math> L_{y} = 5 </math>, <math> L_{x} =2 </math> e <math> L_{y} = 4 </math> e <math> L_{x} =1.5 </math> e <math> L_{y} = 4.5 </math>:


[[Arquivo:erro.png|frameless|400px|Erro obtido para diferentes comprimentos do retangulo]]
[[Arquivo:erro.png|frameless|400px|Erro obtido para diferentes comprimentos do retangulo]]


=== Método da Super Relaxação ===
==== Super Relaxação ====


[[Arquivo:errooverr.png|thumb|400px|Grafico mostra a quantidades de Iterações para convergência do exemplo 1 em função de omega]]
[[Arquivo:errooverr.png|thumb|400px|Grafico mostra a quantidades de Iterações para convergência do exemplo 1 em função de omega]]


Podemos, assim como no caso não estacionário da condução do calor (Método de Crank Nicholson), que realiza uma média entre os valores explícito e Implícito da Equação, o método da Super relaxação é da seguinte forma:
Como o método depende de <math> \omega  </math>, e se <math> \omega = 1 </math> temos a Relaxação normal, então podemos já observar que a convergência será mais demorada se <math> \omega < 1 </math>, já quando maior que 1, a solução sempre converge mais rápido até um <math> \omega_{otimo} </math> após o valor ótimo o valor cresce até 2, valor que diverge para o SOR.


<math> u_{ij}^{k+1} = u_{i j}^{k}(1 - \omega) + \omega u_{ij}^{R} </math>
Usaremos <math> \omega = \frac{2}{1+sen(\pi \Delta x)} </math> que é o valor ótimo.


tal que <math> u_{ij}^{R} </math> é o valor calculado através do método da Relaxação e <math> \omega \in [0:2]</math>.


=== Exemplos ===
=== Exemplos ===
Linha 212: Linha 240:
==== Exemplo 1 ====
==== Exemplo 1 ====


O primeiro problema é descrito pela seguinte expressão, para o domínio <math> \Omega = (0;L) \times (0;L) </math>:
O primeiro problema,que pode ser descrito como a temperatura <math> T = T(x, y)</math> de uma chapa que está em equilíbrio térmico com um dos lados a temperatura <math>T(x, 0) = L</math>, todos os outros lados iguais com temperatura 0 e sem geração. Para o domínio <math> \Omega = (0;L) \times (0;L) </math>, tal problema é descrito pela seguinte expressão:


<math>
<math>
\begin{cases}
\begin{cases}
\Delta u = 0  \forall x \in \Omega \\  
\Delta u = 0  \forall x,y \in \Omega \\  
u(x,0) = L  \big( \forall x \in [0;L] \big)  \\  
u(x,0) = L  \big( \forall x \in [0;L] \big)  \\  
u(x, L) = 0 \big( \forall x \in [0;L] \big) \\
u(x, L) = 0 \big( \forall x \in [0;L] \big) \\
Linha 224: Linha 252:
</math>
</math>


Foram obtidas as soluções mostradas nos gráficos a seguir, através do algoritmo de Gauss-Seidel.
A solução analítica para o problema é:
 
<math> u(x, y) = \frac{4L}{\pi}\sum_{n=0}^{\infty} \frac{sen(\frac{(2n+1) \pi x}{L})}{2n+1} \frac{senh\Big(\frac{(2n+1) \pi (L-y)}{L}\Big)}{ senh((2n+1)\pi)} </math>
 
Foi possível obter o gráfico da função analítica aproximada com os 110 primeiros termos parciais e estimar o erro entre a solução numérica e a solução analítica, que foi alto devido ao truncamento da série que forma a solução analítica.
 
[[Arquivo:solucaoanalitica1.png|450px|thumb|Solução analítica do Problema de Contorno 1]]
 
[[Arquivo:erroanalitico1.png|450px|thumb|Erro entre solução numérica e solução analítica em função do número de iterações]]


[[Arquivo:gaussseidel1.png|500px|Solução do Problema de Contorno através do Método de Relaxação (utilizando o algoritmo de Gauss-Seidel)]]
Foram obtidas as soluções mostradas nos gráficos a seguir, através do algoritmo de Gauss-Seidel, no caso da relaxação foram necessárias 1471 iterações para a solução convergir, já SOR levou 1000.
[[Arquivo:overr.png|500px|Solução do Problema de Contorno através do Método de Super Relaxação (utilizando o algoritmo de Gauss-Seidel)]]


[[Arquivo:relaxacao1k1471.png|450px|Solução do Problema de Contorno através do Método de Relaxação (utilizando o algoritmo de Gauss-Seidel)]]
[[Arquivo:sol1k1001.png|450px|Solução do Problema de Contorno através do Método de Super Relaxação (utilizando o algoritmo de Gauss-Seidel)]]


==== Exemplo 2 ====
==== Exemplo 2 ====
Linha 236: Linha 273:
<math>
<math>
\begin{cases}
\begin{cases}
\Delta u = 0  \forall x \in \Omega \\  
\Delta u = 0  \forall x,y \in \Omega \\  
u(x,0) = Lsin(\frac{\pi x}{L})  \big( \forall x \in [0;L] \big)  \\  
u(x,0) = Lsin(\frac{\pi x}{L})  \big( \forall x \in [0;L] \big)  \\  
u(x, L) = Lcos(\frac{\pi x}{L}) \big( \forall x \in [0;L] \big) \\
u(x, L) = Lcos(\frac{\pi x}{L}) \big( \forall x \in [0;L] \big) \\
Linha 246: Linha 283:
Foram obtidas as seguintes soluções mostradas nos gráficos, através do algoritmo de Gauss-Seidel.
Foram obtidas as seguintes soluções mostradas nos gráficos, através do algoritmo de Gauss-Seidel.


[[Arquivo:gaussseidel2.png|500px|Solução do Problema de Contorno através do Método de Relaxação (utilizando o algoritmo de Gauss-Seidel)]]
[[Arquivo:relaxacao2k2122.png|300px|Solução do Problema de Contorno através do Método de Relaxação (utilizando o algoritmo de Gauss-Seidel)]]
[[Arquivo:overr2.png|500px|Solução do Problema de Contorno através do Método de Super Relaxação (utilizando o algoritmo de Gauss-Seidel)]]
[[Arquivo:relaxacao2k1001.png|300px|Solução do Problema de Contorno através do Método de Super Relaxação (utilizando o algoritmo de Gauss-Seidel)]]
 
==== Exemplo 3 ====
 
O primeiro problema sobre a equação de Poisson , possui contornos fixos e uma gaussiana, pode ser considerado um problema de eletrostática com distribuição de carga <math> \frac{\rho}{\epsilon_{0}} = xy e^{-(x^2+ y^2)} </math> . Tal problema é descrito pela seguinte expressão, para o domínio <math> \Omega = (0;L) \times (0;L) </math>:
 
<math>
\begin{cases}
\Delta u = -g(x, y) = -xye^{-(x^2 + y^2)}  \forall x,y \in \Omega \\
u(x,0) = 0  \big( \forall x \in [0;L] \big)  \\
u(x, L) = 0 \big( \forall x \in [0;L] \big) \\
u(0, y) = 0 \big( \forall y \in [0;L] \big) \\
u(L, y) = 0 \big( \forall y \in [0;L] \big) .
\end{cases}
</math>
 
Foram obtidas as soluções mostradas nos gráficos a seguir, através do algoritmo de Gauss-Seidel.
 
[[Arquivo:poisson0.png|300px|Solução do Problema de Contorno através do Método de Relaxação]]
[[Arquivo:overrP0.png|300px|Solução do Problema de Contorno através do Método de Super Relaxação]]
 
==== Exemplo 4 ====
 
O segundo problema sobre a equação de Poisson, que é descrito pela seguinte expressão, para o domínio <math> \Omega = (0;L) \times (0;L) </math>:
 
<math>
\begin{cases}
\Delta u = -g(x, y) = -xye^{-(x^2 + y^2)}  \forall x,y \in \Omega \\
u(x,0) = 0  \big( \forall x \in [0;L] \big)  \\
u(x, L) = L \big( \forall x \in [0;L] \big) \\
u(0, y) = 0 \big( \forall y \in [0;L] \big) \\
u(L, y) = 0 \big( \forall y \in [0;L] \big) .
\end{cases}
</math>
 
Foram obtidas as soluções mostradas nos gráficos a seguir, através do algoritmo de Gauss-Seidel.
 
[[Arquivo:poisson1r.png|300px|Solução do Problema de Contorno através do Método de Relaxação]]
[[Arquivo:poisson1sor.png|300px|Solução do Problema de Contorno através do Método de Super Relaxação]]
 
==== Exemplo 5 ====
 
O segundo problema sobre a equação de Poisson, que é descrito pela seguinte expressão, para o domínio <math> \Omega = (0;L) \times (0;L) </math>:
 
<math>
\begin{cases}
\Delta u = -g(x, y) = -xye^{-(x^2 + y^2)}  \forall x,y \in \Omega \\
u(x,0) = Lsin(\frac{\pi x}{L})  \big( \forall x \in [0;L] \big)  \\
u(x, L) = Lcos(\frac{\pi x}{L}) \big( \forall x \in [0;L] \big) \\
u(0, y) = \frac{x^{2}}{L} \big( \forall y \in [0;L] \big) \\
u(L, y) = \frac{(L-x)^{2}}{L} \big( \forall y \in [0;L] \big) .
\end{cases}
</math>
 
Foram obtidas as soluções mostradas nos gráficos a seguir, através do algoritmo de Gauss-Seidel.
 
 
[[Arquivo:poisson2r.png|300px|Solução do Problema de Contorno através do Método de Relaxação]]
[[Arquivo:poisson2sor.png|300px|Solução do Problema de Contorno através do Método de Super Relaxação]]
 
O programa utilizado para gerar as soluçoes e erros foi o seguir (ou com pequenas alteraçoes):
 
==== Programa ====
Trechos do programa realizado para os exemplos acima.
 
Programa para o método de Relaxação (Equação de Laplace):
 
<source lang="c">
#include<stdio.h>
#include<math.h>
#define N 1000
#define M 70
#define P 1
#define PI 3.141529
void gaussseidelL(){ 
  double u[M+2][M+2];
  double dx=0, dy=0;
  double L=5., parada=0, erro=0.00001, up=0;
  int i=0, j=0, k=1, a=0;
  dx = L/(M+1);
  dy = L/(M+1);
  for(i=1;i<M+1;i++){
    for(j=1;j<M+1;j++){
      u[i][j] = 1.;
    }
  }
  /* Primeira Solução, GaussSeidel1  */
  for(i=0;i<M+1;i++){
    u[0][i] = L;
    u[M+1][i] = 0.0;
    u[i][0] = 0.0;
    u[i][M+1] = 0.0;
  }
  u[0][M+1] = L;
  u[M+1][M+1] = 0.0;
  /* Segunda Solução, GaussSeidel2
  for(i=0;i<M+1;i++){
    u[0][i] = L*pow(cos(i*dx*PI/L),2);
    u[M+1][i] = L*pow(sin(i*dx*PI/L),2);
    u[i][0] = L -  pow(i*dx,2)/L;
    u[i][M+1] = pow(((M+1)-i)*dx,2)/L;
  }
  u[0][M+1] = L;
  u[M+1][M+1] = 0.0;
  */
  do{
    up = (u[1][1] + u[M/2+1][M/2+1] + u[M][M])/3;
    for(i=1;i<M+1;i++){
      for(j=1;j<M+1;j++){
u[i][j] = (u[i+1][j] + u[i-1][j] + u[i][j+1] + u[i][j-1])/4;
      }
    }
    k++;
    parada = fabs((up-(u[1][1] + u[M/2+1][M/2+1] + u[M][M])/3)/up);
    if(parada < erro && k>N/5){
      a = 1;
    }else{
      a = 0;
    }
  }while(a == 0);
  for(i=0;i<M+2;i++){
    for(j=0;j<M+2;j++){
      printf("%lf %lf %lf \n",i*dx, j*dy, u[i][j]);
    }
  }
  printf("\n\n");
  return;
}
</source>
 
Segundo trecho para método de Relaxação, Equação de Poisson
 
<source lang="c">
 
#include<stdio.h>
#include<math.h>
#define N 1000
#define M 70
#define P 1
#define PI 3.141529
 
 
void gaussseidelP(){
 
  double u[M+2][M+2], F[M+2][M+2];
  double dx=0, dy=0;
  double L=5., parada=0, erro=0.00001, up=0;
  int i=0, j=0, k=1, a=0;
 
  dx = L/(M+1);
  dy = L/(M+1);
 
  for(i=1;i<M+1;i++){
    for(j=1;j<M+1;j++){
      u[i][j] = 0.;
      F[i][j] = i*dx*j*dy*exp(-(pow(i*dx,2) + pow(j*dy,2))/L);
    }
  }
 
/* Solução Zero, Poisson0  */
  for(i=0;i<M+1;i++){
    u[0][i] = 0;
    u[M+1][i] = 0.0;
    u[i][0] = 0.0;
    u[i][M+1] = 0.0;
  }
  u[0][M+1] = 0;
  u[M+1][M+1] = 0.0;
 
 
  /* Primeira Solução, PoissonGS1
  for(i=0;i<M+1;i++){
    u[0][i] = L;
    u[M+1][i] = 0.0;
    u[i][0] = 0.0;
    u[i][M+1] = 0.0;
  }
  u[0][M+1] = L;
  u[M+1][M+1] = 0.0;
  */
 
 
  /* Segunda Solução, PoissonGS2
  for(i=0;i<M+1;i++){
  u[0][i] = L*pow(cos(i*dx*PI/L),2);
  u[M+1][i] = L*pow(sin(i*dx*PI/L),2);
  u[i][0] = L -  pow(i*dx,2)/L;
  u[i][M+1] = pow(((M+1)-i)*dx,2)/L;
  }
  u[0][M+1] = L;
  u[M+1][M+1] = 0.0;
  */
 
 
  do{   
    up = (u[1][1] + u[M/2+1][M/2+1] + u[M][M])/3;
   
    for(i=1;i<M+1;i++){
      for(j=1;j<M+1;j++){
u[i][j] = (u[i+1][j] + u[i-1][j] + u[i][j+1] + u[i][j-1] + 4*dx*dx*F[i][j])/4;
      }
    }
   
    k++;
    parada = fabs((up-(u[1][1] + u[M/2+1][M/2+1] + u[M][M])/3)/up);
    if(parada < erro && k>N/5){
      a = 1;
    }else{
      a = 0;
    }
  }while(a == 0);
 
  for(i=0;i<M+2;i++){
    for(j=0;j<M+2;j++){
      printf("%lf %lf %lf \n",i*dx, j*dy, u[i][j]);
    }
  }
 
  printf("\n\n");
  return;
}
 
</source>
 
Trecho de programa que utiliza o método de Super Relaxação para Equação de Laplace:
 
<source lang="c">
 
#include<stdio.h>
#include<math.h>
#define N 1000
#define M 70
#define P 1
#define PI 3.141529
 
 
void overrelaxationL(){
 
  double u[M+2][M+2], un[M+2][M+2];
  double dx=0, dy=0, omega=1.;
  double L=5., parada=0, erro=0.00005, up=0;
  int i=0, j=0, k=1, a=0;
 
  omega = 2/(1+PI*dx);
 
    k = 0;
   
    dx = L/(M+1);
    dy = L/(M+1);
   
    for(i=1;i<M+1;i++){
      for(j=1;j<M+1;j++){
un[i][j] = 1.;
u[i][j] = un[i][j];
      }
    }
    /*
    for(i=0;i<M+1;i++){
      un[0][i] = 0.0;
      un[M+1][i] = L;
      un[i][0] = 0.0;
      un[i][M+1] = 0.0;
    }
    un[0][M+1] = 0.0;
    un[M+1][M+1] = L;
    */
   
    for(i=0;i<M+1;i++){
      un[0][i] = L*pow(cos(i*dx*PI/L),2);
      un[M+1][i] = L*pow(sin(i*dx*PI/L),2);
      un[i][0] = L -  pow(i*dx,2)/L;
      un[i][M+1] = pow(((M+1)-i)*dx,2)/L;
    }
    un[0][M+1] = L;
    un[M+1][M+1] = 0.0; 
   
    do{
     
      up = (un[1][1] + 4*un[M/2+1][M/2+1] + un[M][M] + un[1][M] + un[M][1])/8;
     
      for(i=1;i<M+1;i++){
for(j=1;j<M+1;j++){
  un[i][j] = (un[i+1][j] + un[i-1][j] + un[i][j+1] + un[i][j-1])/4;
  u[i][j] = u[i][j]*(1 - omega) + omega*un[i][j];
  un[i][j] = u[i][j];
}
      }
      k++;
      parada = fabs((up-(un[1][1] + 4*un[M/2+1][M/2+1] + un[M][M] + un[1][M] + un[M][1])/8)/up);
      if(parada < erro && k>N){
a = 1;
      }else{
a = 0;
      }
    }while(a == 0);
   
    for(i=0;i<M+2;i++){
      for(j=0;j<M+2;j++){
printf("%lf %lf %lf \n",i*dx, j*dy, un[i][j]);
      }
    }
 
  printf("\n\n");
}
 
</source>
 
Trecho de programa do algoritmo de Super Relaxação para Equação de Poisson:
 
 
<source lang="c">
 
#include<stdio.h>
#include<math.h>
#define N 1000
#define M 70
#define P 1
#define PI 3.141529
 
 
void overrelaxationP(){
 
  double u[M+2][M+2], un[M+2][M+2], F[M+2][M+2];
  double dx=0, dy=0, omega=1.;
  double L=5., parada=0, erro=0.00005, up=0;
  int i=0, j=0, k=1, a=0;
 
  omega = 2/(1+PI*dx);
 
    k = 0;
   
    dx = L/(M+1);
    dy = L/(M+1);
   
    for(i=1;i<M+1;i++){
      for(j=1;j<M+1;j++){
un[i][j] = 1.;
u[i][j] = un[i][j];
F[i][j] = i*dx*j*dy*exp(-(pow(i*dx,2) + pow(j*dy,2))/L);
      }
    }
 
     
    for(i=0;i<M+1;i++){
      un[0][i] = 0.0;
      un[M+1][i] = 0.0;
      un[i][0] = 0.0;
      un[i][M+1] = 0.0;
    }
    un[0][M+1] = 0.0;
    un[M+1][M+1] = 0.0;
   
 
    /*
    for(i=0;i<M+1;i++){
      un[0][i] = L*pow(cos(i*dx*PI/L),2);
      un[M+1][i] = L*pow(sin(i*dx*PI/L),2);
      un[i][0] = L -  pow(i*dx,2)/L;
      un[i][M+1] = pow(((M+1)-i)*dx,2)/L;
    }
    un[0][M+1] = L;
    un[M+1][M+1] = 0.0; 
    */
   
    do{
     
      up = (un[1][1] + 4*un[M/2+1][M/2+1] + un[M][M] + un[1][M] + un[M][1])/8;
     
      for(i=1;i<M+1;i++){
for(j=1;j<M+1;j++){
  un[i][j] = (un[i+1][j] + un[i-1][j] + un[i][j+1] + un[i][j-1] + 4*dx*dx*F[i][j])/4;
  u[i][j] = u[i][j]*(1 - omega) + omega*un[i][j];
  un[i][j] = u[i][j];
}
      }
      k++;
      parada = fabs((up-(un[1][1] + 4*un[M/2+1][M/2+1] + un[M][M] + un[1][M] + un[M][1])/8)/up);
      if(parada < erro && k>N){
a = 1;
      }else{
a = 0;
      }
    }while(a == 0);
   
    for(i=0;i<M+2;i++){
      for(j=0;j<M+2;j++){
printf("%lf %lf %lf \n",i*dx, j*dy, un[i][j]);
      }
    }
 
  printf("\n\n");
}
 
</source>
 
== Referências ==
 
#[[http://www.mat.ufmg.br/~rodney/notas_de_aula/metodos_numericos_EDPs_elipticas.pdf]] Biezuner, Rodney Josué,''Notas de Aula Métodos Numéricos para Equações Diferenciais Parciais Elípticas'' (UFMG, 2007)
#[[http://www.mat.ufmg.br/~rodney/notas_de_aula/iedp.pdf]] Biezuner, Rodney Josué,''Introdução às Equações Diferenciais Parciais'' (UFMG, 2007)
#Jianping Zhu, ''Solving Partial Differential Equations on Parallel Computers'' (World Scientific, 1994)
 
Grupo: Pedro Barbisan Widniczck

Edição atual tal como às 13h18min de 2 de dezembro de 2017

A equação de Poisson:

é uma equação do tipo Elíptica que representa fenômenos físicos estácionarios relacionados a Eletrostatica, Dinâmica de Fluídos e Transferência de Calor. Se a equação passa a ser chamada de Equação de Laplace. Os problemas relacionados a equação de Laplace são estudados pela "Teoria do Potencial".

As soluções da Equação de Laplace são denominadas funções Harmônicas. Os problemas mais habituais na vida de um físico, engenheiro ou matemático ao se depararem com uma EDP, são os problemas com Condições de Contorno em um dominío , essencialmente será trabalhada a Condição de Dirichlet, que possui fronteiras () conhecidas, tendo o seguinte formato:

A equação de Poisson possui forma parecida para o Problema de Dirichlet, que fica:


Para tais problemas, estudaremos os métodos de Relaxação e Super-Relaxação para encontrar as soluções da Equação de Laplace na região de Quadrado de Lado .

Dominio Quadrado de Lado

Solução Analítica da Equação de Laplace

Seja o problema em , temos:

sendo

Separamos o problema geral de Dirchlet em 4 problemas "menores", com condições de contorno diferentes de zero em apenas um trecho da fronteira, de modo que obtemos desde:

...

até:


Podemos então utilizar o Método da Separação de Variáveis para resolver os 4 problemas e, como a Equação de Laplace é linear, sua soma será a solução completa do Problema de Dirichlet. O método consiste em supor , para então, ao substituirmos na equação obtermos a seguinte expressão:

Podemos isolar as funções e , de fato ficamos com com duas relações que dependem ou apenas de e portanto para elas serem sempre iguais, é necessário que sejam constantes ():

Assim obtemos 2 EDOs de segunda ordem, que podem ser resolvidas pelo Método dos Coeficientes a Determinar. Como não é objetivo aqui realizar cálculos analíticos (especialmente "na mão") apenas será resolvido o primeiro problema ():

As condições de contorno mostram que , e .

Dividindo o problema, temos a parte de


Supondo uma solução da forma :


Ou seja, temos a solução de sendo

Partindo para a segunda equação ,


supondo solução do tipo temos:

Ou seja, temos solução

Utilizando a primeira C.C. obtemos

ou seja, temos que

Utilizando a segunda C.C. temos

ou seja, existem infinitos tal que é solução.

Voltando a , temos

Finalmente unindo as respostas, temos

sendo

Para os outros problemas, temos soluções parecidas:

sendo

sendo

sendo

A solução completa do problema de Dirichlet no quadrado de Lado é a soma das quatro soluções parciais: .

Método da Relaxação

O Método da Relaxação é um método iterativo utilizado para obter a solução numérica para a equação de Laplace e Poisson. A ideia do método é de, utilizando a vizinhança iterar os pontos da malha até que convirjam para uma solução.

Discretizando a equação temos e para , e a função , nos deparamos com uma matriz quadrada sendo as bordas , , e .

Realizando-se a discretização, podemos tomar as derivadas:

e

Substituindo na Equação, temos

, ou seja:

,

ou mais geralmente (supondo ):

para

Para condição de parada, foi convencionado tomar o erro relativo entre as iterações e , para estimar o erro, se optou por tomar como valores, o ponto médio da malha (já que é o ultimo ponto a ser alcançado nas iterações, portanto, quando sua variação diminuir é sinal de que a solução já está convergindo), para observarmos a evolução em relação a outros pontos que variam desde o inicio, foram utilizados as diagonais interiores, tal que o erro relativo é:

fazendo a média ponderada com peso 4 para o ponto médio:

Método da Super Relaxação (SOR)

Podemos também realizar uma média entre os valores já calculados e os ainda não calculados na iteração, o método da Super Relaxação ou Sobrerrelaxação (SOR) é da seguinte forma:

tal que é o valor calculado através do método da Relaxação e é o fator de relaxamento, se temos a Relaxação normal.

Estabilidade

Relaxação

Erro relativo (método da Relaxação) em função da quantidade de iterações

A relaxação é um método Iterativo sobre os pontos vizinhos que pode ser feita de 2 modos, pelo Algoritmo de Jacobi, e pelo de Gauss-Seidel.

O algoritmo de Jacobi pega valores "antigos" para a iteração e possui convergencia muito lenta, por isso não é muito utilizado. Já o algoritmo de Gauss-Seidel pega os valores "novos" (que ja foram calculados na iteração) e os "antigos" (que não foram calculados na iteração ), possui convergência mais rapida, porém ainda é lenta.

Como a média definida anteriormente foi feita utilizando o ponto médio do domínio, o erro cresce após decair, pois é quando efetivamente ocorrem variações maiores no ponto médio.

Para a relaxação, o algoritmo de Jacobi faz o seguinte cálculo:

,

já o algoritmo de Gauss-Seidel faz:

,

Algoritmos iterativos tendem a convergir para solução unica, se a matriz que as representa for Diagonal Dominante, ou seja:

Sendo e , então e e

Como , a desigualdade vale e o método converge.

De fato, podemos ver que a equação de Laplace respeita tal desigualdade.


Caso façamos um retangulo com , obtemos os erros da imagem a seguir, feito utilizando e , e e e :

Erro obtido para diferentes comprimentos do retangulo

Super Relaxação

Grafico mostra a quantidades de Iterações para convergência do exemplo 1 em função de omega

Como o método depende de , e se temos a Relaxação normal, então podemos já observar que a convergência será mais demorada se , já quando maior que 1, a solução sempre converge mais rápido até um após o valor ótimo o valor cresce até 2, valor que diverge para o SOR.

Usaremos que é o valor ótimo.


Exemplos

Foram realizados 5 exemplos, 2 sobre a equação de Laplace e 3 sobre a equação de Poisson.

Exemplo 1

O primeiro problema,que pode ser descrito como a temperatura de uma chapa que está em equilíbrio térmico com um dos lados a temperatura , todos os outros lados iguais com temperatura 0 e sem geração. Para o domínio , tal problema é descrito pela seguinte expressão:

A solução analítica para o problema é:

Foi possível obter o gráfico da função analítica aproximada com os 110 primeiros termos parciais e estimar o erro entre a solução numérica e a solução analítica, que foi alto devido ao truncamento da série que forma a solução analítica.

Solução analítica do Problema de Contorno 1
Erro entre solução numérica e solução analítica em função do número de iterações

Foram obtidas as soluções mostradas nos gráficos a seguir, através do algoritmo de Gauss-Seidel, no caso da relaxação foram necessárias 1471 iterações para a solução convergir, já SOR levou 1000.

Solução do Problema de Contorno através do Método de Relaxação (utilizando o algoritmo de Gauss-Seidel) Solução do Problema de Contorno através do Método de Super Relaxação (utilizando o algoritmo de Gauss-Seidel)

Exemplo 2

O segundo problema é descrito pela seguinte expressão, para o domínio :

Foram obtidas as seguintes soluções mostradas nos gráficos, através do algoritmo de Gauss-Seidel.

Solução do Problema de Contorno através do Método de Relaxação (utilizando o algoritmo de Gauss-Seidel) Solução do Problema de Contorno através do Método de Super Relaxação (utilizando o algoritmo de Gauss-Seidel)

Exemplo 3

O primeiro problema sobre a equação de Poisson , possui contornos fixos e uma gaussiana, pode ser considerado um problema de eletrostática com distribuição de carga . Tal problema é descrito pela seguinte expressão, para o domínio :

Foram obtidas as soluções mostradas nos gráficos a seguir, através do algoritmo de Gauss-Seidel.

Solução do Problema de Contorno através do Método de Relaxação Solução do Problema de Contorno através do Método de Super Relaxação

Exemplo 4

O segundo problema sobre a equação de Poisson, que é descrito pela seguinte expressão, para o domínio :

Foram obtidas as soluções mostradas nos gráficos a seguir, através do algoritmo de Gauss-Seidel.

Solução do Problema de Contorno através do Método de Relaxação Solução do Problema de Contorno através do Método de Super Relaxação

Exemplo 5

O segundo problema sobre a equação de Poisson, que é descrito pela seguinte expressão, para o domínio :

Foram obtidas as soluções mostradas nos gráficos a seguir, através do algoritmo de Gauss-Seidel.


Solução do Problema de Contorno através do Método de Relaxação Solução do Problema de Contorno através do Método de Super Relaxação

O programa utilizado para gerar as soluçoes e erros foi o seguir (ou com pequenas alteraçoes):

Programa

Trechos do programa realizado para os exemplos acima.

Programa para o método de Relaxação (Equação de Laplace):

 #include<stdio.h>
 #include<math.h>
 #define N 1000
 #define M 70
 #define P 1
 #define PI 3.141529
 void gaussseidelL(){  
  double u[M+2][M+2];
  double dx=0, dy=0;
  double L=5., parada=0, erro=0.00001, up=0;
  int i=0, j=0, k=1, a=0;
  dx = L/(M+1);
  dy = L/(M+1);
  for(i=1;i<M+1;i++){
    for(j=1;j<M+1;j++){
      u[i][j] = 1.;
    }
  }
  /* Primeira Solução, GaussSeidel1  */
  for(i=0;i<M+1;i++){
    u[0][i] = L;
    u[M+1][i] = 0.0;
    u[i][0] = 0.0;
    u[i][M+1] = 0.0;
  }
  u[0][M+1] = L;
  u[M+1][M+1] = 0.0;
  /* Segunda Solução, GaussSeidel2
  for(i=0;i<M+1;i++){
    u[0][i] = L*pow(cos(i*dx*PI/L),2);
    u[M+1][i] = L*pow(sin(i*dx*PI/L),2);
    u[i][0] = L -  pow(i*dx,2)/L;
    u[i][M+1] = pow(((M+1)-i)*dx,2)/L;
  }
  u[0][M+1] = L;
  u[M+1][M+1] = 0.0;
  */
  do{
    up = (u[1][1] + u[M/2+1][M/2+1] + u[M][M])/3;
    for(i=1;i<M+1;i++){
      for(j=1;j<M+1;j++){
	u[i][j] = (u[i+1][j] + u[i-1][j] + u[i][j+1] + u[i][j-1])/4;
      }
    }
    k++;
    parada = fabs((up-(u[1][1] + u[M/2+1][M/2+1] + u[M][M])/3)/up);
    if(parada < erro && k>N/5){
      a = 1;
    }else{
      a = 0;
    }
  }while(a == 0);
  for(i=0;i<M+2;i++){
    for(j=0;j<M+2;j++){
      printf("%lf %lf %lf \n",i*dx, j*dy, u[i][j]);
    }
  }
  printf("\n\n");
  return;
 }

Segundo trecho para método de Relaxação, Equação de Poisson

#include<stdio.h>
#include<math.h>
#define N 1000
#define M 70
#define P 1
#define PI 3.141529


void gaussseidelP(){
  
  double u[M+2][M+2], F[M+2][M+2];
  double dx=0, dy=0;
  double L=5., parada=0, erro=0.00001, up=0;
  int i=0, j=0, k=1, a=0;

  dx = L/(M+1);
  dy = L/(M+1);
  
  for(i=1;i<M+1;i++){
    for(j=1;j<M+1;j++){
      u[i][j] = 0.;
      F[i][j] = i*dx*j*dy*exp(-(pow(i*dx,2) + pow(j*dy,2))/L);
    }
  }

 /* Solução Zero, Poisson0  */
  for(i=0;i<M+1;i++){
    u[0][i] = 0;
    u[M+1][i] = 0.0;
    u[i][0] = 0.0;
    u[i][M+1] = 0.0;
  }
  u[0][M+1] = 0;
  u[M+1][M+1] = 0.0;

  
  /* Primeira Solução, PoissonGS1 
  for(i=0;i<M+1;i++){
    u[0][i] = L;
    u[M+1][i] = 0.0;
    u[i][0] = 0.0;
    u[i][M+1] = 0.0;
  }
  u[0][M+1] = L;
  u[M+1][M+1] = 0.0;
  */
  
  
  /* Segunda Solução, PoissonGS2 
  for(i=0;i<M+1;i++){
  u[0][i] = L*pow(cos(i*dx*PI/L),2);
  u[M+1][i] = L*pow(sin(i*dx*PI/L),2);
  u[i][0] = L -  pow(i*dx,2)/L;
  u[i][M+1] = pow(((M+1)-i)*dx,2)/L;
  }
  u[0][M+1] = L;
  u[M+1][M+1] = 0.0;
  */
  
  
  do{    
    up = (u[1][1] + u[M/2+1][M/2+1] + u[M][M])/3;
    
    for(i=1;i<M+1;i++){
      for(j=1;j<M+1;j++){
	u[i][j] = (u[i+1][j] + u[i-1][j] + u[i][j+1] + u[i][j-1] + 4*dx*dx*F[i][j])/4;
      }
    }
    
    k++;
    parada = fabs((up-(u[1][1] + u[M/2+1][M/2+1] + u[M][M])/3)/up);
    if(parada < erro && k>N/5){
      a = 1;
    }else{
      a = 0;
    }
  }while(a == 0);

  for(i=0;i<M+2;i++){
    for(j=0;j<M+2;j++){
      printf("%lf %lf %lf \n",i*dx, j*dy, u[i][j]);
    }
  }

  printf("\n\n");
  return;
}

Trecho de programa que utiliza o método de Super Relaxação para Equação de Laplace:

#include<stdio.h>
#include<math.h>
#define N 1000
#define M 70
#define P 1
#define PI 3.141529


void overrelaxationL(){
  
  double u[M+2][M+2], un[M+2][M+2];
  double dx=0, dy=0, omega=1.;
  double L=5., parada=0, erro=0.00005, up=0;
  int i=0, j=0, k=1, a=0;
  
  omega = 2/(1+PI*dx);

    k = 0;
    
    dx = L/(M+1);
    dy = L/(M+1);
    
    for(i=1;i<M+1;i++){
      for(j=1;j<M+1;j++){
	un[i][j] = 1.;
	u[i][j] = un[i][j];
      }
    }
    /*
    for(i=0;i<M+1;i++){
      un[0][i] = 0.0;
      un[M+1][i] = L;
      un[i][0] = 0.0;
      un[i][M+1] = 0.0;
    }
    un[0][M+1] = 0.0;
    un[M+1][M+1] = L;
    */
    
    for(i=0;i<M+1;i++){
      un[0][i] = L*pow(cos(i*dx*PI/L),2);
      un[M+1][i] = L*pow(sin(i*dx*PI/L),2);
      un[i][0] = L -  pow(i*dx,2)/L;
      un[i][M+1] = pow(((M+1)-i)*dx,2)/L;
    }
    un[0][M+1] = L;
    un[M+1][M+1] = 0.0;   
    
    do{
      
      up = (un[1][1] + 4*un[M/2+1][M/2+1] + un[M][M] + un[1][M] + un[M][1])/8;
      
      for(i=1;i<M+1;i++){
	for(j=1;j<M+1;j++){
	  un[i][j] = (un[i+1][j] + un[i-1][j] + un[i][j+1] + un[i][j-1])/4;
	  u[i][j] = u[i][j]*(1 - omega) + omega*un[i][j];
	  un[i][j] = u[i][j];
	}
      }
      k++;
      parada = fabs((up-(un[1][1] + 4*un[M/2+1][M/2+1] + un[M][M] + un[1][M] + un[M][1])/8)/up);
      if(parada < erro && k>N){
	a = 1;
      }else{
	a = 0;
      }
    }while(a == 0);
    
    for(i=0;i<M+2;i++){
      for(j=0;j<M+2;j++){
	printf("%lf %lf %lf \n",i*dx, j*dy, un[i][j]);
      }
    }

  printf("\n\n");
}

Trecho de programa do algoritmo de Super Relaxação para Equação de Poisson:


#include<stdio.h>
#include<math.h>
#define N 1000
#define M 70
#define P 1
#define PI 3.141529


void overrelaxationP(){
  
  double u[M+2][M+2], un[M+2][M+2], F[M+2][M+2];
  double dx=0, dy=0, omega=1.;
  double L=5., parada=0, erro=0.00005, up=0;
  int i=0, j=0, k=1, a=0;
  
  omega = 2/(1+PI*dx);

    k = 0;
    
    dx = L/(M+1);
    dy = L/(M+1);
    
    for(i=1;i<M+1;i++){
      for(j=1;j<M+1;j++){
	un[i][j] = 1.;
	u[i][j] = un[i][j];
	F[i][j] = i*dx*j*dy*exp(-(pow(i*dx,2) + pow(j*dy,2))/L);
      }
    }

       
    for(i=0;i<M+1;i++){
      un[0][i] = 0.0;
      un[M+1][i] = 0.0;
      un[i][0] = 0.0;
      un[i][M+1] = 0.0;
    }
    un[0][M+1] = 0.0;
    un[M+1][M+1] = 0.0;
    

    /*
    for(i=0;i<M+1;i++){
      un[0][i] = L*pow(cos(i*dx*PI/L),2);
      un[M+1][i] = L*pow(sin(i*dx*PI/L),2);
      un[i][0] = L -  pow(i*dx,2)/L;
      un[i][M+1] = pow(((M+1)-i)*dx,2)/L;
    }
    un[0][M+1] = L;
    un[M+1][M+1] = 0.0;   
    */
    
    do{
      
      up = (un[1][1] + 4*un[M/2+1][M/2+1] + un[M][M] + un[1][M] + un[M][1])/8;
      
      for(i=1;i<M+1;i++){
	for(j=1;j<M+1;j++){
	  un[i][j] = (un[i+1][j] + un[i-1][j] + un[i][j+1] + un[i][j-1] + 4*dx*dx*F[i][j])/4;
	  u[i][j] = u[i][j]*(1 - omega) + omega*un[i][j];
	  un[i][j] = u[i][j];
	}
      }
      k++;
      parada = fabs((up-(un[1][1] + 4*un[M/2+1][M/2+1] + un[M][M] + un[1][M] + un[M][1])/8)/up);
      if(parada < erro && k>N){
	a = 1;
      }else{
	a = 0;
      }
    }while(a == 0);
    
    for(i=0;i<M+2;i++){
      for(j=0;j<M+2;j++){
	printf("%lf %lf %lf \n",i*dx, j*dy, un[i][j]);
      }
    }

  printf("\n\n");
}

Referências

  1. [[1]] Biezuner, Rodney Josué,Notas de Aula Métodos Numéricos para Equações Diferenciais Parciais Elípticas (UFMG, 2007)
  2. [[2]] Biezuner, Rodney Josué,Introdução às Equações Diferenciais Parciais (UFMG, 2007)
  3. Jianping Zhu, Solving Partial Differential Equations on Parallel Computers (World Scientific, 1994)

Grupo: Pedro Barbisan Widniczck