INTELLIGENT WORK FORUMS
FOR COMPUTER PROFESSIONALS

Log In

Come Join Us!

Are you a
Computer / IT professional?
Join Tek-Tips Forums!
  • Talk With Other Members
  • Be Notified Of Responses
    To Your Posts
  • Keyword Search
  • One-Click Access To Your
    Favorite Forums
  • Automated Signatures
    On Your Posts
  • Best Of All, It's Free!

*Tek-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

Jobs

integral of sin(x)

integral of sin(x)

(OP)
Using Monte Carlos Method, how to make a program to calculate the integral points of sin(x) from 0 to pi and how to calculate the errors?

RE: integral of sin(x)

(OP)
What is wrong and what is missing in this code?

PROGRAM area_sin_monte_carlo
USE nrtype
USE nrutil
USE ran_state, ONLY: ran_seed
USE nr, ONLY: ran1
IMPLICIT none
INTEGER, PARAMETER :: u_file = 1, angle_i = 0.0, angle_f=pi, p= 5000
REAL, PARAMETER :: v_angle = (angle_f - angle_i)/REAL(p)
INTEGER ::n, k
REAL :: x, y, area_mc, area, erro, harvest, angle= angulo_i
OPEN(unit= u_file, file= "p_sen.dat", action = "write")
k = 0
DO n = 1, p
call ran1(harvest)
x = harvest*PI
y = sin(x)
WRITE(u_file, *) x, y
k = k + 1
angle= angle + v_angle
END DO
area_mc = ?????????
! area of circle by traditional method
area= -cos(pi) + cos(0.0)
erro=abs((area - area_mc)/(area))*100 ! the erro
WRITE(*,*) p
WRITE(*,*)area_mc
WRITE(*,*) area
WRITE(*,*)'Erro(%): ', erro
END PROGRAM area_sin_monte_carlo

RE: integral of sin(x)

If you try compiling it, it will tell you what is wrong. If you can't work it out, show us what the compiler is moaning about and someone may be able to tell you how to fix it.

Once you have built it, examine what you are getting against what you are expecting. That will possibly tell you what is missing. Again, tell us what you are getting and what you are expecting and someone may be able to tell you how to fix it.

RE: integral of sin(x)

(OP)
I can't transforma the information in the following link you sugested me --> Link in a fortran code. The second link is in java and I don't know nothing about it.
My code calculate just the slope and not a set of points as I was waiting, by Monte Carlos Method.
IN the program that I wrote before there was a little wront, not it is ok.

PROGRAM area_sin_monte_carlo
USE nrtype
USE nrutil
USE ran_state, ONLY: ran_seed
USE nr, ONLY: ran1
IMPLICIT none
INTEGER, PARAMETER :: u_file = 1, angle_f=pi, p= 5000
REAL, PARAMETER :: angle_i = 0.0, v_angle = (angle_f - angle_i)/REAL(p)
INTEGER ::n, k
REAL :: x, y, area_mc, area, erro, harvest, angle

OPEN(unit= u_file, file= "p_sen.dat", action = "write")
k = 0
angle= angle_i
DO n = 1, p
call ran1(harvest)
x = harvest*PI
y = sin(x)
WRITE(u_file, *) x, y
k = k + 1
angle= angle + v_angle
END DO
CLOSE(unit=u_file)
! area_mc = ?????????
! area of circle by the traditional method
area= -cos(pi) + cos(0.0)
erro=abs((area - area_mc)/(area))*100 ! the erro

WRITE(*,*) p
WRITE(*,*)area_mc
WRITE(*,*) area
WRITE(*,*)'Erro(%): ', erro
END PROGRAM area_sin_monte_carlo

RE: integral of sin(x)

Quote (j0zn)


I can't transforma the information in the following link you sugested me --> Link in a fortran code.

It says, that for a given big N, you should generate from the interval [a, b] N random points: x1, x2, .., xN

Then the integral of the function f(x) on the interval [a, b] should be approximately given by a formula:
I(f, a, b) = (b-a) / N * ( f(x1) + f(x2) + ... + f(xN) )

RE: integral of sin(x)

I spent a week on a business trip without a Fortran compiler. Now I have the compiler again available smile

Earlier I posted a code with some classical integration rules in this thread
http://www.tek-tips.com/viewthread.cfm?qid=1570169

Now I have added to this code the Monte Carlo method and f3(x) as the function sin(x)

Here is the code

monte_carlo_integration.f95

CODE

module integration_methods
  ! This module contains the integration methods  
contains
  real function trapezoid_rule(f, a, b, n)
    ! Approximate the definite integral of f from a to b by the
    ! composite trapezoidal rule, using N subintervals
    ! see: <a href="http://en.wikipedia.org/wiki/Trapezoid_rule">http://en.wikipedia.org/wiki/Trapezoid_rule</a>
    ! function arguments ---------------------
    ! f: real function
    real :: f
    ! [a, b] : the interval of integration
    real, intent(in) ::  a, b
    ! n : number of subintervals used
    integer, intent(in) :: n
    ! ----------------------------------------
    ! temporary variables
    integer :: k
    real :: s
    s = 0
    do k=1, n-1
      s = s + f(a + (b-a)*k/n)
    end do  
    trapezoid_rule = (b-a) * (f(a)/2 + f(b)/2 + s) / n
  end function trapezoid_rule

  real function simpson_rule(f, a, b, n)
    ! Approximate the definite integral of f from a to b by the
    ! composite Simpson's rule, using N subintervals
    ! see: <a href="http://en.wikipedia.org/wiki/Simpson%27s_rule">http://en.wikipedia.org/wiki/Simpson%27s_rule</a>
    ! function arguments ---------------------
    ! f: real function
    real :: f
    ! [a, b] : the interval of integration
    real, intent(in) ::  a, b
    ! n : number of subintervals used
    integer, intent(in) :: n
    ! ----------------------------------------
    ! temporary variables
    integer :: k
    real :: s
    s = 0
    do k=1, n-1
      s = s + 2**(mod(k,2)+1) * f(a + (b-a)*k/n)
    end do  
    simpson_rule = (b-a) * (f(a) + f(b) + s) / (3*n)
  end function simpson_rule

  real function monte_carlo_integration(f, a, b, n)
    ! Approximate the definite integral of f from a to b by the
    ! Monte Carlo Method
    ! function arguments ---------------------
    ! f: real function
    real :: f
    ! [a, b] : the interval of integration
    real, intent(in) ::  a, b
    ! n : number of random points used
    integer, intent(in) :: n
    ! ----------------------------------------
    ! temporary variables
    integer :: k
    real :: s, x
    s = 0
    do k=1, n
      ! generate random number from [0, 1]
      !x = rand()
      call random_number(x)
      ! compute function value in random point of interval [a, b]      
      s = s + f(a + (b-a)*x)
    end do  
    monte_carlo_integration = (b-a) * s / n
  end function monte_carlo_integration
end module integration_methods

module user_functions
  ! Define here your own functions to be integrated 
contains
  real function f1(x)
    implicit none
    real, intent(in) :: x
    f1 = x*x
  end function f1

  real function f2(y)
    implicit none
    real, intent(in) :: y
    f2 = 1/(y-(0.55*y**2 + 0.0165*y + 0.00012375))
  end function f2

  real function f3(x)
    implicit none
    real, intent(in) :: x
    f3 = sin(x)
  end function f3
  
end module user_functions

program integration
  use integration_methods
  use user_functions
  
  implicit none
  real, parameter :: pi=3.141592653589793d0
  real :: integral

  ! integrating he function f1
  integral=trapezoid_rule(f1, 0.0, 1.0, 1000)
  write (*,*) 'Trapezoid rule     = ', integral
  integral=simpson_rule(f1, 0.0, 1.0, 10000)
  write (*,*) "Simpson's rule     = ", integral
  integral=monte_carlo_integration(f1, 0.0, 1.0, 10000)
  write (*,*) "Monte Carlo method = ", integral
  write (*,*)

  ! integrating the function f2
  integral=trapezoid_rule(f2, 0.005, 0.1, 1000)
  write (*,*) 'Trapezoid rule     = ', integral
  integral=simpson_rule(f2, 0.005, 0.1, 1000)
  write (*,*) "Simpson's rule     = ", integral
  integral=monte_carlo_integration(f2, 0.005, 0.1, 10000)
  write (*,*) "Monte Carlo method = ", integral
  write (*,*)

  ! integrating the function f3
  integral=trapezoid_rule(f3, 0.0, pi, 1000)
  write (*,*) 'Trapezoid rule     = ', integral
  integral=simpson_rule(f3, 0.0, pi, 1000)
  write (*,*) "Simpson's rule     = ", integral
  integral=monte_carlo_integration(f3, 0.0, pi, 100000)
  write (*,*) "Monte Carlo method = ", integral 
end program integration 
Here are the results which seems to be reasonable

CODE

$ gfortran monte_carlo_integration.f95 -o monte_carlo_integration

$ monte_carlo_integration
 Trapezoid rule     =   0.333333433
 Simpson's rule     =   0.333333194
 Monte Carlo method =   0.337485731

 Trapezoid rule     =    3.12676668
 Simpson's rule     =    3.12673712
 Monte Carlo method =    3.14724326

 Trapezoid rule     =    1.99999881
 Simpson's rule     =    2.00000048
 Monte Carlo method =    2.00680566 

Note, that this is the probability method, so you cannot be sure to increase the accuracy of the result simply by increasing the number of attempts.
If you are interested see this thread
http://www.tek-tips.com/viewthread.cfm?qid=1730131

RE: integral of sin(x)

(OP)
Thank you mikrom! From this I will be able to solve my problems.

Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

Red Flag Submitted

Thank you for helping keep Tek-Tips Forums free from inappropriate posts.
The Tek-Tips staff will check this out and take appropriate action.

Reply To This Thread

Posting in the Tek-Tips forums is a member-only feature.

Click Here to join Tek-Tips and talk with other members!

Resources

Close Box

Join Tek-Tips® Today!

Join your peers on the Internet's largest technical computer professional community.
It's easy to join and it's free.

Here's Why Members Love Tek-Tips Forums:

Register now while it's still free!

Already a member? Close this window and log in.

Join Us             Close