INTELLIGENT WORK FORUMS
FOR COMPUTER PROFESSIONALS

Are you a
Computer / IT professional?
Join Tek-Tips Forums!
• Talk With Other Members
• Be Notified Of Responses
• Keyword Search
Favorite Forums
• Automated Signatures
• 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.

# 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

Earlier I posted a code with some classical integration rules in this thread

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

### 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.

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:

• Talk To Other Members
• Notification Of Responses To Questions
• Favorite Forums One Click Access
• Keyword Search Of All Posts, And More...

Register now while it's still free!