Exibir código-fonte para Código 1

       MODULE globais
           IMPLICIT NONE
           INTEGER :: N,L 
       END MODULE

       PROGRAM exemplo1_laplace
       USE globais
       IMPLICIT NONE
       INTEGER :: i,j,k,c,t,edge
       REAL :: dV
       REAL,DIMENSION(:,:),ALLOCATABLE :: V,V1,E1,E2,x,y 
     
       N = 10
       edge = 6*N
       L = edge + 1

       ALLOCATE(V(L,L),V1(L,L),E1(L,L),E2(L,L),x(L,L),y(L,L)) 
       
       DO i=1,L
           DO j=1,L
               x(j,i) = -1.0 + (i-1)*2.0/(L-1)
               y(j,i) = -1.0 + (j-1)*2.0/(L-1)
           END DO
       END DO

!    Calcula utilizando o algoritmo de Jacobi 
       CALL inicia_zeros(E1)
       CALL inicia_zeros(E2)
       CALL initialize_V(V)
       
       c = 0
       V1 = V
       
       CALL laplace_jacobi(V,V1,dV,c)
       CALL electric_field(V,E1,E2)
       
       WRITE(*,*) 'Jacobi:',dV,c
       
       CALL dados(x,y,V,E1,E2)

!    Calcula utilizando o algoritmo de Gauss-Siedel
       CALL inicia_zeros(E1)
       CALL inicia_zeros(E2)
       CALL initialize_V(V)
       
       c = 0
       
       CALL laplace_gauss_siedel(V,dV,c)
       CALL electric_field(V,E1,E2)
       
       WRITE(*,*) 'Gauss-Siedel:',dV,c
       
       CALL dados(x,y,V,E1,E2)

       END PROGRAM

       SUBROUTINE inicia_zeros(M)
!        Inicializa uma matriz em zeros
           USE globais
           IMPLICIT NONE
           INTEGER :: i,j
           REAL,DIMENSION(L,L) :: M
           DO i=1,L
               DO j=1,L
                   M(j,i) = 0.0
               END DO
           END DO
       END SUBROUTINE inicia_zeros
       
       SUBROUTINE initialize_V(V0)
!        Inicialização do potencial
           USE globais
           IMPLICIT NONE
           INTEGER :: i
           REAL,DIMENSION(L,L) :: V0
           DO i=1,L
               V0(i,1) = -1.0 
               V0(i,L) = 1.0 
               V0(1,i) = -1.0 + (i-1)*2.0/(L-1) 
               V0(L,i) = -1.0 + (i-1)*2.0/(L-1) 
           END DO
       END SUBROUTINE initialize_V

       SUBROUTINE update_V_jacobi(Vn,Vn1,deltaV)
!        Calcula a atualização uma vez
           USE globais
           IMPLICIT NONE
           INTEGER :: i,j
           REAL :: deltaV 
           REAL,DIMENSION(L,L) :: Vn,Vn1
           deltaV = 0.0 
           DO i=2,L-1
               DO j=2,L-1
                   Vn1(j,i) = (Vn(j+1,i) + Vn(j-1,i) + Vn(j,i+1) + Vn(j &
     ,i-1))/4.0
                   deltaV = deltaV + abs(Vn(j,i) - Vn1(j,i))
               END DO
           END DO
       END SUBROUTINE update_V_jacobi

       SUBROUTINE laplace_jacobi(Vm,Vm1,delta_v,cont)
!        Implementa o método de Jacobi
           USE globais
           IMPLICIT NONE
           INTEGER :: cont,x
           REAL :: delta_v,tol
           REAL,DIMENSION(L,L) :: Vm,Vm1
           tol = 1e-5*L*L 
           delta_v = 10.0
           DO WHILE (delta_v .gt. tol)
               CALL update_V_jacobi(Vm,Vm1,delta_v)
               CALL update_V_jacobi(Vm1,Vm,delta_v)
               cont = cont+1
           END DO
       END SUBROUTINE laplace_jacobi

       SUBROUTINE update_V_gauss_siedel(Vn,deltaV)
!        Calcula a atualização uma vez
           USE globais
           IMPLICIT NONE
           INTEGER :: i,j
           REAL :: deltaV,v1,v2
           REAL,DIMENSION(L,L) :: Vn
           deltaV = 0.0
           DO i=2,L-1
               DO j=2,L-1
                   v1 = Vn(j,i)
                   Vn(j,i) = (Vn(j+1,i) + Vn(j-1,i) + Vn(j,i+1) + Vn(j, &
     i-1))/4.0
                   v2 = Vn(j,i)
                   deltaV = deltaV + abs(v2 - v1)
               END DO
           END DO
       END SUBROUTINE update_V_gauss_siedel
       
       SUBROUTINE laplace_gauss_siedel(Vm,delta_v,cont)
!        Implementa o método de Gauss_Siedel
           USE globais
           IMPLICIT NONE
           INTEGER :: cont
           REAL :: delta_v,tol
           REAL,DIMENSION(L,L) :: Vm
           tol = 1e-5*L*L 
           delta_v = 10.0
           DO WHILE (delta_v .gt. tol)
               CALL update_V_gauss_siedel(Vm,delta_v)
               CALL update_V_gauss_siedel(Vm,delta_v)
               cont = cont+1
           END DO
       END SUBROUTINE laplace_gauss_siedel

       SUBROUTINE electric_field(Vk,Ex,Ey)
!        Calcula o campo elétrico a partir do potencial   
           USE globais
           IMPLICIT NONE
           INTEGER :: i,j
           REAL,DIMENSION(L,L) :: Vk,Ex,Ey
!        Derivada para frente (bordas)         
           DO i=1,L-1
               Ex(1,i) = - (Vk(1,i+1) - Vk(1,i)) 
               Ex(L,i) = - (Vk(L,i+1) - Vk(L,i))               
               Ey(i,1) = - (Vk(i+1,1) - Vk(i,1)) 
               Ey(i,L) = - (Vk(i+1,L) - Vk(i,L)) 
           END DO
!        Derivada para trás (bordas)
           Ex(1,L) = - (Vk(1,L) - Vk(1,L-1))
           Ex(L,L) = - (Vk(L,L) - Vk(L,L-1))
           Ey(L,1) = - (Vk(L,1) - Vk(L-1,1))
           Ey(L,L) = - (Vk(L,L) - Vk(L-1,L))
!        Derivada centrada
           DO i=2,L-1
               DO j=2,L-1 
                   Ex(j,i) = - (Vk(j,i+1) - Vk(j,i-1))/2.0 
                   Ey(j,i) = - (Vk(j+1,i) - Vk(j-1,i))/2.0 
               END DO
           END DO
       END SUBROUTINE electric_field

       SUBROUTINE dados(x1,y1,Vk,Ex,Ey)
!        Imprime os dados em arquivos
           USE globais
           IMPLICIT NONE
           INTEGER :: i,j,k
           REAL,DIMENSION(L,L) :: Vk,Ex,Ey,x1,y1
           CHARACTER(28),DIMENSION(2) :: a
           CHARACTER(27),DIMENSION(2) :: b

           a(1) = 'potencial_laplace_jacob1.dat'
           a(2) = 'potencial_laplace_gauss1.dat'
           b(1) = 'electric_laplace_jacob1.dat'
           b(2) = 'electric_laplace_gauss1.dat'

           DO k=1,2 
           OPEN(10,FILE=a(k),STATUS='unknown',ACTION='write')
               DO i=1,L
                   DO j=1,L
                       WRITE(10,*) x1(j,i),y1(j,i),Vk(j,i)
                   END DO
                   WRITE(10,*) ''
               END DO
           CLOSE(10)

           OPEN(20,FILE=b(k),STATUS='unknown',ACTION='write')
               DO i=1,L
                   DO j=1,L
                       WRITE(20,*) x1(j,i),y1(j,i),Ex(j,i),Ey(j,i)
                   END DO
                   WRITE(20,*) ''
               END DO
           CLOSE(20)
           END DO 
       END SUBROUTINE dados