From 88fc0b666c58334009bc274105ea0d2b90edf75d Mon Sep 17 00:00:00 2001 From: Benjamin Culkin Date: Mon, 9 Apr 2018 15:42:35 -0700 Subject: Initial commit --- CSMath/.classpath | 6 + CSMath/.gitignore | 1 + CSMath/.project | 17 + CSMath/.settings/org.eclipse.jdt.core.prefs | 11 + CSMath/src/CulkinAsssignmentNine.java | 215 +++++++++++ CSMath/src/bisection/Bisection.java | 157 ++++++++ CSMath/src/bisection/Dual.java | 55 +++ CSMath/src/bisection/DualExpr.java | 269 ++++++++++++++ CSMath/src/bisection/NeoBisection.java | 551 ++++++++++++++++++++++++++++ CSMath/src/bisection/PreBisection.java | 74 ++++ 10 files changed, 1356 insertions(+) create mode 100644 CSMath/.classpath create mode 100644 CSMath/.gitignore create mode 100644 CSMath/.project create mode 100644 CSMath/.settings/org.eclipse.jdt.core.prefs create mode 100644 CSMath/src/CulkinAsssignmentNine.java create mode 100644 CSMath/src/bisection/Bisection.java create mode 100644 CSMath/src/bisection/Dual.java create mode 100644 CSMath/src/bisection/DualExpr.java create mode 100644 CSMath/src/bisection/NeoBisection.java create mode 100644 CSMath/src/bisection/PreBisection.java diff --git a/CSMath/.classpath b/CSMath/.classpath new file mode 100644 index 0000000..e461bea --- /dev/null +++ b/CSMath/.classpath @@ -0,0 +1,6 @@ + + + + + + diff --git a/CSMath/.gitignore b/CSMath/.gitignore new file mode 100644 index 0000000..ae3c172 --- /dev/null +++ b/CSMath/.gitignore @@ -0,0 +1 @@ +/bin/ diff --git a/CSMath/.project b/CSMath/.project new file mode 100644 index 0000000..453ea31 --- /dev/null +++ b/CSMath/.project @@ -0,0 +1,17 @@ + + + CSMath + + + + + + org.eclipse.jdt.core.javabuilder + + + + + + org.eclipse.jdt.core.javanature + + diff --git a/CSMath/.settings/org.eclipse.jdt.core.prefs b/CSMath/.settings/org.eclipse.jdt.core.prefs new file mode 100644 index 0000000..bb35fa0 --- /dev/null +++ b/CSMath/.settings/org.eclipse.jdt.core.prefs @@ -0,0 +1,11 @@ +eclipse.preferences.version=1 +org.eclipse.jdt.core.compiler.codegen.inlineJsrBytecode=enabled +org.eclipse.jdt.core.compiler.codegen.targetPlatform=1.8 +org.eclipse.jdt.core.compiler.codegen.unusedLocal=preserve +org.eclipse.jdt.core.compiler.compliance=1.8 +org.eclipse.jdt.core.compiler.debug.lineNumber=generate +org.eclipse.jdt.core.compiler.debug.localVariable=generate +org.eclipse.jdt.core.compiler.debug.sourceFile=generate +org.eclipse.jdt.core.compiler.problem.assertIdentifier=error +org.eclipse.jdt.core.compiler.problem.enumIdentifier=error +org.eclipse.jdt.core.compiler.source=1.8 diff --git a/CSMath/src/CulkinAsssignmentNine.java b/CSMath/src/CulkinAsssignmentNine.java new file mode 100644 index 0000000..da4290d --- /dev/null +++ b/CSMath/src/CulkinAsssignmentNine.java @@ -0,0 +1,215 @@ +import java.awt.BorderLayout; +import java.awt.Component; +import java.awt.GridLayout; +import java.util.Arrays; + +import javax.swing.DefaultListModel; +import javax.swing.JButton; +import javax.swing.JFrame; +import javax.swing.JLabel; +import javax.swing.JList; +import javax.swing.JPanel; +import javax.swing.JScrollPane; +import javax.swing.JSplitPane; +import javax.swing.ListCellRenderer; + +public class CulkinAsssignmentNine { + public static void main(String[] args) { + JFrame fram = new JFrame(); + fram.setLayout(new GridLayout(1, 1)); + + JPanel canvas = new JPanel(); + canvas.setLayout(new GridLayout(1, 1)); + + JPanel points = new JPanel(); + points.setLayout(new BorderLayout()); + + JPanel listPanel = new JPanel(); + listPanel.setLayout(new GridLayout(1, 1)); + + DefaultListModel pointModel = new DefaultListModel<>(); + + JList pointList = new JList<>(pointModel); + pointList.setCellRenderer(new TDPointRenderer()); + + JScrollPane listScroller = new JScrollPane(pointList); + + listPanel.add(listScroller); + + JPanel buttonPanel = new JPanel(); + buttonPanel.setLayout(new GridLayout(1, 2)); + + JButton addPoint = new JButton("Add Control Point"); + + JButton remPoint = new JButton("Remove Control Point"); + + buttonPanel.add(addPoint); + buttonPanel.add(remPoint); + + points.add(listPanel, BorderLayout.CENTER); + points.add(buttonPanel, BorderLayout.PAGE_END); + + JSplitPane main = new JSplitPane(JSplitPane.HORIZONTAL_SPLIT, points, canvas); + + fram.add(main); + + fram.setSize(640, 480); + fram.setDefaultCloseOperation(JFrame.DISPOSE_ON_CLOSE); + fram.setVisible(true); + } +} + +class Bezier { + public final TDPoint[] controls; + + public Bezier(TDPoint... points) { + controls = points; + } + + public TDPoint eval(double t) { + return evalIntern(t, controls.length, 0); + } + + private TDPoint evalIntern(double t, int j, int k) { + if (j == 0) + return controls[k]; + + TDPoint l = evalIntern(t, j - 1, k); + TDPoint r = evalIntern(t, j - 1, k + 1); + + l = l.multiply(1 - t); + r = r.multiply(t); + + return TDPoint.add(l, r); + } + + public Bezier[] decompose(double t) { + Bezier[] curves = new Bezier[2]; + + { + TDPoint[] points = new TDPoint[controls.length]; + for (int i = 0; i < points.length; i++) { + points[i] = evalIntern(t, i, 0); + } + curves[0] = new Bezier(points); + } + + { + TDPoint[] points = new TDPoint[controls.length]; + for (int i = 0; i < points.length; i++) { + points[i] = evalIntern(t, (points.length - 1) - i, i); + } + curves[1] = new Bezier(points); + } + + return curves; + } + + @Override + public int hashCode() { + final int prime = 31; + int result = 1; + result = prime * result + Arrays.hashCode(controls); + return result; + } + + @Override + public boolean equals(Object obj) { + if (this == obj) + return true; + if (obj == null) + return false; + if (getClass() != obj.getClass()) + return false; + Bezier other = (Bezier) obj; + if (!Arrays.equals(controls, other.controls)) + return false; + return true; + } + + @Override + public String toString() { + return "Bezier [controls=" + Arrays.toString(controls) + "]"; + } +} + +class TDPoint { + public final double x; + public final double y; + + public TDPoint(double x, double y) { + this.x = x; + this.y = y; + } + + public static TDPoint p2(double x, double y) { + return new TDPoint(x, y); + } + + public TDPoint multiply(double s) { + return p2(x * s, y * s); + } + + public static TDPoint add(TDPoint p1, TDPoint p2) { + return p2(p1.x + p2.x, p1.y + p2.y); + } + + @Override + public int hashCode() { + final int prime = 31; + int result = 1; + long temp; + temp = Double.doubleToLongBits(x); + result = prime * result + (int) (temp ^ (temp >>> 32)); + temp = Double.doubleToLongBits(y); + result = prime * result + (int) (temp ^ (temp >>> 32)); + return result; + } + + @Override + public boolean equals(Object obj) { + if (this == obj) + return true; + if (obj == null) + return false; + if (getClass() != obj.getClass()) + return false; + TDPoint other = (TDPoint) obj; + if (Double.doubleToLongBits(x) != Double.doubleToLongBits(other.x)) + return false; + if (Double.doubleToLongBits(y) != Double.doubleToLongBits(other.y)) + return false; + return true; + } + + @Override + public String toString() { + return "TDPoint [x=" + x + ", y=" + y + "]"; + } +} + +final class TDPointRenderer extends JLabel implements ListCellRenderer { + private static final long serialVersionUID = 629873168260730449L; + + public TDPointRenderer() { + setOpaque(true); + setHorizontalAlignment(CENTER); + setVerticalAlignment(CENTER); + } + + @Override + public Component getListCellRendererComponent(JList list, TDPoint value, int index, + boolean isSelected, boolean cellHasFocus) { + if (isSelected) { + setBackground(list.getSelectionBackground()); + setForeground(list.getSelectionForeground()); + } else { + setBackground(list.getBackground()); + setForeground(list.getForeground()); + } + + setText(String.format("(%.2d, %.2d)", value.x, value.y)); + + return this; + } +} \ No newline at end of file diff --git a/CSMath/src/bisection/Bisection.java b/CSMath/src/bisection/Bisection.java new file mode 100644 index 0000000..edca743 --- /dev/null +++ b/CSMath/src/bisection/Bisection.java @@ -0,0 +1,157 @@ +package bisection; +/* + * Benjamin Culkin + * 1/16/2018 + * CS 479 + * Bisection + * + * Use the bisection method to find bracketed roots of arbitrary real-valued functions. + */ + +public class Bisection { + /** + * Maximum number of iterations to attempt. + */ + public static final int MAXITR = 500; + + public static void main(String[] args) { + /* 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(DualExpr.ExprType.MULTIPLICATION, varXExpr, varXExpr), + new DualExpr(new Dual(3))); + + /* 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", 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); + } + + /** + * 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 = prevmid - (res.real / res.dual); + + /* 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; + } + + 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 top1 = guess1 * res2.real; + double top2 = guess2 * res1.real; + + double top = top1 - top2; + + double bot = res2.real - res1.real; + + /* Use the secant method to refine our guesses. */ + guess1 = guess2; + guess2 = top / bot; + } + + /* 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; + } +} diff --git a/CSMath/src/bisection/Dual.java b/CSMath/src/bisection/Dual.java new file mode 100644 index 0000000..de5b004 --- /dev/null +++ b/CSMath/src/bisection/Dual.java @@ -0,0 +1,55 @@ +package bisection; + +/** + * Represents a 'dual' number. + * + * Think imaginary numbers, where instead of i, we add a value d such that d^2 = + * 0. + */ +public class Dual { + /** + * The real part of the dual number. + */ + public double real; + /** + * The dual part of the dual number. + */ + public double dual; + + /** + * Create a new dual with both parts zero. + */ + public Dual() { + real = 0; + dual = 0; + } + + /** + * Create a new dual number with a zero dual part. + * + * @param real + * The real part of the number. + */ + public Dual(double real) { + this.real = real; + this.dual = 0; + } + + /** + * Create a new dual number with a specified dual part. + * + * @param real + * The real part of the number. + * @param dual + * The dual part of the number. + */ + public Dual(double real, double dual) { + this.real = real; + this.dual = dual; + } + + @Override + public String toString() { + return String.format("<%f, %f>", real, dual); + } +} \ No newline at end of file diff --git a/CSMath/src/bisection/DualExpr.java b/CSMath/src/bisection/DualExpr.java new file mode 100644 index 0000000..0e16b3a --- /dev/null +++ b/CSMath/src/bisection/DualExpr.java @@ -0,0 +1,269 @@ +package bisection; +/** + * Represents an expression using dual numbers. + * + * Useful for automatically differentiating expressions. + */ +public class DualExpr { + /** + * Represents the various types of dual expressions. + */ + public static enum ExprType { + /** + * A fixed number. + */ + CONSTANT, + /** + * An addition operation. + */ + ADDITION, + /** + * A subtraction operation. + */ + SUBTRACTION, + /** + * A multiplication operation. + */ + MULTIPLICATION, + /** + * A division operation. + */ + DIVISION, + /** + * A sine operation. + */ + SIN, + /** + * A cosine operation. + */ + COS, + /** + * An exponential function. + */ + EXPONENTIAL, + /** + * A logarithm function. + */ + LOGARITHM, + /** + * A power operation. + */ + POWER, + /** + * An absolute value. + */ + ABSOLUTE + } + + /** + * The type of the expression. + */ + public final DualExpr.ExprType type; + + /** + * The dual number value, for constants. + */ + public Dual number; + + /** + * The left (or first) part of the expression. + */ + public DualExpr left; + /** + * The right (or second) part of the expression. + */ + public DualExpr right; + + /** + * The power to use, for power operations. + */ + public int power; + + /** + * Create a new constant dual number. + * + * @param num + * The value of the dual number. + */ + public DualExpr(Dual num) { + this.type = ExprType.CONSTANT; + + number = num; + } + + /** + * Create a new unary dual number. + * + * @param type + * The type of operation to perform. + * @param val + * The parameter to the value. + */ + public DualExpr(DualExpr.ExprType type, DualExpr val) { + this.type = type; + + left = val; + } + + /** + * Create a new binary dual number. + * + * @param type + * The type of operation to perform. + * @param val + * The parameter to the value. + */ + public DualExpr(DualExpr.ExprType type, DualExpr left, DualExpr right) { + this.type = type; + + this.left = left; + this.right = right; + } + + /** + * Create a new power expression. + * + * @param left + * The expression to raise. + * @param power + * The power to raise it by. + */ + public DualExpr(DualExpr left, int power) { + this.type = ExprType.POWER; + + this.left = left; + this.power = power; + } + + /** + * Evaluate an expression to a number. + * + * Uses the rules provided in + * https://en.wikipedia.org/wiki/Automatic_differentiation + * + * @return The evaluated expression. + */ + public Dual evaluate() { + /* The evaluated dual numbers. */ + Dual lval, rval; + + /* Perform the right operation for each type. */ + switch (type) { + case CONSTANT: + return number; + case ADDITION: + lval = left.evaluate(); + rval = right.evaluate(); + + return new Dual(lval.real + rval.real, lval.dual + rval.dual); + case SUBTRACTION: + lval = left.evaluate(); + rval = right.evaluate(); + + return new Dual(lval.real - rval.real, lval.dual - rval.dual); + case MULTIPLICATION: + lval = left.evaluate(); + rval = right.evaluate(); + + { + double lft = lval.dual * rval.real; + double rght = lval.real * rval.dual; + + return new Dual(lval.real * rval.real, lft + rght); + } + case DIVISION: + lval = left.evaluate(); + rval = right.evaluate(); + + { + if (rval.real == 0) { + throw new IllegalArgumentException("ERROR: Attempted to divide by zero."); + } + + double lft = lval.dual * rval.real; + double rght = lval.real * rval.dual; + + double val = (lft - rght) / (rval.real * rval.real); + + return new Dual(lval.real / rval.real, val); + } + case SIN: + lval = left.evaluate(); + + return new Dual(Math.sin(lval.real), lval.dual * Math.cos(lval.real)); + case COS: + lval = left.evaluate(); + + return new Dual(Math.cos(lval.real), -lval.dual * Math.sin(lval.real)); + case EXPONENTIAL: + lval = left.evaluate(); + + { + double val = Math.exp(lval.real); + + return new Dual(val, lval.dual * val); + } + case LOGARITHM: + lval = left.evaluate(); + + if (lval.real <= 0) { + throw new IllegalArgumentException( + "ERROR: Attempted to take non-positive log."); + } + + return new Dual(Math.log(lval.real), lval.dual / lval.real); + case POWER: + lval = left.evaluate(); + + if (lval.real == 0) { + throw new IllegalArgumentException("ERROR: Raising zero to a power."); + } + + { + double rl = Math.pow(lval.real, power); + + double lft = Math.pow(lval.real, power - 1); + + return new Dual(rl, power * lft * lval.dual); + } + case ABSOLUTE: + lval = left.evaluate(); + + return new Dual(Math.abs(lval.real), lval.dual * Math.signum(lval.real)); + default: + String msg = "ERROR: Unknown expression type %s"; + + throw new IllegalArgumentException(String.format(msg, type)); + } + } + + @Override + public String toString() { + switch (type) { + case ABSOLUTE: + return String.format("abs(%s)", left.toString()); + case ADDITION: + return String.format("(%s + %s)", left.toString(), right.toString()); + case CONSTANT: + return String.format("%s", number.toString()); + case COS: + return String.format("cos(%s)", left.toString()); + case DIVISION: + return String.format("(%s / %s)", left.toString(), right.toString()); + case EXPONENTIAL: + return String.format("exp(%s)", left.toString()); + case LOGARITHM: + return String.format("log(%s)", left.toString()); + case MULTIPLICATION: + return String.format("(%s * %s)", left.toString(), right.toString()); + case POWER: + return String.format("(%s ^ %d)", left.toString(), power); + case SIN: + return String.format("sin(%s)", left.toString()); + case SUBTRACTION: + return String.format("(%s - %s)", left.toString(), right.toString()); + default: + return String.format("UNKNOWN_EXPR"); + } + } +} \ No newline at end of file diff --git a/CSMath/src/bisection/NeoBisection.java b/CSMath/src/bisection/NeoBisection.java new file mode 100644 index 0000000..73e3797 --- /dev/null +++ b/CSMath/src/bisection/NeoBisection.java @@ -0,0 +1,551 @@ +package bisection; +import java.util.function.DoubleUnaryOperator; + +/* + * Benjamin Culkin + * 1/16/2018 + * CS 479 + * Bisection + * + * Use the bisection method to find bracketed roots of arbitrary real-valued functions. + */ + +public class NeoBisection { + /** + * Maximum number of iterations to attempt. + */ + public static final int MAXITR = 500; + + 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; + } + + 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; + } + + /** + * Represents a 'dual' number. + * + * Think imaginary numbers, where instead of i, we add a value d such that d^2 = + * 0. + */ + public static class Dual { + /** + * The real part of the dual number. + */ + public double real; + /** + * The dual part of the dual number. + */ + public double dual; + + /** + * Create a new dual with both parts zero. + */ + public Dual() { + real = 0; + dual = 0; + } + + /** + * Create a new dual number with a zero dual part. + * + * @param real + * The real part of the number. + */ + public Dual(double real) { + this.real = real; + this.dual = 0; + } + + /** + * Create a new dual number with a specified dual part. + * + * @param real + * The real part of the number. + * @param dual + * The dual part of the number. + */ + public Dual(double real, double dual) { + this.real = real; + this.dual = dual; + } + + @Override + public String toString() { + return String.format("<%f, %f>", real, dual); + } + } + + /** + * Represents an expression using dual numbers. + * + * Useful for automatically differentiating expressions. + */ + public static class DualExpr { + /** + * Represents the various types of dual expressions. + */ + public static enum ExprType { + /** + * A fixed number. + */ + CONSTANT, + /** + * An addition operation. + */ + ADDITION, + /** + * A subtraction operation. + */ + SUBTRACTION, + /** + * A multiplication operation. + */ + MULTIPLICATION, + /** + * A division operation. + */ + DIVISION, + /** + * A sine operation. + */ + SIN, + /** + * A cosine operation. + */ + COS, + /** + * An exponential function. + */ + EXPONENTIAL, + /** + * A logarithm function. + */ + LOGARITHM, + /** + * A power operation. + */ + POWER, + /** + * An absolute value. + */ + ABSOLUTE + } + + /** + * The type of the expression. + */ + public final ExprType type; + + /** + * The dual number value, for constants. + */ + public Dual number; + + /** + * The left (or first) part of the expression. + */ + public DualExpr left; + /** + * The right (or second) part of the expression. + */ + public DualExpr right; + + /** + * The power to use, for power operations. + */ + public int power; + + /** + * Create a new constant dual number. + * + * @param num + * The value of the dual number. + */ + public DualExpr(Dual num) { + this.type = ExprType.CONSTANT; + + number = num; + } + + /** + * Create a new unary dual number. + * + * @param type + * The type of operation to perform. + * @param val + * The parameter to the value. + */ + public DualExpr(ExprType type, DualExpr val) { + this.type = type; + + left = val; + } + + /** + * Create a new binary dual number. + * + * @param type + * The type of operation to perform. + * @param val + * The parameter to the value. + */ + public DualExpr(ExprType type, DualExpr left, DualExpr right) { + this.type = type; + + this.left = left; + this.right = right; + } + + /** + * Create a new power expression. + * + * @param left + * The expression to raise. + * @param power + * The power to raise it by. + */ + public DualExpr(DualExpr left, int power) { + this.type = ExprType.POWER; + + this.left = left; + this.power = power; + } + + /** + * Evaluate an expression to a number. + * + * Uses the rules provided in + * https://en.wikipedia.org/wiki/Automatic_differentiation + * + * @return The evaluated expression. + */ + public Dual evaluate() { + /* The evaluated dual numbers. */ + Dual lval, rval; + + /* Perform the right operation for each type. */ + switch (type) { + case CONSTANT: + return number; + case ADDITION: + lval = left.evaluate(); + rval = right.evaluate(); + + return new Dual(lval.real + rval.real, lval.dual + rval.dual); + case SUBTRACTION: + lval = left.evaluate(); + rval = right.evaluate(); + + return new Dual(lval.real - rval.real, lval.real - rval.real); + case MULTIPLICATION: + lval = left.evaluate(); + rval = right.evaluate(); + + { + double lft = lval.dual * rval.real; + double rght = lval.real * rval.dual; + + return new Dual(lval.real * rval.real, lft + rght); + } + case DIVISION: + lval = left.evaluate(); + rval = right.evaluate(); + + { + if (rval.real == 0) { + throw new IllegalArgumentException("ERROR: Attempted to divide by zero."); + } + + double lft = lval.dual * rval.real; + double rght = lval.real * rval.dual; + + double val = (lft - rght) / (rval.real * rval.real); + + return new Dual(lval.real / rval.real, val); + } + case SIN: + lval = left.evaluate(); + + return new Dual(Math.sin(lval.real), lval.dual * Math.cos(lval.real)); + case COS: + lval = left.evaluate(); + + return new Dual(Math.cos(lval.real), -lval.dual * Math.sin(lval.real)); + case EXPONENTIAL: + lval = left.evaluate(); + + { + double val = Math.exp(lval.real); + + return new Dual(val, lval.dual * val); + } + case LOGARITHM: + lval = left.evaluate(); + + if (lval.real <= 0) { + throw new IllegalArgumentException( + "ERROR: Attempted to take non-positive log."); + } + + return new Dual(Math.log(lval.real), lval.dual / lval.real); + case POWER: + lval = left.evaluate(); + + if (lval.real == 0) { + throw new IllegalArgumentException("ERROR: Raising zero to a power."); + } + + { + double rl = Math.pow(lval.real, power); + + double lft = Math.pow(lval.real, power - 1); + + return new Dual(rl, power * lft * lval.dual); + } + case ABSOLUTE: + lval = left.evaluate(); + + return new Dual(Math.abs(lval.real), lval.dual * Math.signum(lval.real)); + default: + String msg = "ERROR: Unknown expression type %s"; + + throw new IllegalArgumentException(String.format(msg, type)); + } + } + + @Override + public String toString() { + switch (type) { + case ABSOLUTE: + return String.format("abs(%s)", left.toString()); + case ADDITION: + return String.format("(%s + %s)", left.toString(), right.toString()); + case CONSTANT: + return String.format("%s", number.toString()); + case COS: + return String.format("cos(%s)", left.toString()); + case DIVISION: + return String.format("(%s / %s)", left.toString(), right.toString()); + case EXPONENTIAL: + return String.format("exp(%s)", left.toString()); + case LOGARITHM: + return String.format("log(%s)", left.toString()); + case MULTIPLICATION: + return String.format("(%s * %s)", left.toString(), right.toString()); + case POWER: + return String.format("(%s ^ %d)", left.toString(), power); + case SIN: + return String.format("sin(%s)", left.toString()); + case SUBTRACTION: + return String.format("(%s - %s)", left.toString(), right.toString()); + default: + return String.format("UNKNOWN_EXPR"); + } + } + } +} diff --git a/CSMath/src/bisection/PreBisection.java b/CSMath/src/bisection/PreBisection.java new file mode 100644 index 0000000..39bbd60 --- /dev/null +++ b/CSMath/src/bisection/PreBisection.java @@ -0,0 +1,74 @@ +package bisection; +/* + * Benjamin Culkin + * 1/16/2018 + * CS 479 + * Bisection + * + * Use the bisection method to find bracketed roots of arbitrary real-valued functions. + */ +import java.util.function.DoubleUnaryOperator; + +public class PreBisection { + public static void main(String[] args) { + // 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("x = cos(x) => %.4f\n", answerA); + System.out.printf("x^5 - 3x^4 + 4x^2 - 1 => %.4f\n", answerB); + System.out.printf("x^2 - 3 => %.7f\n", answerC); + } + + // 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; + } +} -- cgit v1.2.3