Grupo1 - Dif em 2D

De Física Computacional
Revisão de 09h10min de 24 de outubro de 2017 por Pedrow7 (discussão | contribs)
Ir para navegação Ir para pesquisar

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

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 no Quadrado de Lado .

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" tal que obtemos os problemas desde:

...

até:


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 , 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 de suas variações, 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

Utilizando a C.C. vemos que .

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

Algoritmo de Relaxaçã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

Como condição de parada, foi convencionado tomar o Erro Relativo entre as iterações e , para estimar o erro se faz:


Estabilidade

Erro relativo ao método de 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) e os "antigos" (que não foram calculados), possui convergencia mais rapida, porém ainda é lenta.

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

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


Caso façamos um retangulo com , obtemos o erro da imagem a seguir:

Erro obtido para diferentes comprimentos do retangulo

Método da Super Relaxação

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:

tal que é o valor calculado através do método da Relaxação e .

Usaremos que é um otimizador.

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 é 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 (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 não poderia ser diferente, 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 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/(M+1));

    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/(M+1));

    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");
}