Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

I need help with Fortran 90 : Simpson's rule

  1. Apr 29, 2008 #1

    fluidistic

    User Avatar
    Gold Member

    I must calculate [tex]\int_0^1 e^{-x} dx[/tex] using the composite Simpson's rule, i.e. the common Simpson's rule but applied on many intervals between 0 and 1. This is not all : I must divide the interval [0,1] in 100 subintervals and then in 200, to compare the value obtained of the integral. And then I must calculate the coefficient of precision, that is Q=absolute value of [tex]\frac{s-I}{r-I}[/tex], where I is the (correct) value of the integral we are searching, s its approximation by using the Simpson's rule on 100 subintervals and r its approximation by using the S' rule on 200 subintervals. And the answer, Q, is 16. While I find almost 2. I've thought a lot about my error, and I know it's in the subroutine when computing the integral, but can't figure where exactly. The bad thing is that my result of the integral is quite similar to the real one... Should be a little mistake then.
    I forgot to say, my program may be a little hard to understand. When you execute it, first it will ask the bounds of the integral, just put 0,1. Then chose 100 for n and 200 for m and we're done.
    Program Simpson
    implicit none

    Real(8) :: x,s,x_i,x_f,r,q
    Integer :: k,n,m

    Write(*,*)'Enter x_i y x_f'
    Read(*,*)x_i,x_f
    Write(*,*)'Enter n'
    Read(*,*)n
    Write(*,*)'Enter m'
    Read(*,*)m

    Call Simp(x_i,x_f,s,n)
    Write(*,*)'When n=100, s is worth',s

    Call Simp2(x_i,x_f,r,m)
    Write(*,*)'When n=200, s is worth',r


    Q=(s-(1-exp(-1.)))/(r-(1-exp(-1.)))
    Write(*,*)'The quotient of precision is',Q



    Contains
    Subroutine Simp(x_i,x_f,s,n)
    implicit none
    Real(8), intent(in) :: x_i,x_f
    Real(8), Intent(out) :: s
    Real(8) :: dx
    Integer :: n

    dx=(x_f-x_i)/n
    s=(dx/3.)/(f(x_i)+f(x_f))

    Do k=2,n-2,2
    x=x_i+k*dx
    s=s+2*f(x)
    end do

    do k=1,n-1,2
    x=x_i+k*dx
    s=s+4*f(x)
    end do
    s=(dx/3.)*s

    end subroutine

    Subroutine Simp2(x_i,x_f,r,m)
    implicit none
    Real(8), Intent(in) :: x_i,x_f
    Real(8), Intent(out) :: r
    Real(8) :: dx
    Integer :: m

    dx=(x_f-x_i)/m
    r=(dx/3.)/(f(x_i)+f(x_f))

    Do k=2,m-2,2
    x=x_i+k*dx
    r=r+2*f(x)
    end do

    do k=1,m-1,2
    x=x_i+k*dx
    r=r+4*f(x)
    end do
    r=(dx/3.)*r
    end subroutine


    Real(8) Function f(x)
    Real(8), Intent(in) :: x
    f=exp(-x)
    end function
    end program
     
  2. jcsd
  3. Apr 29, 2008 #2
  4. Apr 29, 2008 #3

    fluidistic

    User Avatar
    Gold Member

    All works fine, I tested it under both linux and windows, with and without the double precision of digits and on 2 different machines. Still don't know where my subroutines are wrong with the Simpson's rule.
     
  5. Apr 29, 2008 #4

    alphysicist

    User Avatar
    Homework Helper

    Hi fluidistic,

    In the second executable line of your subroutine simp and simp2 you have:

    Code (Text):

    s=(dx/3.)/(f(x_i)+f(x_f))
     
    I believe this should just be

    Code (Text):

    s=(f(x_i)+f(x_f))
     
    I definitely don't think the endpoint values should be in a denominator; and also since you multiply the entire s value by (dx/3.) at the last statement of the subroutine, having it here also in this statement makes it a factor (dx/3.)^2 for the endpoints.
     
  6. Apr 29, 2008 #5
    here's my simpson's rule integration subroutine, the (prec) has to do with the precision, just remove it.


    subroutine simpson ( a,b,n_points, approx ) ! simpsons rule integration

    real(prec), intent(in) :: a,b
    integer, intent(in) :: n_points
    integer :: k,i
    real(prec), intent(out) :: approx
    real(prec) :: del_x, height, area, t,x
    real(prec) :: n_points_1
    n_points_1 = real(n_points,DP)

    del_x = abs(b-a)/n_points_1


    approx = 0
    t = a

    do k=1, n_points-1

    x = t + del_x/two
    area = del_x/six * ( four*functio_n ( x ) + functio_n(t) + functio_n(t+del_x) )
    t = t + del_x

    approx = approx + area

    end do

    end subroutine
     
  7. Apr 29, 2008 #6

    fluidistic

    User Avatar
    Gold Member

    Oh... you're right alphysicist! But even changed that I don't get a quotient of 16, but 1 now.
     
  8. Apr 29, 2008 #7

    fluidistic

    User Avatar
    Gold Member

    Thanks for your answer ice109. But I don't understand all your program. I'm quite new with fortran. For example, when you first define " n_points_1 = real(n_points,DP)" real() means something special? And what is DP?
    I tried without the precision once again in my program and it's the same, I don't get the 16 value I should.
     
  9. Apr 29, 2008 #8
    use mine and declare functio_n to be what ever you want.
     
  10. Apr 29, 2008 #9

    alphysicist

    User Avatar
    Homework Helper

    fluidistic,

    You're losing your precision when you calculate Q. You have:

    Q=(s-(1-exp(-1.)))/(r-(1-exp(-1.)))

    which mixes type 8 real variables ( s and r) with an integer constant (1) and a default type real constant ( -1.). You need to specify that your numeric constants are type 8 reals by using:

    Q=(s-(1._8-exp(-1._8)))/(r-(1._8-exp(-1._8)))
     
  11. Apr 29, 2008 #10

    fluidistic

    User Avatar
    Gold Member

    Alphysicist I can't believe it! It works now... I would never have thought about that. Thanks a lot.
    Ice109, if you're still interested in offering me your program, you should give it to me entirely otherwise I would have to define not only functio_n but everything, and as I said, I don't understand this part : "n_points_1 = real(n_points,DP)". Thanks to all of you that pay some attention and time to my problem.
     
  12. Apr 29, 2008 #11
    change that to simply real(n_points). that's the only other thing that's not portable.
     
  13. Apr 29, 2008 #12

    alphysicist

    User Avatar
    Homework Helper

    Hi ice109,

    I think there is a problem in one line of your subroutine. Where you have the line:

    del_x = abs(b-a)/n_points_1

    I think this needs to be

    del_x = abs(b-a)/(n_points_1-1.)

    (where the 1. needs to match the precision of the variables.) For example, if you were choosing 3 points, you would divide the interval (b-a) by 2 to find the stepsize.


    When I have functio_n(x)=1, for example, for the integral from 0 to 1 with n=50, I get a result 0.9800000000000005 with a 2% error; if I make that change in the denominator, I get the result 1.0000000000000007.
     
  14. Apr 30, 2008 #13

    fluidistic

    User Avatar
    Gold Member

    Nicely seen Alphysicist. You're right, using the program of ice109 I got a close result, but changing it in your way I obtain exactly what I obtain using my program.
    I post his correct program here :

    Program Simpsons
    implicit none

    real, Parameter :: a=0,b=1
    integer :: n_points
    integer :: k,i
    real :: approx
    real :: del_x, height, area, t,x
    real :: n_points_1

    Write(*,*)'Enter number of points'
    Read(*,*)n_points

    n_points_1 = real(n_points)

    del_x = abs(b-a)/(n_points_1-1.)


    approx = 0
    t = a

    do k=1, n_points-1

    x = t + del_x/2
    area = del_x/6 * ( 4*functio_n ( x ) + functio_n(t) + functio_n(t+del_x) )
    t = t + del_x

    approx = approx + area

    end do



    Write(*,*)approx


    Contains
    Real Function functio_n(x)
    Real, Intent(in) :: x
    functio_n=exp(-x)
    end function

    end program
     
  15. May 8, 2008 #14
    can't get this to work

    I tried running this as it's written using the Lahey compilier and i get a load of errors.
    What am I doing wrong?
     
  16. May 8, 2008 #15

    fluidistic

    User Avatar
    Gold Member

    Which program doesn't work for you?
    If it's the first or the last, then it's strange. Maybe try to compile under gfortran, the program I'm using for it.
     
  17. May 8, 2008 #16
    the first program...i get errors along the lines of.."Contains" requires stop or return....also errors that state n & m must be used with "intent". In total its about 17 errors.
    I don't have the compiler on my home pc or i'd post the actual error messages for you.
     
  18. May 8, 2008 #17

    fluidistic

    User Avatar
    Gold Member

    Sadly I'm not a specialist, so I cannot help you... I can try post my program on an attached file.
    It's almost the same I posted here at first (but in Spanish and corrected). Download it and tell us what happens.
     

    Attached Files:

  19. May 8, 2008 #18
    thanks i'll try it...i don't have class again till monday so i'll post my results then.
     
  20. May 9, 2008 #19

    alphysicist

    User Avatar
    Homework Helper

    Hi surfernj,

    I think the Lahey compiler is strict on what it allows. I don't have access to a Lahey compiler but you might try the following things to correct the errors messages you mentioned:

    put the statement

    RETURN

    before every END FUNCTION and END SUBROUTINE (3 places)


    Also put:

    STOP

    before the line with CONTAINS





    In the subroutines, you'll find two lines that say:

    Integer::m

    and

    Integer::n


    Change these to:

    Integer,Intent(in)::m

    and

    Integer,Intent(in)::n




    One of my compilers did not accept kind=8 reals; I did a find and replace to convert every 8 to a 2 and it compiled fine.
     
  21. May 9, 2008 #20
    ok its running now however i'm trying to find the area from the function f(x)= e^(x/2)^2 from 0 to 2
    and i keep getting 2.9.... using 20 or 50 divisions when the answer should be 4.7....

    heres the code:
    11 IMPLICIT NONE
    12 INTEGER, PARAMETER:: DOUBLE=kind(1.0D0)
    13 REAL(KIND=DOUBLE) :: e,x,s,x_i,x_f,r,q
    14 Integer :: k,n,m
    15 !
    16 ! Open Printer
    17 OPEN (UNIT= 2, FILE='iii1')
    18 !
    19 ! Get value of e
    20 e = EXP(1.0)
    21 !
    22 !
    23 Write(UNIT=*,FMT=*)'Enter x_i y x_f'
    24 Read(UNIT=*,FMT=*)x_i,x_f
    25 Write(UNIT=*,FMT=*)'Enter n'
    26 Read(UNIT=*,FMT=*)n
    27 Write(UNIT=*,FMT=*)'Enter m'
    28 Read(UNIT=*,FMT=*)m
    29 !
    30 Call Simp(x_i,x_f,s,n)
    31 Write(UNIT=*,FMT=*)'When n=20, s is worth',s
    32 !
    33 Call Simp2(x_i,x_f,r,m)
    34 Write(UNIT=*,FMT=*)'When n=50, s is worth',r
    35 !
    36 !
    37 Q=(s-(1-exp(-1.)))/(r-(1-exp(-1.)))
    38 Write(UNIT=*,FMT=*)'The quotient of precision is',Q
    39 !
    40 !
    41 stop
    42 Contains
    43 Subroutine Simp(x_i,x_f,s,n)
    44 implicit none
    45 Real(KIND=DOUBLE), intent(in) :: x_i,x_f
    46 Real(KIND=DOUBLE), Intent(out) :: s
    47 Real(KIND=DOUBLE) :: dx
    48 INTEGER, INTENT(IN) :: n
    49 !
    50 !
    51 dx=(x_f-x_i)/n
    52 s=(f(x_i)+f(x_f))
    53 !
    54 !
    55 Do k=2,n-2,2
    56 x=x_i+k*dx
    57 s=s+2*f(x)
    58 end do
    59 !
    60 !
    61 Do k=1,n-1,2
    62 x=x_i+k*dx
    63 s=s+4*f(x)
    64 END DO
    65 !
    66 s=(dx/3.)*s
    67 !
    68 return
    69 END Subroutine Simp
    70 !
    71 !
    72 Subroutine Simp2(x_i,x_f,r,m)
    73 implicit none
    74 Real(KIND=DOUBLE), Intent(in) :: x_i,x_f
    75 Real(KIND=DOUBLE), Intent(out) :: r
    76 Real(KIND=DOUBLE) :: dx
    77 INTEGER, INTENT(IN) :: m
    78 !
    79 !
    80 dx=(x_f-x_i)/m
    81 r=(f(x_i)+f(x_f))
    82 !
    83 !
    84 DO k=2,m-2,2
    85 x=x_i+k*dx
    86 r=r+2*f(x)
    87 END Do
    88 !
    89 !
    90 DO k=1,m-1,2
    91 x=x_i+k*dx
    92 r=r+4*f(x)
    93 END DO
    94 r=(dx/3.)*r
    95 return
    96 end Subroutine Simp2
    97 !
    98 ! Defining the internal function
    99 FUNCTION F(x)
    100 REAL (KIND = double) :: F
    101 REAL (KIND = double), INTENT(IN) :: X
    102 ! Defining the equation
    103 f= e**((x/2)**2)
    104 return
    105 END function F
     
    Last edited: May 9, 2008
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?



Similar Discussions: I need help with Fortran 90 : Simpson's rule
Loading...