 # 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.

×