001 package physics3d;
002
003 public strictfp class ZeroFinder {
004
005 /**
006 * ZeroFinder.Function is an interface which specifies a function whose roots
007 * can be found by this class.
008 */
009 public static interface Function {
010 public abstract double evaluate(double t);
011 }
012
013 /**
014 * @requires t_min <= t_max, stepSize > 0
015 *
016 * Searches for possible roots of <code>function</code> by looking for sign
017 * changes. If a possible root is found, regula falsi is performed to try to
018 * determine the precise value of the root. This function searches for roots
019 * between <code>t_min</code> and <code>t_max</code>. Will return the
020 * smallest zero within the first interval of length stepSize where function
021 * changes sign. If no solution is found returns <code>NaN</code>, else
022 * returns the solution.
023 *
024 * To improve speed, use smaller intervals or increase the step size. However,
025 * increasing the step size also increases the chance that the root you get is
026 * not the smallest. Also increasing the step size could decreases the
027 * probability that a root will be found, as less points are tested.
028 */
029 public static double findRoot(Function function, double t_min, double t_max,
030 int subintervals) {
031 double stepSize = (t_max - t_min) / subintervals;
032
033 if (Vect3.isAbout(function.evaluate(0), 0)) {
034 return 0;
035 }
036 double t = findInterval(function, t_min, t_max, stepSize);
037
038 if (t > t_max) {
039 return Double.NaN;
040 }
041 t = mullersMethod(function, t - stepSize, t);
042
043 return t;
044 }
045
046 // finds subinterval [a, b]of [t_min, t_max] of length stepSize
047 // with the function values at a and b opposite sign
048 private static double findInterval(Function function, double t_min,
049 double t_max, double stepSize) {
050 double f_told = 0;
051 double f_tnew = function.evaluate(t_min);
052 double t = t_min;
053
054 boolean foundChange = false;
055
056 while ((t <= t_max + .00000000000001) && !foundChange) {
057 t = t + stepSize;
058
059 f_told = f_tnew;
060 f_tnew = function.evaluate(t);
061 if (f_told * f_tnew < 0) {
062 foundChange = true;
063 }
064 }
065 return t;
066 }
067
068 // requires function values at a and b opposite
069 // and a < b
070 // finds a zero in the interval a, b of function f
071 private static double regulaFalsi(Function function, double a, double b) {
072 double min = a;
073 double max = b;
074 double middle = (a + b) / 2;
075
076 double fmin = function.evaluate(min);
077 double fmax = function.evaluate(max);
078 double fmiddle = function.evaluate(middle);
079
080 int trials = 0;
081 int maxTrials = 30;
082 boolean leftSide = false;
083 while (Math.abs(min - max) > Math.pow(10, -15) && trials < maxTrials) {
084 trials++;
085 middle = max - fmax / ((fmax - fmin) / (max - min));
086 //System.out.println(middle);
087
088 fmiddle = function.evaluate(middle);
089 if (fmiddle * fmin < 0) {
090 max = middle;
091 fmax = fmiddle;
092 if (leftSide) {
093 fmin = fmin / 2;
094 }
095 leftSide = true;
096 } else {
097 min = middle;
098 fmin = fmiddle;
099 if (!leftSide) {
100 fmax = fmax / 2;
101 }
102 leftSide = false;
103 }
104 }
105 return middle;
106 }
107
108 // requires function values at a and b opposite
109 // and a < b
110 // finds a zero in the interval a, b of function f
111 private static double mullersMethod(Function function, double a, double b) {
112 double p1 = a;
113 double p3 = b;
114 double p2 = (a + b) / 2;
115 double next = 0;
116
117 double fp1 = function.evaluate(p1);
118 double fp2 = function.evaluate(p2);
119 double fp3 = function.evaluate(p3);
120 double fnext = function.evaluate(next);
121 double q = 0;
122 double A = 0;
123 double B = 0;
124 double C = 0;
125 double disc = 0;
126
127 int trials = 0;
128 int maxTrials = 30;
129 while (Math.abs(p3 - p2) > Math.pow(10, -15) && trials < maxTrials) {
130 trials++;
131 q = (p3 - p2) / (p2 - p1);
132 A = q * fp3 - q * (1 + q) * fp2 + q * q * fp1;
133 B = (2 * q + 1) * fp3 - (1 + q) * (1 + q) * fp2 + q * q * fp1;
134 C = (1 + q) * fp3;
135
136 disc = 0;
137 if (B * B - 4 * A * C > 0)
138 disc = Math.sqrt(B * B - 4 * A * C);
139 if (B < 0)
140 disc = -disc;
141
142 if (B + disc == 0)
143 return 0;
144
145 next = p3 - (p3 - p2) * 2 * C / (B + disc);
146 fnext = function.evaluate(next);
147 p1 = p2;
148 fp1 = fp2;
149 p2 = p3;
150 fp2 = fp3;
151 p3 = next;
152 fp3 = fnext;
153 }
154 return next;
155 }
156
157 }