# Successive Over Relaxation Method in FORTRAN

• Fortran
Hey everyone.
I've attached a code I have written in FORTRAN implementing the SOR Method for a 2 D Laplace Equation.
Subroutines have been created, grid has been created and initial and BC's applied appropriately.
While computing the error of each iteration, I am unable to achieve the desired graph. In the graphs I obtained, the error after a few iterations, instantly dives to a very low value and remains there. However, the graph must display a lowest value of error for an 'optimum value' of omega.
I'll be grateful if somebody with FORTRAN experience could look into this and kindly point out where the error is being made. I've tried it a hundred times now and can't find where the error lies. Thanks.

#### Attachments

• 37.2 KB Views: 911
• 5.3 KB Views: 503

DrClaude
Mentor
Next time, please post your code directly, not as an attachment, and use code tags:
Fortran:
PROGRAM ELLIPTIC_SOLVER
USE NRTYPE; USE NRUTIL, ONLY:NRERROR
IMPLICIT NONE

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!*************            VARIABLES AND PARAMETERS            *************!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

INTEGER(I4B)::i,j,k,l
INTEGER(I4B),PARAMETER::N=21,MAXITRS=20,maxw=501
REAL(DP)::DX,DY,w,ERROR
REAL(DP), PARAMETER::deltaw=0.002
REAL(DP),DIMENSION(:,:),ALLOCATABLE::X,Y,S,U0,UOLD,UNEW,R
ALLOCATE(X(N,N),Y(N,N),S(N,N),U0(N,N),UOLD(N,N),UNEW(N,N),R(N,N))

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!*************                   GRID SIZING                  *************!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

CALL GRID(X,Y,DX,DY,N)
WRITE(*,*)'The Matrix X(i,j) is given as:'
DO i=LBOUND(X,1),UBOUND(X,1)
WRITE(*,*) (X(i,j), j=LBOUND(X,2),UBOUND(X,2))
END DO

WRITE(*,*)'The Matrix Y(i,j) is given as:'
DO i=LBOUND(Y,1),UBOUND(Y,1)
WRITE(*,*) (Y(i,j), j=LBOUND(Y,2),UBOUND(Y,2))
END DO

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!*************               ANALYTICAL SOLUTION              *************!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

DO j=1,N
DO i=1,N
U0(i,j)=(X(i,j))**2-(Y(i,j))**2
END DO
END DO

WRITE(*,*)'The Analytical Solution of U0(i,j) is given as:'
DO i=LBOUND(U0,1),UBOUND(U0,1)
WRITE(*,*)(U0(i,j), j=LBOUND(U0,2),UBOUND(U0,2))
END DO

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!*************         INITIAL AND BOUNDARY CONDITIONS        *************!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

DO i=1,N
S(i,1)=U0(i,1)              !Boundary Condition from Analytical Solution
S(i,N)=U0(i,N)              !Boundary Condition from Analytical Solution
END DO

DO j=1,N
S(1,j)=U0(1,j)              !Boundary Condition from Analytical Solution
S(N,j)=U0(N,j)              !Boundary Condition from Analytical Solution
END DO

DO i=2,N-1
DO j=2,N-1
S(i,j)=0.0                               !Initial-Value from Problem
END DO
END DO

DO i=1,N
DO j=1,N
UOLD(i,j)=S(i,j)           !Equating Values from S(i,j) to UOLD(i,j)
END DO
END DO

WRITE(*,*)'The Initial Value Matrix UOLD(i,j) is given as:'
DO i=LBOUND(UOLD,1),UBOUND(UOLD,1)
WRITE(*,*) (UOLD(i,j), j=LBOUND(UOLD,2),UBOUND(UOLD,2))
END DO

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!*************               NUMERICAL SOLUTION               *************!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

DO i=1,N                                !Boundary Conditions Remain Constant
UNEW(i,1) = UOLD(i,1)
UNEW(i,N) = UOLD(i,N)
END DO

DO j=1,N
UNEW(1,j) = UOLD(1,j)
UNEW(N,j) = UOLD(N,j)
END DO

UNEW(i,j)=UOLD(i,j)
ERROR=0.0
w=1.0
DO l=1,maxw
DO k=1,MAXITRS
DO j=2,N-1
DO i=2,N-1
CALL PSOR(UOLD,UNEW)
CALL ERROR_CALCULATION(UOLD,UNEW,ERROR)
CALL RR(UNEW,U0,R)

UOLD(i,j)=UNEW(i,j)
END DO
END DO
WRITE(*,*)'The Initial Value Matrix R(i,j) is given as:'
DO i=LBOUND(R,1),UBOUND(R,1)
WRITE(*,*) (R(i,j), j=LBOUND(R,2),UBOUND(R,2))
END DO
END DO
WRITE(*,*)w,ERROR
ERROR=0.0
w=w+deltaw
END DO

CONTAINS

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!*************                   SUBROUTINES                  *************!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE GRID(X,Y,DX,DY,N)
USE NRTYPE; USE NRUTIL, ONLY:NRERROR
IMPLICIT NONE
INTEGER(I4B)::i,j
REAL(DP), DIMENSION(N,N)::X,Y
INTEGER(I4B),INTENT(IN)::N
REAL(DP),INTENT(OUT)::DX,DY

DX=1.0/(FLOAT(N)-1.0)
DY=1.0/(FLOAT(N)-1.0)

DO j=1,N
DO i=1,N
X(i,j)=(i-1)*DX
Y(i,j)=(j-1)*DY
END DO
END DO
END SUBROUTINE GRID

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE PSOR(UOLD,UNEW)
USE NRTYPE; USE NRUTIL, ONLY:NRERROR
IMPLICIT NONE
REAL(DP)::SUM1
REAL(DP),DIMENSION(N,N),INTENT(IN)::UOLD
REAL(DP),DIMENSION(N,N),INTENT(OUT)::UNEW

SUM1=0.25*(UOLD(i+1,j)+UNEW(i-1,j)+UOLD(i,j+1)+UNEW(i,j-1)-(4.0*UOLD(i,j)))
UNEW(i,j)=UOLD(i,j)+(w*SUM1)
END SUBROUTINE PSOR

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE ERROR_CALCULATION(UOLD,UNEW,ERROR)
USE NRTYPE; USE NRUTIL, ONLY:NRERROR
IMPLICIT NONE
REAL(DP),DIMENSION(N,N),INTENT(IN)::UOLD,UNEW
REAL(DP),INTENT(OUT)::ERROR

ERROR=ERROR+ABS(UNEW(i,j)-UOLD(i,j))
END SUBROUTINE ERROR_CALCULATION

SUBROUTINE RR(UNEW,U0,R)
USE NRTYPE; USE NRUTIL, ONLY:NRERROR
IMPLICIT NONE
REAL(DP),DIMENSION(N,N),INTENT(IN)::U0,UNEW
REAL(DP),DIMENSION(N,N),INTENT(OUT)::R
R(i,j)=U0(i,j)-UNEW(i,j)
END SUBROUTINE RR

END PROGRAM ELLIPTIC_SOLVER

DrClaude
Mentor
Not related to your problem, but the following code is extremely inefficient:
Fortran:
DO i=2,N-1
DO j=2,N-1
S(i,j)=0.0                               !Initial-Value from Problem
END DO
END DO

DO i=1,N
DO j=1,N
UOLD(i,j)=S(i,j)           !Equating Values from S(i,j) to UOLD(i,j)
END DO
END DO
2D arrays are stored in column-major order in Fortran, so the leftmost index such change first. And since you are using Fortran 90, you should simply write
Fortran:
S = 0.0 !Initial-Value from Problem

DO i=1,N
S(i,1)=U0(i,1)              !Boundary Condition from Analytical Solution
S(i,N)=U0(i,N)              !Boundary Condition from Analytical Solution
END DO

DO j=1,N
S(1,j)=U0(1,j)              !Boundary Condition from Analytical Solution
S(N,j)=U0(N,j)              !Boundary Condition from Analytical Solution
END DO

UOLD = S           !Equating Values from S(i,j) to UOLD(i,j)

DrClaude
Mentor
There is a problem with this here:
Fortran:
SUBROUTINE ERROR_CALCULATION(UOLD,UNEW,ERROR)
USE NRTYPE; USE NRUTIL, ONLY:NRERROR
IMPLICIT NONE
REAL(DP),DIMENSION(N,N),INTENT(IN)::UOLD,UNEW
REAL(DP),INTENT(OUT)::ERROR

ERROR=ERROR+ABS(UNEW(i,j)-UOLD(i,j))
END SUBROUTINE ERROR_CALCULATION
ERROR is declared with an intent of out, which means that its value on entry is discarded. Therefore it is executable statement there is equivalent to
Fortran:
    ERROR = ABS(UNEW(i,j)-UOLD(i,j))

However, the results are still the same.
I'm plotting 'omega' versus 'logError'. I'm unable to find wopt(optimum) from the graph.
The blue-lined graph below shows my results.
The red-lined graph is the desired output.

#### Attachments

• 29.1 KB Views: 339
• 15.5 KB Views: 334
DrClaude
Mentor
I can't comment much because I am not familiar with the method you are using. But looking at your graph, the values are in the range 104 to 10-12. That's 16 orders of magnitude, which is basically equal to the machine epsilon for a double precision real.

This might be a case where you can't do any better because you are running out of decimal places.

What you are suggesting seems plausible.
The results however must display a general trend similar to the desired result. The Error must decrease, reach a minima and then begin to increase again. Sadly that is not true in my case.

Hey everyone.
I've attached a code I have written in FORTRAN implementing the SOR Method for a 2 D Laplace Equation.
Subroutines have been created, grid has been created and initial and BC's applied appropriately.
While computing the error of each iteration, I am unable to achieve the desired graph. In the graphs I obtained, the error after a few iterations, instantly dives to a very low value and remains there. However, the graph must display a lowest value of error for an 'optimum value' of omega.
I'll be grateful if somebody with FORTRAN experience could look into this and kindly point out where the error is being made. I've tried it a hundred times now and can't find where the error lies. Thanks.
Hey, I see the following the problems:
1) Your code isn't computing the error of each iterations, as you state. Instead you are summing up the errors of all iterations. That is, ERROR = Error after one iteration + Error after two iterations + ... . In order to solve this you need to put "Error =0.0d0" after the line DO k=1,MAXITRS . Also don't use "Error = 0.0", since "0.0" is in single-precision, so a lot of precision is lost.
2) According to your formula you should divide by (N-2)*(N-2), but I cannot see that it is done in the code.

I've changed the value from single to double precision as you stated. But could you please reiterate what the formula attached below actually means?
In my understanding,
Omega, or w, is being varied from 1 to 2 with a step-size of 0.002. This makes 501 iterations.
For every value of w, I'm carrying out 20 iterations for the calculation of UNEW.
Within these 20 iterations, i and j are varied from 2 to N-1 for the whole grid.
I'm trying to calculate the value of error for every value of w. Therefore, the sum of errors for the 20 iterations has to be summed up right?

#### Attachments

• 18.6 KB Views: 370
Hey, sorry for the delay of my reply. No, in my understanding the formula isn't saying that.
Written out it is ##\mathrm{Error} = \frac{\sum_{ij}|U_{i,j}-S_{i,j}|}{(N-2)*(N-2))}##, where ##U## is the exact solution and ##V## is the approximate one.
In other words, for a fixed number of iterations, the error is computed as the average error over the considered mesh in ##x## and ##y##.
Typically, you should therefore compute the error after in your case 20 iterations. If you want to check convergence with respect to the number of iterations, you keep #w# fixed and compare the errors for let's say 19 and 20 iterations. However, in none of the cases you sum up the the errors.
If you do so, the total error will be dominated by the error coming from the first iterations (where errors are large).