summaryrefslogtreecommitdiff
path: root/CSMath/src/bisection/NeoBisection.java
blob: 530fbeb5a086b9778a73c2fbe323029f8b2b9e18 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
package bisection;

import java.util.function.DoubleUnaryOperator;

import bjc.utils.math.Dual;
import bjc.utils.math.DualExpr;

/*
 * Benjamin Culkin
 * 1/16/2018
 * CS 479
 * Bisection
 * 
 * Use the bisection method to find bracketed roots of arbitrary real-valued functions.
 */

/**
 * Use the bisection method to find bracketed roots of arbitrary real-valued
 * functions.
 * 
 * @author student
 *
 */
public class NeoBisection {
	/**
	 * Maximum number of iterations to attempt.
	 */
	public static final int MAXITR = 500;

	/**
	 * Main method
	 * 
	 * @param args
	 *            Unused CLI args
	 */
	public static void main(String[] args) {
		/* Bisection Method. */

		/* The functions we're approximating a root for. */
		DoubleUnaryOperator functionA = x -> Math.cos(x) - x;

		DoubleUnaryOperator functionB = x -> {
			double valA = Math.pow(x, 5);
			double valB = 3 * Math.pow(x, 4);
			double valC = 4 * Math.pow(x, 3);

			return valA - valB + valC - 1;
		};

		DoubleUnaryOperator functionC = x -> (x * x) - 3;

		/* The approximated roots. */
		double answerA = bisect(functionA, 0, 1, 0.0001);
		double answerB = bisect(functionB, 0, 1, 0.0001);
		double answerC = bisect(functionC, 1, 2, 0.0000001);

		/* Print out the answers + their function. */
		System.out.printf("Bisection Method:\n");
		System.out.printf("\tx   = cos(x)          => %.4f\n", answerA);
		System.out.printf("\tx^5 - 3x^4 + 4x^2 - 1 => %.4f\n", answerB);
		System.out.printf("\tx^2 - 3               => %.7f\n\n", answerC);

		/* Newton's / Secant Method. */
		/* The variable to manipulate in the expressions. */
		Dual varX = new Dual();
		DualExpr varXExpr = new DualExpr(varX);

		/* The functions to find the roots of. */
		DualExpr dualA = new DualExpr(DualExpr.ExprType.SUBTRACTION, new DualExpr(DualExpr.ExprType.COS, varXExpr),
				varXExpr);

		DualExpr dualB;

		{
			/* Construct the second function. */
			DualExpr tempA = new DualExpr(varXExpr, 5);
			DualExpr tempB = new DualExpr(DualExpr.ExprType.MULTIPLICATION, new DualExpr(new Dual(3)),
					new DualExpr(varXExpr, 4));
			DualExpr tempC = new DualExpr(DualExpr.ExprType.MULTIPLICATION, new DualExpr(new Dual(4)),
					new DualExpr(varXExpr, 3));

			dualB = new DualExpr(DualExpr.ExprType.ADDITION, new DualExpr(DualExpr.ExprType.SUBTRACTION, tempA, tempB),
					new DualExpr(DualExpr.ExprType.SUBTRACTION, tempC, new DualExpr(new Dual(1))));
		}

		DualExpr dualC = new DualExpr(DualExpr.ExprType.SUBTRACTION, new DualExpr(varXExpr, 2),
				new DualExpr(new Dual(3)));

		/* Print out dualized expressions. */
		System.out.printf("Expressions:\n");
		System.out.printf("\t%s\n", dualA);
		System.out.printf("\t%s\n", dualB);
		System.out.printf("\t%s\n", dualC);

		/* The approximated roots, using Newton's method. */
		double newtonA = newton(dualA, varX, 0, 1, 0.0001);
		double newtonB = newton(dualB, varX, 0, 1, 0.0001);
		double newtonC = newton(dualC, varX, 1, 2, 0.0000001);

		/* Print out the answers + their function. */
		System.out.printf("Newton's Method:\n");
		System.out.printf("\tx   = cos(x)          => %.4f\n", newtonA);
		System.out.printf("\tx^5 - 3x^4 + 4x^2 - 1 => %.4f\n", newtonB);
		System.out.printf("\tx^2 - 3               => %.7f\n\n", newtonC);

		/* The approximated roots, using the secant method. */
		double secantA = secant(dualA, varX, 0, 1, 0.0001);
		double secantB = secant(dualB, varX, 0, 1, 0.0001);
		double secantC = secant(dualC, varX, 1, 2, 0.0000001);

		/* Print out the answers + their function. */
		System.out.printf("Secant Method:\n");
		System.out.printf("\tx   = cos(x)          => %.4f\n", secantA);
		System.out.printf("\tx^5 - 3x^4 + 4x^2 - 1 => %.4f\n", secantB);
		System.out.printf("\tx^2 - 3               => %.7f\n", secantC);
	}

	/*
	 * Calculate the number of iterations to approximate to a specified tolerance.
	 */
	private static double iterCount(double a, double b, double tol) {
		double inside = Math.log((b - a) / tol);

		return Math.ceil(inside / Math.log(2));
	}

	/* Calculate a root for an equation by the bisection method. */
	private static double bisect(DoubleUnaryOperator func, double a, double b, double tol) {
		/* Calculate the number of iterations. */
		double N = iterCount(a, b, tol);

		double val = 0.0;

		/* Set our values. */
		double newa = a;
		double newb = b;
		double newc = (a + b) / 2;

		/*
		 * Continue onwards until we get the approximation to within the right
		 * tolerance.
		 */
		for (int i = 0; i < N; i++) {
			newc = (newa + newb) / 2;

			val = func.applyAsDouble(newc);

			/* Pick the right direction to bisect in. */
			if (func.applyAsDouble(newa) * val < 0) {
				newb = newc;
			} else {
				newa = newc;
			}
		}

		/* Return the right value. */
		return newc;
	}

	/**
	 * Use Newton's method to find the root of an equation.
	 *
	 * @param func
	 *            The function to find a root for.
	 *
	 * @param var
	 *            The variable in the function.
	 *
	 * @param lo
	 *            The lower bound for the root
	 *
	 * @param hi
	 *            The higher bound for the root
	 *
	 * @param tol
	 *            The tolerance for the answer
	 *
	 * @return The estimated value for the equation.
	 */
	public static double newton(DualExpr func, Dual var, double lo, double hi, double tol) {
		/* Initial guess for root. */
		double newmid = (lo + hi) / 2;

		for (int i = 0; i < MAXITR; i++) {
			/* Previous root guess. */
			double prevmid = newmid;

			/* Set the variable properly. */
			var.real = newmid;
			var.dual = 1;

			/* Evaluate the function and its derivative. */
			Dual res = func.evaluate();

			/* Use Newton's method to refine our solution. */
			newmid = newmid - (res.real / res.dual);
			System.out.printf("\tTRACE: prevmid: %f\t\tnewmid: %f\n", prevmid, newmid);
			/* Hand back the solution if it's good enough. */
			if (Math.abs(newmid - prevmid) < tol) {
				return newmid;
			}
		}

		System.out.println("Newton's method iteration limit reached.");

		/* Give back the solution. */
		return newmid;
	}

	/**
	 * Calculate a root using the secant method.
	 * 
	 * @param func
	 *            The function to calculate a root for.
	 * @param var
	 *            The variable in the function.
	 * @param lo
	 *            The lower bound of the root.
	 * @param hi
	 *            The upper bound of the root.
	 * @param tol
	 *            The tolerance to find the root to.
	 * @return The root, to within the desired tolerance.
	 */
	public static double secant(DualExpr func, Dual var, double lo, double hi, double tol) {
		/* Initial guesses for root. */
		double guess1 = (lo + hi) / 3; // 1/3 into the range
		double guess2 = ((lo + hi) / 3) * 2; // 2/3 into the range

		for (int i = 0; i < MAXITR; i++) {
			var.real = guess1;
			var.dual = 1;
			/* Evaluate the first guess. */
			Dual res1 = func.evaluate();

			var.real = guess2;
			var.dual = 1;
			/* Evaluate the first guess. */
			Dual res2 = func.evaluate();

			double oldGuess1 = guess1;

			/* Use the secant method to refine our guesses. */
			guess1 = guess2;
			guess2 = ((oldGuess1 * res2.real) - (guess1 * res1.real)) / (res1.real - res2.real);

			System.out.printf("\tTRACE: guess1: %f\t\tguess2: %f\n", guess1, guess2);
			/* Hand back the solution if it's good enough. */
			if (Math.abs(guess1 - guess2) < tol) {
				return guess2;
			}
		}

		System.out.println("Secant method iteration limit reached.");

		/* Give back the solution. */
		return guess2;
	}
}