Successive Over Relaxation Method in FORTRAN

In summary: END DOIn summary, the programmer is trying to find the error in the code. He has tried it a hundred times and is stumped. He asks for help from a FORTRAN expert.
  • #1
Aun Muhammad
14
0
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

  • for forum.jpg
    for forum.jpg
    34.2 KB · Views: 1,414
  • Forum.txt
    5.3 KB · Views: 715
Technology news on Phys.org
  • #2
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
 
  • Like
Likes Aun Muhammad
  • #3
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)
 
  • Like
Likes Aun Muhammad
  • #4
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))
 
  • Like
Likes Aun Muhammad
  • #5
I've made the necessary changes.
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

  • getImage.jpg
    getImage.jpg
    15.3 KB · Views: 524
  • forumgraph.jpg
    forumgraph.jpg
    13.4 KB · Views: 503
  • #6
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.
 
  • Like
Likes Aun Muhammad
  • #7
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.
 
  • #8
Aun Muhammad said:
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.
 
  • Like
Likes Aun Muhammad
  • #9
Thank you for your reply.
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

  • for forum - Copy.jpg
    for forum - Copy.jpg
    15.9 KB · Views: 546
  • #10
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).
 

FAQ: Successive Over Relaxation Method in FORTRAN

1. What is the Successive Over Relaxation Method in FORTRAN?

The Successive Over Relaxation (SOR) Method in FORTRAN is a numerical method used to solve linear systems of equations. It is an iterative method that uses relaxation to improve the convergence rate of the solution.

2. How does the SOR Method work?

The SOR Method works by first dividing the system of equations into two parts: a diagonal part and a non-diagonal part. Then, it uses a relaxation factor to update the solution iteratively. The relaxation factor determines the amount of correction applied to the previous solution to generate the updated solution.

3. What are the advantages of using the SOR Method?

The SOR Method has several advantages, including its ability to converge faster than other iterative methods, its versatility in solving a wide range of linear systems, and its ease of implementation in FORTRAN. It is also a stable method, meaning it will always produce a solution as long as the relaxation factor is chosen correctly.

4. How do you choose the optimal relaxation factor for the SOR Method?

The optimal relaxation factor for the SOR Method depends on the system of equations being solved. It can be determined experimentally by trying different values and choosing the one that results in the fastest convergence. Alternatively, there are also theoretical methods for estimating the optimal relaxation factor, such as the successive substitution method.

5. Can the SOR Method be used to solve non-linear systems of equations?

No, the SOR Method is only suitable for solving linear systems of equations. Non-linear systems require different numerical methods, such as the Newton-Raphson method or the bisection method.

Back
Top