gistfile1.f90

 program simpson 
!Summary: Integral from a to b of f(x) = h/3( f(a) + 4*f(a+h) + 2*f(a+2h) + 4*f(a+3h) + ...4*f(b-h) + f(b)) 
implicit none 
integer, parameter :: dp = selected_real_kind(15, 307) 
real, parameter :: pi = 3.1415927 
real (kind=dp) :: h1, h2, part1, part2, exact, approx, diff, FUNC_1, FUNC_2, a, b, c, d, arg, temp 
integer, parameter :: nn=50 
integer	:: N, cntr, coord 
character(len=nn) :: stars=REPEAT('*',nn) 

!----- Initialise ---- 
exact = (2*pi)/sqrt(3.0)	! Given 
a=0.0**(1/3.0) 
b=0.5**(1/3.0)	! 1st function integrates from 0 to 0.5 - change of variable u=t(1/3) 
d=(1.0-0.5)**(3/2.0)	! 2nd function from 0.5 to 1.0 - after change of variable v=(1-t)**(3/2), 
c=(1.0-1.0)**(3/2.0)	! - swap limits and change sign of integral 
!----- header ---- 
print *, " " 
print *, 'Program Simpson starts' 
print *, 'Simpson approximation of Integral for f(t) = (t**(-2/3)) * ((1-t)**(-1/3))' 
print *, 'Exact value = ', exact 
print *, stars 
print *, 'Enter value of N (multiples of 4) (0 to stop)' 
print 20, 'N approximate ','difference ' 
read (*,*) N 
do while ((mod(N,4).ne.0).OR.(n.lt.0)) 
print *, 'N must be positive and a multiple of 4 (or 0) - please re-enter N' 
read *, N 
end do 
!----- end header ---- 
do while (N.ne.0) ! allows repeated trials of diff. values of h 
! Evaluate 1st part 
h1=(b-a)/N	! N intervals of h between a and b 
print 10, ' Evalute first part from a= ', a, ' to b= ', b, ' at intervals h= ', h1 !debug 
print *, " " 
arg=a 
approx=FUNC_1(arg)	!1st term, f0 
Do cntr = 1,N-1	
if (mod(cntr,2)>0) then !set co-ord - either 4 (even) or 2 (odd) 
coord = 2 
else 
coord = 4 
end if 
arg=arg+h1 !co-ord terms, f(a+h) 
temp=FUNC_1(arg) 
approx=approx+(coord*temp) 
! print 40, ' term= ', cntr, ' coord= ', coord, ' arg= ', arg, ' function= ', temp, ' approx= ',approx	!debug 
end do 
part1 = (h1/3.0)*(approx + FUNC_1(b))	!last term + simpson h/3 multiplier 
! Evaluate 2nd part 
h2=(d-c)/N	! N intervals of h between c and d 
print 10, ' Evalute 2nd part from c= ', c, ' to d= ', d, ' at intervals h= ', h2 
arg=c 
approx=FUNC_2(arg)	!1st term 
Do cntr = 1,N-1	
if (mod(cntr,2)>0) then !set co-ord - either 4 (even) or 2 (odd) 
coord = 2 
else 
coord = 4 
end if 
arg=arg+h2 
temp=FUNC_2(arg) 
approx=approx+(coord*temp) 
! print 40, ' term= ', cntr, ' coord= ', coord, ' arg= ', arg, ' function= ', temp, ' approx= ',approx	!debug 
end do 
part2=(h2/3.0)*(approx + FUNC_2(d))	!last term + simpson h/3 multiplier 
print *, part1, ' ', part2 !debug 
approx = part1 + part2	
diff=exact-approx 
print 30, approx, diff 
read (*,*) N 
do while ((mod(N,4).ne.0).OR.(n.lt.0)) 
print *, 'N must be positive and a multiple of 4 (or 0) - please re-enter N' 
read *,