1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Need fortran help! Trapazoid riemann sum!

  1. Oct 5, 2011 #1
    Alright, I cannot seem to get this subroutine to return the correct sums for the trapezoidal rule... Where do I need to fix?

    SUBROUTINE atrap(i)
    USE space_data
    IMPLICIT NONE
    INTEGER :: i, j
    REAL :: f_b1, f_b2, f_x1, f_x2, trap_area
    REAL :: delta_x
    trap_area = 0
    f_b1 = lower
    f_b2 = delta_x + f_b2
    f_x1 = 0
    f_x2 = 0
    trap_area = 0

    DO j = 1, n
    delta_x = (upper - lower)/(2**j)
    f_b2 = delta_x + f_b1
    f_x1 = a*((f_b1)**2) + b*(f_b1) + c
    f_x2 = a*((f_b2)**2) + b*(f_b2) + c
    trap_area = trap_area + .5*delta_x*(f_x1 + f_x2)

    DO i = 2, (2**j)
    f_x1 = a*((f_b1)**2) + b*(f_b1) + c
    f_x2 = a*((f_b2)**2) + b*(f_b2) + c
    f_b1 = f_b2
    f_b2 = delta_x + f_b2
    trap_area = trap_area + .5*delta_x*(f_x1 + f_x2)

    END DO
    WRITE(*,'(A, I2, A, I7, A, A, f16.3)') ' 2^',j,' = ',2**j,' :', ' ', trap_area
    WRITE(*,*) '---------------------------------'
    END DO

    END SUBROUTINE atrap
     
  2. jcsd
  3. Oct 6, 2011 #2

    Mark44

    Staff: Mentor

    In the future, please put a [ code ] tag at the beginning of your code, and a [ /code ] tag after the last line, but leave out the extra spaces.

    I suspect that this code won't even compile, let alone produce any results. You have a variable named lower that hasn't been declared, that I can see. If it actually has been declared outside of this routine, then it would be better to pass it as a parameter to your routine.

    Another problem is that many of your variables are declared but never initialized. Two examples are delta_x and fb_2. A variable that isn't initialized will have a garbage value in it.
     
  4. Oct 6, 2011 #3
    All the variables have been initialized, I'm sorry I forgot to include the code for the module that declares some of the variables that the main program and my subroutines use. This code compiles fine (when used with the main program of course) but the results it produces are off by some strange number. I just figured out how to fix it last night, but I still don't fully understand why the fix works:
    Code (Text):

    SUBROUTINE atrap(i)
      USE space_data
      IMPLICIT NONE
      INTEGER :: i, j
      REAL :: f_b1, f_b2, f_x1, f_x2, trap_area
      REAL :: delta_x  
      DO j = 1, n
        delta_x = (upper - lower)/(2**j)
        trap_area = 0
        DO i = 1, (2**j)
          f_b1 = lower + delta_x*(i-1)
          f_b2 = lower + delta_x*i  
          f_x1 = a*((f_b1)**2) + b*(f_b1) + c  
          f_x2 = a*((f_b2)**2) + b*(f_b2) + c  
          trap_area = trap_area + .5*delta_x*(f_x1 + f_x2)
        END DO
        WRITE(*,'(A, I2, A, I7, A, A, f16.3)') ' 2^',j,' = ',2**j,' |', ' ', trap_area
        WRITE(*,*) '---------------------------------'
      END DO  
      WRITE(*,*) '                  '
    END SUBROUTINE atrap
     
    The fix is in f_b1 and f_b2 and what numbers they initialize at and refresh at.
     
    Last edited by a moderator: Oct 6, 2011
  5. Oct 6, 2011 #4

    Mark44

    Staff: Mentor

    mattmac,
    I took the liberty of reformatting your code to 1) make it easier to understand, and 2) so that the long WRITE statements don't extend so far to the right.

    Your code is difficult to follow, in part because it uses variables that are initialized elsewhere (in MAIN?) and because your routine is passed a parameter (i) that is ignored.

    The variables that are initialized elsewhere are upper, lower, a, b, and c. Without seeing what these values are, it's impossible for someone reading your code to be able to determine the values it produces.

    Your subroutine is passed a value i, but doesn't do anything with it, so why is this being passed?

    Also, it would be better to write a FUNCTION that calculates and returns function values, than to hard-code this into your subroutine, as you have done. Apparently you are using the trapezoid rule to calculate the area between the graph of f(x) = ax^2 + bx + c and the x-axis.
     
  6. Oct 6, 2011 #5
    I apologize for my lack of formatting, I'm not particularly skilled yet at using the formatting statements in Fortran, it's even worse in the main code with all the write statements that had to be added.

    I definitely should have specified that the variables are initialized from input by the user, it was late when I posted and I was in a tremendous rush. The variables are declared in a module before the main program, and initialized in a subroutine that is called by the main program. I made subroutines for different processes in the program because I needed to refer to them at different times, under different conditions and I couldn't figure out a simpler way to write it. I admit that I have absolutely over complicated the program.

    yeah, 'i' was originally a variable declared in the main program, which was where the computations were being made before I wrote a subroutine. 'i' is pretty much a dummy variable that I just forgot to get rid of.... unfortunately, I have many variables like that in the program, in my haste to fix and compile.

    I actually hadn't thought of that, I've written functions for other calculations before, but I was so fixated on the idea that I needed a subroutine to perform both the computation and write statements with the computation that I didn't even figure on building a function for it.

    I probably should have posted the project instructions too, but yes basically the goal was to use 2 different riemann sum methods to calculate the area under a quadratic bounded below by the x-axis.

    Here is the entire code:
    Code (Text):

    MODULE space_data
    REAL :: a, b, c, upper, lower
    INTEGER :: n
    END MODULE

    program riemann
    USE space_data
    IMPLICIT NONE
    INTEGER :: i, l
    CHARACTER :: choice, m, t, q, yn
    REAL :: mid, sum, trap_area

            WRITE(*,*) '          COP 2271 - Project 3'
            WRITE(*,*) '    Integration Approximation Methods       '
            WRITE(*,*) '-----------------------------------------'
            WRITE(*,*) '                                '
    a = 0
    b = 0
    c = 0
    n = 0
    DO
            WRITE(*,*) 'Choose a method:                '
            WRITE(*,*) '     M) Riemann "Middle" Rule   '
            WRITE(*,*) '     T) Riemann "Trapezoidal" Rule'
            WRITE(*,*) '     Q) Quit                        '
            WRITE(*,'(A,\)') 'Choice?   '
            READ(*,*) choice

    IF (choice == 'Q' .OR. choice == 'q') THEN
    WRITE(*,*) '            '
        STOP
    END IF
    IF (choice /= 'M' .AND. choice /= 'm' .AND. choice /= 'T' .AND. choice /= 't' .AND. choice /= 'Q' .AND. choice /= 'q') THEN
                                WRITE(*,*) '===================================='
                                WRITE(*,*) ' Invalid input. Must be M, T, or Q.'
                                WRITE(*,*) '===================================='
                                WRITE(*,*) ' '
                                END IF

    IF (choice == 'M' .OR. choice == 'm' .OR. choice == 'T' .OR. choice == 't') THEN
    IF (n > 0) THEN
      DO
            WRITE(*,*) '            '
            WRITE(*,'(A,\)') ' Would you like to reenter A, B, C, and the interval over which to integrate (y/n)? '
            READ(*,*) yn
            IF ( yn /= 'Y' .AND. yn /= 'y' .AND. yn /= 'N' .AND. yn /= 'n') THEN
                    WRITE(*,*) '                                    '
                    WRITE(*,*) '==================================='
                    WRITE(*,*) ' Invalid input. Must be Y or N. '
                    WRITE(*,*) '==================================='
                    ELSE
                      EXIT
                    END IF
                 END DO
                SELECT CASE(yn)  
                CASE('Y','y')
                CALL choi_mod(choice, m, t, q) 
                CASE('N','n')
                IF (choice == 'm' .OR. choice == 'M') THEN   
                        CALL rite(a, b, c, lower, upper)    
                        CALL midsum(upper, lower, a, b, c, i, n)
                        ELSE IF (choice == 't' .OR. choice == 'T') THEN
                        WRITE(*,*) 'Using Riemann "Trapazoidal" Rule to evaluate'
                        WRITE(*,*) 'y = ',a,'* x^2 +',b,'* x +', c
                        WRITE(*,*) 'from [',lower,',',upper,'].'
                        WRITE(*,*) '                            '
                        WRITE(*,*) '# of intervals  |  Approximate area'
                        WRITE(*,*) '-----------------------------------'
                        CALL atrap(i)
                ELSE
                        STOP
                        END IF
                END SELECT
     
           
        ELSE IF ( a == 0. .AND. b == 0. .AND. c == 0.) THEN
                    CALL choi(choice, m, t, q)                 
        END IF
        END IF
    END DO
    end program riemann

    SUBROUTINE midsum(upper, lower, a, b, c, i, n)
    IMPLICIT NONE
    REAL, INTENT(IN) :: upper, lower, a, b, c
    INTEGER :: i
    INTEGER, INTENT(IN) :: n
    REAL :: f_x, x_par
    REAL :: delta_x, sum
    INTEGER :: j
    DO j = 1, n
    sum = 0
    x_par = 0
    f_x = 0
        delta_x = (upper - lower)/(2**j)       
         x_par = lower + .5*delta_x
         f_x = a*((x_par)**2) + b*(x_par) + c
         sum = delta_x*f_x
    DO i = 2, (2**j)
        x_par = delta_x + x_par
        f_x = a*((x_par)**2) + b*(x_par) + c
        sum = sum + (delta_x*f_x)
    END DO
                WRITE(*,'(A, I2, A, I7, A, A, f16.3)') ' 2^',j,' = ',2**j,' |', ' ', sum
                WRITE(*,*) '---------------------------------'
    END DO
                WRITE(*,*) '                    '
    END SUBROUTINE midsum

    SUBROUTINE rite(a, b, c, lower, upper)
    IMPLICIT NONE
    REAL, INTENT(IN) :: a, b, c, lower, upper
                    WRITE(*,*) 'Using Riemann "Middle" Rule to evaluate'
                    WRITE(*,*) 'y =',a,'* x^2 +',b,'* x +', c
                    WRITE(*,*) 'from [',lower,',',upper,'].'
                    WRITE(*,*) '                            '
                    WRITE(*,*) '# of intervals | Approximate area'
                    WRITE(*,*) '---------------------------------'
    END SUBROUTINE RITE

    SUBROUTINE data()
    USE space_data
    IMPLICIT NONE
                    WRITE(*,*) '                            '
                    WRITE(*,*) 'Enter A, B, and C for the quadratic formula to integrate in the form of:'
                    WRITE(*,*) 'y = Ax^2 + Bx + C'
                    WRITE(*,*) '                    '
                    WRITE(*,'(A,\)') ' Enter A: '
                    READ(*,*) a
                    WRITE(*,'(A,\)') ' Enter B: '
                    READ(*,*) b
                    WRITE(*,'(A,\)') ' Enter C: '
                    READ(*,*) c
                    WRITE(*,*) '            '
    DO
                    WRITE(*,'(A,\)') ' Enter minimum of integration range: '
                    READ (*,*) lower
                    WRITE(*,'(A,\)') ' Enter maximum of integration range: '
                    READ (*,*) upper
                    WRITE(*,*) '                                          '
                IF (lower < upper) THEN
                    EXIT
                ELSE
                    WRITE(*,*) '==================================='
                    WRITE(*,*) ' Invalid input. Max must be > Min. '
                    WRITE(*,*) '==================================='
                    WRITE(*,*) '                                   '
                END IF
                END DO
                DO
                    WRITE(*,*) 'Enter number of approximations to examine'
                    WRITE(*,*) '(ex., entering 3 calculates the integral three times using'
                    WRITE(*,'(A,\)') ' 2^1, 2^2 and 2^3 intervals respectively): '
                    READ(*,*) n
                    IF (n <= 0.) THEN
                    WRITE(*,*) '==================================='
                    WRITE(*,*) ' Invalid input. Number must be > 0. '
                    WRITE(*,*) '==================================='
                    ELSE
                      EXIT
               END IF
               END DO
                               
    END SUBROUTINE data

    SUBROUTINE choi(choice, m, t, q)
    USE space_data
    IMPLICIT NONE
    CHARACTER :: choice, m, t, q
    INTEGER :: i
    SELECT CASE (choice)
                            CASE('M','m')
                                CALL data()
                                CALL rite(a, b, c, lower, upper)
                                CALL midsum(upper, lower, a, b, c, i, n)       
                            CASE('T','t')
                                CALL data()
                                WRITE(*,*) 'Using Riemann "Trapazoidal" Rule to evaluate'
                                WRITE(*,*) 'y = ',a,'* x^2 +',b,'* x +', c
                                WRITE(*,*) 'from [',lower,',',upper,'].'
                                WRITE(*,*) '# of intervals | Approximate area'
                                WRITE(*,*) '---------------------------------'
                                CALL atrap(i)
                            CASE('Q','q')
                            WRITE(*,*) '                    '
                                STOP
    END SELECT
    END SUBROUTINE choi
     
    SUBROUTINE atrap(i)
    USE space_data
    IMPLICIT NONE
    INTEGER :: i, j
    REAL :: f_b1, f_b2, f_x1, f_x2, trap_area
    REAL :: delta_x  
    DO j = 1, n
        delta_x = (upper - lower)/(2**j)
        trap_area = 0
    DO i = 1, (2**j)
        f_b1 = lower + delta_x*(i-1)
        f_b2 = lower + delta_x*i  
        f_x1 = a*((f_b1)**2) + b*(f_b1) + c  
        f_x2 = a*((f_b2)**2) + b*(f_b2) + c  
        trap_area = trap_area + .5*delta_x*(f_x1 + f_x2)

    END DO
                WRITE(*,'(A, I2, A, I7, A, A, f16.3)') ' 2^',j,' = ',2**j,' |', ' ', trap_area
                WRITE(*,*) '---------------------------------'
    END DO  
                WRITE(*,*) '                    '
    END SUBROUTINE atrap





    SUBROUTINE choi_mod(choice, m, t, q)
    USE space_data
    IMPLICIT NONE
    CHARACTER :: choice, m, t, q
    INTEGER :: i
    SELECT CASE (choice)
                            CASE('M','m')
                                CALL data_mod()
                                CALL rite(a, b, c, lower, upper)
                                CALL midsum(upper, lower, a, b, c, i, n)       
                            CASE('T','t')
                                CALL data_mod()
                                WRITE(*,*) 'Using Riemann "Trapazoidal" Rule to evaluate'
                                WRITE(*,*) 'y = ',a,'* x^2 +',b,'* x +', c
                                WRITE(*,*) 'from [',lower,',',upper,'].'
                                WRITE(*,*) '# of intervals | Approximate area'
                                WRITE(*,*) '---------------------------------'
                                CALL atrap(i)
                            CASE('Q','q')
                                                            WRITE(*,*) '                    '
                                STOP
    END SELECT
    END SUBROUTINE choi_mod

    SUBROUTINE data_mod()
    USE space_data
    IMPLICIT NONE
                    WRITE(*,*) '                    '
                    WRITE(*,'(A,\)') ' Enter A: '
                    READ(*,*) a
                    WRITE(*,'(A,\)') ' Enter B: '
                    READ(*,*) b
                    WRITE(*,'(A,\)') ' Enter C: '
                    READ(*,*) c
                    WRITE(*,*) '            '
    DO
                    WRITE(*,'(A,\)') ' Enter minimum of integration range: '
                    READ (*,*) lower
                    WRITE(*,'(A,\)') ' Enter maximum of integration range: '
                    READ (*,*) upper
                    WRITE(*,*) '                                          '
                IF (lower < upper) THEN
                    EXIT
                ELSE
                    WRITE(*,*) '==================================='
                    WRITE(*,*) ' Invalid input. Max must be > Min. '
                    WRITE(*,*) '==================================='
                    WRITE(*,*) '                                   '
                END IF
                END DO
                DO
                    WRITE(*,*) 'Enter number of approximations to examine'
                    WRITE(*,*) '(ex., entering 3 calculates the integral three times using'
                    WRITE(*,'(A,\)') ' 2^1, 2^2 and 2^3 intervals respectively): '
                    READ(*,*) n
                    IF (n <= 0.) THEN
                                    WRITE(*,*) '==================================='
                    WRITE(*,*) ' Invalid input. Number must be > 0. '
                    WRITE(*,*) '==================================='
                    ELSE
                      EXIT
               END IF
               END DO
                               
    END SUBROUTINE data_mod
     
    NOTE: If you see a bunch of extra write statements, among other things, the only reason those are there is to match the output of the executable they give us to compare to. if everything doesn't match, they deduct massive points.
     
  7. Oct 6, 2011 #6

    Mark44

    Staff: Mentor

    The formatting I was talking about has nothing to do with the FORMAT statement. It has to do with how you line program elements up in your code.

    Your formatting is so random and scattered across the screen that I have to scroll off to the right to see what all is in a particular line. The random formatting also takes a lot of work to try to understand what's going on. If I were your instructor, I would definitely take off points for it.

    For example, here is a section of your code in the PROGRAM block:
    Code (Text):
    DO
            WRITE(*,*) '            '
            WRITE(*,'(A,\)') ' Would you like to reenter A, B, C, and the interval over which to integrate (y/n)? '
            READ(*,*) yn
            IF ( yn /= 'Y' .AND. yn /= 'y' .AND. yn /= 'N' .AND. yn /= 'n') THEN
                    WRITE(*,*) '                                    '
                    WRITE(*,*) '==================================='
                    WRITE(*,*) ' Invalid input. Must be Y or N. '
                    WRITE(*,*) '==================================='
                    ELSE
                      EXIT
                    END IF
                 END DO
                SELECT CASE(yn)  
                CASE('Y','y')
                CALL choi_mod(choice, m, t, q) 
                CASE('N','n')
                IF (choice == 'm' .OR. choice == 'M') THEN   
                        CALL rite(a, b, c, lower, upper)    
                        CALL midsum(upper, lower, a, b, c, i, n)
                        ELSE IF (choice == 't' .OR. choice == 'T') THEN
                        WRITE(*,*) 'Using Riemann "Trapazoidal" Rule to evaluate'
                        WRITE(*,*) 'y = ',a,'* x^2 +',b,'* x +', c
                        WRITE(*,*) 'from [',lower,',',upper,'].'
                        WRITE(*,*) '                            '
                        WRITE(*,*) '# of intervals  |  Approximate area'
                        WRITE(*,*) '-----------------------------------'
                        CALL atrap(i)
                ELSE
                        STOP
                        END IF
                END SELECT

     
    Here is how I would format the same block of code.
    Code (Text):

      DO
        WRITE(*,*) '            '
        WRITE(*,'(A,\)') ' Would you like to reenter A, B, C, and the interval over which to integrate (y/n)? '
        READ(*,*) yn
        IF ( yn /= 'Y' .AND. yn /= 'y' .AND. yn /= 'N' .AND. yn /= 'n') THEN
          WRITE(*,*) '                              '
          WRITE(*,*) '==================================='
          WRITE(*,*) ' Invalid input. Must be Y or N. '
          WRITE(*,*) '==================================='
        ELSE
           EXIT
        END IF
      END DO

      SELECT CASE(yn)  
        CASE('Y','y')
          CALL choi_mod(choice, m, t, q)   

        CASE('N','n')
          IF (choice == 'm' .OR. choice == 'M') THEN     
            CALL rite(a, b, c, lower, upper)    
            CALL midsum(upper, lower, a, b, c, i, n)
          ELSE IF (choice == 't' .OR. choice == 'T') THEN
            WRITE(*,*) 'Using Riemann "Trapazoidal" Rule to evaluate'
            WRITE(*,*) 'y = ',a,'* x^2 +',b,'* x +', c
            WRITE(*,*) 'from [',lower,',',upper,'].'
            WRITE(*,*) '                            '
            WRITE(*,*) '# of intervals  |  Approximate area'
            WRITE(*,*) '-----------------------------------'
            CALL atrap(i)
          ELSE
            STOP
          END IF
       END SELECT
    The basic idea is that statements that execute sequentially are aligned. Statements in which the flow of execution is not sequential (DO loops, IF..ELSE blocks, SELECT blocks, etc.) have their bodies indented.

    Back to your problems. Your atrap routine should have reasonable parameters passed to it, with a call looking like this: CALL atrap(upper, lower, n)
    This is similar to but not identical to what you are doing in your midsum routine.

    If you have a FUNCTION to return the values of your quadratic, the values of the coefficients a, b, and c can be hard-coded their, so that they don't need to be passed to atrap or midsum.

    For debugging pick an easy function, say y = x2, between x = 0 and x = 1. Calculate a trapezoid approximation to the integral, by hand, using 2 subintervals. Then use your program to see if you get the same result.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Need fortran help! Trapazoid riemann sum!
  1. Help with fortran (Replies: 1)

  2. Fortran Help! (Replies: 4)

  3. Fortran help (Replies: 16)

Loading...