1.

program trapez
! integrate f(x) using the trapezoid rule
! the integral function of int (sin x)**2 is x/2 - (sin 2x)/4
! the exact value is then pi/4
  implicit none
  real, parameter :: pi = 3.141592654
  integer :: i
  real :: integral, exact, error
  
  exact= pi/4 ;
  do i=1,10
    integral = trapezint(0.0, pi/2, i)
    error = abs((integral-exact)/exact)
    write(*,'(I3,2F10.6)') i,integral, error
  end do

contains

real function f(x)
  real, intent(in) :: x
  f = sin(x)**2
end function

real function trapezint(a, b, n)
! integrate f over [a,b] yli
! using the trapezoid rule and n subdivisions
  real, intent(in) :: a, b
  integer, intent(in) :: n
  real :: sum, h, f0, f1
  integer :: i

  sum=0
  h=(b-a)/n
  f0 = f(a)
  do i=1, n
     f1=f(a+i*h)
     sum = sum + h*(f0+f1)/2
     f0=f1
  end do
  trapezint=sum
end function
end program


./a.out
  1  0.785398  0.000000
  2  0.785398  0.000000
  3  0.785398  0.000000
  4  0.785398  0.000000
  5  0.785398  0.000000
  6  0.785398  0.000000
  7  0.785398  0.000000
  8  0.785398  0.000000
  9  0.785398  0.000000
 10  0.785398  0.000000



2.


program gaussint
! evaluate an integral using the Gaussian quadrature
  implicit none
  real, parameter :: pi = 3.141592654
  integer :: i
  real :: int, exact, error
  
  exact=pi/4
  int = gaussint5(0.0, pi/2)
  error = abs((int-exact)/exact)
  write(*,'(F10.6,E12.3)') int, error

contains
real function f(x)
  real, intent(in) :: x
  f = sin(x)**2
end function

real function gaussint5(a, b)
! lasketaan funktion f integraali valin [a,b] yli
! 5 pisteen Gaussin menetelmalla
  real, intent(in) :: a, b
  real, dimension(5) :: &
    x = (/ &
        -0.906179845938664, &
        -0.538469310105684, &
         0.000000000000000, & 
         0.538469310105684, &
         0.906179845938664 /), & 
    w = (/ &
         0.236926885056189, &
         0.478628670499366, &
         0.568888888888889, &
         0.478628670499366, &
         0.236926885056189 /)
  real :: sum
  integer :: i

  sum=0
  do i=1, 5
     sum = sum + w(i) * f(x(i)*(b-a)/2 + (b+a)/2)
  end do
  gaussint5=sum*(b-a)/2
end function
end program

./a.out
  0.785398   0.000E+00



3.


program simpson2
! evaluate a 2D integral using the Simpson's rule
  implicit none
  integer :: i
  do i=4,10,2
    write(*,*) simpsonint2(0.0, 1.0, i, 0.0, 1.0, i)
  end do
contains
  
real function f(x, y) 
  ! the integrand
  real, intent(in) :: x, y
  f=exp(x*y)
end function

real function simpsonint1(a, b, n, y)
  ! integratef(x, y) over x=a..b
  ! using n subdivisions
  real, intent(in) :: y, a, b
  integer, intent(in) :: n

  real :: h, &    ! step length
          s2, s4  ! partial sums
  integer i

  h = (b-a)/n
  s2=0.0
  s4=0.0
  do i=1,n-1,2
     s4=s4+f(a+i*h, y)
  end do
  do i=2,n-2,2
     s2=s2+f(a+i*h, y)
  end do
  simpsonint1 = (h/3)*(f(a, y)+f(b, y)+2*s2+4*s4)
end function

real function simpsonint2(a, b, nx, c, d, ny)
  ! double integral of f(x,y) ovar x=a..b, y=c..d 
  ! use nx subdivisions ib the x direction and 
  ! ny subdivisions in the  y direction
  real, intent(in) :: a, b, c, d
  integer, intent(in) :: nx, ny

  real :: h, &    ! askelpituus
          s2, s4  ! osasummat
  integer i

  if (2*(nx/2)/= nx .or. 2*(ny/2)/=ny) then
     write(*,*) 'nr of subdivisons must be even:',nx,ny
     simpsonint2=0.0
     return
  end if

  h = (d-c)/ny
  s2=0.0
  s4=0.0
  do i=1,ny-1,2
     s4=s4+simpsonint1(a, b, nx, c+i*h)
  end do
  do i=2,ny-2,2
     s2=s2+simpsonint1(a, b, nx, c+i*h)
  end do
  simpsonint2 = (h/3)*(simpsonint1(a, b, nx, c) &
                      +simpsonint1(a, b, nx, d) &
                      +2*s2+4*s4)
end function

end program


./a.out
   1.317916    
   1.317905    
   1.317903    
   1.317903    


4.


program sphere
! find the volume of an n dimensional sphere
  implicit none
  call montecarlo(3)
  call montecarlo(4)
  call montecarlo(5)

contains

subroutine montecarlo(dim)
! volume of an n dimensional sphere by MC integration
  integer, intent(in) :: dim
  real, dimension(dim) :: r
  real dist, vol, exact
  integer i, j, iv, total
  
  exact=volume(1.0, dim)
  write(*,'("n=",I2," exact=",F8.3)') dim, exact

  total=0;   iv=0
  do i=1,10
  do j=1,100000
     call random_number(r)
     r=-1.0+2.0*r
     dist=sum(r**2)
     if (dist < 1.0) iv=iv+1
     total=total+1
  end do
  vol = real(iv)/total*(2**dim)
  write(*,'(I8,2F8.3)') total,vol,abs((vol-exact)/exact)
  end do
end subroutine

real function volume(radius, n)
! exact volume of an n dimensional sphere
  real, intent(in) :: radius
  integer, intent(in) :: n
  real, parameter :: pi=3.141592654
  real t
  integer i, f, m

  f=1
  if (mod(n,2)==0) then
    do i=1,n/2
       f=f*i
    end do
    t=real(f)    
  else
    m=(n-1)/2
    do i=m+2, 2*m+2
       f=f*i
    end do
    t=sqrt(pi)*f/(4.0**(m+1))
  end if

  volume=radius**n * pi**(n/2.0) / t
end function
end program


n= 3 exact=   4.189
  100000   4.208   0.005
  200000   4.197   0.002
  300000   4.199   0.002
  400000   4.195   0.001
  500000   4.193   0.001
  600000   4.191   0.000
  700000   4.191   0.000
  800000   4.190   0.000
  900000   4.192   0.001
 1000000   4.191   0.000
n= 4 exact=   4.935
  100000   4.919   0.003
  200000   4.911   0.005
  300000   4.931   0.001
  400000   4.925   0.002
  500000   4.930   0.001
  600000   4.930   0.001
  700000   4.932   0.001
  800000   4.929   0.001
  900000   4.931   0.001
 1000000   4.929   0.001
n= 5 exact=   5.264
  100000   5.372   0.021
  200000   5.312   0.009
  300000   5.300   0.007
  400000   5.266   0.000
  500000   5.258   0.001
  600000   5.253   0.002
  700000   5.255   0.002
  800000   5.259   0.001
  900000   5.253   0.002
 1000000   5.247   0.003