A_04 GUASSIAN QUADRATURE METHOD, Fortran

Language: Layout:

Program Guassian_Quadrature
 call GausiianQuadSingleIntegral(3)
 call GausiianQuadSingleIntegral(4)
 call GausiianQuadDoubleIntegral(3,4)
 call GausiianQuadDoubleIntegral(4,3)
End Program
Subroutine GausiianQuadSingleIntegral(n)
 Implicit None
 Real::c(4),t(4),f,s
 Integer::i,j
 Integer,intent(in)::n
 if (n==4)then
 c(1) = 0.3478548451
 c(2) = 0.6521451549
 c(3) = c(2)
 c(4) = c(1)
 t(1) = 0.8611363116
 t(2) = 0.3399810436
 t(3) = -1.0 * t(2)
 t(4) = -1.0 * t(1)
 elseif(n==3)then
 c(1) = 0.5555555556
 c(2) = 0.8888888889
 c(3) = c(1)
 t(1) = 0.7745966692
 t(2) = 0
 t(3) = -1.0 * t(1)
 end if
 s=0
 Do i=1,n
 s = s+ c(i)*F(t(i))
 End Do
 Write(*,'(5x,A,I1,A,F20.18)')"Single Integration using Guassian Quadrature for n=",n," : ",s
End Subroutine

Subroutine GausiianQuadDoubleIntegral(n,m)
 Implicit None
Real::cn(4),tn(4),cm(4),tm(4),f,s,h1,h2,k1&
&,k2,a,b,c1,d1,sx,x,y,fxy
 Integer::i,j
 Integer,intent(in)::n,m
 if (n==4)then
 cn(1) = 0.3478548451
 cn(2) = 0.6521451549
 cn(3) = cn(2)
 cn(4) = cn(1)
 tn(1) = 0.8611363116
 tn(2) = 0.3399810436
 tn(3) = -1.0 * tn(2)
 tn(4) = -1.0 * tn(1)
 elseif(n==3)then
 cn(1) = 0.5555555556
 cn(2) = 0.8888888889
 cn(3) = cn(1)
 tn(1) = 0.7745966692
 tn(2) = 0
 tn(3) = -1.0 * tn(1)
 end if
 if (m==4)then
 cm(1) = 0.3478548451
 cm(2) = 0.6521451549
cm(3) = cm(2)
 cm(4) = cm(1)
 tm(1) = 0.8611363116
 tm(2) = 0.3399810436
 tm(3) = -1.0 * tm(2)
 tm(4) = -1.0 * tm(1)
 elseif(m==3)then
 cm(1) = 0.5555555556
 cm(2) = 0.8888888889
 cm(3) = cm(1)
 tm(1) = 0.7745966692
 tm(2) = 0
 tm(3) = -1.0 * tm(1)
 end if
 a=1.4
 b=2.0
 c1 = 1.0
 d1 = 1.5
 h1=(b-a)/2
 h2 = (b+a)/2
 s=0
 Do i = 1,m
 sx =0
 x = h1*tm(i) + h2
 k1 = (d1- c1)/2.0
 k2 = (d1+c1)/2.0
 Do j=1,n
 y = k1*tn(j)+k2
 sx = sx + cn(j) * fxy(x,y)
 End Do
s = s + cm(i)*k1*sx
 End Do
 s = s*h1
 Write(*,'(5x,A,I1,A,I1,A,F20.18)')"Single& &Integration using Guassian Quadrature for& &n=",n,", m= ",m," : ",s
End Subroutine

Function f(t)
 Implicit None
 Real::t,f
 f = 0.5*(((1.0+t)/2.0)**2.)*exp((1.0+t)/2.0)
End Function

function fxy(x,y)
 Implicit None
 Real::x,fxy,y
 fxy = Log(x+2.0*y)
end function

Show compiler warnings [ + ] Show input
Compilation time: 0.12 sec, absolute running time: 0.15 sec, cpu time: 0.01 sec, memory peak: 6 Mb, absolute service time: 0,37 sec