As part of my course I have written program to evaluate a complex integral numerically using Simpsons rule. To avoid singularities from the limits, I split the integral into parts, each with its own equivalent limits. I coded each of these parts as a function. I believe my maths is right, but the program is not converging as it should - I have spent days checking both. I suspect that due to my newness to fortran, I have a subtle error in my code and would appreciate some help finding it, I've reached the stage I have looked at it too often and a fresh pair of more experienced eyes would be a good idea. (sorry the code is not indented, I don't know how to preserve the indentation when copying)
program simpson !Summary: Integral from a to b of f(x) = h/3( f(a) + 4f(a+h) + 2f(a+2h) + 4f(a+3h) + ...4f(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+(coordtemp)
! 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 *, N
end do
end do
Print *, 'User selected end'
print *, stars
!---------- Admin --------------
10 format (3(A,F7.5))
20 format (A,8X,A)
30 format (4X, 2(F15.10, 5X))
40 format (2(A,I3),3(A,F8.6))
end program simpson
!-----------Functions ----------
Function FUNC_1(u) !calc 3* [(1-(u3))(-1/3)]
integer, parameter :: dp = selected_real_kind(15, 307)
real(kind=dp) :: u, FUNC_1, temp, out
temp=(1.0-(u3.0))
out=(temp)(-1/3.0)
out=3*out
! print 10, ' u= ',u, ' temp= ', temp, ' out= ', out !debug
Func_1= out
10 format (3(A,F10.7))
End Function
Function FUNC_2(v) !calc (3/2)* [(1-(v**(3/2))**(-2/3)]
integer, parameter :: dp = selected_real_kind(15, 307)
real(kind=dp) :: v, FUNC_2, out
temp=(1.0-(v**(3/2.0)))
out=(1/temp)**(2/3.0)
out=(3)*out
Func_2= out
10 format (2(A,F10.7))
End Function