Comp Sci Need fortran help Trapazoid riemann sum

[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

 

[Mark44]

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.

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

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.

 

[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:

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.

 

[Mark44]

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.

 

[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:

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.

 

[Mark44]

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:

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.

  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 = x², 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.