1 /* 2 Copyright 2008,2009 3 Matthias Ehmann, 4 Michael Gerhaeuser, 5 Carsten Miller, 6 Bianca Valentin, 7 Alfred Wassermann, 8 Peter Wilfahrt 9 10 This file is part of JSXGraph. 11 12 JSXGraph is free software: you can redistribute it and/or modify 13 it under the terms of the GNU Lesser General Public License as published by 14 the Free Software Foundation, either version 3 of the License, or 15 (at your option) any later version. 16 17 JSXGraph is distributed in the hope that it will be useful, 18 but WITHOUT ANY WARRANTY; without even the implied warranty of 19 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 20 GNU Lesser General Public License for more details. 21 22 You should have received a copy of the GNU Lesser General Public License 23 along with JSXGraph. If not, see <http://www.gnu.org/licenses/>. 24 */ 25 26 /** 27 * @fileoverview In this file the namespace Math.Numerics is defined, which holds numerical 28 * algorithms for solving linear equations etc. 29 * @author graphjs 30 */ 31 32 33 /** 34 * The JXG.Math.Numerics namespace holds numerical algorithms, constants, and variables. 35 * @namespace 36 */ 37 JXG.Math.Numerics = (function(JXG, Math) { 38 39 // Predefined butcher tableaus for the common Runge-Kutta method (fourth order), Heun method (second order), and Euler method (first order). 40 var predefinedButcher = { 41 rk4: { 42 s: 4, 43 A: [ 44 [ 0, 0, 0, 0], 45 [0.5, 0, 0, 0], 46 [ 0, 0.5, 0, 0], 47 [ 0, 0, 1, 0] 48 ], 49 b: [1. / 6., 1. / 3., 1. / 3., 1. / 6.], 50 c: [0, 0.5, 0.5, 1] 51 }, 52 heun: { 53 s: 2, 54 A: [ 55 [0, 0], 56 [1, 0] 57 ], 58 b: [0.5, 0.5], 59 c: [0, 1] 60 }, 61 euler: { 62 s: 1, 63 A: [ 64 [0] 65 ], 66 b: [1], 67 c: [0] 68 } 69 }; 70 71 72 /** @lends JXG.Math.Numerics */ 73 return { 74 75 /** 76 * Solves a system of linear equations given by A and b using the Gauss-Jordan-elimination. 77 * The algorithm runs in-place. I.e. the entries of A and b are changed. 78 * @param {Array} A Square matrix represented by an array of rows, containing the coefficients of the lineare equation system. 79 * @param {Array} b A vector containing the linear equation system's right hand side. 80 * @throws {Error} If a non-square-matrix is given or if b has not the right length or A's rank is not full. 81 * @returns {Array} A vector that solves the linear equation system. 82 */ 83 Gauss: function(A, b) { 84 var eps = JXG.Math.eps, 85 // number of columns of A 86 n = A.length > 0 ? A[0].length : 0, 87 // copy the matrix to prevent changes in the original 88 Acopy, 89 // solution vector, to prevent changing b 90 x, 91 i, j, k, 92 // little helper to swap array elements 93 swap = function(i, j) { 94 var temp = this[i]; 95 96 this[i] = this[j]; 97 this[j] = temp; 98 }; 99 100 if ((n !== b.length) || (n !== A.length)) 101 throw new Error("JXG.Math.Numerics.Gauss: Dimensions don't match. A must be a square matrix and b must be of the same length as A."); 102 103 // initialize solution vector 104 Acopy = new Array(n); 105 x = b.slice(0, n); 106 for (i = 0; i < n; i++) { 107 Acopy[i] = A[i].slice(0, n); 108 } 109 110 // Gauss-Jordan-elimination 111 for (j = 0; j < n; j++) { 112 for (i = n - 1; i > j; i--) { 113 // Is the element which is to eliminate greater than zero? 114 if (Math.abs(Acopy[i][j]) > eps) { 115 // Equals pivot element zero? 116 if (Math.abs(Acopy[j][j]) < eps) { 117 // At least numerically, so we have to exchange the rows 118 swap.apply(Acopy, [i, j]); 119 swap.apply(x, [i, j]); 120 } else { 121 // Saves the L matrix of the LR-decomposition. unnecessary. 122 Acopy[i][j] /= Acopy[j][j]; 123 // Transform right-hand-side b 124 x[i] -= Acopy[i][j] * x[j]; 125 // subtract the multiple of A[i][j] / A[j][j] of the j-th row from the i-th. 126 for (k = j + 1; k < n; k ++) { 127 Acopy[i][k] -= Acopy[i][j] * Acopy[j][k]; 128 } 129 } 130 } 131 } 132 if (Math.abs(Acopy[j][j]) < eps) { // The absolute values of all coefficients below the j-th row in the j-th column are smaller than JXG.Math.eps. 133 throw new Error("JXG.Math.Numerics.Gauss(): The given matrix seems to be singular."); 134 } 135 } 136 this.backwardSolve(Acopy, x, true); // return Array 137 138 return x; 139 }, 140 141 /** 142 * Solves a system of linear equations given by the right triangular matrix R and vector b. 143 * @param {Array} R Right triangular matrix represented by an array of rows. All entries a_(i,j) with i < j are ignored. 144 * @param {Array} b Right hand side of the linear equation system. 145 * @param {Boolean} [canModify=false] If true, the right hand side vector is allowed to be changed by this method. 146 * @returns {Array} An array representing a vector that solves the system of linear equations. 147 */ 148 backwardSolve: function(R, b, canModify) { 149 var x, m, n, i, j; 150 151 if (canModify) { 152 x = b; 153 } else { 154 x = b.slice(0, b.length); 155 } 156 157 // m: number of rows of R 158 // n: number of columns of R 159 m = R.length; 160 n = R.length > 0 ? R[0].length : 0; 161 for (i = m - 1; i >= 0; i--) { 162 for (j = n - 1; j > i; j--) { 163 x[i] -= R[i][j] * x[j]; 164 } 165 x[i] /= R[i][i]; 166 } 167 168 return x; 169 }, 170 171 /** 172 * @private 173 * Gauss-Bareiss algorithm to compute the 174 * determinant of matrix without fractions. 175 * See H. Cohen, "A course in computational 176 * algebraic number thoery". 177 */ 178 gaussBareiss: function(mat) { 179 var k, c, s, i, j, p, n, M, t, 180 eps = JXG.Math.eps; 181 182 n = mat.length; 183 if (n<=0) { return 0; } 184 if (mat[0].length<n) { n = mat[0].length; } 185 186 // Copy the input matrix to M 187 M = new Array(n); 188 for (i = 0; i < n; i++) { 189 M[i] = mat[i].slice(0, n); 190 } 191 192 c = 1; 193 s = 1; 194 for (k=0;k<n-1;k++) { 195 p = M[k][k]; 196 if (Math.abs(p)<eps) { // Pivot step 197 for (i=0;i<n;i++) { 198 if (Math.abs(M[i][k])>=eps) break 199 } 200 if (i==n) { // No nonzero entry found in column k -> det(M) = 0 201 return 0.0; 202 } 203 for (j=k;j<n;j++) { // swap row i and k partially 204 t = M[i][j]; M[i][j] = M[k][j]; M[k][j] = t; 205 } 206 s = -s; 207 p = M[k][k]; 208 } 209 210 // Main step 211 for (i=k+1;i<n;i++) { 212 for (j=k+1;j<n;j++) { 213 t = p*M[i][j] - M[i][k]*M[k][j]; 214 M[i][j] = t/c; 215 } 216 } 217 c = p; 218 } 219 return s*M[n-1][n-1]; 220 221 }, 222 223 /** 224 * Computes the determinant of a square nxn matrix with the 225 * Gauss-Bareiss algorithm. 226 * @param {Array} mat Matrix. 227 * @returns {Number} The determinant pf the matrix mat. 228 The empty matrix returns 0. 229 */ 230 det: function(mat) { 231 return this.gaussBareiss(mat); 232 }, 233 234 /** 235 * Compute the Eigenvalues and Eigenvectors of a symmetric 3x3 matrix with the Jacobi method 236 * Adaption of a FORTRAN program by Ed Wilson, Dec. 25, 1990 237 * @param {Array} Ain A symmetric 3x3 matrix. 238 * @returns {Array} [A,V] the matrices A and V. The diagonal of A contains the Eigenvalues, V contains the Eigenvectors. 239 */ 240 Jacobi: function(Ain) { 241 var i,j,k,aa,si,co,tt,eps = JXG.Math.eps, 242 sum = 0.0, 243 ssum, amax, 244 n = Ain.length, 245 V = [ 246 [0,0,0], 247 [0,0,0], 248 [0,0,0] 249 ], 250 A = [ 251 [0,0,0], 252 [0,0,0], 253 [0,0,0] 254 ], nloops=0; 255 256 // Initialization. Set initial Eigenvectors. 257 for (i = 0; i < n; i++) { 258 for (j = 0; j < n; j++) { 259 V[i][j] = 0.0; 260 A[i][j] = Ain[i][j]; 261 sum += Math.abs(A[i][j]); 262 } 263 V[i][i] = 1.0; 264 } 265 // Trivial problems 266 if (n == 1) { 267 return [A,V]; 268 } 269 if (sum <= 0.0) { 270 return [A,V]; 271 } 272 273 sum /= (n * n); 274 275 // Reduce matrix to diagonal 276 do { 277 ssum = 0.0; 278 amax = 0.0; 279 for (j = 1; j < n; j++) { 280 for (i = 0; i < j; i++) { 281 // Check if A[i][j] is to be reduced 282 aa = Math.abs(A[i][j]); 283 if (aa > amax) { 284 amax = aa; 285 } 286 ssum += aa; 287 if (aa >= eps) { 288 // calculate rotation angle 289 aa = Math.atan2(2.0 * A[i][j], A[i][i] - A[j][j]) * 0.5; 290 si = Math.sin(aa); 291 co = Math.cos(aa); 292 // Modify 'i' and 'j' columns 293 for (k = 0; k < n; k++) { 294 tt = A[k][i]; 295 A[k][i] = co * tt + si * A[k][j]; 296 A[k][j] = -si * tt + co * A[k][j]; 297 tt = V[k][i]; 298 V[k][i] = co * tt + si * V[k][j]; 299 V[k][j] = -si * tt + co * V[k][j]; 300 } 301 // Modify diagonal terms 302 A[i][i] = co * A[i][i] + si * A[j][i]; 303 A[j][j] = -si * A[i][j] + co * A[j][j]; 304 A[i][j] = 0.0; 305 // Make 'A' matrix symmetrical 306 for (k = 0; k < n; k++) { 307 A[i][k] = A[k][i]; 308 A[j][k] = A[k][j]; 309 } 310 // A[i][j] made zero by rotation 311 } 312 } 313 } 314 nloops++; 315 } while (Math.abs(ssum) / sum > eps && nloops<2000); 316 return [A,V]; 317 }, 318 319 /** 320 * Calculates the integral of function f over interval using Newton-Cotes-algorithm. 321 * @param {Array} interval The integration interval, e.g. [0, 3]. 322 * @param {function} f A function which takes one argument of type number and returns a number. 323 * @param {object} [config={number_of_nodes:28,integration_type:'milne'}] The algorithm setup. Accepted properties are number_of_nodes of type number and integration_type 324 * with value being either 'trapez', 'simpson', or 'milne'. 325 * @returns {Number} Integral value of f over interval 326 * @throws {Error} If config.number_of_nodes doesn't match config.integration_type an exception is thrown. If you want to use 327 * simpson rule respectively milne rule config.number_of_nodes must be dividable by 2 respectively 4. 328 * @example 329 * function f(x) { 330 * return x*x; 331 * } 332 * 333 * // calculates integral of <tt>f</tt> from 0 to 2. 334 * var area1 = JXG.Math.Numerics.NewtonCotes([0, 2], f); 335 * 336 * // the same with an anonymous function 337 * var area2 = JXG.Math.Numerics.NewtonCotes([0, 2], function (x) { return x*x; }); 338 * 339 * // use trapez rule with 16 nodes 340 * var area3 = JXG.Math.Numerics.NewtonCotes([0, 2], f, 341 * {number_of_nodes: 16, intergration_type: 'trapez'}); 342 */ 343 NewtonCotes: function(interval, f, config) { 344 var integral_value = 0.0, 345 number_of_nodes = config && typeof config.number_of_nodes === 'number' ? config.number_of_nodes : 28, 346 available_types = {trapez: true, simpson: true, milne: true}, 347 integration_type = config && config.integration_type && available_types.hasOwnProperty(config.integration_type) && available_types[config.integration_type] ? config.integration_type : 'milne', 348 step_size = (interval[1] - interval[0]) / number_of_nodes, 349 evaluation_point, i, number_of_intervals; 350 351 switch (integration_type) { 352 case 'trapez': 353 integral_value = (f(interval[0]) + f(interval[1])) * 0.5; 354 355 evaluation_point = interval[0]; 356 for (i = 0; i < number_of_nodes - 1; i++) { 357 evaluation_point += step_size; 358 integral_value += f(evaluation_point); 359 } 360 integral_value *= step_size; 361 362 break; 363 case 'simpson': 364 if (number_of_nodes % 2 > 0) { 365 throw new Error("JSXGraph: INT_SIMPSON requires config.number_of_nodes dividable by 2."); 366 } 367 number_of_intervals = number_of_nodes / 2.0; 368 integral_value = f(interval[0]) + f(interval[1]); 369 evaluation_point = interval[0]; 370 for (i = 0; i < number_of_intervals - 1; i++) { 371 evaluation_point += 2.0 * step_size; 372 integral_value += 2.0 * f(evaluation_point); 373 } 374 evaluation_point = interval[0] - step_size; 375 for (i = 0; i < number_of_intervals; i++) { 376 evaluation_point += 2.0 * step_size; 377 integral_value += 4.0 * f(evaluation_point); 378 } 379 integral_value *= step_size / 3.0; 380 break; 381 default: 382 if (number_of_nodes % 4 > 0) { 383 throw new Error("JSXGraph: Error in INT_MILNE: config.number_of_nodes must be a multiple of 4"); 384 } 385 number_of_intervals = number_of_nodes * 0.25; 386 integral_value = 7.0 * (f(interval[0]) + f(interval[1])); 387 evaluation_point = interval[0]; 388 for (i = 0; i < number_of_intervals - 1; i++) { 389 evaluation_point += 4.0 * step_size; 390 integral_value += 14.0 * f(evaluation_point); 391 } 392 evaluation_point = interval[0] - 3.0 * step_size; 393 for (i = 0; i < number_of_intervals; i++) { 394 evaluation_point += 4.0 * step_size; 395 integral_value += 32.0 * (f(evaluation_point) + f(evaluation_point + 2 * step_size)); 396 } 397 evaluation_point = interval[0] - 2.0 * step_size; 398 for (i = 0; i < number_of_intervals; i++) { 399 evaluation_point += 4.0 * step_size; 400 integral_value += 12.0 * f(evaluation_point); 401 } 402 integral_value *= 2.0 * step_size / 45.0; 403 } 404 return integral_value; 405 }, 406 407 /** 408 * Integral of function f over interval. 409 * @param {Array} interval The integration interval, e.g. [0, 3]. 410 * @param {function} f A function which takes one argument of type number and returns a number. 411 * @returns {Number} The value of the integral of f over interval 412 * @see JXG.Math.Numerics.NewtonCotes 413 */ 414 I: function(interval, f) { 415 return this.NewtonCotes(interval, f, {number_of_nodes: 16, integration_type: 'milne'}); 416 }, 417 418 /** 419 * Newton's method to find roots of a funtion in one variable. 420 * @param {function} f We search for a solution of f(x)=0. 421 * @param {Number} x initial guess for the root, i.e. start value. 422 * @param {object} object optional object that is treated as "this" in the function body. This is useful, if the function is a 423 * method of an object and contains a reference to its parent object via "this". 424 * @returns {Number} A root of the function f. 425 */ 426 Newton: function(f, x, object) { 427 var i = 0, 428 h = JXG.Math.eps, 429 newf = f.apply(object, [x]), // set "this" to "object" in f 430 nfev = 1, 431 df; 432 if (JXG.isArray(x)) { // For compatibility 433 x = x[0]; 434 } 435 while (i < 50 && Math.abs(newf) > h) { 436 df = this.D(f, object)(x); nfev += 2; 437 if (Math.abs(df) > h) { 438 x -= newf / df; 439 } else { 440 x += (Math.random() * 0.2 - 1.0); 441 } 442 newf = f.apply(object, [x]); nfev++; 443 i++; 444 } 445 //JXG.debug("N nfev="+nfev); 446 return x; 447 }, 448 449 /** 450 * Abstract method to find roots of univariate functions. 451 * @param {function} f We search for a solution of f(x)=0. 452 * @param {Number} x initial guess for the root, i.e. staring value. 453 * @param {object} object optional object that is treated as "this" in the function body. This is useful, if the function is a 454 * method of an object and contains a reference to its parent object via "this". 455 * @returns {Number} A root of the function f. 456 */ 457 root: function(f, x, object) { 458 //return this.Newton(f, x, object); 459 return this.fzero(f, x, object); 460 }, 461 462 /** 463 * Returns the Lagrange polynomials for curves with equidistant nodes, see 464 * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation, 465 * SIAM Review, Vol 46, No 3, (2004) 501-517. 466 * The graph of the parametric curve [x(t),y(t)] runs through the given points. 467 * @param {Array} p Array of JXG.Points 468 * @returns {Array} An array consisting of two functions x(t), y(t) which define a parametric curve 469 * f(t) = (x(t), y(t)) and two numbers x1 and x2 defining the curve's domain. x1 always equals zero. 470 */ 471 Neville: function(p) { 472 var w = [], 473 makeFct = function(fun) { 474 return function(t, suspendedUpdate) { 475 var i, d, s, 476 bin = JXG.Math.binomial, 477 len = p.length, 478 len1 = len - 1, 479 num = 0.0, 480 denom = 0.0; 481 482 if (!suspendedUpdate) { 483 s = 1; 484 for (i = 0; i < len; i++) { 485 w[i] = bin(len1, i) * s; 486 s *= (-1); 487 } 488 } 489 490 d = t; 491 for (i = 0; i < len; i++) { 492 if (d === 0) { 493 return p[i][fun](); 494 } else { 495 s = w[i] / d; 496 d--; 497 num += p[i][fun]() * s; 498 denom += s; 499 } 500 } 501 return num / denom; 502 }; 503 }, 504 505 xfct = makeFct('X'), 506 yfct = makeFct('Y'); 507 508 return [xfct, yfct, 0, function() { 509 return p.length - 1; 510 }]; 511 }, 512 513 /** 514 * Calculates second derivatives at the given knots. 515 * @param {Array} x x values of knots 516 * @param {Array} y y values of knots 517 * @returns {Array} Second derivatives of the interpolated function at the knots. 518 * @see #splineEval 519 */ 520 splineDef: function(x, y) { 521 var n = Math.min(x.length, y.length), 522 pair, i, l, 523 diag = [], 524 z = [], 525 data = [], 526 dx = [], 527 delta = [], 528 F = []; 529 530 if (n === 2) { 531 return [0, 0]; 532 } 533 534 for (i = 0; i < n; i++) { 535 pair = {X: x[i], Y: y[i]}; 536 data.push(pair); 537 } 538 data.sort(function (a, b) { 539 return a.X - b.X; 540 }); 541 for (i = 0; i < n; i++) { 542 x[i] = data[i].X; 543 y[i] = data[i].Y; 544 } 545 546 for (i = 0; i < n - 1; i++) { 547 dx.push(x[i + 1] - x[i]); 548 } 549 for (i = 0; i < n - 2; i++) { 550 delta.push(6 * (y[i + 2] - y[i + 1]) / (dx[i + 1]) - 6 * (y[i + 1] - y[i]) / (dx[i])); 551 } 552 553 // ForwardSolve 554 diag.push(2 * (dx[0] + dx[1])); 555 z.push(delta[0]); 556 557 for (i = 0; i < n - 3; i++) { 558 l = dx[i + 1] / diag[i]; 559 diag.push(2 * (dx[i + 1] + dx[i + 2]) - l * dx[i + 1]); 560 z.push(delta[i + 1] - l * z[i]); 561 } 562 563 // BackwardSolve 564 F[n - 3] = z[n - 3] / diag[n - 3]; 565 for (i = n - 4; i >= 0; i--) { 566 F[i] = (z[i] - (dx[i + 1] * F[i + 1])) / diag[i]; 567 } 568 569 // Generate f''-Vector 570 for (i = n - 3; i >= 0; i--) { 571 F[i + 1] = F[i]; 572 } 573 574 // natural cubic spline 575 F[0] = 0; 576 F[n - 1] = 0; 577 return F; 578 }, 579 580 /** 581 * Evaluate points on spline. 582 * @param {Number,Array} x0 A single float value or an array of values to evaluate 583 * @param {Array} x x values of knots 584 * @param {Array} y y values of knots 585 * @param {Array} F Second derivatives at knots, calculated by {@link #splineDef} 586 * @see #splineDef 587 * @returns {Number,Array} A single value or an array, depending on what is given as x0. 588 */ 589 splineEval: function(x0, x, y, F) { 590 var n = Math.min(x.length, y.length), 591 l = 1, 592 asArray = false, 593 y0 = [], 594 i, j, a, b, c, d, x_; 595 596 // number of points to be evaluated 597 if (JXG.isArray(x0)) { 598 l = x0.length; 599 asArray = true; 600 } else 601 x0 = [x0]; 602 603 for (i = 0; i < l; i++) { 604 // is x0 in defining interval? 605 if ((x0[i] < x[0]) || (x[i] > x[n - 1])) 606 return NaN; 607 608 // determine part of spline in which x0 lies 609 for (j = 1; j < n; j++) { 610 if (x0[i] <= x[j]) 611 break; 612 } 613 j--; 614 615 // we're now in the j-th partial interval, i.e. x[j] < x0[i] <= x[j+1]; 616 // determine the coefficients of the polynomial in this interval 617 a = y[j]; 618 b = (y[j + 1] - y[j]) / (x[j + 1] - x[j]) - (x[j + 1] - x[j]) / 6 * (F[j + 1] + 2 * F[j]); 619 c = F[j] / 2; 620 d = (F[j + 1] - F[j]) / (6 * (x[j + 1] - x[j])); 621 // evaluate x0[i] 622 x_ = x0[i] - x[j]; 623 //y0.push(a + b*x_ + c*x_*x_ + d*x_*x_*x_); 624 y0.push(a + (b + (c + d * x_) * x_) * x_); 625 } 626 627 if (asArray) 628 return y0; 629 else 630 return y0[0]; 631 }, 632 633 /** 634 * Generate a string containing the function term of a polynomial. 635 * @param {Array} coeffs Coefficients of the polynomial. The position i belongs to x^i. 636 * @param {Number} deg Degree of the polynomial 637 * @param {String} varname Name of the variable (usually 'x') 638 * @param {Number} prec Precision 639 * @returns {String} A string containg the function term of the polynomial. 640 */ 641 generatePolynomialTerm: function(coeffs, deg, varname, prec) { 642 var t = [], i; 643 for (i = deg; i >= 0; i--) { 644 t = t.concat(['(', coeffs[i].toPrecision(prec), ')']); 645 if (i > 1) { 646 t = t.concat(['*', varname, '<sup>', i, '<', '/sup> + ']); 647 } else if (i === 1) { 648 t = t.concat(['*', varname, ' + ']); 649 } 650 } 651 return t.join(''); 652 }, 653 654 /** 655 * Computes the polynomial through a given set of coordinates in Lagrange form. 656 * Returns the Lagrange polynomials, see 657 * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation, 658 * SIAM Review, Vol 46, No 3, (2004) 501-517. 659 * @param {Array} p Array of JXG.Points 660 * @returns {function} A function of one parameter which returns the value of the polynomial, whose graph runs through the given points. 661 */ 662 lagrangePolynomial: function(p) { 663 var w = [], 664 fct = function(x, suspendedUpdate) { 665 var i, k, len, xi, s, 666 num = 0, denom = 0, 667 M, j; 668 669 len = p.length; 670 if (!suspendedUpdate) { 671 for (i = 0; i < len; i++) { 672 w[i] = 1.0; 673 xi = p[i].X(); 674 for (k = 0; k < len; k++) if (k != i) { 675 w[i] *= (xi - p[k].X()); 676 } 677 w[i] = 1 / w[i]; 678 } 679 680 M = []; 681 for (j = 0; j < len; j++) { 682 M.push([1]); 683 } 684 } 685 686 for (i = 0; i < len; i++) { 687 xi = p[i].X(); 688 if (x === xi) { 689 return p[i].Y(); 690 } else { 691 s = w[i] / (x - xi); 692 denom += s; 693 num += s * p[i].Y(); 694 } 695 } 696 return num / denom; 697 }; 698 699 fct.getTerm = function() { 700 return ''; 701 }; 702 703 return fct; 704 }, 705 706 /** 707 * Computes the regression polynomial of a given degree through a given set of coordinates. 708 * Returns the regression polynomial function. 709 * @param {Number} degree number, function or slider. 710 * Either 711 * @param {Array} dataX array containing the x-coordinates of the data set 712 * @param {Array} dataY array containing the y-coordinates of the data set, 713 * or 714 * @param {Array} dataX Array consisting of JXG.Points. 715 * @param {} dataY Ignored 716 * @returns {function} A function of one parameter which returns the value of the regression polynomial of the given degree. 717 * It possesses the method getTerm() which returns the string containing the function term of the polynomial. 718 */ 719 regressionPolynomial: function(degree, dataX, dataY) { 720 var coeffs, 721 deg, dX, dY, 722 inputType, 723 fct, 724 term = ''; 725 726 if (JXG.isPoint(degree) && typeof degree.Value == 'function') { // Slider 727 deg = function() {return degree.Value();}; 728 } else if (JXG.isFunction(degree)) { 729 deg = degree; 730 } else if (JXG.isNumber(degree)) { 731 deg = function() {return degree;}; 732 } else { 733 throw new Error("JSXGraph: Can't create regressionPolynomial from degree of type'" + (typeof degree) + "'."); 734 } 735 736 if (arguments.length == 3 && JXG.isArray(dataX) && JXG.isArray(dataY)) { // Parameters degree, dataX, dataY 737 inputType = 0; 738 } else if (arguments.length == 2 && JXG.isArray(dataX) && dataX.length > 0 && JXG.isPoint(dataX[0])) { // Parameters degree, point array 739 inputType = 1; 740 } else { 741 throw new Error("JSXGraph: Can't create regressionPolynomial. Wrong parameters."); 742 } 743 744 fct = function(x, suspendedUpdate) { 745 var i, j, M, MT, y, B, c, s, 746 d, len = dataX.length; // input data 747 748 d = Math.floor(deg()); // input data 749 if (!suspendedUpdate) { 750 if (inputType === 1) { // point list as input 751 dX = []; 752 dY = []; 753 for (i = 0; i < len; i++) { 754 dX[i] = dataX[i].X(); 755 dY[i] = dataX[i].Y(); 756 } 757 } 758 759 if (inputType === 0) { // check for functions 760 dX = []; 761 dY = []; 762 for (i = 0; i < len; i++) { 763 if (JXG.isFunction(dataX[i])) 764 dX.push(dataX[i]()); 765 else 766 dX.push(dataX[i]); 767 if (JXG.isFunction(dataY[i])) 768 dY.push(dataY[i]()); 769 else 770 dY.push(dataY[i]); 771 } 772 } 773 774 M = []; 775 for (j = 0; j < len; j++) { 776 M.push([1]); 777 } 778 for (i = 1; i <= d; i++) { 779 for (j = 0; j < len; j++) { 780 M[j][i] = M[j][i - 1] * dX[j]; // input data 781 } 782 } 783 784 y = dY; // input data 785 MT = JXG.Math.transpose(M); 786 B = JXG.Math.matMatMult(MT, M); 787 c = JXG.Math.matVecMult(MT, y); 788 coeffs = JXG.Math.Numerics.Gauss(B, c); 789 term = JXG.Math.Numerics.generatePolynomialTerm(coeffs, d, 'x', 3); 790 } 791 792 // Horner's scheme to evaluate polynomial 793 s = coeffs[d]; 794 for (i = d - 1; i >= 0; i--) { 795 s = (s * x + coeffs[i]); 796 } 797 return s; 798 }; 799 800 fct.getTerm = function() { 801 return term; 802 }; 803 804 return fct; 805 }, 806 807 /** 808 * Computes the cubic Bezier curve through a given set of points. 809 * @param {Array} points Array consisting of 3*k+1 JXG.Points. 810 * The points at position k with k mod 3 = 0 are the data points, 811 * points at position k with k mod 3 = 1 or 2 are the control points. 812 * @returns {Array} An array consisting of two functions of one parameter t which return the 813 * x resp. y coordinates of the Bezier curve in t, one zero value, and a third function accepting 814 * no parameters and returning one third of the length of the points. 815 */ 816 bezier: function(points) { 817 var len, 818 makeFct = function(which) { 819 return function(t, suspendedUpdate) { 820 var z = Math.floor(t) * 3, 821 t0 = t % 1, 822 t1 = 1 - t0; 823 824 if (!suspendedUpdate) { 825 len = Math.floor(points.length / 3); 826 } 827 828 if (t < 0) { 829 return points[0][which](); 830 } 831 if (t >= len) { 832 return points[points.length - 1][which](); 833 } 834 if (isNaN(t)) { 835 return NaN; 836 } 837 return t1 * t1 * (t1 * points[z][which]() + 3 * t0 * points[z + 1][which]()) + (3 * t1 * points[z + 2][which]() + t0 * points[z + 3][which]()) * t0 * t0; 838 }; 839 }; 840 841 return [ 842 makeFct('X'), 843 makeFct('Y'), 844 0, 845 function() { 846 return Math.floor(points.length / 3); 847 }]; 848 }, 849 850 /** 851 * Computes the B-spline curve of order k (order = degree+1) through a given set of points. 852 * @param {Array} points Array consisting of JXG.Points. 853 * @param {Number} order Order of the B-spline curve. 854 * @todo closed B-spline curves 855 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 856 * which return the x resp. y coordinates of the B-spline curve in t, a zero value, and a function simply 857 * returning the length of the points array minus one. 858 */ 859 bspline: function(points, order) { 860 var knots, N = [], 861 _knotVector = function(n, k) { 862 var j, kn = []; 863 for (j = 0; j < n + k + 1; j++) { 864 if (j < k) { 865 kn[j] = 0.0; 866 } else if (j <= n) { 867 kn[j] = j - k + 1; 868 } else { 869 kn[j] = n - k + 2; 870 } 871 } 872 return kn; 873 }, 874 875 _evalBasisFuncs = function(t, kn, n, k, s) { 876 var i,j,a,b,den, 877 N = []; 878 879 if (kn[s] <= t && t < kn[s + 1]) { 880 N[s] = 1; 881 } else { 882 N[s] = 0; 883 } 884 for (i = 2; i <= k; i++) { 885 for (j = s - i + 1; j <= s; j++) { 886 if (j <= s - i + 1 || j < 0) { 887 a = 0.0; 888 } else { 889 a = N[j]; 890 } 891 if (j >= s) { 892 b = 0.0; 893 } else { 894 b = N[j + 1]; 895 } 896 den = kn[j + i - 1] - kn[j]; 897 if (den == 0) { 898 N[j] = 0; 899 } else { 900 N[j] = (t - kn[j]) / den * a; 901 } 902 den = kn[j + i] - kn[j + 1]; 903 if (den != 0) { 904 N[j] += (kn[j + i] - t) / den * b; 905 } 906 } 907 } 908 return N; 909 }, 910 makeFct = function(which) { 911 return function(t, suspendedUpdate) { 912 var len = points.length, 913 y, j, s, 914 n = len - 1, 915 k = order; 916 917 if (n <= 0) return NaN; 918 if (n + 2 <= k) k = n + 1; 919 if (t <= 0) return points[0][which](); 920 if (t >= n - k + 2) return points[n][which](); 921 922 knots = _knotVector(n, k); 923 s = Math.floor(t) + k - 1; 924 N = _evalBasisFuncs(t, knots, n, k, s); 925 926 y = 0.0; 927 for (j = s - k + 1; j <= s; j++) { 928 if (j < len && j >= 0) y += points[j][which]() * N[j]; 929 } 930 return y; 931 }; 932 }; 933 934 return [makeFct('X'), 935 makeFct('Y'), 936 0, 937 function() { 938 return points.length - 1; 939 } 940 ]; 941 }, 942 943 /** 944 * Numerical (symmetric) approximation of derivative. suspendUpdate is piped through, see {@link JXG.Curve#updateCurve} 945 * and {@link JXG.Curve#hasPoint}. 946 * @param {function} f Function in one variable to be differentiated. 947 * @param {object} [obj] Optional object that is treated as "this" in the function body. This is useful, if the function is a 948 * method of an object and contains a reference to its parent object via "this". 949 * @returns {function} Derivative function of a given function f. 950 */ 951 D: function(f, obj) { 952 var h = 0.00001, 953 h2 = 1.0 / (h * 2.0); 954 955 if (arguments.length == 1 || (arguments.length > 1 && !JXG.exists(arguments[1]))) { 956 return function(x, suspendUpdate) { 957 return (f(x + h, suspendUpdate) - f(x - h, suspendUpdate)) * h2; 958 }; 959 } else { // set "this" to "obj" in f 960 return function(x, suspendUpdate) { 961 return (f.apply(obj, [x + h,suspendUpdate]) - f.apply(obj, [x - h,suspendUpdate])) * h2; 962 }; 963 } 964 }, 965 966 /** 967 * Helper function to create curve which displays Riemann sums. 968 * Compute coordinates for the rectangles showing the Riemann sum. 969 * @param {function} f Function, whose integral is approximated by the Riemann sum. 970 * @param {Number} n number of rectangles. 971 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', or 'trapezodial'. 972 * @param {Number} start Left border of the approximation interval 973 * @param {Number} end Right border of the approximation interval 974 * @returns {Array} An array of two arrays containing the x and y coordinates for the rectangles showing the Riemann sum. This 975 * array may be used as parent array of a {@link JXG.Curve}. 976 */ 977 riemann: function(f, n, type, start, end) { 978 var xarr = [], 979 yarr = [], 980 i, j = 0, 981 delta, 982 x = start, y, 983 x1, y1, delta1; 984 985 n = Math.floor(n); 986 xarr[j] = x; 987 yarr[j] = 0.0; 988 989 if (n > 0) { 990 delta = (end - start) / n; 991 delta1 = delta * 0.01; // for 'lower' and 'upper' 992 993 for (i = 0; i < n; i++) { 994 if (type === 'right') { 995 y = f(x + delta); 996 } else if (type === 'middle') { 997 y = f(x + delta * 0.5); 998 } else if ((type === 'left') || (type === 'trapezodial')) { 999 y = f(x); 1000 } else if (type === 'lower') { 1001 y = f(x); 1002 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1003 y1 = f(x1); 1004 if (y1 < y) { 1005 y = y1; 1006 } 1007 } 1008 } else { // (type=='upper') 1009 y = f(x); 1010 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1011 y1 = f(x1); 1012 if (y1 > y) { 1013 y = y1; 1014 } 1015 } 1016 } 1017 1018 j++; 1019 xarr[j] = x; 1020 yarr[j] = y; 1021 j++; 1022 x += delta; 1023 if (type === 'trapezodial') { 1024 y = f(x); 1025 } 1026 xarr[j] = x; 1027 yarr[j] = y; 1028 j++; 1029 xarr[j] = x; 1030 yarr[j] = 0.0; 1031 } 1032 } 1033 return [xarr,yarr]; 1034 }, 1035 1036 /** 1037 * Approximate the integral by Riemann sums. 1038 * Compute the area described by the riemann sum rectangles. 1039 * @param {function} f Function, whose integral is approximated by the Riemann sum. 1040 * @param {Number} n number of rectangles. 1041 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', or 'trapezodial'. 1042 * @param {Number} start Left border of the approximation interval 1043 * @param {Number} end Right border of the approximation interval 1044 * @returns {Number} The sum of the areas of the rectangles. 1045 */ 1046 riemannsum: function(f, n, type, start, end) { 1047 var sum = .0, 1048 i, delta, 1049 x = start, y, 1050 x1, y1, delta1; 1051 1052 n = Math.floor(n); 1053 if (n > 0) { 1054 delta = (end - start) / n; 1055 delta1 = delta * 0.01; // for 'lower' and 'upper' 1056 for (i = 0; i < n; i++) { 1057 if (type === 'right') { 1058 y = f(x + delta); 1059 } else if (type === 'middle') { 1060 y = f(x + delta * 0.5); 1061 } else if (type === 'trapezodial') { 1062 y = 0.5 * (f(x + delta) + f(x)); 1063 } else if (type === 'left') { 1064 y = f(x); 1065 } else if (type === 'lower') { 1066 y = f(x); 1067 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1068 y1 = f(x1); 1069 if (y1 < y) { 1070 y = y1; 1071 } 1072 } 1073 } else { // (type=='upper') 1074 y = f(x); 1075 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1076 y1 = f(x1); 1077 if (y1 > y) { 1078 y = y1; 1079 } 1080 } 1081 } 1082 sum += delta * y; 1083 x += delta; 1084 } 1085 } 1086 return sum; 1087 }, 1088 1089 /** 1090 * Solve initial value problems numerically using Runge-Kutta-methods. 1091 * See {@link http://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm. 1092 * @param {object,String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing 1093 * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure 1094 * <pre> 1095 * { 1096 * s: <Number>, 1097 * A: <matrix>, 1098 * b: <Array>, 1099 * c: <Array> 1100 * } 1101 * </pre> 1102 * which corresponds to the Butcher tableau structure shown here: http://en.wikipedia.org/w/index.php?title=List_of_Runge%E2%80%93Kutta_methods&oldid=357796696 1103 * @param {Array} x0 Initial value vector. If the problem is of one-dimensional, the initial value also has to be given in an array. 1104 * @param {Array} I Interval on which to integrate. 1105 * @param {Number} N Number of evaluation points. 1106 * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode 1107 * is given by the equation <pre>dx/dt = f(t, x(t)).</pre> So f has to take two parameters, a number <tt>t</tt> and a 1108 * vector <tt>x</tt>, and has to return a vector of the same dimension as <tt>x</tt> has. 1109 * @returns {Array} An array of vectors describing the solution of the ode on the given interval I. 1110 * @example 1111 * // A very simple autonomous system dx(t)/dt = x(t); 1112 * function f(t, x) { 1113 * return x; 1114 * } 1115 * 1116 * // Solve it with initial value x(0) = 1 on the interval [0, 2] 1117 * // with 20 evaluation points. 1118 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 1119 * 1120 * // Prepare data for plotting the solution of the ode using a curve. 1121 * var dataX = []; 1122 * var dataY = []; 1123 * var h = 0.1; // (I[1] - I[0])/N = (2-0)/20 1124 * for(var i=0; i<data.length; i++) { 1125 * dataX[i] = i*h; 1126 * dataY[i] = data[i][0]; 1127 * } 1128 * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'}); 1129 * </pre><div id="d2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div> 1130 * <script type="text/javascript"> 1131 * var board = JXG.JSXGraph.initBoard('d2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false}); 1132 * function f(t, x) { 1133 * // we have to copy the value. 1134 * // return x; would just return the reference. 1135 * return [x[0]]; 1136 * } 1137 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 1138 * var dataX = []; 1139 * var dataY = []; 1140 * var h = 0.1; 1141 * for(var i=0; i<data.length; i++) { 1142 * dataX[i] = i*h; 1143 * dataY[i] = data[i][0]; 1144 * } 1145 * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'}); 1146 * </script><pre> 1147 */ 1148 rungeKutta: function(butcher, x0, I, N, f) { 1149 var x = [], 1150 y = [], 1151 h = (I[1] - I[0]) / N, 1152 t = I[0], 1153 e, i, j, 1154 k, l, 1155 dim = x0.length, 1156 s, 1157 result = [], 1158 r = 0; 1159 1160 if (JXG.isString(butcher)) { 1161 butcher = predefinedButcher[butcher] || predefinedButcher.euler; 1162 } 1163 s = butcher.s; 1164 1165 // don't change x0, so copy it 1166 for (e = 0; e < dim; e++) 1167 x[e] = x0[e]; 1168 1169 for (i = 0; i < N; i++) { 1170 // Optimization doesn't work for ODEs plotted using time 1171 // if((i % quotient == 0) || (i == N-1)) { 1172 result[r] = []; 1173 for (e = 0; e < dim; e++) 1174 result[r][e] = x[e]; 1175 r++; 1176 // } 1177 // init k 1178 k = []; 1179 for (j = 0; j < s; j++) { 1180 // init y = 0 1181 for (e = 0; e < dim; e++) 1182 y[e] = 0.; 1183 1184 // Calculate linear combination of former k's and save it in y 1185 for (l = 0; l < j; l++) { 1186 for (e = 0; e < dim; e++) { 1187 y[e] += (butcher.A[j][l]) * h * k[l][e]; 1188 } 1189 } 1190 1191 // add x(t) to y 1192 for (e = 0; e < dim; e++) { 1193 y[e] += x[e]; 1194 } 1195 1196 // calculate new k and add it to the k matrix 1197 k.push(f(t + butcher.c[j] * h, y)); 1198 } 1199 1200 // init y = 0 1201 for (e = 0; e < dim; e++) 1202 y[e] = 0.; 1203 1204 for (l = 0; l < s; l++) { 1205 for (e = 0; e < dim; e++) 1206 y[e] += butcher.b[l] * k[l][e]; 1207 } 1208 1209 for (e = 0; e < dim; e++) { 1210 x[e] = x[e] + h * y[e]; 1211 } 1212 1213 t += h; 1214 } 1215 1216 return result; 1217 }, 1218 1219 1220 /** 1221 * 1222 * Find zero of an univariate function f. 1223 * @param {function} f Function, whose root is to be found 1224 * @param {array or number} x0 Start value or start interval enclosing the root 1225 * @param {object} object Parent object in case f is method of it 1226 * @return {number} the approximation of the root 1227 * Algorithm: 1228 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 1229 * computations. M., Mir, 1980, p.180 of the Russian edition 1230 * 1231 * if x0 is an array containing lower and upper bound for the zero 1232 * algorithm 748 is applied. Otherwise, if x0 is a number, 1233 * the algorithm tries to bracket a zero of f starting from x0. 1234 * If this fails, we fall back to Newton's method. 1235 * 1236 **/ 1237 fzero: function(f, x0, object) { 1238 var tol = JXG.Math.eps, 1239 maxiter = 50, niter = 0, 1240 nfev = 0, 1241 eps = tol, 1242 a,b,c, 1243 fa, fb, fc, 1244 aa, blist, i, len, u, fu, 1245 prev_step, 1246 tol_act, // Actual tolerance 1247 p, // Interpolation step is calcu- 1248 q, // lated in the form p/q; divi- 1249 // sion operations is delayed 1250 // until the last moment 1251 new_step, // Step at this iteration 1252 t1, cb, t2; 1253 1254 if (JXG.isArray(x0)) { 1255 if (x0.length<2) 1256 throw new Error("JXG.Math.Numerics.fzero: length of array x0 has to be at least two."); 1257 a = x0[0]; fa = f.apply(object,[a]); nfev++; 1258 b = x0[1]; fb = f.apply(object,[b]); nfev++; 1259 } else { 1260 a = x0; fa = f.apply(object,[a]); nfev++; 1261 // Try to get b. 1262 if (a == 0) { 1263 aa = 1; 1264 } else { 1265 aa = a; 1266 } 1267 blist = [0.9*aa, 1.1*aa, aa-1, aa+1, 0.5*aa, 1.5*aa, -aa, 2*aa, -10*aa, 10*aa]; 1268 len = blist.length; 1269 for (i=0;i<len;i++) { 1270 b = blist[i]; 1271 fb = f.apply(object,[b]); nfev++; 1272 if (fa*fb<=0) { 1273 break; 1274 } 1275 } 1276 if (b < a) { 1277 u = a; a = b; b = u; 1278 fu = fa; fa = fb; fb = fu; 1279 } 1280 } 1281 1282 if (fa*fb > 0) { 1283 // Bracketing not successful, fall back to Newton's method or to fminbr 1284 if (JXG.isArray(x0)) { 1285 //JXG.debug("fzero falls back to fminbr"); 1286 return this.fminbr(f, [a,b], object); 1287 } else { 1288 //JXG.debug("fzero falls back to Newton"); 1289 return this.Newton(f, a, object); 1290 } 1291 } 1292 1293 // OK, we have enclosed a zero of f. 1294 // Now we can start Brent's method 1295 1296 c = a; 1297 fc = fa; 1298 while (niter<maxiter) { // Main iteration loop 1299 prev_step = b-a; // Distance from the last but one 1300 // to the last approximation 1301 1302 if ( Math.abs(fc) < Math.abs(fb) ) { 1303 // Swap data for b to be the 1304 a = b; b = c; c = a; // best approximation 1305 fa=fb; fb=fc; fc=fa; 1306 } 1307 tol_act = 2*eps*Math.abs(b) + tol*0.5; 1308 new_step = (c-b)*0.5; 1309 1310 if ( Math.abs(new_step) <= tol_act || Math.abs(fb) <= eps ) { 1311 //JXG.debug("nfev="+nfev); 1312 return b; // Acceptable approx. is found 1313 } 1314 1315 // Decide if the interpolation can be tried 1316 if ( Math.abs(prev_step) >= tol_act // If prev_step was large enough 1317 && Math.abs(fa) > Math.abs(fb) ) { // and was in true direction, 1318 // Interpolatiom may be tried 1319 cb = c-b; 1320 if ( a==c ) { // If we have only two distinct 1321 // points linear interpolation 1322 t1 = fb/fa; // can only be applied 1323 p = cb*t1; 1324 q = 1.0 - t1; 1325 } else { // Quadric inverse interpolation 1326 q = fa/fc; t1 = fb/fc; t2 = fb/fa; 1327 p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) ); 1328 q = (q-1.0) * (t1-1.0) * (t2-1.0); 1329 } 1330 if ( p>0 ) { // p was calculated with the op- 1331 q = -q; // posite sign; make p positive 1332 } else { // and assign possible minus to 1333 p = -p; // q 1334 } 1335 1336 if( p < (0.75*cb*q-Math.abs(tol_act*q)*0.5) // If b+p/q falls in [b,c] 1337 && p < Math.abs(prev_step*q*0.5) ) { // and isn't too large 1338 new_step = p/q; // it is accepted 1339 } // If p/q is too large then the 1340 // bissection procedure can 1341 // reduce [b,c] range to more 1342 // extent 1343 } 1344 1345 if ( Math.abs(new_step) < tol_act ) { // Adjust the step to be not less 1346 if ( new_step > 0 ) { // than tolerance 1347 new_step = tol_act; 1348 } else { 1349 new_step = -tol_act; 1350 } 1351 } 1352 1353 a = b; fa = fb; // Save the previous approx. 1354 b += new_step; 1355 fb = f.apply(object,[b]); nfev++; // Do step to a new approxim. 1356 if ( (fb>0 && fc>0) || (fb<0 && fc<0) ) { 1357 // Adjust c for it to have a sign 1358 c = a; fc = fa; // opposite to that of b 1359 } 1360 niter++; 1361 } // End while 1362 1363 JXG.debug("fzero: maxiter="+maxiter+" reached."); 1364 return b; 1365 }, 1366 1367 /** 1368 * 1369 * Find minimum of an univariate function f. 1370 * @param {function} f Function, whose minimum is to be found 1371 * @param {array} x0 Start interval enclosing the minimum 1372 * @param {object} object Parent object in case f is method of it 1373 * @return {number} the approximation of the minimum 1374 * Algorithm: 1375 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 1376 * computations. M., Mir, 1980, p.180 of the Russian edition 1377 * x0 1378 **/ 1379 1380 fminbr: function(f, x0, object) { // An estimate to the min location 1381 var a, b, x, v, w, 1382 fx, fv, fw, 1383 r = (3.-Math.sqrt(5.0))*0.5, // Golden section ratio 1384 tol = JXG.Math.eps, 1385 sqrteps = Math.sqrt(JXG.Math.eps), 1386 maxiter = 50, niter = 0, 1387 range, middle_range, tol_act, new_step, 1388 p, q, t, ft, 1389 nfev = 0; 1390 1391 if (!JXG.isArray(x0) || x0.length<2) { 1392 throw new Error("JXG.Math.Numerics.fminbr: length of array x0 has to be at least two."); 1393 } 1394 1395 a = x0[0]; 1396 b = x0[1]; 1397 v = a + r*(b-a); 1398 fv = f.apply(object,[v]); nfev++; // First step - always gold section 1399 x = v; w = v; 1400 fx=fv; fw=fv; 1401 1402 while (niter<maxiter) { 1403 range = b-a; // Range over which the minimum 1404 // is seeked for 1405 middle_range = (a+b)*0.5; 1406 tol_act = sqrteps*Math.abs(x) + tol/3; // Actual tolerance 1407 if( Math.abs(x-middle_range) + range*0.5 <= 2*tol_act ) { 1408 //JXG.debug(nfev); 1409 return x; // Acceptable approx. is found 1410 } 1411 // Obtain the golden section step 1412 new_step = r * ( x<middle_range ? b-x : a-x ); 1413 // Decide if the interpolation can be tried 1414 if ( Math.abs(x-w) >= tol_act ) { // If x and w are distinct 1415 // interpolatiom may be tried 1416 // Interpolation step is calculated as p/q; 1417 // division operation is delayed until last moment 1418 t = (x-w) * (fx-fv); 1419 q = (x-v) * (fx-fw); 1420 p = (x-v)*q - (x-w)*t; 1421 q = 2*(q-t); 1422 1423 if ( q>0 ) { // q was calculated with the op- 1424 p = -p; // posite sign; make q positive 1425 } else { // and assign possible minus to 1426 q = -q; // p 1427 } 1428 if ( Math.abs(p) < Math.abs(new_step*q) && // If x+p/q falls in [a,b] 1429 p > q*(a-x+2*tol_act) && // not too close to a and 1430 p < q*(b-x-2*tol_act) ) { // b, and isn't too large */ 1431 new_step = p/q; // it is accepted 1432 } 1433 // If p/q is too large then the 1434 // golden section procedure can 1435 // reduce [a,b] range to more 1436 // extent 1437 } 1438 1439 if ( Math.abs(new_step) < tol_act ) { // Adjust the step to be not less 1440 if( new_step > 0 ) { // than tolerance 1441 new_step = tol_act; 1442 } else { 1443 new_step = -tol_act; 1444 } 1445 } 1446 1447 // Obtain the next approximation to min 1448 // and reduce the enveloping range 1449 t = x + new_step; // Tentative point for the min 1450 ft = f.apply(object,[t]); nfev++; 1451 if ( ft <= fx ) { // t is a better approximation 1452 if ( t < x ) { // Reduce the range so that 1453 b = x; // t would fall within it 1454 } else { 1455 a = x; 1456 } 1457 v = w; w = x; x = t; // Assign the best approx to x 1458 fv=fw; fw=fx; fx=ft; 1459 } else { // x remains the better approx 1460 if ( t < x ) { // Reduce the range enclosing x 1461 a = t; 1462 } else { 1463 b = t; 1464 } 1465 if ( ft <= fw || w==x ) { 1466 v = w; w = t; 1467 fv=fw; fw=ft; 1468 } else if ( ft<=fv || v==x || v==w ) { 1469 v = t; 1470 fv=ft; 1471 } 1472 } 1473 niter++; 1474 } 1475 JXG.debug("fminbr: maxiter="+maxiter+" reached."); 1476 return x; 1477 }, 1478 1479 /** 1480 * Helper function to create curve which displays Reuleaux polygons. 1481 * @param {array} points Array of points which should be the vertices of the Reuleaux polygon. Typically, 1482 * these point list is the array vrtices of a regular polygon. 1483 * @param {number} nr Number of vertices 1484 * @returns {array} An array containing the two functions defining the Reuleaux polygon and the two values 1485 * for the start and the end of the paramtric curve. 1486 * array may be used as parent array of a {@link JXG.Curve}. 1487 * 1488 * @example 1489 * var A = brd.create('point',[-2,-2]); 1490 * var B = brd.create('point',[0,1]); 1491 * var pol = brd.create('regularpolygon',[A,B,3], {withLines:false, fillColor:'none', highlightFillColor:'none', fillOpacity:0.0}); 1492 * var reuleauxTriangle = brd.create('curve', JXG.Math.Numerics.reuleauxPolygon(pol.vertices, 3), 1493 * {strokeWidth:6, strokeColor:'#d66d55', fillColor:'#ad5544', highlightFillColor:'#ad5544'}); 1494 * 1495 * </pre><div id="2543a843-46a9-4372-abc1-94d9ad2db7ac" style="width: 300px; height: 300px;"></div> 1496 * <script type="text/javascript"> 1497 * var brd = JXG.JSXGraph.initBoard('2543a843-46a9-4372-abc1-94d9ad2db7ac', {boundingbox: [-5, 5, 5, -5], axis: true, showcopyright:false, shownavigation: false}); 1498 * var A = brd.create('point',[-2,-2]); 1499 * var B = brd.create('point',[0,1]); 1500 * var pol = brd.create('regularpolygon',[A,B,3], {withLines:false, fillColor:'none', highlightFillColor:'none', fillOpacity:0.0}); 1501 * var reuleauxTriangle = brd.create('curve', JXG.Math.Numerics.reuleauxPolygon(pol.vertices, 3), 1502 * {strokeWidth:6, strokeColor:'#d66d55', fillColor:'#ad5544', highlightFillColor:'#ad5544'}); 1503 * </script><pre> 1504 */ 1505 reuleauxPolygon: function(points, nr) { 1506 var pi2 = Math.PI*2, 1507 pi2_n = pi2/nr, 1508 diag = (nr-1)/2, 1509 beta, d = 0, 1510 makeFct = function(which, trig) { 1511 return function(t, suspendUpdate) { 1512 if (!suspendUpdate) { 1513 d = points[0].Dist(points[diag]); 1514 beta = JXG.Math.Geometry.rad([points[0].X()+1,points[0].Y()],points[0],points[(diag)%nr]); 1515 } 1516 var t1 = (t%pi2 + pi2) % pi2; 1517 var j = Math.floor(t1 / pi2_n)%nr; 1518 if (isNaN(j)) return j; 1519 //t1 = (t1-j*pi2_n)*0.5 + beta+j*pi2_n; 1520 t1 = t1*0.5+j*pi2_n*0.5 + beta; 1521 return points[j][which]()+d*Math[trig](t1); 1522 }; 1523 }; 1524 return [ 1525 makeFct('X','cos'), 1526 makeFct('Y','sin'), 1527 0, 1528 Math.PI*2 1529 ]; 1530 } 1531 1532 } 1533 })(JXG, Math); 1534 1535