001/** 002 * Roots -- A collection of static methods for finding the roots of functions. 003 * 004 * Copyright (C) 1999-2015, by Joseph A. Huwaldt. All rights reserved. 005 * 006 * This library is free software; you can redistribute it and/or modify it under the terms 007 * of the GNU Lesser General Public License as published by the Free Software Foundation; 008 * either version 2.1 of the License, or (at your option) any later version. 009 * 010 * This library is distributed in the hope that it will be useful, but WITHOUT ANY 011 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 012 * PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. 013 * 014 * You should have received a copy of the GNU Lesser General Public License along with 015 * this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - 016 * Suite 330, Boston, MA 02111-1307, USA. Or visit: http://www.gnu.org/licenses/lgpl.html 017 */ 018package jahuwaldt.tools.math; 019 020import static java.lang.Math.*; 021 022/** 023 * A collection of static routines to find the roots of functions or sets of functions. 024 * 025 * <p> Modified by: Joseph A. Huwaldt </p> 026 * 027 * @author Joseph A. Huwaldt, Date: October 8, 1997 028 * @version September 11, 2015 029 */ 030public class Roots { 031 032 // Debug flag. 033 private static final boolean DEBUG = false; 034 035 /** 036 * Machine floating point precision. 037 */ 038 private static final double EPS = MathTools.EPS; 039 040 /** 041 * Maximum number of allowed iterations. 042 */ 043 private static final int MAXIT = 1000; 044 045 // Constants used by newtonND. 046 private static final double TOLF = 1.0e-8, TOLMIN = 1.0e-12, STPMX = 100.0; 047 private static final double TINY = 1.0e-20; 048 private static final double ALF = 1.0e-4; 049 050 //----------------------------------------------------------------------------------- 051 /** 052 * Given a function, <code>func</code>, defined on the interval <code>x1</code> to 053 * <code>x2</code>, this routine subdivides the interval into <code>n</code> equally 054 * spaced segments, and searches for zero crossings of the function. Brackets around 055 * any zero crossings found are returned. 056 * 057 * @param func The function that is being search for zero crossings. 058 * @param x1 The start of the interval to be searched. 059 * @param x2 The end of the interval to be searched. 060 * @param n The number segments to divide the interval into. 061 * @param xb1 Lower value of each bracketing pair found (number of pairs is returned 062 * by function). The size of the array determines the maximum number of 063 * bracketing pairs that will be found. 064 * @param xb2 Upper value of each bracketing pair found (number of pairs is returned 065 * by function). The size of this array must match the size of 066 * <code>xb1</code>. 067 * @return The number of bracketing pairs found. 068 * @throws jahuwaldt.tools.math.RootException if the supplied function can not accept 069 * the interval provided. 070 */ 071 public static int bracket(Evaluatable1D func, double x1, double x2, int n, double[] xb1, double[] xb2) throws RootException { 072 int nb = xb1.length; 073 int nbb = 0; 074 double dx = (x2 - x1) / n; 075 if (dx == 0) 076 return 0; 077 double x = x1; 078 double fp = func.function(x); 079 for (int i = 0; i < n; ++i) { 080 x += dx; 081 if (x > x2) 082 x = x2; 083 double fc = func.function(x); 084 if (fc * fp <= 0.0) { 085 xb1[nbb] = x - dx; 086 xb2[nbb++] = x; 087 if (nb == nbb) 088 return nb; 089 } 090 fp = fc; 091 } 092 return nbb; 093 } 094 095 /** 096 * Find the root of a general 1D function <code>f(x) = 0</code> known to lie between 097 * <code>x1</code> and <code>x2</code>. The root will be refined until it's accuracy 098 * is <code>tol</code>. 099 * <p> 100 * Have your Evaluatable1D <code>derivative()</code> function return 101 * <code>Double.NaN</code> when you want this routine to use Brent's method. 102 * Newton-Raphson will be used if a derivative function is provided that returns 103 * something other than <code>Double.NaN</code>. The <code>function()</code> method 104 * will always be called before the <code>derivative()</code> method. The value 105 * returned from the 1st call to the <code>function()</code> method is ignored. 106 * </p> 107 * 108 * @param eval An evaluatable 1D function that returns the function value at x and 109 * optionally returns the function derivative value (d(fx)/dx) at x. It 110 * is guaranteed that "function()" will always be called before 111 * "derivative()". 112 * @param bracket The upper and lower bracket surrounding the root (assumed that root 113 * lies inside the bracket). 114 * @param tol The root will be refined until it's accuracy is better than this. 115 * @return The value of x that solves the equation f(x) = 0. 116 * @exception RootException if unable to find a root of the function. 117 */ 118 public static double findRoot1D(Evaluatable1D eval, BracketRoot1D bracket, double tol) throws RootException { 119 return findRoot1D(eval, bracket.x1, bracket.x2, tol); 120 } 121 122 /** 123 * Find the root of a general 1D function <code>f(x) = 0</code> known to lie between 124 * <code>x1</code> and <code>x2</code>. The root will be refined until it's accuracy 125 * is <code>tol</code>. 126 * <p> 127 * Have your Evaluatable1D <code>derivative()</code> function return 128 * <code>Double.NaN</code> when you want this routine to use Brent's method. 129 * Newton-Raphson will be used if a derivative function is provided that returns 130 * something other than <code>Double.NaN</code>. The <code>function()</code> method 131 * will always be called before the <code>derivative()</code> method. The value 132 * returned from the 1st call to the <code>function()</code> method is ignored. 133 * </p> 134 * 135 * @param eval An evaluatable 1D function that returns the function value at x and 136 * optionally returns the function derivative value (d(fx)/dx) at x. It is 137 * guaranteed that "function()" will always be called before 138 * "derivative()". 139 * @param x1 The lower bracket surrounding the root (assumed that root lies between 140 * x1 and x2). 141 * @param x2 The upper bracket surrounding the root (assumed that root lies between 142 * x1 and x2). 143 * @param tol The root will be refined until it's accuracy is better than this. 144 * @return The value of x that solves the equation f(x) = 0. 145 * @exception RootException Unable to find a root of the function. 146 */ 147 public static double findRoot1D(Evaluatable1D eval, double x1, 148 double x2, double tol) throws RootException { 149 if (tol < EPS) 150 tol = EPS; 151 eval.function(x1); 152 double value = eval.derivative(x1); 153 154 if (Double.isNaN(value)) 155 // No derivative information available, use Brent's method. 156 value = zeroin(eval, x1, x2, tol); 157 158 else 159 // Derivative info is available, use Newton-Raphson. 160 value = newton(eval, x1, x2, tol); 161 162 return value; 163 } 164 165 /** 166 * Find the root of a function <code>f(x) = 0</code> known to lie between 167 * <code>x1</code> and <code>x2</code>. The root will be refined until it's accuracy 168 * is <code>tol</code>. This is the method of choice to find a bracketed root of a 169 * general 1D function when you can not easily compute the function's derivative. 170 * <p> 171 * Uses Van Wijngaarden-Dekker-Brent method (Brent's method). Algorithm: 172 * <pre> 173 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 174 * computations. M., Mir, 1980, p.180 of the Russian edition. 175 * </pre></p> 176 * <p> 177 * Ported to Java from Netlib C version by: Joseph A. Huwaldt, July 10, 2000 </p> 178 * <p> 179 * This function makes use of the bisection procedure combined with linear or inverse 180 * quadric interpolation. At every step the program operates on three abscissae - a, 181 * b, and c.<pre> 182 * b - the last and the best approximation to the root, 183 * a - the last but one approximation, 184 * c - the last but one or even earlier approximation than a that 185 * 1) |f(b)| ≤ |f(c)|, 186 * 2) f(b) and f(c) have opposite signs, i.e. b and c confine 187 * the root. 188 * </pre> 189 * At every step zeroin() selects one of the two new approximations, the former 190 * being obtained by the bisection procedure and the latter resulting in the 191 * interpolation (if a,b, and c are all different the quadric interpolation is 192 * utilized, otherwise the linear one). If the latter (i.e. obtained by the 193 * interpolation) point is reasonable (i.e. lies within the current interval [b,c] not 194 * being too close to the boundaries) it is accepted. The bisection result is used 195 * otherwise. Therefore, the range of uncertainty is ensured to be reduced at least by 196 * the factor 1.6 on each iteration. 197 * </p> 198 * 199 * @param eval An Evaluatable1D function that returns the function value at x. 200 * @param x1, x2 The bracket surrounding the root (assumed that root lies between x1 201 * and x2). 202 * @param tol The root will be refined until it's accuracy is better than this. 203 * @return The value x that solves the equation f(x) = 0. 204 * @exception RootException if unable to find a root of the function. 205 */ 206 private static double zeroin(Evaluatable1D eval, double x1, double x2, double tol) 207 throws RootException { 208 double a = x1, b = x2, c = x1; // Abscissae, descr. see above. 209 double fa = eval.function(a); // f(a) 210 double fb = eval.function(b); // f(b) 211 double fc = fa; // f(c) 212 213 if (DEBUG) 214 System.out.println("Using zeroin()..."); 215 216 // Check for bracketing. 217 if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 & fb < 0.0)) 218 throw new RootException("Root must be bracketed"); 219 220 // Main iteration loop 221 for (int iter = 0; iter < MAXIT; ++iter) { 222 if (DEBUG) 223 System.out.println("Iteration #" + iter); 224 225 double prev_step = b - a; // Distance from the last but one to the last approx. 226 double tol_act; // Actual tolerance 227 double p; // Interp step is calculated in the form p/q; 228 double q; // division is delayed until the last moment. 229 double new_step; // Step at this iteration 230 231 if (abs(fc) < abs(fb)) { 232 a = b; // Swap data for b to be the best approximation. 233 b = c; 234 c = a; 235 fa = fb; 236 fb = fc; 237 fc = fa; 238 } 239 240 // Convergence check. 241 tol_act = 2. * EPS * abs(b) + tol * 0.5; 242 new_step = 0.5 * (c - b); 243 244 if (abs(new_step) <= tol_act || fb == 0.0) 245 return b; // An acceptable approximation was found. 246 247 // Decide if the interpolation can be tried. 248 if (abs(prev_step) >= tol_act && abs(fa) > abs(fb)) { 249 // If prev_step was large enough and was in true direction, 250 // inverse quadratic interpolation may be tried. 251 double s = fb / fa; 252 double cb = c - b; 253 if (a == c) { // If we have only two distinct points then only 254 p = cb * s; // linear interpolation can be applied. 255 q = 1.0 - s; 256 } else { 257 // Quadric inverse interpolation. 258 double t = fa / fc; 259 double r = fb / fc; 260 p = s * (cb * t * (t - r) - (b - a) * (r - 1.0)); 261 q = (t - 1.0) * (r - 1.0) * (s - 1.0); 262 } 263 // Check wether in bounds 264 if (p > 0.0) // p was calculated with the opposite sign; 265 q = -q; // make p positive and assign possible minus to q. 266 else 267 p = -p; 268 269 double min1 = (0.75 * cb * q - abs(tol_act * q * 0.5)); 270 double min2 = abs(prev_step * q * 0.5); 271 if (p < min1 && p < min2) { 272 /* If b+p/q falls in [b,c] and isn't too large it is accepted. 273 If p/q is too large then use the bisection procedure which 274 will reduce the [b,c] range to a greater extent. */ 275 new_step = p / q; 276 } 277 } 278 279 if (abs(new_step) < tol_act) // Adjust the step to be not less than tolerance. 280 if (new_step > 0.0) 281 new_step = tol_act; 282 else 283 new_step = -tol_act; 284 285 a = b; // Save the previous approximation. 286 fa = fb; 287 b += new_step; 288 fb = eval.function(b); // Do step to a new approxim. 289 if ((fb > 0. && fc > 0.) || (fb < 0. && fc < 0.)) { 290 c = a; // Adjust c for it to have a sign opp. to that of b. 291 fc = fa; 292 } 293 } 294 295 // Max iterations exceeded. 296 throw new RootException("Maximun number of iterations exceeded"); 297 298 } 299 300 /** 301 * Find the root of a function <code>f(x) = 0</code> known to lie between 302 * <code>x1</code> and <code>x2</code>. The root will be refined until it's accuracy 303 * is known within <code>tol</code>. This is the method of choice to find a bracketed 304 * root of a general 1D function when you are able to supply the derivative of the 305 * function with x. 306 * <p> 307 * Uses a combination of Newton-Raphson and bisection. 308 * </p> 309 * <p> 310 * Original FORTRAN program, obtained from Netlib, bore this message: 311 * <pre> 312 * NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994. 313 * NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992, 314 * Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. 315 * This free software is complements of the author. 316 * Algorithm 2.5 (Newton-Raphson Iteration). 317 * Section 2.4, Newton-Raphson and Secant Methods, Page 84. 318 * </pre></p> 319 * <p> 320 * Ported from FORTRAN to Java by: Joseph A. Huwaldt, July 11, 2000. Root bracketing 321 * and bisection added to avoid pathologies as suggested by: Numerical Recipes in C, 322 * 2nd Edition, pg. 366. 323 * </p> 324 * 325 * @param eval An Evaluatable1D function that returns the function value at x and 326 * optionally returns the function derivative value (d(fx)/dx) at x. It is 327 * guaranteed that "function()" will always be called before 328 * "derivative()". 329 * @param x1, x2 The bracket surrounding the root (assumed that root lies between x1 330 * and x2). 331 * @param tol The root will be refined until it's accuracy is better than this. 332 * @return The value x that solves the equation f(x) = 0. 333 * @exception RootException if unable to find a root of the function. 334 */ 335 private static double newton(Evaluatable1D eval, double x1, 336 double x2, double tol) throws RootException { 337 338 if (DEBUG) 339 System.out.println("Using newton()..."); 340 341 // Evaluate the function low and high. 342 double fl = eval.function(x1); 343 double fh = eval.function(x2); 344 345 // Check for bracketing. 346 if ((fl > 0.0 && fh > 0.0) || (fl < 0.0 && fh < 0.0)) 347 throw new RootException("Root must be bracketed"); 348 349 // If either end of the bracket is the root, we are done. 350 if (fl == 0.0) 351 return x1; 352 if (fh == 0.0) 353 return x2; 354 355 double xh, xl; 356 if (fl < 0.0) { 357 // Orient the search so that fn(x1) < 0. 358 xl = x1; 359 xh = x2; 360 } else { 361 xh = x1; 362 xl = x2; 363 } 364 365 double p0 = 0.5 * (x1 + x2); // Initialize the guess for root. 366 double y0 = eval.function(p0); // Do initial function evaluation. 367 368 double dp_old = abs(x2 - x1); // Step size before last. 369 double dp = dp_old; // Last step size. 370 371 // Loop over allowed iterations. 372 for (int k = 0; k < MAXIT; ++k) { 373 if (DEBUG) 374 System.out.println("Iteration #" + k); 375 376 double df = eval.derivative(p0); // Evaluate the derivative. 377 378 if ((((p0 - xh) * df - y0) * ((p0 - xl) * df - y0) >= 0.0) 379 || (abs(2.0 * y0) > abs(dp_old * df))) { 380 // Use bisection if Newton is out of range or not decreasing fast enough. 381 dp_old = dp; 382 dp = 0.5 * (xh - xl); 383 p0 = xl + dp; 384 if (xl == p0) 385 return p0; // Change in root is negligable. 386 387 } else { 388 // Do a regular Newton step. 389 dp_old = dp; 390 dp = y0 / df; 391 double p1 = p0; 392 p0 = p0 - dp; 393 if (p1 == p0) 394 return p0; // Change in root is negligable. 395 } 396 397 // Do a convergence check. 398 if (abs(dp) < tol) 399 return p0; 400 401 // Do a function evaluation. 402 y0 = eval.function(p0); 403 if (y0 < 0.0) // Maintain the bracket on the root. 404 xl = p0; 405 else 406 xh = p0; 407 } 408 409 throw new RootException("Maximun number of iterations exceeded"); 410 411 } 412 413 /** 414 * Find a root of a 1D equation using the Aitken method. 415 * 416 * @param eval An evaluatable 1D function that returns the function value at x and 417 * optionally returns the function derivative value (d(fx)/dx) at x. This 418 * method never calls "derivative()". 419 * @param x A starting guess at the root value. 420 * @param tol The root will be refined until it's accuracy is better than this. 421 * @return The value of x that solves the equation f(x) = 0 to within the specified 422 * tolerance. 423 * @exception RootException if unable to find a root of the function. 424 */ 425 public static double aitken(Evaluatable1D eval, double x, double tol) throws RootException { 426 // Reference: http://www.tm.bi.ruhr-uni-bochum.de/profil/mitarbeiter/meyers/Aitken.html 427 if (tol < EPS) 428 tol = EPS; 429 double fx0 = 1e99; 430 double fac = 1.0; 431 int m = MAXIT; 432 int n = 0; 433 while (m > n) { 434 double fx1 = eval.function(x); 435 ++n; 436 double d = fx0 - fx1; 437 if (abs(d) < tol) { 438 x += fx1; 439 return x; 440 } 441 fac *= fx0 / d; 442 x += fac * fx1; 443 fx0 = fx1; 444 } 445 throw new RootException("Maximun number of iterations exceeded"); 446 } 447 448 /** 449 * Find the roots of a set of <code>N</code> non-linear equations in <code>N</code> 450 * variables. 451 * <p> 452 * Uses a globally convergent Newton's method with either a Jacobian supplied by the 453 * VectorFunction or one computed by differencing. 454 * </p> 455 * <p> Reference: Numerical Recipes in C, 2nd Edition, pg 386. </p> 456 * 457 * @param vecfunc A VectorFunction that computes the values for each equation using 458 * the input values in x[] and optionally the Jacobian. 459 * @param x A pre-existing array that contains the initial guesses at the x 460 * values. On completion, it will contain the roots of the equations. 461 * @param n The number of equations to be solved (and the number of elements in 462 * x. 463 * @return <code>false</code> if the routine completed normally or <code>true</code> 464 * if the routine has converged to a local minimum of the "fmin" function. In 465 * this case, try restarting from a different initial guess. 466 * @exception RootException Unable to find the roots of the functions. 467 */ 468 public static boolean findRootsND(VectorFunction vecfunc, double[] x, int n) throws RootException { 469 return newtonND(vecfunc, x, n); 470 } 471 472 private static boolean newtonND(VectorFunction vecfunc, double[] x, int n) throws RootException { 473 474 int[] indx = new int[n]; 475 double[] g = new double[n]; 476 double[] p = new double[n]; 477 double[] xold = new double[n]; 478 double[][] fjac = new double[n][n]; 479 double[] fvec = new double[n]; 480 481 double f = fmin(x, fvec, n, vecfunc); 482 double test = 0.0; 483 for (int i = 0; i < n; i++) 484 if (abs(fvec[i]) > test) 485 test = abs(fvec[i]); 486 if (test < 0.01 * TOLF) { 487 return false; 488 } 489 double sum = 0.0; 490 for (int i = 0; i < n; i++) { 491 double xi = x[i]; 492 sum += xi * xi; 493 } 494 double stpmax = STPMX * max(sqrt(sum), n); 495 for (int its = 0; its < MAXIT; its++) { 496 if (!vecfunc.jacobian(n, x, fjac)) 497 fdjac(x, fvec, fjac, n, vecfunc); 498 for (int i = 0; i < n; i++) { 499 sum = 0.0; 500 for (int j = 0; j < n; j++) 501 sum += fjac[j][i] * fvec[j]; 502 g[i] = sum; 503 } 504 System.arraycopy(x, 0, xold, 0, n); // for (int i=0;i < n;i++) xold[i]=x[i]; 505 double fold = f; 506 for (int i = 0; i < n; i++) 507 p[i] = -fvec[i]; 508 double d = ludcmp(fjac, indx, n); 509 lubksb(fjac, indx, p, n); 510 boolean check = lnsrch(xold, n, fold, g, p, x, f, stpmax, fvec, vecfunc); 511 test = 0.0; 512 for (int i = 0; i < n; i++) 513 if (abs(fvec[i]) > test) 514 test = abs(fvec[i]); 515 if (test < TOLF) { 516 return false; 517 } 518 if (check) { 519 test = 0.0; 520 double den = max(f, 0.5 * n); 521 for (int i = 0; i < n; i++) { 522 double temp = abs(g[i]) * max(abs(x[i]), 1.0) / den; 523 if (temp > test) 524 test = temp; 525 } 526 return (test < TOLMIN); 527 } 528 test = 0.0; 529 for (int i = 0; i < n; i++) { 530 double temp = (abs(x[i] - xold[i])) / max(abs(x[i]), 1.0); 531 if (temp > test) 532 test = temp; 533 } 534 if (test < EPS) { 535 return check; 536 } 537 } 538 539 throw new RootException("MAXIT exceeded in newtonND"); 540 } 541 542 private static boolean lnsrch(double[] xold, int n, double fold, double[] g, double[] p, 543 double[] x, double f, double stpmax, double[] fvec, VectorFunction vecfunc) throws RootException { 544 545 double sum = 0; 546 for (int i = 0; i < n; i++) { 547 double pi = p[i]; 548 sum += pi * pi; 549 } 550 sum = sqrt(sum); 551 if (sum > stpmax) 552 for (int i = 0; i < n; i++) 553 p[i] *= stpmax / sum; 554 double slope = 0.0; 555 for (int i = 0; i < n; i++) 556 slope += g[i] * p[i]; 557 if (slope >= 0.0) 558 throw new RootException("Roundoff problem in lnsrch."); 559 double test = 0; 560 for (int i = 0; i < n; i++) { 561 double temp = abs(p[i]) / max(abs(xold[i]), 1.0); 562 if (temp > test) 563 test = temp; 564 } 565 double alamin = EPS / test; 566 double alam = 1.0; 567 double tmplam, alam2 = 0, f2 = 0; 568 while (true) { 569 for (int i = 0; i < n; i++) 570 x[i] = xold[i] + alam * p[i]; 571 f = fmin(x, fvec, n, vecfunc); 572 if (alam < alamin) { 573 System.arraycopy(xold, 0, x, 0, n); // for (int i=0;i < n;i++) x[i]=xold[i]; 574 return true; 575 } else if (f <= fold + ALF * alam * slope) { 576 return false; 577 } else { 578 if (alam == 1.0) 579 tmplam = -slope / (2.0 * (f - fold - slope)); 580 else { 581 double rhs1 = f - fold - alam * slope; 582 double rhs2 = f2 - fold - alam2 * slope; 583 double a = (rhs1 / (alam * alam) - rhs2 / (alam2 * alam2)) / (alam - alam2); 584 double b = (-alam2 * rhs1 / (alam * alam) + alam * rhs2 / (alam2 * alam2)) / (alam - alam2); 585 if (a == 0.0) 586 tmplam = -slope / (2.0 * b); 587 else { 588 double disc = b * b - 3.0 * a * slope; 589 if (disc < 0.0) 590 tmplam = 0.5 * alam; 591 else if (b <= 0.0) 592 tmplam = (-b + sqrt(disc)) / (3.0 * a); 593 else 594 tmplam = -slope / (b + sqrt(disc)); 595 } 596 if (tmplam > 0.5 * alam) 597 tmplam = 0.5 * alam; 598 } 599 600 } 601 alam2 = alam; 602 f2 = f; 603 alam = max(tmplam, 0.1 * alam); 604 } 605 } 606 607 private static void lubksb(double[][] a, int[] indx, double[] b, int n) throws RootException { 608 int ii = 0; 609 for (int i = 0; i < n; i++) { 610 int ip = indx[i]; 611 double sum = b[ip]; 612 b[ip] = b[i]; 613 if (ii != 0) 614 for (int j = ii - 1; j < i; j++) 615 sum -= a[i][j] * b[j]; 616 else if (sum != 0.0) 617 ii = i + 1; 618 b[i] = sum; 619 } 620 for (int i = n - 1; i >= 0; i--) { 621 double sum = b[i]; 622 for (int j = i + 1; j < n; j++) 623 sum -= a[i][j] * b[j]; 624 b[i] = sum / a[i][i]; 625 } 626 } 627 628 private static double ludcmp(double[][] a, int[] indx, int n) throws RootException { 629 630 double[] vv = new double[n]; 631 double d = 1.0; 632 for (int i = 0; i < n; i++) { 633 double big = 0.0; 634 double temp; 635 for (int j = 0; j < n; j++) 636 if ((temp = abs(a[i][j])) > big) 637 big = temp; 638 if (big == 0.0) 639 throw new RootException("Singular matrix in routine ludcmp"); 640 vv[i] = 1.0 / big; 641 } 642 int imax = 0; 643 for (int j = 0; j < n; j++) { 644 for (int i = 0; i < j; i++) { 645 double sum = a[i][j]; 646 for (int k = 0; k < i; k++) 647 sum -= a[i][k] * a[k][j]; 648 a[i][j] = sum; 649 } 650 double big = 0.0; 651 for (int i = j; i < n; i++) { 652 double sum = a[i][j]; 653 for (int k = 0; k < j; k++) 654 sum -= a[i][k] * a[k][j]; 655 a[i][j] = sum; 656 double dum; 657 if ((dum = vv[i] * abs(sum)) >= big) { 658 big = dum; 659 imax = i; 660 } 661 } 662 if (j != imax) { 663 for (int k = 0; k < n; k++) { 664 double dum = a[imax][k]; 665 a[imax][k] = a[j][k]; 666 a[j][k] = dum; 667 } 668 d = -d; 669 vv[imax] = vv[j]; 670 } 671 indx[j] = imax; 672 if (a[j][j] == 0.0) 673 a[j][j] = TINY; 674 if (j != n - 1) { 675 double dum = 1.0 / (a[j][j]); 676 for (int i = j + 1; i < n; i++) 677 a[i][j] *= dum; 678 } 679 } 680 681 return d; 682 } 683 684 private static void fdjac(double[] x, double[] fvec, double[][] df, int n, VectorFunction vecfunc) throws RootException { 685 double[] f = new double[n]; 686 for (int j = 0; j < n; j++) { 687 double temp = x[j]; 688 double h = TOLF * abs(temp); 689 if (h == 0.0) 690 h = TOLF; 691 x[j] = temp + h; 692 h = x[j] - temp; 693 vecfunc.function(n, x, f); 694 x[j] = temp; 695 for (int i = 0; i < n; i++) 696 df[i][j] = (f[i] - fvec[i]) / h; 697 } 698 } 699 700 private static double fmin(double[] x, double[] fvec, int n, VectorFunction vecfunc) throws RootException { 701 vecfunc.function(n, x, fvec); 702 double sum = 0; 703 for (int i = 0; i < n; i++) { 704 double fveci = fvec[i]; 705 sum += fveci * fveci; 706 } 707 return 0.5 * sum; 708 } 709 710 /** 711 * Used to test out the methods in this class. 712 * 713 * @param args Command line arguments (not used). 714 */ 715 public static void main(String args[]) { 716 717 System.out.println(); 718 System.out.println("Testing Roots..."); 719 720 try { 721 // Find a root of f(x) = x^3 - 3*x + 2. 722 723 // First create an instance of our evaluatable function. 724 Evaluatable1D function = new DemoFunction1D(); 725 726 // Then call the root finder with brackets on the root and a tolerance. 727 // The brackets can be used to isolate which root you want (for instance). 728 double root = findRoot1D(function, -10, 0, 0.0001); 729// double root = aitken( function, 0, 0.0001 ); 730 731 System.out.println(" A root of f(x) = x^3 - 3*x + 2 is " + root); 732 733 // Find the roots to: x0^2 + x1^2 - 2; and e^(x0 - 1) + x1^3 - 2 734 VectorFunction funcND = new DemoFunctionND(); 735 double[] f = new double[2]; 736 double[] x = new double[2]; 737 x[0] = 2.0; 738 x[1] = 0.5; 739 boolean check = findRootsND(funcND, x, 2); 740 funcND.function(2, x, f); 741 if (check) 742 System.out.println("Convergence problems."); 743 System.out.println("Roots of x0^2 + x1^2 - 2; and e^(x0 - 1) + x1^3 - 2 are:"); 744 System.out.printf("%7s %3s %12s\n", "Index", "x", "f"); 745 for (int i = 0; i < 2; i++) 746 System.out.printf("%5d %12.6f %12.6f\n", i + 1, x[i], f[i]); 747 748 } catch (Exception e) { 749 e.printStackTrace(); 750 } 751 752 } 753 754 /** 755 * A simple demonstration class that shows how to use the Evaluatable1D class to find 756 * the roots of a 1D equation. Since this class returns a value for derivative(), we 757 * are signaling the root finder that we want it to use the Newton-Raphson technique 758 * to find the root of the equation. 759 */ 760 private static class DemoFunction1D implements Evaluatable1D { 761 762 /** 763 * User supplied method that calculates the function y = fn(x). This demonstration 764 * returns the value of y = x^3 - 3*x + 2. 765 * 766 * @param x Independent parameter to the function. 767 * @return The x^3 - 3*x + 2 evaluated at x. 768 */ 769 @Override 770 public double function(double x) { 771 return (x * x * x - 3. * x + 2.); 772 } 773 774 /** 775 * Calculates the derivative of the function dy/dx = d fn(x)/dx. This 776 * demonstration returns dy/dx = 3*x^2 - 3. 777 * 778 * @param x Independent parameter to the function. 779 * @return The 3*x^2 - 3 evaluated at x. 780 */ 781 @Override 782 public double derivative(double x) { 783 return (3. * x * x - 3.); 784 } 785 786 } 787 788 /** 789 * A demonstration function for use with the newtonND multi-dimensional root finder. 790 */ 791 private static class DemoFunctionND implements VectorFunction { 792 793 /** 794 * User supplied method that calculates the function y[] = fn(x[]). 795 * 796 * @param n The number of variables in the x & y arrays. 797 * @param x Independent parameters to the function, passed in as input. 798 * @param y An existing array that is filled in with the outputs of the function 799 */ 800 @Override 801 public void function(int n, double x[], double[] y) throws RootException { 802 double x0 = x[0]; 803 double x1 = x[1]; 804 double x12 = x1 * x1; 805 y[0] = x0 * x0 + x12 - 2.0; 806 y[1] = exp(x0 - 1.0) + x1 * x12 - 2.0; 807 } 808 809 /** 810 * User supplied method that calculates the Jacobian of the function. 811 * 812 * @param n The number of rows and columns in the Jacobian. 813 * @param jac The Jacobian array. 814 * @return True if the Jacobian was computed by this method, false if it was not. 815 */ 816 @Override 817 public boolean jacobian(int n, double[] x, double[][] jac) { 818 return false; 819 } 820 } 821}