I with Fortran 90 : Simpson's rule

  • Context: Fortran 
  • Thread starter Thread starter fluidistic
  • Start date Start date
  • Tags Tags
    Fortran
Click For Summary

Discussion Overview

The discussion revolves around implementing the composite Simpson's rule in Fortran 90 to calculate the integral of the function \( e^{-x} \) over the interval [0,1]. Participants are addressing issues related to the implementation, debugging, and precision of the calculations, as well as comparing results obtained with different numbers of subintervals.

Discussion Character

  • Technical explanation
  • Debate/contested
  • Mathematical reasoning

Main Points Raised

  • One participant describes their implementation of the composite Simpson's rule and notes discrepancies in the expected and calculated precision coefficient.
  • Another participant suggests testing the program with simpler integrals and using print statements for debugging.
  • A participant reports successful testing of the program on multiple platforms but remains uncertain about the source of errors in the subroutines.
  • Concerns are raised about the formulation of the initial conditions in the subroutine, particularly regarding the denominator in the calculation of \( s \) and \( r \).
  • One participant provides an alternative subroutine for Simpson's rule integration, indicating that the precision aspect should be removed.
  • Another participant points out potential issues with mixing data types in the calculation of the precision coefficient \( Q \), suggesting the use of explicit type specifications for constants.
  • Corrections are proposed regarding the calculation of the step size \( del_x \) in relation to the number of points used in the integration.
  • Participants share modified versions of the original program, highlighting changes that lead to improved accuracy in results.
  • Some participants express difficulty in compiling the provided code, citing errors related to the structure and requirements of the Fortran language.

Areas of Agreement / Disagreement

Participants generally agree on the need for precision in calculations and the importance of correctly implementing the Simpson's rule. However, there are multiple competing views regarding the specific implementation details and corrections needed, leading to an unresolved discussion on the best approach.

Contextual Notes

Some participants mention issues with the Fortran compiler and specific requirements for variable declarations. There are also unresolved questions about the handling of precision and the structure of the subroutines.

Who May Find This Useful

This discussion may be useful for individuals interested in numerical methods for integration, particularly those using Fortran for scientific computing, as well as those troubleshooting similar programming issues.

  • #31
Is the program in this case a good one to use for simpson's composite rule because it gives a value with a very small difference when entering an ood number compared to an even number when simpson's composite needs only an even number of intervals
 
Technology news on Phys.org
  • #32
bigev234 said:
Is the program in this case a good one to use for simpson's composite rule because it gives a value with a very small difference when entering an ood number compared to an even number when simpson's composite needs only an even number of intervals

I think I would have to disagree with that. The difference is not small at all.

The value of the integral (integrating exp(-x) from 0 to 1) should be 0.63212 to five digits.

When using n=100 intervals, you get that value; if you use n=101 intervals, you get 0.628448. That's about a half a percent difference, which I would say is a huge error for this integral.


For that matter, you can just replace the function given in the program with 1; that is, calculate the integral of 1 from 0 to 1. Of course it should be give 1, right? But the results are:

n=100 ===> integral=1.
n=101 ===> integral=0.990099


The error formula predicts zero error for this type of integral. But for an odd number of intervals it is very close to 1% error--much too large for the integral of a constant. I would say you need to use an even number of subintervals.
 
  • #33
hellloo..how can i make my compoosite simpson fortran code general.i mean...take the function as user's input and solve according..heres the codde..what needs to be changed?.

PROGRAM SIMPSON
C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994
C To accompany the text:
C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
C This free software is complements of the author.
C
C Algorithm 7.2 (Composite Simpson Rule).
C Section 7.2, Composite Trapezoidal and Simpson's Rule, Page 365
C
INTEGER M
REAL A,B,Srule
CHARACTER*60 ANS*1,DFUN,FUN
EXTERNAL F
10 CALL INPUTS(A,B,M,DFUN,FUN)
CALL SIMPRU(F,A,B,M,Srule)
CALL RESULT(A,B,M,Srule,DFUN,FUN)
WRITE(9,*)'WANT TO TRY ANOTHER INTERVAL ? <Y/N> '
READ (9,'(A)') ANS
IF (ANS.EQ.'Y' .OR. ANS.EQ.'y') GOTO 10
END

REAL FUNCTION F(X)
REAL X
F=X/(1+X*X)
RETURN
END

SUBROUTINE PRINTFUN(DFUN,FUN)
CHARACTER*(*) DFUN,FUN
FUN ='X/(1+X*X)'
DFUN='X/(1+X*X) DX'
RETURN
END

SUBROUTINE SIMPRU(F,A,B,M,Srule)
INTEGER K,M
REAL A,B,H,Sum,SumEven,SumOdd,Srule,X
EXTERNAL F
H=(B-A)/(2*M)
SumEven=0
DO K=1,(M-1)
X=A+H*2*K
SumEven=SumEven+F(X)
ENDDO
SumOdd=0
DO K=1,M
X=A+H*(2*K-1)
SumOdd=SumOdd+F(X)
ENDDO
Sum=H*(F(A)+F(B)+2*SumEven+4*SumOdd)/3
Srule=Sum
RETURN
END

SUBROUTINE XSIMPRU(F,A,B,M,Srule)
C This subroutine uses labeled DO loop(s).
INTEGER K,M
REAL A,B,H,Sum,SumEven,SumOdd,Srule,X
EXTERNAL F
H=(B-A)/(2*M)
SumEven=0
DO 10 K=1,(M-1)
X=A+H*2*K
SumEven=SumEven+F(X)
10 CONTINUE
SumOdd=0
DO 20 K=1,M
X=A+H*(2*K-1)
SumOdd=SumOdd+F(X)
20 CONTINUE
Sum=H*(F(A)+F(B)+2*SumEven+4*SumOdd)/3
Srule=Sum
RETURN
END

SUBROUTINE INPUTS(A,B,M,DFUN,FUN)
INTEGER I,M
REAL A,B
CHARACTER*(*) DFUN,FUN
CALL PRINTFUN(DFUN,FUN)
DO 10 I=1,18
WRITE(9,*)' '
10 CONTINUE
WRITE(9,*)' SIMPSON`S RULE IS USED TO COMPUTE AN APPROXIMATION'
WRITE(9,*)' '
WRITE(9,*)'FOR THE VALUE OF THE DEFINITE INTEGRAL:'
WRITE(9,*)' '
WRITE(9,*)' '
WRITE(9,*)' B'
WRITE(9,*)' /'
WRITE(9,*)' | ',DFUN
WRITE(9,*)' /'
WRITE(9,*)' A'
WRITE(9,*)' '
WRITE(9,*)' '
WRITE(9,*)'ENTER THE LEFT ENDPOINT A = '
READ(9,*) A
WRITE(9,*)' '
WRITE(9,*)'ENTER THE RIGHT ENDPOINT B = '
READ(9,*) B
WRITE(9,*)' '
WRITE(9,*)'THE NUMBER OF SUBINTERVALS USED IS 2*M'
WRITE(9,*)' ENTER THE NUMBER M = '
READ(9,*) M
WRITE(9,*)' '
RETURN
END

SUBROUTINE RESULT(A,B,M,Srule,DFUN,FUN)
INTEGER I,M
REAL A,B,Srule
CHARACTER*(*) DFUN,FUN
CALL PRINTFUN(DFUN,FUN)
DO 10 I=1,18
WRITE(9,*)' '
10 CONTINUE
WRITE(9,*)' ',B
WRITE(9,*)' /'
WRITE(9,*)' |'
WRITE(9,*)Srule,' ~ | ',DFUN
WRITE(9,*)' |'
WRITE(9,*)' /'
WRITE(9,*)' ',A
WRITE(9,*)' '
WRITE(9,*)' '
WRITE(9,*)'AN APPROXIMATE VALUE FOR THE DEFINITE INTEGRAL OF'
WRITE(9,*)' '
WRITE(9,*)' '
WRITE(9,*)'F(X) = ',FUN
WRITE(9,*)' '
WRITE(9,*)' '
WRITE(9,*)'TAKEN OVER [',A,' ,',B,' ] WAS FOUND.'
WRITE(9,*)' '
WRITE(9,*)'SINCE M = ',M,', THERE WERE ',2*M,' SUBINTERVALS.'
WRITE(9,*)' '
WRITE(9,*)'THE SIMPSON RULE APPROXIMATION IS ',Srule
WRITE(9,*)' '
RETURN
END
 

Similar threads

Fortran Errors in Fortran code for solving Laplace's equation in 3D
  • · Replies 4 ·
Replies
4
Views
3K
Fortran Spring Pendulum system in FORTRAN
  • · Replies 8 ·
Replies
8
Views
2K
Fortran Need help with Jacobi relaxation method for Dirichlet boundary conditions
  • · Replies 1 ·
Replies
1
Views
2K
Fortran Fortran jump command for fortran 95
  • · Replies 12 ·
Replies
12
Views
2K
Fortran How do I create a pointer to an element of an array with a variable index?
  • · Replies 20 ·
Replies
20
Views
3K
Fortran Dimension of arrays (RESHAPE) in Fortran 90
  • · Replies 4 ·
Replies
4
Views
2K
Fortran Fortran subroutine solving for the velocity and displacement of a spring
  • · Replies 8 ·
Replies
8
Views
2K
Fortran Getting my method of bisection to output the correct value
  • · Replies 11 ·
Replies
11
Views
2K
Fortran I can't find a small bug in this Fortran code
  • · Replies 2 ·
Replies
2
Views
1K
Fortran Assign 2D Array with Reals in Fortran - Help Needed
  • · Replies 4 ·
Replies
4
Views
2K