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