Professor Diomar Cesar Lobao

Universidade Federal Fluminense-Volta Redonda, RJ, Brasil

Diomar Cesar


Dept. Ciências Exatas - Exact Science Dept.

Search

trapz2.f

c<html>
c<head><title>trapz2.f</title></head>
c<body>
c<pre>
      module commdata
         integer ntrys, nint1,npoints,np,dbl
         parameter (ntrys=8,nint1=10,npoints=nint1*2**(ntrys-1)+1 )
         parameter (dbl = selected_real_kind(13,200))
         real(dbl) x(npoints),y(npoints)
         real(dbl) ansreal,ansnum
c<a name="em"><font color="FF0000">
      end module commdata
c</a></font>
      program trapezoid
c
c    The use statements must be the first thing found after the PROGRAM
c    statement
c
      use commdata
c
c    Test Trapezoidal Rule of Integration
c
c   John Mahaffy 5/10/96
c
      implicit none
c   PARAMETERS
c    ntrys - Number of integrations, each with half the x step of the last
c    nint1 - Number of steps between lowest and highest x values
c    npoints - Largest number of points evaluated for use in the integration
c
c   Variables
c    x   -   x values at which the function is evaluated
c    y   -   values of the function for corresponding values of x
c    np  -   number of points currently used in x and y
c    ansreal - Actual Value of the integral
c    ansnum  - Value of the integral obtained from a Trapezoidal Rule
c
      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    Use statements must be the first thing found after the SUBROUTINE
c    statement
c
      use commdata
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 5/10/96
c
      implicit none
c
      integer i
      real(dbl) dx,xlow,xhigh
      data xlow,xhigh/0.0,3.0/
c
      dx=(xhigh-xlow)/(np-1)
      x(1)=xlow
      y(1)=sin(x(1))
      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
c
      subroutine integrate
c<a name="use"><font color="FF0000">
      use commdata
c</font></a>
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 5/10/96
c
      implicit none
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(dbl) sum2, esterr
c<a name=1><font color=FF0000>
      ansnum=0.5*dot_product(y(1:np-1)+y(2:np),x(2:np)-x(1:np-1))
c</font>
c
c   The following only works if np is an odd integer
c
      sum2 = 0.5*dot_product(y(1:np-1:2)+y(3:np:2),
     &                       x(3:np:2)-x(1:np-1:2))
      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
c</pre>
c</body>
c</html>
      subroutine compare
c
c    The following USE only picks up information about ANSNUM and
c    ANSREAL from the MODULE named COMMDATA, ignores other information
c    If I used the variable names "x" or "y" here they would be treated
c    as local variables, unrelated to the "x" and "y" defined in
c    COMMDATA
c
c   John Mahaffy 5/10/96
c
      use commdata , only : ansnum, ansreal, dbl
c
c  Compare results of the numerical integration with the actual result
c
      implicit none
      real (dbl) error
      error=ansnum-ansreal
      write(6,2001) error
 2001 format(' Actual Error = ',1p,e10.2)
      return
      end
c</pre>
c</body>
c</html>
Skip to content