SpoonTune

Untitled

Mar 4th, 2021
974
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1.       ! Michael Wiedner
  2.       ! Dr Fuse
  3.       ! CMS 495 - Numerical Methods
  4.       ! 3/5/2021
  5.       ! Exam 1 Question 2 Code - Bisection Search and False-Positive Search
  6.      
  7.       PROGRAM Problem 2
  8.    
  9.       integer :: i, degree, iterat, itera2
  10.       real :: a,b,tol,cnstnt,m,ya,yb,ym,atemp,btemp,tollo,
  11.      &tollo2,a2,b2,z
  12.      
  13.       real, dimension(5) :: c ! array of coefficients
  14.      
  15.       iterat = 1
  16.       itera2 = 1
  17.       do i = 1, 5
  18.           c(i) = 0
  19.       end do
  20.      
  21.       write(*,*) "Input the following information in order:"
  22.       write(*,*) "Min of range, max of range, tolerance, degree of pol."
  23.       write(*,*) "The first 3 values must be doubles and the last int."
  24.       read(*,*) a
  25.       read(*,*) b
  26.       if (a.ge.b) then
  27.           write(*,*) "Minimum is larger than maximum. Abort!"
  28.           stop
  29.       end if
  30.       read(*,*) tol
  31.       if (tol.LE.0.0) then
  32.           write(*,*) "Tolerance must be greater than zero. Abort!"
  33.           stop
  34.       end if
  35.       read(*,*) degree
  36.       if (degree.gt.5.or.degree.lt.1) then
  37.           write(*,*) "The degree must be 1, 2, 3, 4, or 5. Abort!"
  38.           stop
  39.       end if
  40.      
  41.       write(*,*) "Enter the coefficients of your polynomial from
  42.     &right to left (the coefficient in front of the x^1 first):"
  43.       do i = 1, degree
  44.           read(*,*) c(i)
  45.       end do
  46.      
  47.       write(*,*) "Enter the constant of your polynomial:"
  48.       read(*,*) cnstnt
  49.      
  50.       write(*,*)
  51.       write(*,*) "Your polynomial is being processed as the following:"
  52.       do i = 1, degree
  53.           write(*,*) c(degree - i + 1), "x to the ", (degree - i + 1)
  54.           write(*,*) "plus"
  55.       end do
  56.      
  57.       write(*,*) cnstnt
  58.       write(*,*) ! Empty line for neat printing
  59.      
  60.       ya = fvalue(a, c(1), c(2), c(3), c(4), c(5), cnstnt) ! y value at left bound
  61.       ya = fvalue(b, c(1), c(2), c(3), c(4), c(5), cnstnt) ! y value at right bound
  62.      
  63.       if (ya.eq.0.0) then ! Test if the boundaries are zeroes
  64.           write(*,*) "There exists a zero at the boundary", a
  65.           write(*,*) "This solution was found with no iterations."
  66.       else if (yb.eq.0.0) then
  67.           write(*,*) "There exists a zero at the boundary", b
  68.           write(*,*) "This solution was found with no iterations."
  69.       end if
  70.      
  71.       if ((a1*b1).gt.0.0) then ! Test if there exists one zero within the bounds
  72.           write(*,*) "The zero does not exist. Goodbye"
  73.           stop
  74.       end if
  75.      
  76.       tollo = tmptol(a,b,c(1),c(2),c(3),c(4),c(5),cnstnt) ! Variable for the Absolute Relative Error
  77.       tollo2 = tollo ! Store the user-inputs as copied variables so that the two different algorithms don't interfere
  78.       a2 = a
  79.       b2 = b
  80.      
  81.       do while(tollo.GE.tol) ! Main while loop for BISECTION
  82.           tollo = tmptol(a,b,c(1),c(2),c(3),c(4),c(5),cnstnt)  
  83.           m = (a+b)/2.0 ! define the midpoint as the sum of the bounds divided by 2
  84.          
  85.           atemp = a ! Placeholding variables for the values of a and b. No idea why but this fixed an error
  86.           btemp = b
  87.          
  88.           ya = fvalue(a, c(1), c(2), c(3), c(4), c(5), cnstnt) ! y value at left bound
  89.           yb = fvalue(b, c(1), c(2), c(3), c(4), c(5), cnstnt) ! y value at right bound
  90.           ym = fvalue(m, c(1), c(2), c(3), c(4), c(5), cnstnt) ! y value at midpoint
  91.          
  92.           if((ya*ym).lt.0.0) then ! if f of left bound times f of midpoint is negative, then the 0 is in it
  93.              ! Then 0 is in first half
  94.              ! Left bound stays the same
  95.             b=m ! Right bound becomes the m
  96.             a = atemp
  97.           end if
  98.           if((yb*ym).lt.0.0) then ! Else 0 is in the second half
  99.             a=m ! Left bound becomes the m
  100.              ! Right bound stays the same
  101.             b = btemp
  102.           end if
  103.          
  104.        iterat = iterat + 1
  105.        
  106.       end do ! End main while loop for BISECTION
  107.      
  108.      
  109.       write(*,*) "Using the bisection method,
  110.     &the zero has been found after", iterat, "iterations."
  111.      
  112.       if (ya.lt.0.0) then
  113.           ya = ya * (-1.0)
  114.       end if
  115.       if (yb.lt.0.0) then
  116.           yb = yb * (-1.0)
  117.       end if
  118.      
  119.       if (yb.GT.ya) then
  120.               write(*,*) "The zero is", a
  121.           else
  122.               write(*,*) "The zero is", b
  123.           end if
  124.      
  125.           write(*,*) ! Empty line for formatting
  126.          
  127.           ! Since the bisection algorithm has finished running by this point, we can replace the necessary variables with the ones
  128.           ! we stored previously.
  129.           tollo = tollo2
  130.           a = a2
  131.           b = b2
  132.       do while(tollo.GE.tol) ! Main while loop for FALSE-POSITION
  133.           tollo = tmptol(a,b,c(1),c(2),c(3),c(4),c(5),cnstnt)
  134.           z = seca(a,b,c(1),c(2),c(3),c(4),c(5),cnstnt) ! z is the x value where secant line of a and b intersects the x axis
  135.          
  136.           atemp = a ! Placeholding variables for the values of a and b. No idea why but this fixed an error
  137.           btemp = b
  138.          
  139.           ya = fvalue(a, c(1), c(2), c(3), c(4), c(5), cnstnt) ! y value at left bound
  140.           yb = fvalue(b, c(1), c(2), c(3), c(4), c(5), cnstnt) ! y value at right bound
  141.           yz = fvalue(z, c(1), c(2), c(3), c(4), c(5), cnstnt) ! y value at z
  142.          
  143.          
  144.           if((ya*yz).gt.0.0) then ! if f of left bound times f of c is negative, then the 0 is in it
  145.           a = z
  146.           b = btemp
  147.           end if
  148.           if((yb*yz).gt.0.0) then ! Else 0 is in the second half
  149.             b = z ! Left bound becomes the c
  150.              ! Right bound stays the same
  151.             a = atemp
  152.           end if
  153.          
  154.        itera2 = itera2 + 1
  155.        
  156.        write(*,*) "Bounds:", a, b
  157.        write(*,*) "C:", z
  158.        write(*,*) "Tolerance:", tollo
  159.        !if (itera2.gt.10) then
  160.        !stop
  161.        !end if
  162.        
  163.        
  164.       end do ! End main while loop for FALSE POSITION  
  165.      
  166.       write(*,*) "Using the false-position method,
  167.     &the zero has been found after", itera2, "iterations."
  168.      
  169.       if (ya.lt.0.0) then
  170.           ya = ya * (-1.0)
  171.       end if
  172.       if (yb.lt.0.0) then
  173.           yb = yb * (-1.0)
  174.       end if
  175.      
  176.       if (yb.GT.ya) then
  177.               write(*,*) "The zero is", a
  178.           else
  179.               write(*,*) "The zero is", b
  180.           end if
  181.      
  182.       STOP
  183.       END ! End the program
  184.      
  185.      
  186.       FUNCTION fvalue (x, c1, c2, c3, c4, c5, con)
  187.           real :: fvalue
  188.          
  189.           real :: x, c1, c2, c3, c4, c5, con
  190.          
  191.       fvalue = c5*(x**5) + c4*(x**4) + c3*(x**3) + c2*(x**2)
  192.      &+ c1*(x) + con
  193.       END FUNCTION
  194.      
  195.      
  196.      
  197.       FUNCTION tmptol (p1, p2, c1, c2, c3, c4, c5, con)
  198.           real :: tmptol
  199.          
  200.       real :: p1,p2,temp,c1,c2,c3,c4,c5,con,yp1,yp2,bst,oth,ty1,ty2
  201.          
  202.           yp1 = fvalue(p1, c1, c2, c3, c4, c5, con) ! y value of point one
  203.           yp2 = fvalue(p2, c1, c2, c3, c4, c5, con) ! y value of point two
  204.          
  205.           ! Find which of the two points is the best
  206.           ! Start by getting the absolute value of the y values of the points to see which is furthest from zero
  207.           if(yp1.lt.0.0) then
  208.               ty1 = yp1 * (-1.0)
  209.           else
  210.               ty1 = yp1
  211.           end if
  212.          
  213.           if (yp2.lt.0.0) then
  214.               ty2 = yp2 * (-1.0)
  215.           else
  216.               ty2 = yp2
  217.           end if
  218.          
  219.          
  220.           if (ty1.GT.ty2) then
  221.               bst = p2
  222.               oth = p1
  223.           else
  224.               bst = p1
  225.               oth = p2
  226.           end if
  227.          
  228.           temp = ((bst - oth) / bst)
  229.           if (temp.LT.0.0) then
  230.               temp = temp * (-1.0)
  231.           end if
  232.           tmptol = temp
  233.       END FUNCTION
  234.      
  235.      
  236.      
  237.       FUNCTION seca(p1, p2, c1, c2, c3, c4, c5, con)
  238.       real :: seca
  239.          
  240.       real :: p1,p2,c1,c2,c3,c4,c5,con,yp1,yp2
  241.      
  242.       yp1 = fvalue(p1, c1, c2, c3, c4, c5, con) ! y value of point one
  243.       yp2 = fvalue(p2, c1, c2, c3, c4, c5, con) ! y value of point two
  244.          
  245.       seca = (p2 - (yp2*((p2-p1)/(yp2-yp1)))) ! temp is the x value of the C point
  246.       END FUNCTION
RAW Paste Data

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×