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 cubic Catmull-Rom spline curve through a given set of points. The curve 708 * is uniformly parametrized. 709 * Two artificial control points at the beginning and the end are added. 710 * @param {Array} points Array consisting of JXG.Points. 711 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 712 * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, and a function simply 713 * returning the length of the points array minus three. 714 */ 715 CatmullRomSpline: function(points) { 716 var coeffs = [], 717 p, first = {}, last = {}, // control point at the beginning and at the end 718 719 makeFct = function(which) { 720 return function(t, suspendedUpdate) { 721 var len = points.length, 722 s, c; 723 724 if (len < 2 ) { return NaN; } 725 if (!suspendedUpdate) { 726 first[which] = function() {return 2*points[0][which]()-points[1][which]();} 727 last[which] = function() {return 2*points[len-1][which]()-points[len-2][which]();} 728 p = [first].concat(points, [last]); 729 coeffs[which] = []; 730 for (s=0; s < len-1; s++) { 731 coeffs[which][s] = [ 732 2*p[s+1][which](), 733 -p[s][which]() + p[s+2][which](), 734 2*p[s][which]() - 5*p[s+1][which]() + 4*p[s+2][which]() - p[s+3][which](), 735 -p[s][which]() + 3*p[s+1][which]() - 3*p[s+2][which]() + p[s+3][which]() 736 ]; 737 } 738 } 739 len += 2; // add the two control points 740 if (isNaN(t)) { return NaN; } 741 // This is necessay for our advanced plotting algorithm: 742 if (t < 0) { 743 return p[1][which](); 744 } else if (t >= len - 3) { 745 return p[len-2][which](); 746 } 747 748 s = Math.floor(t); 749 if (s==t) { 750 return p[s][which](); 751 } 752 t -= s; 753 c = coeffs[which][s]; 754 return 0.5*(((c[3]*t + c[2])*t + c[1])*t + c[0]); 755 }; 756 }; 757 758 return [makeFct('X'), 759 makeFct('Y'), 760 0, 761 function() { 762 return points.length - 1; 763 } 764 ]; 765 }, 766 767 768 /** 769 * Computes the regression polynomial of a given degree through a given set of coordinates. 770 * Returns the regression polynomial function. 771 * @param {Number} degree number, function or slider. 772 * Either 773 * @param {Array} dataX array containing the x-coordinates of the data set 774 * @param {Array} dataY array containing the y-coordinates of the data set, 775 * or 776 * @param {Array} dataX Array consisting of JXG.Points. 777 * @param {} dataY Ignored 778 * @returns {function} A function of one parameter which returns the value of the regression polynomial of the given degree. 779 * It possesses the method getTerm() which returns the string containing the function term of the polynomial. 780 */ 781 regressionPolynomial: function(degree, dataX, dataY) { 782 var coeffs, 783 deg, dX, dY, 784 inputType, 785 fct, 786 term = ''; 787 788 if (JXG.isPoint(degree) && typeof degree.Value == 'function') { // Slider 789 deg = function() {return degree.Value();}; 790 } else if (JXG.isFunction(degree)) { 791 deg = degree; 792 } else if (JXG.isNumber(degree)) { 793 deg = function() {return degree;}; 794 } else { 795 throw new Error("JSXGraph: Can't create regressionPolynomial from degree of type'" + (typeof degree) + "'."); 796 } 797 798 if (arguments.length == 3 && JXG.isArray(dataX) && JXG.isArray(dataY)) { // Parameters degree, dataX, dataY 799 inputType = 0; 800 } else if (arguments.length == 2 && JXG.isArray(dataX) && dataX.length > 0 && JXG.isPoint(dataX[0])) { // Parameters degree, point array 801 inputType = 1; 802 } else { 803 throw new Error("JSXGraph: Can't create regressionPolynomial. Wrong parameters."); 804 } 805 806 fct = function(x, suspendedUpdate) { 807 var i, j, M, MT, y, B, c, s, 808 d, len = dataX.length; // input data 809 810 d = Math.floor(deg()); // input data 811 if (!suspendedUpdate) { 812 if (inputType === 1) { // point list as input 813 dX = []; 814 dY = []; 815 for (i = 0; i < len; i++) { 816 dX[i] = dataX[i].X(); 817 dY[i] = dataX[i].Y(); 818 } 819 } 820 821 if (inputType === 0) { // check for functions 822 dX = []; 823 dY = []; 824 for (i = 0; i < len; i++) { 825 if (JXG.isFunction(dataX[i])) 826 dX.push(dataX[i]()); 827 else 828 dX.push(dataX[i]); 829 if (JXG.isFunction(dataY[i])) 830 dY.push(dataY[i]()); 831 else 832 dY.push(dataY[i]); 833 } 834 } 835 836 M = []; 837 for (j = 0; j < len; j++) { 838 M.push([1]); 839 } 840 for (i = 1; i <= d; i++) { 841 for (j = 0; j < len; j++) { 842 M[j][i] = M[j][i - 1] * dX[j]; // input data 843 } 844 } 845 846 y = dY; // input data 847 MT = JXG.Math.transpose(M); 848 B = JXG.Math.matMatMult(MT, M); 849 c = JXG.Math.matVecMult(MT, y); 850 coeffs = JXG.Math.Numerics.Gauss(B, c); 851 term = JXG.Math.Numerics.generatePolynomialTerm(coeffs, d, 'x', 3); 852 } 853 854 // Horner's scheme to evaluate polynomial 855 s = coeffs[d]; 856 for (i = d - 1; i >= 0; i--) { 857 s = (s * x + coeffs[i]); 858 } 859 return s; 860 }; 861 862 fct.getTerm = function() { 863 return term; 864 }; 865 866 return fct; 867 }, 868 869 /** 870 * Computes the cubic Bezier curve through a given set of points. 871 * @param {Array} points Array consisting of 3*k+1 JXG.Points. 872 * The points at position k with k mod 3 = 0 are the data points, 873 * points at position k with k mod 3 = 1 or 2 are the control points. 874 * @returns {Array} An array consisting of two functions of one parameter t which return the 875 * x resp. y coordinates of the Bezier curve in t, one zero value, and a third function accepting 876 * no parameters and returning one third of the length of the points. 877 */ 878 bezier: function(points) { 879 var len, 880 makeFct = function(which) { 881 return function(t, suspendedUpdate) { 882 var z = Math.floor(t) * 3, 883 t0 = t % 1, 884 t1 = 1 - t0; 885 886 if (!suspendedUpdate) { 887 len = Math.floor(points.length / 3); 888 } 889 890 if (t < 0) { 891 return points[0][which](); 892 } 893 if (t >= len) { 894 return points[points.length - 1][which](); 895 } 896 if (isNaN(t)) { 897 return NaN; 898 } 899 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; 900 }; 901 }; 902 903 return [ 904 makeFct('X'), 905 makeFct('Y'), 906 0, 907 function() { 908 return Math.floor(points.length / 3); 909 }]; 910 }, 911 912 /** 913 * Computes the B-spline curve of order k (order = degree+1) through a given set of points. 914 * @param {Array} points Array consisting of JXG.Points. 915 * @param {Number} order Order of the B-spline curve. 916 * @todo closed B-spline curves 917 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 918 * which return the x resp. y coordinates of the B-spline curve in t, a zero value, and a function simply 919 * returning the length of the points array minus one. 920 */ 921 bspline: function(points, order) { 922 var knots, N = [], 923 _knotVector = function(n, k) { 924 var j, kn = []; 925 for (j = 0; j < n + k + 1; j++) { 926 if (j < k) { 927 kn[j] = 0.0; 928 } else if (j <= n) { 929 kn[j] = j - k + 1; 930 } else { 931 kn[j] = n - k + 2; 932 } 933 } 934 return kn; 935 }, 936 937 _evalBasisFuncs = function(t, kn, n, k, s) { 938 var i,j,a,b,den, 939 N = []; 940 941 if (kn[s] <= t && t < kn[s + 1]) { 942 N[s] = 1; 943 } else { 944 N[s] = 0; 945 } 946 for (i = 2; i <= k; i++) { 947 for (j = s - i + 1; j <= s; j++) { 948 if (j <= s - i + 1 || j < 0) { 949 a = 0.0; 950 } else { 951 a = N[j]; 952 } 953 if (j >= s) { 954 b = 0.0; 955 } else { 956 b = N[j + 1]; 957 } 958 den = kn[j + i - 1] - kn[j]; 959 if (den == 0) { 960 N[j] = 0; 961 } else { 962 N[j] = (t - kn[j]) / den * a; 963 } 964 den = kn[j + i] - kn[j + 1]; 965 if (den != 0) { 966 N[j] += (kn[j + i] - t) / den * b; 967 } 968 } 969 } 970 return N; 971 }, 972 makeFct = function(which) { 973 return function(t, suspendedUpdate) { 974 var len = points.length, 975 y, j, s, 976 n = len - 1, 977 k = order; 978 979 if (n <= 0) return NaN; 980 if (n + 2 <= k) k = n + 1; 981 if (t <= 0) return points[0][which](); 982 if (t >= n - k + 2) return points[n][which](); 983 984 knots = _knotVector(n, k); 985 s = Math.floor(t) + k - 1; 986 N = _evalBasisFuncs(t, knots, n, k, s); 987 988 y = 0.0; 989 for (j = s - k + 1; j <= s; j++) { 990 if (j < len && j >= 0) y += points[j][which]() * N[j]; 991 } 992 return y; 993 }; 994 }; 995 996 return [makeFct('X'), 997 makeFct('Y'), 998 0, 999 function() { 1000 return points.length - 1; 1001 } 1002 ]; 1003 }, 1004 1005 /** 1006 * Numerical (symmetric) approximation of derivative. suspendUpdate is piped through, see {@link JXG.Curve#updateCurve} 1007 * and {@link JXG.Curve#hasPoint}. 1008 * @param {function} f Function in one variable to be differentiated. 1009 * @param {object} [obj] Optional object that is treated as "this" in the function body. This is useful, if the function is a 1010 * method of an object and contains a reference to its parent object via "this". 1011 * @returns {function} Derivative function of a given function f. 1012 */ 1013 D: function(f, obj) { 1014 var h = 0.00001, 1015 h2 = 1.0 / (h * 2.0); 1016 1017 if (arguments.length == 1 || (arguments.length > 1 && !JXG.exists(arguments[1]))) { 1018 return function(x, suspendUpdate) { 1019 return (f(x + h, suspendUpdate) - f(x - h, suspendUpdate)) * h2; 1020 }; 1021 } else { // set "this" to "obj" in f 1022 return function(x, suspendUpdate) { 1023 return (f.apply(obj, [x + h,suspendUpdate]) - f.apply(obj, [x - h,suspendUpdate])) * h2; 1024 }; 1025 } 1026 }, 1027 1028 /** 1029 * Helper function to create curve which displays Riemann sums. 1030 * Compute coordinates for the rectangles showing the Riemann sum. 1031 * @param {function} f Function, whose integral is approximated by the Riemann sum. 1032 * @param {Number} n number of rectangles. 1033 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', or 'trapezodial'. 1034 * @param {Number} start Left border of the approximation interval 1035 * @param {Number} end Right border of the approximation interval 1036 * @returns {Array} An array of two arrays containing the x and y coordinates for the rectangles showing the Riemann sum. This 1037 * array may be used as parent array of a {@link JXG.Curve}. 1038 */ 1039 riemann: function(f, n, type, start, end) { 1040 var xarr = [], 1041 yarr = [], 1042 i, j = 0, 1043 delta, 1044 x = start, y, 1045 x1, y1, delta1; 1046 1047 n = Math.round(n); 1048 1049 xarr[j] = x; 1050 yarr[j] = 0.0; 1051 1052 if (n > 0) { 1053 delta = (end - start) / n; 1054 delta1 = delta * 0.01; // for 'lower' and 'upper' 1055 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 === 'left') || (type === 'trapezodial')) { 1062 y = f(x); 1063 } else if (type === 'lower') { 1064 y = f(x); 1065 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1066 y1 = f(x1); 1067 if (y1 < y) { 1068 y = y1; 1069 } 1070 } 1071 } else { // (type=='upper') 1072 y = f(x); 1073 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1074 y1 = f(x1); 1075 if (y1 > y) { 1076 y = y1; 1077 } 1078 } 1079 } 1080 1081 j++; 1082 xarr[j] = x; 1083 yarr[j] = y; 1084 j++; 1085 x += delta; 1086 if (type === 'trapezodial') { 1087 y = f(x); 1088 } 1089 xarr[j] = x; 1090 yarr[j] = y; 1091 j++; 1092 xarr[j] = x; 1093 yarr[j] = 0.0; 1094 } 1095 } 1096 return [xarr,yarr]; 1097 }, 1098 1099 /** 1100 * Approximate the integral by Riemann sums. 1101 * Compute the area described by the riemann sum rectangles. 1102 * @param {function} f Function, whose integral is approximated by the Riemann sum. 1103 * @param {Number} n number of rectangles. 1104 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', or 'trapezodial'. 1105 * @param {Number} start Left border of the approximation interval 1106 * @param {Number} end Right border of the approximation interval 1107 * @returns {Number} The sum of the areas of the rectangles. 1108 */ 1109 riemannsum: function(f, n, type, start, end) { 1110 var sum = .0, 1111 i, delta, 1112 x = start, y, 1113 x1, y1, delta1; 1114 1115 n = Math.floor(n); 1116 if (n > 0) { 1117 delta = (end - start) / n; 1118 delta1 = delta * 0.01; // for 'lower' and 'upper' 1119 for (i = 0; i < n; i++) { 1120 if (type === 'right') { 1121 y = f(x + delta); 1122 } else if (type === 'middle') { 1123 y = f(x + delta * 0.5); 1124 } else if (type === 'trapezodial') { 1125 y = 0.5 * (f(x + delta) + f(x)); 1126 } else if (type === 'left') { 1127 y = f(x); 1128 } else if (type === 'lower') { 1129 y = f(x); 1130 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1131 y1 = f(x1); 1132 if (y1 < y) { 1133 y = y1; 1134 } 1135 } 1136 } else { // (type=='upper') 1137 y = f(x); 1138 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1139 y1 = f(x1); 1140 if (y1 > y) { 1141 y = y1; 1142 } 1143 } 1144 } 1145 sum += delta * y; 1146 x += delta; 1147 } 1148 } 1149 return sum; 1150 }, 1151 1152 /** 1153 * Solve initial value problems numerically using Runge-Kutta-methods. 1154 * See {@link http://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm. 1155 * @param {object,String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing 1156 * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure 1157 * <pre> 1158 * { 1159 * s: <Number>, 1160 * A: <matrix>, 1161 * b: <Array>, 1162 * c: <Array> 1163 * } 1164 * </pre> 1165 * 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 1166 * @param {Array} x0 Initial value vector. If the problem is of one-dimensional, the initial value also has to be given in an array. 1167 * @param {Array} I Interval on which to integrate. 1168 * @param {Number} N Number of evaluation points. 1169 * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode 1170 * 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 1171 * vector <tt>x</tt>, and has to return a vector of the same dimension as <tt>x</tt> has. 1172 * @returns {Array} An array of vectors describing the solution of the ode on the given interval I. 1173 * @example 1174 * // A very simple autonomous system dx(t)/dt = x(t); 1175 * function f(t, x) { 1176 * return x; 1177 * } 1178 * 1179 * // Solve it with initial value x(0) = 1 on the interval [0, 2] 1180 * // with 20 evaluation points. 1181 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 1182 * 1183 * // Prepare data for plotting the solution of the ode using a curve. 1184 * var dataX = []; 1185 * var dataY = []; 1186 * var h = 0.1; // (I[1] - I[0])/N = (2-0)/20 1187 * for(var i=0; i<data.length; i++) { 1188 * dataX[i] = i*h; 1189 * dataY[i] = data[i][0]; 1190 * } 1191 * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'}); 1192 * </pre><div id="d2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div> 1193 * <script type="text/javascript"> 1194 * var board = JXG.JSXGraph.initBoard('d2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false}); 1195 * function f(t, x) { 1196 * // we have to copy the value. 1197 * // return x; would just return the reference. 1198 * return [x[0]]; 1199 * } 1200 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 1201 * var dataX = []; 1202 * var dataY = []; 1203 * var h = 0.1; 1204 * for(var i=0; i<data.length; i++) { 1205 * dataX[i] = i*h; 1206 * dataY[i] = data[i][0]; 1207 * } 1208 * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'}); 1209 * </script><pre> 1210 */ 1211 rungeKutta: function(butcher, x0, I, N, f) { 1212 var x = [], 1213 y = [], 1214 h = (I[1] - I[0]) / N, 1215 t = I[0], 1216 e, i, j, 1217 k, l, 1218 dim = x0.length, 1219 s, 1220 result = [], 1221 r = 0; 1222 1223 if (JXG.isString(butcher)) { 1224 butcher = predefinedButcher[butcher] || predefinedButcher.euler; 1225 } 1226 s = butcher.s; 1227 1228 // don't change x0, so copy it 1229 for (e = 0; e < dim; e++) 1230 x[e] = x0[e]; 1231 1232 for (i = 0; i < N; i++) { 1233 // Optimization doesn't work for ODEs plotted using time 1234 // if((i % quotient == 0) || (i == N-1)) { 1235 result[r] = []; 1236 for (e = 0; e < dim; e++) 1237 result[r][e] = x[e]; 1238 r++; 1239 // } 1240 // init k 1241 k = []; 1242 for (j = 0; j < s; j++) { 1243 // init y = 0 1244 for (e = 0; e < dim; e++) 1245 y[e] = 0.; 1246 1247 // Calculate linear combination of former k's and save it in y 1248 for (l = 0; l < j; l++) { 1249 for (e = 0; e < dim; e++) { 1250 y[e] += (butcher.A[j][l]) * h * k[l][e]; 1251 } 1252 } 1253 1254 // add x(t) to y 1255 for (e = 0; e < dim; e++) { 1256 y[e] += x[e]; 1257 } 1258 1259 // calculate new k and add it to the k matrix 1260 k.push(f(t + butcher.c[j] * h, y)); 1261 } 1262 1263 // init y = 0 1264 for (e = 0; e < dim; e++) 1265 y[e] = 0.; 1266 1267 for (l = 0; l < s; l++) { 1268 for (e = 0; e < dim; e++) 1269 y[e] += butcher.b[l] * k[l][e]; 1270 } 1271 1272 for (e = 0; e < dim; e++) { 1273 x[e] = x[e] + h * y[e]; 1274 } 1275 1276 t += h; 1277 } 1278 1279 return result; 1280 }, 1281 1282 1283 /* 1284 * Maximum number of iterations in @see #fzero 1285 */ 1286 maxIterationsRoot: 80, 1287 1288 /* 1289 * Maximum number of iterations in @see #fminbr 1290 */ 1291 maxIterationsMinimize: 500, 1292 1293 /** 1294 * 1295 * Find zero of an univariate function f. 1296 * @param {function} f Function, whose root is to be found 1297 * @param {array or number} x0 Start value or start interval enclosing the root 1298 * @param {object} object Parent object in case f is method of it 1299 * @return {number} the approximation of the root 1300 * Algorithm: 1301 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 1302 * computations. M., Mir, 1980, p.180 of the Russian edition 1303 * 1304 * if x0 is an array containing lower and upper bound for the zero 1305 * algorithm 748 is applied. Otherwise, if x0 is a number, 1306 * the algorithm tries to bracket a zero of f starting from x0. 1307 * If this fails, we fall back to Newton's method. 1308 * 1309 **/ 1310 fzero: function(f, x0, object) { 1311 var tol = JXG.Math.eps, 1312 maxiter = this.maxIterationsRoot, niter = 0, 1313 nfev = 0, 1314 eps = tol, 1315 a,b,c, 1316 fa, fb, fc, 1317 aa, blist, i, len, u, fu, 1318 prev_step, 1319 tol_act, // Actual tolerance 1320 p, // Interpolation step is calcu- 1321 q, // lated in the form p/q; divi- 1322 // sion operations is delayed 1323 // until the last moment 1324 new_step, // Step at this iteration 1325 t1, cb, t2; 1326 1327 if (JXG.isArray(x0)) { 1328 if (x0.length<2) 1329 throw new Error("JXG.Math.Numerics.fzero: length of array x0 has to be at least two."); 1330 a = x0[0]; fa = f.apply(object,[a]); nfev++; 1331 b = x0[1]; fb = f.apply(object,[b]); nfev++; 1332 } else { 1333 a = x0; fa = f.apply(object,[a]); nfev++; 1334 // Try to get b. 1335 if (a == 0) { 1336 aa = 1; 1337 } else { 1338 aa = a; 1339 } 1340 blist = [0.9*aa, 1.1*aa, aa-1, aa+1, 0.5*aa, 1.5*aa, -aa, 2*aa, -10*aa, 10*aa]; 1341 len = blist.length; 1342 for (i=0;i<len;i++) { 1343 b = blist[i]; 1344 fb = f.apply(object,[b]); nfev++; 1345 if (fa*fb<=0) { 1346 break; 1347 } 1348 } 1349 if (b < a) { 1350 u = a; a = b; b = u; 1351 fu = fa; fa = fb; fb = fu; 1352 } 1353 } 1354 1355 if (fa*fb > 0) { 1356 // Bracketing not successful, fall back to Newton's method or to fminbr 1357 if (JXG.isArray(x0)) { 1358 //JXG.debug("fzero falls back to fminbr"); 1359 return this.fminbr(f, [a,b], object); 1360 } else { 1361 //JXG.debug("fzero falls back to Newton"); 1362 return this.Newton(f, a, object); 1363 } 1364 } 1365 1366 // OK, we have enclosed a zero of f. 1367 // Now we can start Brent's method 1368 1369 c = a; 1370 fc = fa; 1371 while (niter<maxiter) { // Main iteration loop 1372 prev_step = b-a; // Distance from the last but one 1373 // to the last approximation 1374 1375 if ( Math.abs(fc) < Math.abs(fb) ) { 1376 // Swap data for b to be the 1377 a = b; b = c; c = a; // best approximation 1378 fa=fb; fb=fc; fc=fa; 1379 } 1380 tol_act = 2*eps*Math.abs(b) + tol*0.5; 1381 new_step = (c-b)*0.5; 1382 1383 if ( Math.abs(new_step) <= tol_act && Math.abs(fb) <= eps ) { 1384 //JXG.debug("nfev="+nfev); 1385 return b; // Acceptable approx. is found 1386 } 1387 1388 // Decide if the interpolation can be tried 1389 if ( Math.abs(prev_step) >= tol_act // If prev_step was large enough 1390 && Math.abs(fa) > Math.abs(fb) ) { // and was in true direction, 1391 // Interpolatiom may be tried 1392 cb = c-b; 1393 if ( a==c ) { // If we have only two distinct 1394 // points linear interpolation 1395 t1 = fb/fa; // can only be applied 1396 p = cb*t1; 1397 q = 1.0 - t1; 1398 } else { // Quadric inverse interpolation 1399 q = fa/fc; t1 = fb/fc; t2 = fb/fa; 1400 p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) ); 1401 q = (q-1.0) * (t1-1.0) * (t2-1.0); 1402 } 1403 if ( p>0 ) { // p was calculated with the op- 1404 q = -q; // posite sign; make p positive 1405 } else { // and assign possible minus to 1406 p = -p; // q 1407 } 1408 1409 if( p < (0.75*cb*q-Math.abs(tol_act*q)*0.5) // If b+p/q falls in [b,c] 1410 && p < Math.abs(prev_step*q*0.5) ) { // and isn't too large 1411 new_step = p/q; // it is accepted 1412 } // If p/q is too large then the 1413 // bissection procedure can 1414 // reduce [b,c] range to more 1415 // extent 1416 } 1417 1418 if ( Math.abs(new_step) < tol_act ) { // Adjust the step to be not less 1419 if ( new_step > 0 ) { // than tolerance 1420 new_step = tol_act; 1421 } else { 1422 new_step = -tol_act; 1423 } 1424 } 1425 1426 a = b; fa = fb; // Save the previous approx. 1427 b += new_step; 1428 fb = f.apply(object,[b]); nfev++; // Do step to a new approxim. 1429 if ( (fb>0 && fc>0) || (fb<0 && fc<0) ) { 1430 // Adjust c for it to have a sign 1431 c = a; fc = fa; // opposite to that of b 1432 } 1433 niter++; 1434 } // End while 1435 1436 //JXG.debug("fzero: maxiter="+maxiter+" reached."); 1437 return b; 1438 }, 1439 1440 1441 /** 1442 * 1443 * Find minimum of an univariate function f. 1444 * @param {function} f Function, whose minimum is to be found 1445 * @param {array} x0 Start interval enclosing the minimum 1446 * @param {object} object Parent object in case f is method of it 1447 * @return {number} the approximation of the minimum 1448 * Algorithm: 1449 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 1450 * computations. M., Mir, 1980, p.180 of the Russian edition 1451 * x0 1452 **/ 1453 fminbr: function(f, x0, object) { // An estimate to the min location 1454 var a, b, x, v, w, 1455 fx, fv, fw, 1456 r = (3.-Math.sqrt(5.0))*0.5, // Golden section ratio 1457 tol = JXG.Math.eps, 1458 sqrteps = Math.sqrt(JXG.Math.eps), 1459 maxiter = this.maxIterationsMinimize, 1460 niter = 0, 1461 range, middle_range, tol_act, new_step, 1462 p, q, t, ft, 1463 nfev = 0; 1464 1465 if (!JXG.isArray(x0) || x0.length<2) { 1466 throw new Error("JXG.Math.Numerics.fminbr: length of array x0 has to be at least two."); 1467 } 1468 a = x0[0]; 1469 b = x0[1]; 1470 v = a + r*(b-a); 1471 fv = f.apply(object,[v]); nfev++; // First step - always gold section 1472 x = v; w = v; 1473 fx=fv; fw=fv; 1474 1475 while (niter<maxiter) { 1476 range = b-a; // Range over which the minimum 1477 // is seeked for 1478 middle_range = (a+b)*0.5; 1479 tol_act = sqrteps*Math.abs(x) + tol/3; // Actual tolerance 1480 if( Math.abs(x-middle_range) + range*0.5 <= 2*tol_act ) { 1481 //JXG.debug(nfev); 1482 return x; // Acceptable approx. is found 1483 } 1484 // Obtain the golden section step 1485 new_step = r * ( x<middle_range ? b-x : a-x ); 1486 // Decide if the interpolation can be tried 1487 if ( Math.abs(x-w) >= tol_act ) { // If x and w are distinct 1488 // interpolatiom may be tried 1489 // Interpolation step is calculated as p/q; 1490 // division operation is delayed until last moment 1491 t = (x-w) * (fx-fv); 1492 q = (x-v) * (fx-fw); 1493 p = (x-v)*q - (x-w)*t; 1494 q = 2*(q-t); 1495 1496 if ( q>0 ) { // q was calculated with the op- 1497 p = -p; // posite sign; make q positive 1498 } else { // and assign possible minus to 1499 q = -q; // p 1500 } 1501 if ( Math.abs(p) < Math.abs(new_step*q) && // If x+p/q falls in [a,b] 1502 p > q*(a-x+2*tol_act) && // not too close to a and 1503 p < q*(b-x-2*tol_act) ) { // b, and isn't too large */ 1504 new_step = p/q; // it is accepted 1505 } 1506 // If p/q is too large then the 1507 // golden section procedure can 1508 // reduce [a,b] range to more 1509 // extent 1510 } 1511 1512 if ( Math.abs(new_step) < tol_act ) { // Adjust the step to be not less 1513 if( new_step > 0 ) { // than tolerance 1514 new_step = tol_act; 1515 } else { 1516 new_step = -tol_act; 1517 } 1518 } 1519 1520 // Obtain the next approximation to min 1521 // and reduce the enveloping range 1522 t = x + new_step; // Tentative point for the min 1523 ft = f.apply(object,[t]); nfev++; 1524 if ( ft <= fx ) { // t is a better approximation 1525 if ( t < x ) { // Reduce the range so that 1526 b = x; // t would fall within it 1527 } else { 1528 a = x; 1529 } 1530 v = w; w = x; x = t; // Assign the best approx to x 1531 fv=fw; fw=fx; fx=ft; 1532 } else { // x remains the better approx 1533 if ( t < x ) { // Reduce the range enclosing x 1534 a = t; 1535 } else { 1536 b = t; 1537 } 1538 if ( ft <= fw || w==x ) { 1539 v = w; w = t; 1540 fv=fw; fw=ft; 1541 } else if ( ft<=fv || v==x || v==w ) { 1542 v = t; 1543 fv=ft; 1544 } 1545 } 1546 niter++; 1547 } 1548 //JXG.debug("fminbr: maxiter="+maxiter+" reached."); 1549 return x; 1550 }, 1551 1552 /** 1553 * Helper function to create curve which displays Reuleaux polygons. 1554 * @param {array} points Array of points which should be the vertices of the Reuleaux polygon. Typically, 1555 * these point list is the array vrtices of a regular polygon. 1556 * @param {number} nr Number of vertices 1557 * @returns {array} An array containing the two functions defining the Reuleaux polygon and the two values 1558 * for the start and the end of the paramtric curve. 1559 * array may be used as parent array of a {@link JXG.Curve}. 1560 * 1561 * @example 1562 * var A = brd.create('point',[-2,-2]); 1563 * var B = brd.create('point',[0,1]); 1564 * var pol = brd.create('regularpolygon',[A,B,3], {withLines:false, fillColor:'none', highlightFillColor:'none', fillOpacity:0.0}); 1565 * var reuleauxTriangle = brd.create('curve', JXG.Math.Numerics.reuleauxPolygon(pol.vertices, 3), 1566 * {strokeWidth:6, strokeColor:'#d66d55', fillColor:'#ad5544', highlightFillColor:'#ad5544'}); 1567 * 1568 * </pre><div id="2543a843-46a9-4372-abc1-94d9ad2db7ac" style="width: 300px; height: 300px;"></div> 1569 * <script type="text/javascript"> 1570 * var brd = JXG.JSXGraph.initBoard('2543a843-46a9-4372-abc1-94d9ad2db7ac', {boundingbox: [-5, 5, 5, -5], axis: true, showcopyright:false, shownavigation: false}); 1571 * var A = brd.create('point',[-2,-2]); 1572 * var B = brd.create('point',[0,1]); 1573 * var pol = brd.create('regularpolygon',[A,B,3], {withLines:false, fillColor:'none', highlightFillColor:'none', fillOpacity:0.0}); 1574 * var reuleauxTriangle = brd.create('curve', JXG.Math.Numerics.reuleauxPolygon(pol.vertices, 3), 1575 * {strokeWidth:6, strokeColor:'#d66d55', fillColor:'#ad5544', highlightFillColor:'#ad5544'}); 1576 * </script><pre> 1577 */ 1578 reuleauxPolygon: function(points, nr) { 1579 var pi2 = Math.PI*2, 1580 pi2_n = pi2/nr, 1581 diag = (nr-1)/2, 1582 beta, d = 0, 1583 makeFct = function(which, trig) { 1584 return function(t, suspendUpdate) { 1585 if (!suspendUpdate) { 1586 d = points[0].Dist(points[diag]); 1587 beta = JXG.Math.Geometry.rad([points[0].X()+1,points[0].Y()],points[0],points[(diag)%nr]); 1588 } 1589 var t1 = (t%pi2 + pi2) % pi2; 1590 var j = Math.floor(t1 / pi2_n)%nr; 1591 if (isNaN(j)) return j; 1592 //t1 = (t1-j*pi2_n)*0.5 + beta+j*pi2_n; 1593 t1 = t1*0.5+j*pi2_n*0.5 + beta; 1594 return points[j][which]()+d*Math[trig](t1); 1595 }; 1596 }; 1597 return [ 1598 makeFct('X','cos'), 1599 makeFct('Y','sin'), 1600 0, 1601 Math.PI*2 1602 ]; 1603 }, 1604 1605 /** 1606 * Implements the Ramer-Douglas-Peuker algorithm. 1607 * It discards points which are not necessary from the polygonal line defined by the point array 1608 * pts. The computation is done in screen coordinates. 1609 * Average runtime is O(nlog(n)), worst case runtime is O(n^2), where n is the number of points. 1610 * @param {Array} pts Array of {@link JXG.Coords} 1611 * @param {Number} eps If the absolute value of a given number <tt>x</tt> is smaller than <tt>eps</tt> it is considered to be equal <tt>0</tt>. 1612 * @returns {Array} An array containing points which represent an apparently identical curve as the points of pts do, but contains fewer points. 1613 */ 1614 RamerDouglasPeuker: function(pts, eps) { 1615 var newPts = [], i, k, len, 1616 /** 1617 * RDP() is a private subroutine of {@link JXG.Math.Numerics#RamerDouglasPeuker}. 1618 * It runs recursively through the point set and searches the 1619 * point which has the largest distance from the line between the first point and 1620 * the last point. If the distance from the line is greater than eps, this point is 1621 * included in our new point set otherwise it is discarded. 1622 * If it is taken, we recursively apply the subroutine to the point set before 1623 * and after the chosen point. 1624 * @param {Array} pts Array of {@link JXG.Coords} 1625 * @param {Number} i Index of an element of pts 1626 * @param {Number} j Index of an element of pts 1627 * @param {Number} eps If the absolute value of a given number <tt>x</tt> is smaller than <tt>eps</tt> it is considered to be equal <tt>0</tt>. 1628 * @param {Array} newPts Array of {@link JXG.Coords} 1629 * @private 1630 */ 1631 RDP = function(pts, i, j, eps, newPts) { 1632 var result = findSplit(pts, i, j); 1633 1634 if (result[0] > eps) { 1635 RDP(pts, i, result[1], eps, newPts); 1636 RDP(pts, result[1], j, eps, newPts); 1637 } else { 1638 newPts.push(pts[j]); 1639 } 1640 }, 1641 /** 1642 * findSplit() is a subroutine of {@link JXG.Math.Numerics#RamerDouglasPeuker}. 1643 * It searches for the point between index i and j which 1644 * has the largest distance from the line between the points i and j. 1645 * @param {Array} pts Array of {@link JXG.Coords} 1646 * @param {Number} i Index of a point in pts 1647 * @param {Number} j Index of a point in pts 1648 **/ 1649 findSplit = function(pts, i, j) { 1650 var dist = 0, 1651 f = i, 1652 d, k, ci, cj, ck, 1653 x0, y0, x1, y1, 1654 den, lbda; 1655 1656 if (j - i < 2) return [-1.0,0]; 1657 1658 ci = pts[i].scrCoords; 1659 cj = pts[j].scrCoords; 1660 if (isNaN(ci[1] + ci[2] + cj[1] + cj[2])) return [NaN,j]; 1661 1662 for (k = i + 1; k < j; k++) { 1663 ck = pts[k].scrCoords; 1664 x0 = ck[1] - ci[1]; 1665 y0 = ck[2] - ci[2]; 1666 x1 = cj[1] - ci[1]; 1667 y1 = cj[2] - ci[2]; 1668 den = x1 * x1 + y1 * y1; 1669 /* 1670 if (den >= JXG.Math.eps) { 1671 lbda = (x0 * x1 + y0 * y1) / den; 1672 d = x0 * x0 + y0 * y0 - lbda * (x0 * x1 + y0 * y1); 1673 } else { 1674 lbda = 0.0; 1675 d = x0 * x0 + y0 * y0; 1676 } 1677 if (lbda < 0.0) { 1678 d = x0 * x0 + y0 * y0; 1679 } else if (lbda > 1.0) { 1680 x0 = ck[1] - cj[1]; 1681 y0 = ck[2] - cj[2]; 1682 d = x0 * x0 + y0 * y0; 1683 } 1684 */ 1685 if (den >= JXG.Math.eps) { 1686 lbda = (x0 * x1 + y0 * y1) / den; 1687 //d = x0 * x0 + y0 * y0 - lbda * (x0 * x1 + y0 * y1); 1688 if (lbda<0.0) { 1689 lbda = 0.0; 1690 } else if (lbda>1.0) { 1691 lbda = 1.0; 1692 } 1693 x0 = x0-lbda*x1; 1694 y0 = y0-lbda*y1; 1695 d = x0*x0+y0*y0; 1696 } else { 1697 lbda = 0.0; 1698 d = x0 * x0 + y0 * y0; 1699 } 1700 1701 if (d > dist) { 1702 dist = d; 1703 f = k; 1704 } 1705 } 1706 return [Math.sqrt(dist),f]; 1707 }; 1708 1709 len = pts.length; 1710 1711 // Search for the left most point woithout NaN coordinates 1712 i = 0; 1713 while (i < len && isNaN(pts[i].scrCoords[1] + pts[i].scrCoords[2])) { 1714 i++; 1715 } 1716 // Search for the right most point woithout NaN coordinates 1717 k = len - 1; 1718 while (k > i && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])) { 1719 k--; 1720 } 1721 1722 // Only proceed if something is left 1723 if (!(i > k || i == len)) { 1724 newPts[0] = pts[i]; 1725 RDP(pts, i, k, eps, newPts); 1726 } 1727 1728 return newPts; 1729 } 1730 1731 } 1732 })(JXG, Math); 1733 1734