Need fortran help Trapazoid riemann sum

In summary: I'm not very skilled at it yet...I'm aware that the variable i is not being used, but that's because I haven't figured out what to do with it yet. I'm a bit confused by why it's being passed as a parameter, I didn't write it that way, it's just the way it is. I thought it was necessary for the do loops, but it seems to work without it.I don't know how to write a function yet, so I'll have to look into that. Thank you very much for your input, I greatly appreciate it.In summary, the code provided is a subroutine meant to calculate the area under a curve using the trapezoidal rule. It uses variables that are initialized
  • #1
mattmac.nuke
22
0
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
 
Physics news on Phys.org
  • #2
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.

mattmac.nuke said:
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.
 
  • #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: 
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:
  • #4
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
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: 
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.
 
  • #6
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: 
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: 
  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.
 

1. What is the Need fortran help Trapazoid riemann sum method?

The Need fortran help Trapazoid riemann sum method is a numerical integration technique used to approximate the area under a curve. It involves dividing the area into trapezoids and summing their areas to estimate the total area.

2. How does the Trapazoid riemann sum method work?

The Trapazoid riemann sum method works by dividing the interval of the function into smaller subintervals and approximating the curve within each subinterval with a linear function. The area under each linear function is then calculated and summed to estimate the total area under the curve.

3. What is the difference between the Trapazoid riemann sum and other integration methods?

The Trapazoid riemann sum method is a specific type of numerical integration method that uses trapezoids to approximate the area under a curve. Other methods, such as Simpson's rule or the midpoint rule, use different shapes to approximate the area.

4. When should I use the Need fortran help Trapazoid riemann sum method?

The Trapazoid riemann sum method is best used when the function being integrated is difficult or impossible to integrate analytically. It is a useful tool for numerical integration in these cases, but may not always provide the most accurate approximation.

5. How can I improve the accuracy of the Trapazoid riemann sum method?

To improve the accuracy of the Trapazoid riemann sum method, you can increase the number of subintervals used to approximate the curve. This will create smaller trapezoids and provide a more accurate estimation of the area under the curve.

Similar threads

Comp Sci Can‘t call a subroutine name in a subroutine in Fortran90
    • Engineering and Comp Sci Homework Help
Replies
7
Views
1K
How to read a column of data into Fortran without arrays? (Fortran 77)
    • Engineering and Comp Sci Homework Help
Replies
7
Views
1K
Comp Sci Multiplication of two 2x2 matrices in Fortran
    • Engineering and Comp Sci Homework Help
Replies
2
Views
5K
Comp Sci 5x5 Arrays - Sum & Difference - (Fortran)
    • Engineering and Comp Sci Homework Help
Replies
13
Views
2K
I Original definition of Riemann Integral and Darboux Sums
    • Topology and Analysis
Replies
14
Views
2K
Fortran Errors in Fortran code for solving Laplace's equation in 3D
    • Programming and Computer Science
Replies
4
Views
602
Comp Sci Need help with Fortran homework -- Error: Invalid form of PROGRAM statement
    • Engineering and Comp Sci Homework Help
Replies
4
Views
2K
Comp Sci FORTRAN Help: Functions and Subroutines
    • Engineering and Comp Sci Homework Help
Replies
1
Views
3K
Is the function v(r) = (e^(-mr))/r a solution to the Yukawa PDE equation?
    • Calculus and Beyond Homework Help
Replies
1
Views
619
Need Help with Fortran compiler errors
    • Engineering and Comp Sci Homework Help
Replies
12
Views
3K

Hot Threads

Recent Insights

Back
Top