I am trying to implement a piece of code to get the area under a curve by getting the integral using Simpson's Rule
.
!(file:///D:/1-%20TUD/Semester%201/Numerical%20Methods%20BIWO-04/Lectures/simpson's%20rule.JPG)
I have already tried it using MatchCAD and I got correct results
function: f(x)= x**5+(x-2)*sin(x)+(x-1)
program simpsons
implicit none
real a, b, h
real integ, fa, fb
integer i, m
write(*,*) 'enter the lower boundary'
read(*,*) a
write(*,*) 'enter the upper boundary'
read(*,*) b
do while(a.ge.b)
write(*,*) 'reenter the lower boundary'
read(*,*) a
write(*,*) 'reenter the upper boundary'
read(*,*) b
enddo
write(*,*) 'enter the intervals number'
read(*,*) m
h=(b-a)/2.0
fa=a**5.0+(a-2.0)*sin(a)+(a-1.0)
fb=b**5.0+(b-2.0)*sin(b)+(b-1.0)
integ=0
do i=1, m/2
integ=integ+4*((a+(2*i-1)*h)**5.0+((a+(2*i-1)*h)-2.0)*sin((a+(2*i-1)*h))+ ((a+(2*i-1)*h)-1.0))
if (i.le.((m/2)-1)) then
integ=integ+2*((a+2*i*h)**5.0+((a+2*i*h)-2.0)*sin((a+2*i*h))+((a+2*i*h)-1.0))
endif
enddo
integ=(fa+fb+integ)*(h/3.0)
write(*,*) 'integration = ', integ
end
when I input a=-1
and b=1
and m=20
I should get Integration=-1.398
when I input a=-1
and b=1
and m=40
I should get Integration=-1.398
but somehow I am getting the integration =7015869.0
I was using this code some time ago. As pointed out by @VladimirF, it is better to separate the function f(x)
from the algorithm simpson(f,a,b,integral,n)
.
program main
implicit none
double precision a, b, integral
integer n
a = -1; b = 1; n = 20
call simpson(f,a,b,integral,n)
write (*,*) 'integration = ', integral
contains
function f(x)
implicit none
double precision f, x
f = x**5 + (x-2)*sin(x) + (x-1)
end
Subroutine simpson(f,a,b,integral,n)
implicit none
double precision f, a, b, integral, s
double precision h, x
integer n, i
! if n is odd we add +1 to make it even
if((n/2)*2.ne.n) n = n+1
! loop over n (number of intervals)
s = 0.0
h = (b-a)/dble(n)
do i = 2, n-2, 2
x = a + i*h
s = s + 2.0*f(x) + 4.0*f(x+h)
end do
integral = (s + f(a) + f(b) + 4.0*f(a+h))*h/3.0
end subroutine simpson
end program main
with the correct output:
integration = -1.39766605364858