Advertisement
am1x

lab3b

Mar 14th, 2024
1,309
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Fortran 2.76 KB | Fixit | 0 0
  1. module functions
  2.     implicit none
  3.     contains
  4.    
  5.     function f(x) result(y)
  6.         real(8), intent(in) :: x
  7.         real(8) :: y
  8.         y = ((((x * x) - 3.5d0) * x + 2.5d0) * x - 7d0) * x - 6.4d0
  9.     end function
  10.    
  11.    
  12.     function F1(x) result(y)
  13.         real(8), intent(in) :: x
  14.         real(8) :: y
  15.  
  16. !        y = (1d0/6d0)*(x**6) - (3.5d0/4d0)*(x**4) + (2.5d0/3d0)*(x**3) - (7d0/2d0)*(x**2) - 6.4d0*x
  17.         y = (((((1d0/6d0) * x * x -  3.5d0/4d0) * x + 2.5d0/3d0) * x - 7d0/2d0) * x - 6.4d0) * x
  18.     end function
  19.    
  20.    
  21.     function trap(grid, N, h) result(int1)
  22.         integer, intent(in) :: N
  23.         integer :: i
  24.         real(8), dimension(:) :: grid(N+1)
  25.         real(8) :: h, int1
  26.         int1 = 0d0
  27.         do i = 2, N
  28.             int1 = int1 + grid(i)
  29.         end do
  30.         int1 = (2d0*int1 + grid(1) + grid(N+1))*h/2d0
  31.     end function
  32.    
  33.    
  34.     subroutine intrect(a, b, eps, err, h, N)
  35.         integer :: N, i
  36.         real(8) :: a, b, eps, err, h
  37.         real(8) :: int1, prevint, k, r
  38.         prevint = 0d0
  39.         int1 = 100d0
  40.         N = 1
  41.         do while (dabs(int1 - prevint)/3d0 > (eps))
  42.             N = N * 2
  43.             h = (b - a) / N
  44.             prevint = int1
  45.             int1 = 0
  46.             do i = 0, N - 1
  47.                 int1 = int1 + f(a + (i + 0.5d0) * h)
  48.             end do
  49.             int1 = int1 * h
  50.         end do
  51.         err = dabs((F1(b) - F1(a)) - int1)
  52.         write(*,*) N, eps, err
  53.     end subroutine
  54.  
  55.     subroutine inttrap(a, b, eps, err, h, N)
  56.         integer :: N, i
  57.         real(8) :: a, b, eps, err, h
  58.         real(8) :: int1, prevint
  59.         prevint = 0d0
  60.         int1 = 100d0
  61.         N = 1
  62.         do while (dabs(int1 - prevint) / 3 > (eps))
  63.             N = N * 2
  64.             h = (b - a) / N
  65.             prevint = int1
  66.             int1 = 0.5d0 * (f(a) + f(b))
  67.             do i = 1, N - 1
  68.                 int1 = int1 + f(a + i * h)
  69.             end do
  70.             int1 = int1 * h
  71.         end do
  72.         err = dabs((F1(b) - F1(a)) - int1)
  73.         write(*,*) N, eps, err
  74.     end subroutine
  75.  
  76.     subroutine intsimps(a, b, eps, err, h, N)
  77.         integer :: N, i
  78.         real(8) :: a, b, eps, err, h
  79.         real(8) :: int1, inta, intb, prevint
  80.         prevint = 0d0
  81.         int1 = 100d0
  82.         N = 1
  83.         do while (dabs(int1 - prevint) / 15 > eps)
  84.             N = N * 2
  85.             h = (b - a) / N
  86.  
  87.             inta = 0.5d0 * (f(a) + f(b))
  88.             do i = 1, N - 1
  89.                 inta = inta + f(a + i * h)
  90.             end do
  91.             intb = 0
  92.             do i = 0, N - 1
  93.                 intb = intb + f(a + (i + 0.5d0) * h)
  94.             end do
  95.  
  96.             prevint = int1
  97.             int1 = ( inta + 2 * intb) * h / 3.0d0
  98.  
  99.         end do
  100.  
  101.         err = dabs((F1(b) - F1(a)) - int1)
  102.         write(*,*) N, eps, err
  103.     end subroutine
  104.  
  105.  
  106.  
  107. end module
  108.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement