P131.F90

Demonstrate various Root-Finding techniques


!
! =====> Program - P123.F90
!
        PROGRAM DifDmo
!  Program to demonstrate numerical solutions
!  to ordinary, first order differential equations.

        INTEGER    I, N
        REAL       X, Xf, h,        &
                   Y0, Y_Eul, Y_RK4

        PRINT *, 'This is Program >> P123 = Differential Equations'
!
!     Tell program where data for  READ *  is coming from
      OPEN(UNIT=5, FILE='P123.DAT')      ! UNIT=5 is the default input
!
        PRINT *, 'First order ODEs   F = -y*y  for example'
        PRINT *, 'Enter final (high) x value: '
        READ *, Xf
        Print *, Xf

        PRINT *, 'Number of steps to take: '
        READ *, N
        Print *, N

!       Step size, initial X, Y (given)
        h = Xf/N
        X = 0.0
        Y0 = 1

        PRINT *, 'Step size: ', h
        PRINT 101, 'Step', 'x value', 'Euler', 'RK-4'

!       Compare Euler and Runge-Kutta 4'th order
!       at each step.
        Y_Eul = Y0
        Y_RK4 = Y0
L1:     DO I = 1, N
          Y_Eul = Euler(X, Y_Eul, h)
          Y_RK4 = RngKt4(X, Y_RK4, h)

          X = X + h

          PRINT 103, I, X, Y_Eul, Y_RK4
        END DO L1

        STOP
  101   FORMAT(//1X, A5, 6A12)
  103   FORMAT(1X, I5, 3F12.8)
        END


        FUNCTION F(x, y)
!       The function to be integrated.
!       I.e. y' = dy/dx = F(x, y)

        F = x     ! Just to keep FTN90 quiet
        F = -y*y
        RETURN
        END


! ALGORITHMS

        FUNCTION Euler(x, y, h)
!       Uses F(x, y) to advance the solution
!       from x to x+h.

        REAL            x, y, h

        Euler = y + h*F(x, y)

        RETURN
        END



        FUNCTION RngKt4(x, y, h)
!       Uses F(x, y) to advance the solution
!       from x to x+h.

        REAL            x, y, h,  &
                        h_hf,     &
                        x_hf,     &
                        y_hf_1,   &
                        y_hf_2,   &
                        y_1

        h_hf  = h/2
        x_hf  = x + h_hf

!       Form various estimates of the solution
!       at x+h/2 and x+h, then combine them
!       using the formula below.
        y_hf_1  = Euler(x, y, h_hf)
        y_hf_2  = Euler(x_hf, y_hf_1, h_hf)
        y_1     = Euler(x_hf, y_hf_2, h)

        RngKt4 = y + h/3 *                  &
          ((f(x, y) + f(x+h, y_1))/2 +      &
          f(x_hf, y_hf_1) + f(x_hf, y_hf_2))

      RETURN
      END

DATA:
2.0
10

OUTPUT:
Program entered
 This is Program >> P123 = Differential Equations
 First order ODEs   F = -y*y  for example
 Enter final (high) x value: 
   2.0000000
 Number of steps to take: 
 10
 Step size:    0.2000000


  Step     x value       Euler        RK-4
     1  0.20000000  0.80000001  0.85231536
     2  0.40000000  0.67199999  0.74049305
     3  0.60000002  0.58168321  0.65333771
     4  0.80000001  0.51401215  0.58374416
     5  1.00000000  0.46117046  0.52703171
     6  1.20000005  0.41863483  0.48001244
     7  1.40000010  0.38358381  0.44045072
     8  1.60000014  0.35415649  0.40673780
     9  1.80000019  0.32907113  0.37768960
    10  2.00000024  0.30741357  0.35241693
Fortran-90 STOP

Come back to the previous page

Page builder: Charles Boivin

Last modified: 11/07/95