# Need fortran help! Trapazoid riemann sum!

1. Oct 5, 2011

### mattmac.nuke

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. Oct 6, 2011

### 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.

3. Oct 6, 2011

### mattmac.nuke

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
4. Oct 6, 2011

### 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.

5. Oct 6, 2011

### mattmac.nuke

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?   '

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)? '
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: '
WRITE(*,'(A,\)') ' Enter B: '
WRITE(*,'(A,\)') ' Enter C: '
WRITE(*,*) '            '
DO
WRITE(*,'(A,\)') ' Enter minimum of integration range: '
WRITE(*,'(A,\)') ' Enter maximum of integration range: '
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): '
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: '
WRITE(*,'(A,\)') ' Enter B: '
WRITE(*,'(A,\)') ' Enter C: '
WRITE(*,*) '            '
DO
WRITE(*,'(A,\)') ' Enter minimum of integration range: '
WRITE(*,'(A,\)') ' Enter maximum of integration range: '
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): '
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.

6. Oct 6, 2011

### 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)? '
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)? '
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.