Intégration de gauss

! programme principal

!-------------------------------------------------------------------

      !implicit none

      external quad3d,qtrapx,qtrapy,qtrapz,y1,y2,z1,z2,func

      double precision qtrapx,qtrapy,qtrapz,y1,y2,z1,z2,func

	!double precision x1,x2,ss,pk

 
	x1=0.0

      x2=100.0

 
      call quad3d(x1,x2,ss)

 
      pk=ss

 
      write(6,*)pk

 
      end

 
! subroutine qui calcul une integrale triple

!-------------------------------------------------------------------

      subroutine quad3d(x1,x2,ss)

      real ss,x1,x2,h

      external h

      call qgausx(h,x1,x2,ss)

      return

      end

      function f(zz)

      real f,zz,func,x,y,z

      common /xyz/ x,y,z

      z=zz

      f=func(x,y,z)

      return

      end

 
      function g(yy)

      real g,yy,x,y,z,f,z1,z2

	external f

      common /xyz/ x,y,z

      y=yy

      call qgausz(f,z1(x,y),z2(x,y),ss)

      g=ss

      return

      end

 
      function h(xx)

      real h,xx,g,y1,y2,x,y,z

      external g

      common /xyz/ x,y,z

      real ss

      x=xx

      call qgausy(g,y1(x),y2(x),ss)

      h=ss

      return

      end

 
 
! LA fonction à integrer

!---------------------------------------------------------------------

 
      FUNCTION func(x,y,z)

 
      !func= x*(x+3*y)

	func=1.0

 
      RETURN

 
      END

 
! LA fonction z1(x)

!---------------------------------------------------------------------

 
      FUNCTION z1(x,y)

 
      z1= 0.0

 
      RETURN

 
      END

 
! LA fonction z2(x)

!---------------------------------------------------------------------

 
      FUNCTION z2(x,y)

 
      z2=100.0-x-y

 
      RETURN

 
       END

 
 
! LA fonction y1(x)

!---------------------------------------------------------------------

 
      FUNCTION y1(x)

 
      y1= 0.0

 
      RETURN

 
      END

 
! LA fonction y2(x)

!---------------------------------------------------------------------

 
      FUNCTION y2(x)

 
      y2= 100.0-x

 
      RETURN

 
       END

 
! INTEGRATION PAR LA METHODE DE GAUSS

!-----------------------------------------------------------------

      SUBROUTINE qgausx(func,a,b,ss)

      REAL a,b,ss,func

      EXTERNAL func

      INTEGER j

      REAL dx,xm,xr,w(5),x(5)

      SAVE W,x

      DATA w/.2955242247,.2692667193,.2190863625,.149451349,.0666713449/

      DATA x/.1488743389,.433395394,.6794095682,.8650633666,.9739065285/

      xm=0.5*(b+a)

      xr=0.5*(b-a)

      ss=0

      do 11 j=1,5

      dx=xr*x(j)

      ss=ss+w(j)*(func(xm+dx)+func(xm-dx))

11    continue

      ss=xr*ss

      return

      END

 
 
! INTEGRATION PAR LA METHODE DE GAUSS

 
      SUBROUTINE qgausy(func,a,b,ss)

      REAL a,b,ss,func

      EXTERNAL func

      INTEGER j

      REAL dx,xm,xr,w(5),x(5)

      SAVE W,x

      DATA w/.2955242247,.2692667193,.2190863625,.149451349,.0666713449/

      DATA x/.1488743389,.433395394,.6794095682,.8650633666,.9739065285/

      xm=0.5*(b+a)

      xr=0.5*(b-a)

      ss=0

      do 11 j=1,5

      dx=xr*x(j)

      ss=ss+w(j)*(func(xm+dx)+func(xm-dx))

11    continue

      ss=xr*ss

      return

      END

 
! INTEGRATION PAR LA METHODE DE GAUSS

 
      SUBROUTINE qgausz(func,a,b,ss)

      REAL a,b,ss,func

      EXTERNAL func

      INTEGER j

      REAL dx,xm,xr,w(5),x(5)

      SAVE W,x

      DATA w/.2955242247,.2692667193,.2190863625,.149451349,.0666713449/

      DATA x/.1488743389,.433395394,.6794095682,.8650633666,.9739065285/

      xm=0.5*(b+a)

      xr=0.5*(b-a)

      ss=0

      do 11 j=1,5

      dx=xr*x(j)

      ss=ss+w(j)*(func(xm+dx)+func(xm-dx))

11    continue

      ss=xr*ss

      return

      END