Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

Register Log in
  • Congratulations SkipVought on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!

integral of sin(x)

Status
Not open for further replies.
j0zn

j0zn

Programmer
Jul 19, 2013
36
0
0
BR
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?
 
Sort by date Sort by votes
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
 
Upvote 0 Downvote
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.
 
Upvote 0 Downvote
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
 
Upvote 0 Downvote
j0zn said:
I can't transforma the information in the following link you sugested me --> Link in a fortran code.
Click to expand...

It says, that for a given big N, you should generate from the interval [a, b] N random points: x[sub]1[/sub], x[sub]2[/sub], .., x[sub]N[/sub]

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(x[sub]1[/sub]) + f(x[sub]2[/sub]) + ... + f(x[sub]N[/sub]) )
 
Upvote 0 Downvote
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: 
[COLOR=#a020f0]module[/color] integration_methods
  [COLOR=#0000ff]! This module contains the integration methods  [/color]
[COLOR=#a020f0]contains[/color]
[COLOR=#2e8b57][b]  real[/b][/color] [COLOR=#a020f0]function[/color] trapezoid_rule(f, a, b, n)
    [COLOR=#0000ff]! Approximate the definite integral of f from a to b by the[/color]
    [COLOR=#0000ff]! composite trapezoidal rule, using N subintervals[/color]
    [COLOR=#0000ff]! see: <a href="[URL unfurl="true"]http://en.wikipedia.org/wiki/Trapezoid_rule">http://en.wikipedia.org/wiki/Trapezoid_rule</a>[/URL][/color]
    [COLOR=#0000ff]! function arguments ---------------------[/color]
    [COLOR=#0000ff]! f: real function[/color]
[COLOR=#2e8b57][b]    real[/b][/color] :: f
    [COLOR=#0000ff]! [a, b] : the interval of integration[/color]
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) ::  a, b
    [COLOR=#0000ff]! n : number of subintervals used[/color]
    [COLOR=#2e8b57][b]integer[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: n
    [COLOR=#0000ff]! ----------------------------------------[/color]
    [COLOR=#0000ff]! temporary variables[/color]
    [COLOR=#2e8b57][b]integer[/b][/color] :: k
[COLOR=#2e8b57][b]    real[/b][/color] :: s
    s [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]
    [COLOR=#804040][b]do[/b][/color] k[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], n[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]
      s [COLOR=#804040][b]=[/b][/color] s [COLOR=#804040][b]+[/b][/color] f(a [COLOR=#804040][b]+[/b][/color] (b[COLOR=#804040][b]-[/b][/color]a)[COLOR=#804040][b]*[/b][/color]k[COLOR=#804040][b]/[/b][/color]n)
    [COLOR=#804040][b]end do[/b][/color]  
    trapezoid_rule [COLOR=#804040][b]=[/b][/color] (b[COLOR=#804040][b]-[/b][/color]a) [COLOR=#804040][b]*[/b][/color] (f(a)[COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]2[/color] [COLOR=#804040][b]+[/b][/color] f(b)[COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]2[/color] [COLOR=#804040][b]+[/b][/color] s) [COLOR=#804040][b]/[/b][/color] n
  [COLOR=#a020f0]end function[/color] trapezoid_rule

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

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

[COLOR=#a020f0]module[/color] user_functions
  [COLOR=#0000ff]! Define here your own functions to be integrated [/color]
[COLOR=#a020f0]contains[/color]
[COLOR=#2e8b57][b]  real[/b][/color] [COLOR=#a020f0]function[/color] f1(x)
    [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: x
    f1 [COLOR=#804040][b]=[/b][/color] x[COLOR=#804040][b]*[/b][/color]x
  [COLOR=#a020f0]end function[/color] f1

[COLOR=#2e8b57][b]  real[/b][/color] [COLOR=#a020f0]function[/color] f2(y)
    [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: y
    f2 [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1[/color][COLOR=#804040][b]/[/b][/color](y[COLOR=#804040][b]-[/b][/color]([COLOR=#ff00ff]0.55[/color][COLOR=#804040][b]*[/b][/color]y[COLOR=#804040][b]**[/b][/color][COLOR=#ff00ff]2[/color] [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]0.0165[/color][COLOR=#804040][b]*[/b][/color]y [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]0.00012375[/color]))
  [COLOR=#a020f0]end function[/color] f2

[COLOR=#2e8b57][b]  real[/b][/color] [COLOR=#a020f0]function[/color] f3(x)
    [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: x
    f3 [COLOR=#804040][b]=[/b][/color] [COLOR=#008080]sin[/color](x)
  [COLOR=#a020f0]end function[/color] f3
  
[COLOR=#a020f0]end module[/color] user_functions

[COLOR=#a020f0]program[/color] integration
  [COLOR=#a020f0]use[/color] integration_methods
  [COLOR=#a020f0]use[/color] user_functions
  
  [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
[COLOR=#2e8b57][b]  real[/b][/color], [COLOR=#2e8b57][b]parameter[/b][/color] :: pi[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]3.141592653589793d0[/color]
[COLOR=#2e8b57][b]  real[/b][/color] :: integral

  [COLOR=#0000ff]! integrating he function f1[/color]
  integral[COLOR=#804040][b]=[/b][/color]trapezoid_rule(f1, [COLOR=#ff00ff]0.0[/color], [COLOR=#ff00ff]1.0[/color], [COLOR=#ff00ff]1000[/color])
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Trapezoid rule     = '[/color], integral
  integral[COLOR=#804040][b]=[/b][/color]simpson_rule(f1, [COLOR=#ff00ff]0.0[/color], [COLOR=#ff00ff]1.0[/color], [COLOR=#ff00ff]10000[/color])
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]"Simpson's rule     = "[/color], integral
  integral[COLOR=#804040][b]=[/b][/color]monte_carlo_integration(f1, [COLOR=#ff00ff]0.0[/color], [COLOR=#ff00ff]1.0[/color], [COLOR=#ff00ff]10000[/color])
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]"Monte Carlo method = "[/color], integral
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color])

  [COLOR=#0000ff]! integrating the function f2[/color]
  integral[COLOR=#804040][b]=[/b][/color]trapezoid_rule(f2, [COLOR=#ff00ff]0.005[/color], [COLOR=#ff00ff]0.1[/color], [COLOR=#ff00ff]1000[/color])
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Trapezoid rule     = '[/color], integral
  integral[COLOR=#804040][b]=[/b][/color]simpson_rule(f2, [COLOR=#ff00ff]0.005[/color], [COLOR=#ff00ff]0.1[/color], [COLOR=#ff00ff]1000[/color])
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]"Simpson's rule     = "[/color], integral
  integral[COLOR=#804040][b]=[/b][/color]monte_carlo_integration(f2, [COLOR=#ff00ff]0.005[/color], [COLOR=#ff00ff]0.1[/color], [COLOR=#ff00ff]10000[/color])
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]"Monte Carlo method = "[/color], integral
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color])

  [COLOR=#0000ff]! integrating the function f3[/color]
  integral[COLOR=#804040][b]=[/b][/color]trapezoid_rule(f3, [COLOR=#ff00ff]0.0[/color], pi, [COLOR=#ff00ff]1000[/color])
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Trapezoid rule     = '[/color], integral
  integral[COLOR=#804040][b]=[/b][/color]simpson_rule(f3, [COLOR=#ff00ff]0.0[/color], pi, [COLOR=#ff00ff]1000[/color])
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]"Simpson's rule     = "[/color], integral
  integral[COLOR=#804040][b]=[/b][/color]monte_carlo_integration(f3, [COLOR=#ff00ff]0.0[/color], pi, [COLOR=#ff00ff]100000[/color])
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]"Monte Carlo method = "[/color], integral 
[COLOR=#a020f0]end program[/color] 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
 
Upvote 0 Downvote
Thank you mikrom! From this I will be able to solve my problems.
 
Upvote 0 Downvote
Status
Not open for further replies.

Similar threads

Vahid9
  • Locked
  • Question
Preventing the addition of C-style comments in gfortran compilation 1
Replies
2
Views
141
Vahid9
Vahid9
Benali.2
  • Locked
  • Question
Use a txt file (X, Y, Z), with random values for X and Y, to get a new file X and Y are predefined 1
Replies
6
Views
93
Benali.2
Benali.2
sadeghiut
  • Locked
  • Question
integral programming in fortran
Replies
2
Views
29
sadeghiut
sadeghiut
garynewport2
  • Locked
  • Question
I have a file of almost 5000 lines 1
Replies
5
Views
109
xwb
xwb
AlanJT
  • Locked
  • Question
re-use of unit number
Replies
1
Views
57
mikrom
mikrom

Part and Inventory Search

Sponsor

Back
Top