Programming Simpsons/Trapezium/Newton Raphson method 2

[Fragrancuu]

Hi does anyone know a good method of coding for the solutions of integrals of Simpsons/Trapezoidal functions or the Newton Raphson method?

thanks

 

[mikrom]

http://www.tek-tips.com/viewthread.cfm?qid=1570169&page=2

 

[Fragrancuu]

many thanks, this would be very useful.

 

[mikrom]

Ok, you have the first two metohds and the third, i.e. Newton–Raphson method could be your homework.

 

[Fragrancuu]

I have started a code and need to have each integration method (trapezium and simpsons) in a seperate subroutine, with a module set up up so that global variable can be used for each method. My goal is to end with the integrated equation as the end function. I have had a look at the above link and its good but a little confusing.
Here's my coding

Program Integration

Module integration_methods
!This module contains global data for the integration between the initial set of limits

Implicit none
Save
Double Precision x,f(x)
Integer N,h
real x,f(x),a,b

! a,b=limits
real intent(in) :: a,b

!N,Number of Intervals
real intent(in) :: N

! define limits
a==0
b==1

!the below formula is the general formula to be used by each method
f(x)=x**4+9.31*x**2+x-0.1423

!Temporory values for Trapezoidal function
integer :: i
real :: q
q=0
do i= 1, n-1

its half done but I'm finding it hard to translate the maths into code. I don't know to define the trapeziodal composite function in loops, can anyone help? I believe it would be from 1 to n-1 but I can think of a good formula for input. By the way sorry if it looks messy at the moment, I don't really know what I'm doing.

 

[mikrom]

Hi Fragrancuu,

To integrate your own function, with the program which you will find on the link given above, define the function in the source of the module user_functions.

I can understand that maybe the program is a little bit complicated for a beginner, but I have written it so, because I'm trying to do things universally and I'm trying to use the newer Fortran features.

The code you posted is full of errors and will not compile in any Fortran compiler.
I'm sorry to say, that we cannot help you here with basic understanding of numerical methods and programming. It's on you to study these topics from literature and/or in school.
I wish you good luck!

 

[Fragrancuu]

Hi mikrom, I have re-studied your code and its brilliant however I wanted to ask for your help on one particular aspect of your 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's_rule">http://en.wikipedia.org/wiki/Simpson's_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
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
implicit none
real, intent(in) :: y
f2 = 1/(y-(0.55*y**2 + 0.0165*y + 0.00012375))
end function f2
end module user_functions

program integration
use integration_methods
use user_functions

implicit none
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, 1000)
write (*,*) "Simpson's rule = ", integral

! 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
end program integration

could you help me recode the trapezium rule as a subroutine rathar than a function? I'm just a little slightly confused with the call and subroutine commands. I have tried to look them up on many different websites on the net but I don't fully understand them. Much appreciated.

 

[mikrom]

Hi Fragrancuu,
Ok, I done it for you although it's trivial and you could do it self:

[COLOR=#a020f0]module[/color] integration_methods
  [COLOR=#0000ff]! This module contains the integration methods  [/color]
[COLOR=#a020f0]contains[/color]
  [COLOR=#a020f0]subroutine[/color] trapezoid_rule(f, a, b, n, integral)
[COLOR=#2e8b57][b]    real[/b][/color] :: f
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) ::  a, b
    [COLOR=#2e8b57][b]integer[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: n
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]out[/b][/color]) ::  integral
    [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]  
    integral [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 subroutine[/color] trapezoid_rule  
[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=#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] :: integral

  [COLOR=#0000ff]! integrating he function f1[/color]
  [COLOR=#a020f0]call[/color] trapezoid_rule(f1, [COLOR=#ff00ff]0.0[/color], [COLOR=#ff00ff]1.0[/color], [COLOR=#ff00ff]1000[/color], integral)
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Trapezoid rule = '[/color], integral

  [COLOR=#0000ff]! integrating the function f2[/color]
  [COLOR=#a020f0]call[/color] trapezoid_rule(f2, [COLOR=#ff00ff]0.005[/color], [COLOR=#ff00ff]0.1[/color], [COLOR=#ff00ff]1000[/color], integral)
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Trapezoid rule = '[/color], integral
[COLOR=#a020f0]end program[/color] integration

Compilation and running

$ g95 fragrancuu.f95 -o fragrancuu

$ fragrancuu
 Trapezoid rule =  0.33333343
 Trapezoid rule =  3.1267667

 

[Fragrancuu]

thank you very much.