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 }