Integral calculation using trapezoid or simpsons rule

Question

Basically the assignment says "find the solution to f(x) = 0....". where f(x)= integrate (1/sqrt(2*pi))exp(-t^2/2) dt from x=0 to x - 0.45. Sorry I could not find out how to put the integral in here.

my prof. gave us a little hint how to start.

start x
do i = 1,2,3...
    fp = 1/sqrt(2*pi)exp(-x^2/2)
    f  = use trap,or simpson's rule to find the integration than subtract 0.45
    x = x - (f/fp)
end do 

here is what I did

program main
implicit none

integer :: n, k, i
double precision :: h, a, fp, f, x1, x2, pi, blub, integ, e, dx, j, m

a = 0
n = 25
x1 = 0.5
pi = 4.0d0 * atan(1.0d0)

do i = 1, n
    e = exp((-x1)**2.0d0/2.0d0)
    print*, e

    fp = (1.0d0/sqrt(2.0d0 * pi))* e

    !start of the trapezoid
    dx = (x1-a)/n  !find dx by subtracting starting point from endpoint and divide by number of steps taken
    m = (blub(a) + blub(x1))/2 !find the mean value of the integral


    j = 0
        do k=1, n-1
            h = i
            j = j + blub(h) !calculate the each step of the integral from one to n and add all together. 
        end do

    integ = dx*(m+j) !result of the trapezoid method

    print*, "integ: ", integ
    f = integ - 0.45
    a = x1
    x1 = a - (f/fp)
    print*, "x1: ",x1
    print*, "a: ",a
    print*, "f: ",f
    print*, "fp: ", fp
    print*, (x1-a)/n
end do

stop
end 

double precision function blub (x)
double precision :: x

blub(x) = (1.0d0/sqrt(2.0d0 * (4.0d0 * atan(1.0d0))))*exp((-x)**2.0d0/2.0d0)

return
end

I did all the printing to find out where my mistake might be. I realized that my dx is becoming small based on the x1 being smaller than n. Because of that my integ becomes a number with e^-310 which has almost no effect on the 0.45 and therefore my f stays the same and that messes x1 up...

It would be awesome if somebody could explain to me how I can fix that mistake.

Edit: I'm sure I made a mistake with the trapezoid part, I just don't know where. The f is supposed to change with every loop but dx being so small prevents that.

Answer 1

One of the errors is this. I don't claim it is the only one:

double precision function blub (x)
double precision :: x

blub(x) = (1.0d0/sqrt(2.0d0 * (4.0d0 * atan(1.0d0))))*exp((-x)**2.0d0/2.0d0)

return
end

If you compile it with gfortran -Wall -std=f95 -c fun.f90 you get

fun.f90:5.76:

lub(x) = (1.0d0/sqrt(2.0d0 * (4.0d0 * atan(1.0d0))))*exp((-x)**2.0d0/2.0d0)
                                                                           1
Warning: Obsolescent feature: Statement function at (1)
fun.f90:1.30:

double precision function blub (x)
                              1
Warning: Return value of function 'blub' at (1) not set

This could alert you something strange is going on. You have actually created a new statement function (an obsolete Fortran feature) inside the original function instead of giving it the right value. You should use:

double precision function blub (x)
  implicit none
  double precision, intent(in) :: x

  blub = (1 / sqrt(8 * atan(1.0d0))) * exp((-x)**2 / 2)
end function

You can forget the return end stop you put before the end. There is no reason to use them.

I would recommend you to use even more debugging compiler flags, such a -g -fbacktrace -fcheck=all and put all functions and subroutines into a module or make them internal (using contains).