c<html>
c<head><title>trapz1.f</title></head>
c<body>
c<pre>
program trapezoid
c
c Test Trapezoidal Rule of Integration
c
c John Mahaffy 4/9/95
c
c Include the contents of file "trapcom.h" in this program
c
include 'trapcom.h'
integer i
do 10 i=1,ntrys
np=nint1*2**(i-1)+1
call setcurve
call integrate
call compare
10 continue
stop
end
subroutine setcurve
c
c Set up a Curve and range for integration
c Here we will integrate sin(x) from x=0 to x=3.0
c
c John Mahaffy 4/9/95
c
include 'trapcom.h'
integer i
real dx,xlow,xhigh
data xlow,xhigh/0.0,3.0/
c
dx=(xhigh-xlow)/(np-1)
x(1)=xlow
c<a name=1><font color=FF0000>
y(1)=sin(x(1))
c</font>
do 10 i=2,np
x(i)=x(i-1)+dx
y(i)=sin(x(i))
10 continue
ansreal=cos(xlow)-cos(xhigh)
return
end
subroutine integrate
c
c Use Trapezoidal Rule to Integrate the area under the curve
c specified by the data points in x and y
c
c John Mahaffy 4/9/95
c
include 'trapcom.h'
c
c Local Variables
c
c esterr - estimated error for Trapizoidal integration
c sum2 - Trapizoidal integration using every other available point
c
integer i
real sum2, esterr
ansnum=0.0
do 10 i=1,np-1
10 ansnum=ansnum+.5*(y(i)+y(i+1))*(x(i+1)-x(i))
c
c The following only works if np is an odd integer
c
sum2=0.0
do 20 i=1,np-1,2
20 sum2=sum2+.5*(y(i)+y(i+2))*(x(i+2)-x(i))
esterr=(sum2-ansnum)*.333333
write(6,2000) np,ansnum,esterr
2000 format (' For ',i4,' points the integral is ',1p,e12.5,
& ', Estimated Error=',e10.2)
return
end
subroutine compare
c
c Compare results of the numerical integration with the actual result
c
c John Mahaffy 4/9/95
c
implicit none
c
c ansreal - real answer from analytic integration
c ansnum - answer resulting from the numerical integration
c
include 'trapcom1.h'
real error
error=ansnum-ansreal
write(6,2001) error
2001 format(' Actual Error = ',1p,e10.2)
return
end
c</pre>
c</body>
c</html>
Copyright ©2017 - STI - Todos os direitos reservados