Programs in Fortran 77

Factorial of an integer

===================================================================================
Program to find the factorial of an integer                 
LANGUAGE :: FORTRAN 77  
Compiler :: GNU Fortran
Program By:: G.R. Mohanty , www.numericalmethods.in
===================================================================================
      PROGRAM factorial
	  IMPLICIT NONE
      INTEGER i,n,fac
      PRINT*,'Enter the number'
      READ*,n
      IF (n.EQ.0) THEN
            PRINT*,'0! = 1'
            STOP
      ELSE IF (n.LT.0) THEN
            PRINT*,n,'!= INFINITE'
            STOP
      ELSE
            fac=1
            DO i=1,n
                  fac=fac*i
            END DO
            PRINT *,n,'!=',fac
      END IF
      STOP
      END

---------------OUTPUT 1--------------
 Enter the number
 -1
 -1 != INFINITE
 ---------------OUTPUT 2--------------
 Enter the number
 0
 0! = 1

---------------OUTPUT 3--------------
 Enter the number
 10
 10 !=     3628800

---------------OUTPUT 4--------------
 Enter the number
 5
 5!=   120
 

Adding Matrices

===================================================================================
Program to add two 3X3 matrices                
LANGUAGE :: FORTRAN 77  
Compiler :: GNU Fortran
Program By:: G.R. Mohanty , www.numericalmethods.in
===================================================================================
      PROGRAM matrixAddition
	  IMPLICIT NONE
      REAL a(3,3),b(3,3),c(3,3)
      INTEGER i,j
      PRINT*,'Enter the elements of the first matrix'
      DO i=1,3
            READ*,(a(i,j),j=1,3)
      END DO
      PRINT*,'Enter the elements of the second matrix'
      DO i=1,3
            READ*,(b(i,j),j=1,3)
      END DO
      DO i=1,3
            DO j=1,3
                  c(i,j)=a(i,j)+b(i,j)
            END DO
      END DO
      PRINT*,'Sum of the two matrices is:'
      DO i=1,3
            PRINT*,(C(i,j),j=1,3)
      END DO
      STOP
      END

====================OUTPUT=========================
   Enter the elements of the first matrix
   34 56 8
   5 6 78
   4 5 7
   Enter the elements of the second matrix
   4 5 6
   6 7 8
   7 8 9
   Sum of the two matrices is:
   38.0000000       61.0000000       14.0000000
   11.0000000       13.0000000       86.0000000
   11.0000000       13.0000000       16.0000000

____________________Sample Output___________________
   Enter the elements of the first matrix
   45 56 78
   4.5 5.1 66.75
   1.1 2 50
   Enter the elements of the second matrix
   1 1 1
   2 2 2
   3 3 3
   Sum of the two matrices is:
   46.0000000       57.0000000       79.0000000
   6.50000000       7.09999990       68.7500000
   4.09999990       5.00000000       53.0000000
 

Eigenvalue by Power method

===================================================================================
Program to find the largest eigenvalue of a matrix by power method
LANGUAGE :: FORTRAN 77
Compiler :: GNU Fortran
Program By:: G.R. Mohanty , www.numericalmethods.in
===================================================================================
      program eigenvalue
        implicit none
        integer i,j,n
        parameter (n=3)
        real a(n,n),x0(n),x1(n),lambda, tolerance
        logical approx

        print*,"Enter the ",n,"X",n," matrix"
        do i=1,n
            read*,(a(i,j),j=1,n)
        end do
        print*,"enter the trial eigenvector"
        read*,(x0(i),i=1,n)
        tolerance=0.0001
        x1=0
        
        do while (.not.approx)
            approx=.true.
            x1=0
            
!       multiplying the matrix with the trial vector
            do i=1,n
                do j=1,n
                    x1(i)=x1(i)+a(i,j)*x0(j)
                end do
            end do

            lambda=x1(n)
            print*,lambda


            do i=1,n
                x1(i)=x1(i)/lambda
            end do
            print*,"eigen vector:",(x1(i),i=1,n)
            do i=1,n
                if (abs(x1(i)-x0(i)).gt.tolerance)then
                    approx=.false.
                    exit
                end if
            end do
            do i=1,n
                x0(i)=x1(i)
            end do

        end do
        print*,"the eigenvalue is",lambda
        print*,"eigen vector:",(x1(i),i=1,n)
        stop
        end

====================OUTPUT=========================
 Enter the 3X3 matrix
 1 3 -1
 3 2  4
-1 4 10

 enter the trial eigenvector
0 0 1

 The approximate eigenvalue is:   10.0000000
 The approximate eigenvector is: -0.100000001      0.400000006       1.00000000

 The approximate eigenvalue is:   11.6999998
 The approximate eigenvector is:   8.54701083E-03  0.384615391       1.00000000

 The approximate eigenvalue is:   11.5299149
 The approximate eigenvector is:   1.40845114E-02  0.415863603       1.00000000

 The approximate eigenvalue is:   11.6493702
 The approximate eigenvector is:   2.24626176E-02  0.418390036       1.00000000

 The approximate eigenvalue is:   11.6510973
 The approximate eigenvector is:   2.38288902E-02  0.420919001       1.00000000

 The approximate eigenvalue is:   11.6598473
 The approximate eigenvector is:   2.45788749E-02  0.421388447       1.00000000

 The approximate eigenvalue is:   11.6609745
 The approximate eigenvector is:   2.47615855E-02  0.421621144       1.00000000

 The approximate eigenvalue is:   11.6617231
 The approximate eigenvector is:   2.48355269E-02  0.421681017       1.00000000

 The eigenvalue is:   11.6617231
 The eigen vector:   2.48355269E-02  0.421681017       1.00000000

 

Finding Prime Numbers

===================================================================================
Program to Find Prime Numbers                
LANGUAGE :: FORTRAN 77  
Compiler :: GNU Fortran
Program By:: G.R. Mohanty , www.numericalmethods.in
===================================================================================
      PROGRAM prime
      IMPLICIT NONE
      INTEGER i,j,n,N_div
!     N_div corresponds to the number of divisors of an integer
      PRINT*,'Program to find all prime numbers below N'
      PRINT*,'Enter N'
      READ*, n
      PRINT *,'ALL PRIME NUMBERS BELOW',n,'ARE :'
      DO i=1,n-1
            N_div=0
            DO j=1,i
                  IF(MOD(i,j).EQ.0) N_div=N_div+1
            END DO
            IF (N_div.EQ.2) PRINT*,i
      END DO
      STOP
      END

====================OUTPUT=========================
 Program to find all prime numbers below N
 Enter N
 100
 
 ALL PRIME NUMBERS BELOW         100 ARE :
           2
           3
           5
           7
          11
          13
          17
          19
          23
          29
          31
          37
          41
          43
          47
          53
          59
          61
          67
          71
          73
          79
          83
          89
          97

 

Trapezoidal Rule

===================================================================================
Program to Find Integral by Trapezoidal Rule                
LANGUAGE :: FORTRAN 77  
Compiler :: GNU Fortran
Program By:: G.R. Mohanty , www.numericalmethods.in
===================================================================================
      PROGRAM trapezoidal
      IMPLICIT NONE
      INTEGER j,n
      REAL x,upper,lower,I,summ,h, func
      PRINT*,"INTEGRAL OF A FUNCTION BY TRAPEZOIDAL METHOD"
      PRINT*,('-',j=1,44)
      PRINT*,'Enter the lower limit'
      READ*,lower
      PRINT*,'Enter the upper limit'
      READ*,upper
      PRINT*,'Enter the no of intervals'
      READ*,n
      h=(upper-lower)/n
      x=lower
      summ=0
      WRITE(*,10)
10    FORMAT (8X,"x",18x,"f(x)")
      PRINT*,('-',j=1,32)
      DO j=0,n
          IF ((j.EQ.0).OR.(j.EQ.n)) THEN
             summ=summ+func(x)
          ELSE
             summ=summ+(2*func(x))
          END IF
          PRINT*,x,'  ',func(x)
          x=x+h
      END DO
      I=h*summ/2
      PRINT*,('-',j=1,32)
      PRINT*,'Integral is : ',I
      STOP
      END PROGRAM
!----------------------FUNCTION SUBPROGRAM---------------------
      REAL FUNCTION func(x1)
      func=1/(1+x1)
      RETURN
      END

====================OUTPUT=========================
   INTEGRAL OF A FUNCTION BY TRAPEZOIDAL METHOD
 --------------------------------------------
 Enter the lower limit
0
 Enter the upper limit
1
 Enter the no of intervals
10
        x                  f(x)
 --------------------------------
   0.00000000          1.00000000
  0.100000001         0.909090936
  0.200000003         0.833333313
  0.300000012         0.769230783
  0.400000006         0.714285731
  0.500000000         0.666666687
  0.600000024         0.625000000
  0.700000048         0.588235259
  0.800000072         0.555555522
  0.900000095         0.526315749
   1.00000012         0.499999970
 --------------------------------
 Integral is :   0.693771362
 

Simpsons's 1/3 Rule

===================================================================================
Program to Find Integral by Simpsons's 1/3 Rule                
LANGUAGE :: FORTRAN 77  
Compiler :: GNU Fortran
Program By:: G.R. Mohanty , www.numericalmethods.in
===================================================================================
      PROGRAM simpson1_3rd
      INTEGER j,n
      REAL x,upper,lower,I,summ,h
      PRINT*,'INTEGRAL BY SIMPSON 1/3 RULE'
      PRINT*,('-',j=1,28)
      PRINT*,'Enter the lower limit'
      READ*,lower
      PRINT*,'Enter the upper limit'
      READ*,upper
      PRINT*,'Enter the no of intervals'
      READ*,n
      h=(upper-lower)/n
      x=lower
      summ=0
      PRINT*,('-',j=1,32)
      WRITE(*,10)
10    FORMAT (8X,"x",18x,"f(x)")
      PRINT*,('-',j=1,32)
      DO j=0,n
          IF ((j.EQ.0).OR.(j.EQ.n)) THEN
             summ=summ+func(x)
          ELSE IF (MOD(j,2).NE.0) THEN
             summ=summ+(4*func(x))
          ELSE IF (MOD(j,2).EQ.0) THEN
             summ=summ+(2*func(x))
          END IF
          PRINT*,x,'  ',func(x)
          x=x+h
      END DO
      I=h*summ/3
      PRINT*,('-',j=1,32)
      PRINT*,'Integral is : ',I
      STOP
      END PROGRAM
!-----------------FUNCTION SUBPROGRAM-----------------
      REAL FUNCTION func(x1)
      func=1/(1+x1)
      RETURN
      END

====================OUTPUT=========================
 INTEGRAL BY SIMPSON 1/3 RULE
 ----------------------------
 Enter the lower limit
 0
 Enter the upper limit
 1
 Enter the no of intervals
 10
 --------------------------------
        x                  f(x)
 --------------------------------
   0.00000000          1.00000000
  0.100000001         0.909090936
  0.200000003         0.833333313
  0.300000012         0.769230783
  0.400000006         0.714285731
  0.500000000         0.666666687
  0.600000024         0.625000000
  0.700000048         0.588235259
  0.800000072         0.555555522
  0.900000095         0.526315749
   1.00000012         0.499999970
 --------------------------------
 Integral is :   0.693150222
 

Sine Series

===================================================================================
Program to Find Sum of Sine Series                
LANGUAGE :: FORTRAN 77  
Compiler :: GNU Fortran
Program By:: G.R. Mohanty , www.numericalmethods.in
===================================================================================
      PROGRAM sine
	  IMPLICIT NONE
      INTEGER i,n
      REAL adeg,arad,PI,summ,fac		


      PRINT*,'Enter the angle in degrees'
      READ*,adeg
      PRINT*,'Enter the number of terms up to which you want to sum'
      READ*,n
      arad=adeg*PI/180
      summ=0
      DO i=0,n-1
          summ=summ+((-1)**i)*(arad**((2*i)+1))/fac((2*i)+1)
      END DO
      PRINT 10,adeg,summ
10    FORMAT('Sine(',F5.1,')=',1X,F10.8)
      STOP
      END

! For Cosine series replace line 05 with
! summ=summ+((-1)**i)*(arad**(2*i))/fac(2*i)
!----------FUNCTION to find factorial---------

      FUNCTION fac(n)
      INTEGER i,n
      IF (n.EQ.0) THEN
        fac=1
        RETURN
      ELSE
        fac=1
        DO i=1,n
           fac=fac*i
        END DO
      END IF
      RETURN
      END

====================OUTPUT=========================
 Enter the angle in degrees
 90
 Enter the number of terms up to which you want to sum
 10
 Sine( 90.0)= 1.00000000

 

Forward Difference Table

===================================================================================
Program to Generate forward difference table
LANGUAGE :: FORTRAN 77
Compiler :: GNU Fortran
===================================================================================
      program forwardDtable
        implicit none
        integer i,j,num
        real del(10,10)
        print*,'Enter the size of the data (max 10)'
        read*,num
        print*,'Enter the data'
        print*,'X,Y'
        print*,('-',i=1,8*(num+1))

!       First two columns of the forward difference table contain data
        do i=1,num
            read*,(del(i,j),j=1,2)
        end do

!       Third column onwards contain the forward diffrences
        do j=3,num+1
            do i=1,num-(j-2)
                del(i,j)= del(i+1,j-1)-del(i,j-1)
            end do
        end do

!       printing the table heading---------------------------
        print*,('-',i=1,8*(num+1))
        print*,"The Forward Difference Table:"
        print*,('-',i=1,8*(num+1))
        print 11,'x ','y ',('DEL',j,j=1,num-1)
        print*,('-',i=1,8*(num+1))

!       Printing the table
        do i=1,num
            print 10,(del(i,j),j=1,2),(del(i,j),j=3,num+2-i)
        end do

        print*,('-',i=1,8*(num+1))

10      format(10f8.1)
11      format (2A8,10(4x,a3,i1))
        stop
       end

====================OUTPUT=========================
  Enter the size of the data (max 10)
8
 Enter the data
 X,Y
 ------------------------------------------------------------------------
1,3010
2,3424
3,3802
4,4105
5,4472
6,4771
7,5051
8,5315
 ------------------------------------------------------------------------
 The Forward Difference Table:
 ------------------------------------------------------------------------
      x       y     DEL1    DEL2    DEL3    DEL4    DEL5    DEL6    DEL7
 ------------------------------------------------------------------------
     1.0  3010.0   414.0   -36.0   -39.0   178.0  -449.0   901.0 -1580.0
     2.0  3424.0   378.0   -75.0   139.0  -271.0   452.0  -679.0
     3.0  3802.0   303.0    64.0  -132.0   181.0  -227.0
     4.0  4105.0   367.0   -68.0    49.0   -46.0
     5.0  4472.0   299.0   -19.0     3.0
     6.0  4771.0   280.0   -16.0
     7.0  5051.0   264.0
     8.0  5315.0
 ------------------------------------------------------------------------

 

Newton's Forward Interpolation

===================================================================================
Program to Interpolate using Newtons Forward difference Interpolation Formula for
evenly spaced points
LANGUAGE :: FORTRAN 77
Compiler :: GNU Fortran
===================================================================================
      program newtonFormula
        implicit none
        integer i,j,num,fac
        real del(10,10)
        real x,p,h,y,numerator
        print*,'Enter the size of the data (max 10)'
        read*,num
        print*,'Enter the data'
        print*,'X,Y'
        print*,('-',i=1,8*(num+1))

!       First two columns of the forward difference table contain data
        do i=1,num
            read*,(del(i,j),j=1,2)
        end do

!       Third column onwards contain the forward differences
        do j=3,num+1
            do i=1,num-(j-2)
                del(i,j)= del(i+1,j-1)-del(i,j-1)
                if (i.gt.num-(j-2)) del(i,j)=0
            end do
        end do

!       printing the table heading---------------------------
        print*,('-',i=1,8*(num+1))
        print*,"The Forward Difference Table:"
        print*,('-',i=1,8*(num+1))
        print 11,'x ','y ',('DEL',j,j=1,num-1)
        print*,('-',i=1,8*(num+1))

!       Printing the table
        do i=1,num
            print 10,(del(i,j),j=1,2),(del(i,j),j=3,num+2-i)
        end do

        print*,('-',i=1,8*(num+1))

10      format(10f8.2)
11      format (2A8,10(4x,a3,i1))
!       ------------------------------------------------------

        print*,'Enter x'
        read*,x

        h=del(2,1)-del(1,1)
        p=(x-del(1,1))/h
        y=del(1,2)
        numerator=1

!       Using Newton's Interpolation formula
        do i=1,num-1
            numerator=numerator*(p-i+1)
            y=y+numerator*del(1,i+2)/fac(i)
        end do

        print*,'The value of y(x) at',x,'is:',y

      stop
      end program

      INTEGER FUNCTION fac(n)
      INTEGER i,n
      IF (n.EQ.0) THEN
        fac=1
        RETURN
      ELSE
        fac=1
        DO i=1,n
           fac=fac*i
        END DO
      END IF
      RETURN
      END

====================OUTPUT=========================
 Enter the size of the data (max 10)
7
 Enter the data
 X,Y
 ----------------------------------------------------------------
0,176
1,185
2,194
3,203
4,212
5,220
6,229
 ----------------------------------------------------------------
 The Forward Difference Table:
 ----------------------------------------------------------------
      x       y     DEL1    DEL2    DEL3    DEL4    DEL5    DEL6
 ----------------------------------------------------------------
    0.00  176.00    9.00    0.00    0.00    0.00   -1.00    5.00
    1.00  185.00    9.00    0.00    0.00   -1.00    4.00
    2.00  194.00    9.00    0.00   -1.00    3.00
    3.00  203.00    9.00   -1.00    2.00
    4.00  212.00    8.00    1.00
    5.00  220.00    9.00
    6.00  229.00
 ----------------------------------------------------------------
 Enter x
0.2
 The value of y(x) at  0.200000003     is:   177.672318
 

Lagrange's Interpolation Formula

===================================================================================
Program to Interpolate using Lagrange's Interpolation Formula for
unevenly spaced points
LANGUAGE :: FORTRAN 77
Compiler :: GNU Fortran
By: Gyana Ranjan Mohanty Dt.17.02.2023
===================================================================================

      program lagrange
        implicit none
        integer i,j,n
        real x,y,numerator, denominator,xy(10,2)
        print*,'Enter the number of data points (max 10)'
        read*,n
        print*,'Enter the data'
        print*,'X,Y'
        print*,('----')

        do i=1,n
            read*,(xy(i,j),j=1,2)
        end do

        print*,'Enter the value of x to get y(x)'
        read*,x

        y=0
        do i=1,n
            numerator=1
            denominator=1
            do j=1,n
                if (j.eq.i) cycle
                    numerator = numerator*(x - xy(j,1))
                    denominator = denominator*(xy(i,1)-xy(j,1))
            end do
            y = y + xy(i,2)*numerator/denominator
        end do
        print 10,x,y
10      format('y(',f5.2,')=',f8.3)
        stop
        end

Enter the number of data points (max 10)
4
 Enter the data
 X,Y
 ----
3,168
7,120
9,72
10,63
 Enter the value of x to get y(x)
6
y( 6.00)= 147.000
 

Runge-Kutta Method

===================================================================================
Program to solve differential equations of type y'=f(x,y) using 
Fourth Order Runge-Kutta Method given the value of y at x0
LANGUAGE :: FORTRAN 77
Compiler :: GNU Fortran
By: Gyana Ranjan Mohanty Dt.20.02.2023
===================================================================================

      PROGRAM rkmethod
        IMPLICIT NONE
        REAL k1,k2,k3,k4,h
        REAL x(100),y(100),f
        INTEGER i,n
        PRINT*,'Enter x0'
        READ*,x(1)
        PRINT*,'Enter y(x0)'
        READ*,y(1)
        PRINT*,'Enter step size (h)'
        READ*,h
        PRINT*,'Enter the number of steps'
        READ*,n
        DO i=1,n-1
            k1 = h*f(x(i),y(i))
            k2 = h*f(x(i)+h/2,y(i)+k1/2)
            k3 = h*f(x(i)+h/2,y(i)+k2/2)
            k4 = h*f(x(i)+h,y(i)+k3)
            y(i+1)=y(i)+(k1+2*k2+2*k3+k4)/6
            x(i+1)=x(i)+h
        END DO

        PRINT*,'x  y'
        DO i=1,n
            PRINT*,x(i),y(i)
        END DO
        STOP
      END PROGRAM

      REAL FUNCTION f(x,y)
        REAL x,y
        f=x*y+y*y
        RETURN
      END

---------OUTPUT--------
 Enter x0
0
 Enter y(x0)
1
 Enter step size (h)
0.1
 Enter the number of steps
10
       x            y(x)
 -------------------------
    0.0000         1.0000
    0.1000         1.1169
    0.2000         1.2774
    0.3000         1.5041
    0.4000         1.8389
    0.5000         2.3687
    0.6000         3.3069
    0.7000         5.3539
    0.8000        12.7922
    0.9000       461.3828
 -------------------------