# 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."
26.       if (a.ge.b) then
27.           write(*,*) "Minimum is larger than maximum. Abort!"
28.           stop
29.       end if
31.       if (tol.LE.0.0) then
32.           write(*,*) "Tolerance must be greater than zero. Abort!"
33.           stop
34.       end if
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
45.       end do
46.
47.       write(*,*) "Enter the constant of your polynomial:"
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