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 var n = mat.length; 232 if (n==2 && mat[0].length==2) { 233 return mat[0][0]*mat[1][1] - mat[1][0]*mat[0][1]; 234 } else { 235 return this.gaussBareiss(mat); 236 } 237 }, 238 239 /** 240 * Compute the Eigenvalues and Eigenvectors of a symmetric 3x3 matrix with the Jacobi method 241 * Adaption of a FORTRAN program by Ed Wilson, Dec. 25, 1990 242 * @param {Array} Ain A symmetric 3x3 matrix. 243 * @returns {Array} [A,V] the matrices A and V. The diagonal of A contains the Eigenvalues, V contains the Eigenvectors. 244 */ 245 Jacobi: function(Ain) { 246 var i,j,k,aa,si,co,tt,eps = JXG.Math.eps, 247 sum = 0.0, 248 ssum, amax, 249 n = Ain.length, 250 V = [ 251 [0,0,0], 252 [0,0,0], 253 [0,0,0] 254 ], 255 A = [ 256 [0,0,0], 257 [0,0,0], 258 [0,0,0] 259 ], nloops=0; 260 261 // Initialization. Set initial Eigenvectors. 262 for (i = 0; i < n; i++) { 263 for (j = 0; j < n; j++) { 264 V[i][j] = 0.0; 265 A[i][j] = Ain[i][j]; 266 sum += Math.abs(A[i][j]); 267 } 268 V[i][i] = 1.0; 269 } 270 // Trivial problems 271 if (n == 1) { 272 return [A,V]; 273 } 274 if (sum <= 0.0) { 275 return [A,V]; 276 } 277 278 sum /= (n * n); 279 280 // Reduce matrix to diagonal 281 do { 282 ssum = 0.0; 283 amax = 0.0; 284 for (j = 1; j < n; j++) { 285 for (i = 0; i < j; i++) { 286 // Check if A[i][j] is to be reduced 287 aa = Math.abs(A[i][j]); 288 if (aa > amax) { 289 amax = aa; 290 } 291 ssum += aa; 292 if (aa >= eps) { 293 // calculate rotation angle 294 aa = Math.atan2(2.0 * A[i][j], A[i][i] - A[j][j]) * 0.5; 295 si = Math.sin(aa); 296 co = Math.cos(aa); 297 // Modify 'i' and 'j' columns 298 for (k = 0; k < n; k++) { 299 tt = A[k][i]; 300 A[k][i] = co * tt + si * A[k][j]; 301 A[k][j] = -si * tt + co * A[k][j]; 302 tt = V[k][i]; 303 V[k][i] = co * tt + si * V[k][j]; 304 V[k][j] = -si * tt + co * V[k][j]; 305 } 306 // Modify diagonal terms 307 A[i][i] = co * A[i][i] + si * A[j][i]; 308 A[j][j] = -si * A[i][j] + co * A[j][j]; 309 A[i][j] = 0.0; 310 // Make 'A' matrix symmetrical 311 for (k = 0; k < n; k++) { 312 A[i][k] = A[k][i]; 313 A[j][k] = A[k][j]; 314 } 315 // A[i][j] made zero by rotation 316 } 317 } 318 } 319 nloops++; 320 } while (Math.abs(ssum) / sum > eps && nloops<2000); 321 return [A,V]; 322 }, 323 324 /** 325 * Calculates the integral of function f over interval using Newton-Cotes-algorithm. 326 * @param {Array} interval The integration interval, e.g. [0, 3]. 327 * @param {function} f A function which takes one argument of type number and returns a number. 328 * @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 329 * with value being either 'trapez', 'simpson', or 'milne'. 330 * @returns {Number} Integral value of f over interval 331 * @throws {Error} If config.number_of_nodes doesn't match config.integration_type an exception is thrown. If you want to use 332 * simpson rule respectively milne rule config.number_of_nodes must be dividable by 2 respectively 4. 333 * @example 334 * function f(x) { 335 * return x*x; 336 * } 337 * 338 * // calculates integral of <tt>f</tt> from 0 to 2. 339 * var area1 = JXG.Math.Numerics.NewtonCotes([0, 2], f); 340 * 341 * // the same with an anonymous function 342 * var area2 = JXG.Math.Numerics.NewtonCotes([0, 2], function (x) { return x*x; }); 343 * 344 * // use trapez rule with 16 nodes 345 * var area3 = JXG.Math.Numerics.NewtonCotes([0, 2], f, 346 * {number_of_nodes: 16, intergration_type: 'trapez'}); 347 */ 348 NewtonCotes: function(interval, f, config) { 349 var integral_value = 0.0, 350 number_of_nodes = config && typeof config.number_of_nodes === 'number' ? config.number_of_nodes : 28, 351 available_types = {trapez: true, simpson: true, milne: true}, 352 integration_type = config && config.integration_type && available_types.hasOwnProperty(config.integration_type) && available_types[config.integration_type] ? config.integration_type : 'milne', 353 step_size = (interval[1] - interval[0]) / number_of_nodes, 354 evaluation_point, i, number_of_intervals; 355 356 switch (integration_type) { 357 case 'trapez': 358 integral_value = (f(interval[0]) + f(interval[1])) * 0.5; 359 360 evaluation_point = interval[0]; 361 for (i = 0; i < number_of_nodes - 1; i++) { 362 evaluation_point += step_size; 363 integral_value += f(evaluation_point); 364 } 365 integral_value *= step_size; 366 367 break; 368 case 'simpson': 369 if (number_of_nodes % 2 > 0) { 370 throw new Error("JSXGraph: INT_SIMPSON requires config.number_of_nodes dividable by 2."); 371 } 372 number_of_intervals = number_of_nodes / 2.0; 373 integral_value = f(interval[0]) + f(interval[1]); 374 evaluation_point = interval[0]; 375 for (i = 0; i < number_of_intervals - 1; i++) { 376 evaluation_point += 2.0 * step_size; 377 integral_value += 2.0 * f(evaluation_point); 378 } 379 evaluation_point = interval[0] - step_size; 380 for (i = 0; i < number_of_intervals; i++) { 381 evaluation_point += 2.0 * step_size; 382 integral_value += 4.0 * f(evaluation_point); 383 } 384 integral_value *= step_size / 3.0; 385 break; 386 default: 387 if (number_of_nodes % 4 > 0) { 388 throw new Error("JSXGraph: Error in INT_MILNE: config.number_of_nodes must be a multiple of 4"); 389 } 390 number_of_intervals = number_of_nodes * 0.25; 391 integral_value = 7.0 * (f(interval[0]) + f(interval[1])); 392 evaluation_point = interval[0]; 393 for (i = 0; i < number_of_intervals - 1; i++) { 394 evaluation_point += 4.0 * step_size; 395 integral_value += 14.0 * f(evaluation_point); 396 } 397 evaluation_point = interval[0] - 3.0 * step_size; 398 for (i = 0; i < number_of_intervals; i++) { 399 evaluation_point += 4.0 * step_size; 400 integral_value += 32.0 * (f(evaluation_point) + f(evaluation_point + 2 * step_size)); 401 } 402 evaluation_point = interval[0] - 2.0 * step_size; 403 for (i = 0; i < number_of_intervals; i++) { 404 evaluation_point += 4.0 * step_size; 405 integral_value += 12.0 * f(evaluation_point); 406 } 407 integral_value *= 2.0 * step_size / 45.0; 408 } 409 return integral_value; 410 }, 411 412 /** 413 * Integral of function f over interval. 414 * @param {Array} interval The integration interval, e.g. [0, 3]. 415 * @param {function} f A function which takes one argument of type number and returns a number. 416 * @returns {Number} The value of the integral of f over interval 417 * @see JXG.Math.Numerics.NewtonCotes 418 */ 419 I: function(interval, f) { 420 return this.NewtonCotes(interval, f, {number_of_nodes: 16, integration_type: 'milne'}); 421 }, 422 423 /** 424 * Newton's method to find roots of a funtion in one variable. 425 * @param {function} f We search for a solution of f(x)=0. 426 * @param {Number} x initial guess for the root, i.e. start value. 427 * @param {object} object optional object that is treated as "this" in the function body. This is useful, if the function is a 428 * method of an object and contains a reference to its parent object via "this". 429 * @returns {Number} A root of the function f. 430 */ 431 Newton: function(f, x, object) { 432 var i = 0, 433 h = JXG.Math.eps, 434 newf = f.apply(object, [x]), // set "this" to "object" in f 435 nfev = 1, 436 df; 437 if (JXG.isArray(x)) { // For compatibility 438 x = x[0]; 439 } 440 while (i < 50 && Math.abs(newf) > h) { 441 df = this.D(f, object)(x); nfev += 2; 442 if (Math.abs(df) > h) { 443 x -= newf / df; 444 } else { 445 x += (Math.random() * 0.2 - 1.0); 446 } 447 newf = f.apply(object, [x]); nfev++; 448 i++; 449 } 450 //JXG.debug("N nfev="+nfev); 451 return x; 452 }, 453 454 /** 455 * Abstract method to find roots of univariate functions. 456 * @param {function} f We search for a solution of f(x)=0. 457 * @param {Number} x initial guess for the root, i.e. staring value. 458 * @param {object} object optional object that is treated as "this" in the function body. This is useful, if the function is a 459 * method of an object and contains a reference to its parent object via "this". 460 * @returns {Number} A root of the function f. 461 */ 462 root: function(f, x, object) { 463 //return this.Newton(f, x, object); 464 return this.fzero(f, x, object); 465 }, 466 467 /** 468 * Returns the Lagrange polynomials for curves with equidistant nodes, see 469 * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation, 470 * SIAM Review, Vol 46, No 3, (2004) 501-517. 471 * The graph of the parametric curve [x(t),y(t)] runs through the given points. 472 * @param {Array} p Array of JXG.Points 473 * @returns {Array} An array consisting of two functions x(t), y(t) which define a parametric curve 474 * f(t) = (x(t), y(t)) and two numbers x1 and x2 defining the curve's domain. x1 always equals zero. 475 */ 476 Neville: function(p) { 477 var w = [], 478 makeFct = function(fun) { 479 return function(t, suspendedUpdate) { 480 var i, d, s, 481 bin = JXG.Math.binomial, 482 len = p.length, 483 len1 = len - 1, 484 num = 0.0, 485 denom = 0.0; 486 487 if (!suspendedUpdate) { 488 s = 1; 489 for (i = 0; i < len; i++) { 490 w[i] = bin(len1, i) * s; 491 s *= (-1); 492 } 493 } 494 495 d = t; 496 for (i = 0; i < len; i++) { 497 if (d === 0) { 498 return p[i][fun](); 499 } else { 500 s = w[i] / d; 501 d--; 502 num += p[i][fun]() * s; 503 denom += s; 504 } 505 } 506 return num / denom; 507 }; 508 }, 509 510 xfct = makeFct('X'), 511 yfct = makeFct('Y'); 512 513 return [xfct, yfct, 0, function() { 514 return p.length - 1; 515 }]; 516 }, 517 518 /** 519 * Calculates second derivatives at the given knots. 520 * @param {Array} x x values of knots 521 * @param {Array} y y values of knots 522 * @returns {Array} Second derivatives of the interpolated function at the knots. 523 * @see #splineEval 524 */ 525 splineDef: function(x, y) { 526 var n = Math.min(x.length, y.length), 527 pair, i, l, 528 diag = [], 529 z = [], 530 data = [], 531 dx = [], 532 delta = [], 533 F = []; 534 535 if (n === 2) { 536 return [0, 0]; 537 } 538 539 for (i = 0; i < n; i++) { 540 pair = {X: x[i], Y: y[i]}; 541 data.push(pair); 542 } 543 data.sort(function (a, b) { 544 return a.X - b.X; 545 }); 546 for (i = 0; i < n; i++) { 547 x[i] = data[i].X; 548 y[i] = data[i].Y; 549 } 550 551 for (i = 0; i < n - 1; i++) { 552 dx.push(x[i + 1] - x[i]); 553 } 554 for (i = 0; i < n - 2; i++) { 555 delta.push(6 * (y[i + 2] - y[i + 1]) / (dx[i + 1]) - 6 * (y[i + 1] - y[i]) / (dx[i])); 556 } 557 558 // ForwardSolve 559 diag.push(2 * (dx[0] + dx[1])); 560 z.push(delta[0]); 561 562 for (i = 0; i < n - 3; i++) { 563 l = dx[i + 1] / diag[i]; 564 diag.push(2 * (dx[i + 1] + dx[i + 2]) - l * dx[i + 1]); 565 z.push(delta[i + 1] - l * z[i]); 566 } 567 568 // BackwardSolve 569 F[n - 3] = z[n - 3] / diag[n - 3]; 570 for (i = n - 4; i >= 0; i--) { 571 F[i] = (z[i] - (dx[i + 1] * F[i + 1])) / diag[i]; 572 } 573 574 // Generate f''-Vector 575 for (i = n - 3; i >= 0; i--) { 576 F[i + 1] = F[i]; 577 } 578 579 // natural cubic spline 580 F[0] = 0; 581 F[n - 1] = 0; 582 return F; 583 }, 584 585 /** 586 * Evaluate points on spline. 587 * @param {Number,Array} x0 A single float value or an array of values to evaluate 588 * @param {Array} x x values of knots 589 * @param {Array} y y values of knots 590 * @param {Array} F Second derivatives at knots, calculated by {@link #splineDef} 591 * @see #splineDef 592 * @returns {Number,Array} A single value or an array, depending on what is given as x0. 593 */ 594 splineEval: function(x0, x, y, F) { 595 var n = Math.min(x.length, y.length), 596 l = 1, 597 asArray = false, 598 y0 = [], 599 i, j, a, b, c, d, x_; 600 601 // number of points to be evaluated 602 if (JXG.isArray(x0)) { 603 l = x0.length; 604 asArray = true; 605 } else 606 x0 = [x0]; 607 608 for (i = 0; i < l; i++) { 609 // is x0 in defining interval? 610 if ((x0[i] < x[0]) || (x[i] > x[n - 1])) 611 return NaN; 612 613 // determine part of spline in which x0 lies 614 for (j = 1; j < n; j++) { 615 if (x0[i] <= x[j]) 616 break; 617 } 618 j--; 619 620 // we're now in the j-th partial interval, i.e. x[j] < x0[i] <= x[j+1]; 621 // determine the coefficients of the polynomial in this interval 622 a = y[j]; 623 b = (y[j + 1] - y[j]) / (x[j + 1] - x[j]) - (x[j + 1] - x[j]) / 6 * (F[j + 1] + 2 * F[j]); 624 c = F[j] / 2; 625 d = (F[j + 1] - F[j]) / (6 * (x[j + 1] - x[j])); 626 // evaluate x0[i] 627 x_ = x0[i] - x[j]; 628 //y0.push(a + b*x_ + c*x_*x_ + d*x_*x_*x_); 629 y0.push(a + (b + (c + d * x_) * x_) * x_); 630 } 631 632 if (asArray) 633 return y0; 634 else 635 return y0[0]; 636 }, 637 638 /** 639 * Generate a string containing the function term of a polynomial. 640 * @param {Array} coeffs Coefficients of the polynomial. The position i belongs to x^i. 641 * @param {Number} deg Degree of the polynomial 642 * @param {String} varname Name of the variable (usually 'x') 643 * @param {Number} prec Precision 644 * @returns {String} A string containg the function term of the polynomial. 645 */ 646 generatePolynomialTerm: function(coeffs, deg, varname, prec) { 647 var t = [], i; 648 for (i = deg; i >= 0; i--) { 649 t = t.concat(['(', coeffs[i].toPrecision(prec), ')']); 650 if (i > 1) { 651 t = t.concat(['*', varname, '<sup>', i, '<', '/sup> + ']); 652 } else if (i === 1) { 653 t = t.concat(['*', varname, ' + ']); 654 } 655 } 656 return t.join(''); 657 }, 658 659 /** 660 * Computes the polynomial through a given set of coordinates in Lagrange form. 661 * Returns the Lagrange polynomials, see 662 * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation, 663 * SIAM Review, Vol 46, No 3, (2004) 501-517. 664 * @param {Array} p Array of JXG.Points 665 * @returns {function} A function of one parameter which returns the value of the polynomial, whose graph runs through the given points. 666 */ 667 lagrangePolynomial: function(p) { 668 var w = [], 669 fct = function(x, suspendedUpdate) { 670 var i, k, len, xi, s, 671 num = 0, denom = 0, 672 M, j; 673 674 len = p.length; 675 if (!suspendedUpdate) { 676 for (i = 0; i < len; i++) { 677 w[i] = 1.0; 678 xi = p[i].X(); 679 for (k = 0; k < len; k++) if (k != i) { 680 w[i] *= (xi - p[k].X()); 681 } 682 w[i] = 1 / w[i]; 683 } 684 685 M = []; 686 for (j = 0; j < len; j++) { 687 M.push([1]); 688 } 689 } 690 691 for (i = 0; i < len; i++) { 692 xi = p[i].X(); 693 if (x === xi) { 694 return p[i].Y(); 695 } else { 696 s = w[i] / (x - xi); 697 denom += s; 698 num += s * p[i].Y(); 699 } 700 } 701 return num / denom; 702 }; 703 704 fct.getTerm = function() { 705 return ''; 706 }; 707 708 return fct; 709 }, 710 711 /** 712 * Computes the cubic Catmull-Rom spline curve through a given set of points. The curve 713 * is uniformly parametrized. 714 * Two artificial control points at the beginning and the end are added. 715 * @param {Array} points Array consisting of JXG.Points. 716 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 717 * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, and a function simply 718 * returning the length of the points array minus three. 719 */ 720 CatmullRomSpline: function(points) { 721 var coeffs = [], 722 p, first = {}, last = {}, // control point at the beginning and at the end 723 724 makeFct = function(which) { 725 return function(t, suspendedUpdate) { 726 var len = points.length, 727 s, c; 728 729 if (len < 2 ) { return NaN; } 730 if (!suspendedUpdate) { 731 first[which] = function() {return 2*points[0][which]()-points[1][which]();} 732 last[which] = function() {return 2*points[len-1][which]()-points[len-2][which]();} 733 p = [first].concat(points, [last]); 734 coeffs[which] = []; 735 for (s=0; s < len-1; s++) { 736 coeffs[which][s] = [ 737 2*p[s+1][which](), 738 -p[s][which]() + p[s+2][which](), 739 2*p[s][which]() - 5*p[s+1][which]() + 4*p[s+2][which]() - p[s+3][which](), 740 -p[s][which]() + 3*p[s+1][which]() - 3*p[s+2][which]() + p[s+3][which]() 741 ]; 742 } 743 } 744 len += 2; // add the two control points 745 if (isNaN(t)) { return NaN; } 746 // This is necessay for our advanced plotting algorithm: 747 if (t < 0) { 748 return p[1][which](); 749 } else if (t >= len - 3) { 750 return p[len-2][which](); 751 } 752 753 s = Math.floor(t); 754 if (s==t) { 755 return p[s][which](); 756 } 757 t -= s; 758 c = coeffs[which][s]; 759 return 0.5*(((c[3]*t + c[2])*t + c[1])*t + c[0]); 760 }; 761 }; 762 763 return [makeFct('X'), 764 makeFct('Y'), 765 0, 766 function() { 767 return points.length - 1; 768 } 769 ]; 770 }, 771 772 773 /** 774 * Computes the regression polynomial of a given degree through a given set of coordinates. 775 * Returns the regression polynomial function. 776 * @param {Number} degree number, function or slider. 777 * Either 778 * @param {Array} dataX array containing the x-coordinates of the data set 779 * @param {Array} dataY array containing the y-coordinates of the data set, 780 * or 781 * @param {Array} dataX Array consisting of JXG.Points. 782 * @param {} dataY Ignored 783 * @returns {function} A function of one parameter which returns the value of the regression polynomial of the given degree. 784 * It possesses the method getTerm() which returns the string containing the function term of the polynomial. 785 */ 786 regressionPolynomial: function(degree, dataX, dataY) { 787 var coeffs, 788 deg, dX, dY, 789 inputType, 790 fct, 791 term = ''; 792 793 if (JXG.isPoint(degree) && typeof degree.Value == 'function') { // Slider 794 deg = function() {return degree.Value();}; 795 } else if (JXG.isFunction(degree)) { 796 deg = degree; 797 } else if (JXG.isNumber(degree)) { 798 deg = function() {return degree;}; 799 } else { 800 throw new Error("JSXGraph: Can't create regressionPolynomial from degree of type'" + (typeof degree) + "'."); 801 } 802 803 if (arguments.length == 3 && JXG.isArray(dataX) && JXG.isArray(dataY)) { // Parameters degree, dataX, dataY 804 inputType = 0; 805 } else if (arguments.length == 2 && JXG.isArray(dataX) && dataX.length > 0 && JXG.isPoint(dataX[0])) { // Parameters degree, point array 806 inputType = 1; 807 } else { 808 throw new Error("JSXGraph: Can't create regressionPolynomial. Wrong parameters."); 809 } 810 811 fct = function(x, suspendedUpdate) { 812 var i, j, M, MT, y, B, c, s, 813 d, len = dataX.length; // input data 814 815 d = Math.floor(deg()); // input data 816 if (!suspendedUpdate) { 817 if (inputType === 1) { // point list as input 818 dX = []; 819 dY = []; 820 for (i = 0; i < len; i++) { 821 dX[i] = dataX[i].X(); 822 dY[i] = dataX[i].Y(); 823 } 824 } 825 826 if (inputType === 0) { // check for functions 827 dX = []; 828 dY = []; 829 for (i = 0; i < len; i++) { 830 if (JXG.isFunction(dataX[i])) 831 dX.push(dataX[i]()); 832 else 833 dX.push(dataX[i]); 834 if (JXG.isFunction(dataY[i])) 835 dY.push(dataY[i]()); 836 else 837 dY.push(dataY[i]); 838 } 839 } 840 841 M = []; 842 for (j = 0; j < len; j++) { 843 M.push([1]); 844 } 845 for (i = 1; i <= d; i++) { 846 for (j = 0; j < len; j++) { 847 M[j][i] = M[j][i - 1] * dX[j]; // input data 848 } 849 } 850 851 y = dY; // input data 852 MT = JXG.Math.transpose(M); 853 B = JXG.Math.matMatMult(MT, M); 854 c = JXG.Math.matVecMult(MT, y); 855 coeffs = JXG.Math.Numerics.Gauss(B, c); 856 term = JXG.Math.Numerics.generatePolynomialTerm(coeffs, d, 'x', 3); 857 } 858 859 // Horner's scheme to evaluate polynomial 860 s = coeffs[d]; 861 for (i = d - 1; i >= 0; i--) { 862 s = (s * x + coeffs[i]); 863 } 864 return s; 865 }; 866 867 fct.getTerm = function() { 868 return term; 869 }; 870 871 return fct; 872 }, 873 874 /** 875 * Computes the cubic Bezier curve through a given set of points. 876 * @param {Array} points Array consisting of 3*k+1 JXG.Points. 877 * The points at position k with k mod 3 = 0 are the data points, 878 * points at position k with k mod 3 = 1 or 2 are the control points. 879 * @returns {Array} An array consisting of two functions of one parameter t which return the 880 * x resp. y coordinates of the Bezier curve in t, one zero value, and a third function accepting 881 * no parameters and returning one third of the length of the points. 882 */ 883 bezier: function(points) { 884 var len, 885 makeFct = function(which) { 886 return function(t, suspendedUpdate) { 887 var z = Math.floor(t) * 3, 888 t0 = t % 1, 889 t1 = 1 - t0; 890 891 if (!suspendedUpdate) { 892 len = Math.floor(points.length / 3); 893 } 894 895 if (t < 0) { 896 return points[0][which](); 897 } 898 if (t >= len) { 899 return points[points.length - 1][which](); 900 } 901 if (isNaN(t)) { 902 return NaN; 903 } 904 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; 905 }; 906 }; 907 908 return [ 909 makeFct('X'), 910 makeFct('Y'), 911 0, 912 function() { 913 return Math.floor(points.length / 3); 914 }]; 915 }, 916 917 /** 918 * Computes the B-spline curve of order k (order = degree+1) through a given set of points. 919 * @param {Array} points Array consisting of JXG.Points. 920 * @param {Number} order Order of the B-spline curve. 921 * @todo closed B-spline curves 922 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 923 * which return the x resp. y coordinates of the B-spline curve in t, a zero value, and a function simply 924 * returning the length of the points array minus one. 925 */ 926 bspline: function(points, order) { 927 var knots, N = [], 928 _knotVector = function(n, k) { 929 var j, kn = []; 930 for (j = 0; j < n + k + 1; j++) { 931 if (j < k) { 932 kn[j] = 0.0; 933 } else if (j <= n) { 934 kn[j] = j - k + 1; 935 } else { 936 kn[j] = n - k + 2; 937 } 938 } 939 return kn; 940 }, 941 942 _evalBasisFuncs = function(t, kn, n, k, s) { 943 var i,j,a,b,den, 944 N = []; 945 946 if (kn[s] <= t && t < kn[s + 1]) { 947 N[s] = 1; 948 } else { 949 N[s] = 0; 950 } 951 for (i = 2; i <= k; i++) { 952 for (j = s - i + 1; j <= s; j++) { 953 if (j <= s - i + 1 || j < 0) { 954 a = 0.0; 955 } else { 956 a = N[j]; 957 } 958 if (j >= s) { 959 b = 0.0; 960 } else { 961 b = N[j + 1]; 962 } 963 den = kn[j + i - 1] - kn[j]; 964 if (den == 0) { 965 N[j] = 0; 966 } else { 967 N[j] = (t - kn[j]) / den * a; 968 } 969 den = kn[j + i] - kn[j + 1]; 970 if (den != 0) { 971 N[j] += (kn[j + i] - t) / den * b; 972 } 973 } 974 } 975 return N; 976 }, 977 makeFct = function(which) { 978 return function(t, suspendedUpdate) { 979 var len = points.length, 980 y, j, s, 981 n = len - 1, 982 k = order; 983 984 if (n <= 0) return NaN; 985 if (n + 2 <= k) k = n + 1; 986 if (t <= 0) return points[0][which](); 987 if (t >= n - k + 2) return points[n][which](); 988 989 knots = _knotVector(n, k); 990 s = Math.floor(t) + k - 1; 991 N = _evalBasisFuncs(t, knots, n, k, s); 992 993 y = 0.0; 994 for (j = s - k + 1; j <= s; j++) { 995 if (j < len && j >= 0) y += points[j][which]() * N[j]; 996 } 997 return y; 998 }; 999 }; 1000 1001 return [makeFct('X'), 1002 makeFct('Y'), 1003 0, 1004 function() { 1005 return points.length - 1; 1006 } 1007 ]; 1008 }, 1009 1010 /** 1011 * Numerical (symmetric) approximation of derivative. suspendUpdate is piped through, see {@link JXG.Curve#updateCurve} 1012 * and {@link JXG.Curve#hasPoint}. 1013 * @param {function} f Function in one variable to be differentiated. 1014 * @param {object} [obj] Optional object that is treated as "this" in the function body. This is useful, if the function is a 1015 * method of an object and contains a reference to its parent object via "this". 1016 * @returns {function} Derivative function of a given function f. 1017 */ 1018 D: function(f, obj) { 1019 var h = 0.00001, 1020 h2 = 1.0 / (h * 2.0); 1021 1022 if (arguments.length == 1 || (arguments.length > 1 && !JXG.exists(arguments[1]))) { 1023 return function(x, suspendUpdate) { 1024 return (f(x + h, suspendUpdate) - f(x - h, suspendUpdate)) * h2; 1025 }; 1026 } else { // set "this" to "obj" in f 1027 return function(x, suspendUpdate) { 1028 return (f.apply(obj, [x + h,suspendUpdate]) - f.apply(obj, [x - h,suspendUpdate])) * h2; 1029 }; 1030 } 1031 }, 1032 1033 /** 1034 * Helper function to create curve which displays Riemann sums. 1035 * Compute coordinates for the rectangles showing the Riemann sum. 1036 * @param {function} f Function, whose integral is approximated by the Riemann sum. 1037 * @param {Number} n number of rectangles. 1038 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', or 'trapezodial'. 1039 * @param {Number} start Left border of the approximation interval 1040 * @param {Number} end Right border of the approximation interval 1041 * @returns {Array} An array of two arrays containing the x and y coordinates for the rectangles showing the Riemann sum. This 1042 * array may be used as parent array of a {@link JXG.Curve}. The third parameteris the riemann sum, i.e. the sum of the volumes of all 1043 * rectangles. 1044 */ 1045 riemann: function(f, n, type, start, end) { 1046 var xarr = [], 1047 yarr = [], 1048 i, j = 0, 1049 delta, 1050 x = start, y, 1051 x1, y1, delta1, 1052 sum = 0; 1053 1054 n = Math.round(n); 1055 1056 xarr[j] = x; 1057 yarr[j] = 0.0; 1058 1059 if (n > 0) { 1060 delta = (end - start) / n; 1061 delta1 = delta * 0.01; // for 'lower' and 'upper' 1062 1063 for (i = 0; i < n; i++) { 1064 if (type === 'right') { 1065 y = f(x + delta); 1066 } else if (type === 'middle') { 1067 y = f(x + delta * 0.5); 1068 } else if ((type === 'left') || (type === 'trapezodial')) { 1069 y = f(x); 1070 } else if (type === 'lower') { 1071 y = f(x); 1072 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1073 y1 = f(x1); 1074 if (y1 < y) { 1075 y = y1; 1076 } 1077 } 1078 } else if (type === 'upper') { 1079 y = f(x); 1080 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1081 y1 = f(x1); 1082 if (y1 > y) { 1083 y = y1; 1084 } 1085 } 1086 } else { // if (type === 'random') 1087 y = f(x + delta * Math.random()); 1088 } 1089 1090 j++; 1091 xarr[j] = x; 1092 yarr[j] = y; 1093 j++; 1094 x += delta; 1095 if (type === 'trapezodial') { 1096 y = f(x); 1097 } 1098 xarr[j] = x; 1099 yarr[j] = y; 1100 j++; 1101 xarr[j] = x; 1102 yarr[j] = 0.0; 1103 sum += y*delta; 1104 } 1105 } 1106 return [xarr,yarr,sum]; 1107 }, 1108 1109 /** 1110 * Approximate the integral by Riemann sums. 1111 * Compute the area described by the riemann sum rectangles. 1112 * @deprecated Replaced by JXG.Curve.Value(), @see JXG.Curve#riemannsum 1113 * @param {function} f Function, whose integral is approximated by the Riemann sum. 1114 * @param {Number} n number of rectangles. 1115 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random'. or 'trapezodial'. 1116 * @param {Number} start Left border of the approximation interval 1117 * @param {Number} end Right border of the approximation interval 1118 * @returns {Number} The sum of the areas of the rectangles. 1119 */ 1120 riemannsum: function(f, n, type, start, end) { 1121 var sum = .0, 1122 i, delta, 1123 x = start, y, 1124 x1, y1, delta1; 1125 1126 n = Math.floor(n); 1127 if (n > 0) { 1128 delta = (end - start) / n; 1129 delta1 = delta * 0.01; // for 'lower' and 'upper' 1130 for (i = 0; i < n; i++) { 1131 if (type === 'right') { 1132 y = f(x + delta); 1133 } else if (type === 'middle') { 1134 y = f(x + delta * 0.5); 1135 } else if (type === 'trapezodial') { 1136 y = 0.5 * (f(x + delta) + f(x)); 1137 } else if (type === 'left') { 1138 y = f(x); 1139 } else if (type === 'lower') { 1140 y = f(x); 1141 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1142 y1 = f(x1); 1143 if (y1 < y) { 1144 y = y1; 1145 } 1146 } 1147 } else if (type === 'upper') { 1148 y = f(x); 1149 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1150 y1 = f(x1); 1151 if (y1 > y) { 1152 y = y1; 1153 } 1154 } 1155 } else { // if (type === 'random') 1156 y = f(x + delta * Math.random()); 1157 } 1158 sum += delta * y; 1159 x += delta; 1160 } 1161 } 1162 return sum; 1163 }, 1164 1165 /** 1166 * Solve initial value problems numerically using Runge-Kutta-methods. 1167 * See {@link http://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm. 1168 * @param {object,String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing 1169 * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure 1170 * <pre> 1171 * { 1172 * s: <Number>, 1173 * A: <matrix>, 1174 * b: <Array>, 1175 * c: <Array> 1176 * } 1177 * </pre> 1178 * 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 1179 * @param {Array} x0 Initial value vector. If the problem is of one-dimensional, the initial value also has to be given in an array. 1180 * @param {Array} I Interval on which to integrate. 1181 * @param {Number} N Number of evaluation points. 1182 * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode 1183 * 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 1184 * vector <tt>x</tt>, and has to return a vector of the same dimension as <tt>x</tt> has. 1185 * @returns {Array} An array of vectors describing the solution of the ode on the given interval I. 1186 * @example 1187 * // A very simple autonomous system dx(t)/dt = x(t); 1188 * function f(t, x) { 1189 * return x; 1190 * } 1191 * 1192 * // Solve it with initial value x(0) = 1 on the interval [0, 2] 1193 * // with 20 evaluation points. 1194 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 1195 * 1196 * // Prepare data for plotting the solution of the ode using a curve. 1197 * var dataX = []; 1198 * var dataY = []; 1199 * var h = 0.1; // (I[1] - I[0])/N = (2-0)/20 1200 * for(var i=0; i<data.length; i++) { 1201 * dataX[i] = i*h; 1202 * dataY[i] = data[i][0]; 1203 * } 1204 * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'}); 1205 * </pre><div id="d2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div> 1206 * <script type="text/javascript"> 1207 * var board = JXG.JSXGraph.initBoard('d2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false}); 1208 * function f(t, x) { 1209 * // we have to copy the value. 1210 * // return x; would just return the reference. 1211 * return [x[0]]; 1212 * } 1213 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 1214 * var dataX = []; 1215 * var dataY = []; 1216 * var h = 0.1; 1217 * for(var i=0; i<data.length; i++) { 1218 * dataX[i] = i*h; 1219 * dataY[i] = data[i][0]; 1220 * } 1221 * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'}); 1222 * </script><pre> 1223 */ 1224 rungeKutta: function(butcher, x0, I, N, f) { 1225 var x = [], 1226 y = [], 1227 h = (I[1] - I[0]) / N, 1228 t = I[0], 1229 e, i, j, 1230 k, l, 1231 dim = x0.length, 1232 s, 1233 result = [], 1234 r = 0; 1235 1236 if (JXG.isString(butcher)) { 1237 butcher = predefinedButcher[butcher] || predefinedButcher.euler; 1238 } 1239 s = butcher.s; 1240 1241 // don't change x0, so copy it 1242 for (e = 0; e < dim; e++) 1243 x[e] = x0[e]; 1244 1245 for (i = 0; i < N; i++) { 1246 // Optimization doesn't work for ODEs plotted using time 1247 // if((i % quotient == 0) || (i == N-1)) { 1248 result[r] = []; 1249 for (e = 0; e < dim; e++) 1250 result[r][e] = x[e]; 1251 r++; 1252 // } 1253 // init k 1254 k = []; 1255 for (j = 0; j < s; j++) { 1256 // init y = 0 1257 for (e = 0; e < dim; e++) 1258 y[e] = 0.; 1259 1260 // Calculate linear combination of former k's and save it in y 1261 for (l = 0; l < j; l++) { 1262 for (e = 0; e < dim; e++) { 1263 y[e] += (butcher.A[j][l]) * h * k[l][e]; 1264 } 1265 } 1266 1267 // add x(t) to y 1268 for (e = 0; e < dim; e++) { 1269 y[e] += x[e]; 1270 } 1271 1272 // calculate new k and add it to the k matrix 1273 k.push(f(t + butcher.c[j] * h, y)); 1274 } 1275 1276 // init y = 0 1277 for (e = 0; e < dim; e++) 1278 y[e] = 0.; 1279 1280 for (l = 0; l < s; l++) { 1281 for (e = 0; e < dim; e++) 1282 y[e] += butcher.b[l] * k[l][e]; 1283 } 1284 1285 for (e = 0; e < dim; e++) { 1286 x[e] = x[e] + h * y[e]; 1287 } 1288 1289 t += h; 1290 } 1291 1292 return result; 1293 }, 1294 1295 1296 /* 1297 * Maximum number of iterations in @see #fzero 1298 */ 1299 maxIterationsRoot: 80, 1300 1301 /* 1302 * Maximum number of iterations in @see #fminbr 1303 */ 1304 maxIterationsMinimize: 500, 1305 1306 /** 1307 * 1308 * Find zero of an univariate function f. 1309 * @param {function} f Function, whose root is to be found 1310 * @param {array or number} x0 Start value or start interval enclosing the root 1311 * @param {object} object Parent object in case f is method of it 1312 * @return {number} the approximation of the root 1313 * Algorithm: 1314 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 1315 * computations. M., Mir, 1980, p.180 of the Russian edition 1316 * 1317 * if x0 is an array containing lower and upper bound for the zero 1318 * algorithm 748 is applied. Otherwise, if x0 is a number, 1319 * the algorithm tries to bracket a zero of f starting from x0. 1320 * If this fails, we fall back to Newton's method. 1321 * 1322 **/ 1323 fzero: function(f, x0, object) { 1324 var tol = JXG.Math.eps, 1325 maxiter = this.maxIterationsRoot, niter = 0, 1326 nfev = 0, 1327 eps = tol, 1328 a,b,c, 1329 fa, fb, fc, 1330 aa, blist, i, len, u, fu, 1331 prev_step, 1332 tol_act, // Actual tolerance 1333 p, // Interpolation step is calcu- 1334 q, // lated in the form p/q; divi- 1335 // sion operations is delayed 1336 // until the last moment 1337 new_step, // Step at this iteration 1338 t1, cb, t2; 1339 1340 if (JXG.isArray(x0)) { 1341 if (x0.length<2) 1342 throw new Error("JXG.Math.Numerics.fzero: length of array x0 has to be at least two."); 1343 a = x0[0]; fa = f.apply(object,[a]); nfev++; 1344 b = x0[1]; fb = f.apply(object,[b]); nfev++; 1345 } else { 1346 a = x0; fa = f.apply(object,[a]); nfev++; 1347 // Try to get b. 1348 if (a == 0) { 1349 aa = 1; 1350 } else { 1351 aa = a; 1352 } 1353 blist = [0.9*aa, 1.1*aa, aa-1, aa+1, 0.5*aa, 1.5*aa, -aa, 2*aa, -10*aa, 10*aa]; 1354 len = blist.length; 1355 for (i=0;i<len;i++) { 1356 b = blist[i]; 1357 fb = f.apply(object,[b]); nfev++; 1358 if (fa*fb<=0) { 1359 break; 1360 } 1361 } 1362 if (b < a) { 1363 u = a; a = b; b = u; 1364 fu = fa; fa = fb; fb = fu; 1365 } 1366 } 1367 1368 if (fa*fb > 0) { 1369 // Bracketing not successful, fall back to Newton's method or to fminbr 1370 if (JXG.isArray(x0)) { 1371 //JXG.debug("fzero falls back to fminbr"); 1372 return this.fminbr(f, [a,b], object); 1373 } else { 1374 //JXG.debug("fzero falls back to Newton"); 1375 return this.Newton(f, a, object); 1376 } 1377 } 1378 1379 // OK, we have enclosed a zero of f. 1380 // Now we can start Brent's method 1381 1382 c = a; 1383 fc = fa; 1384 while (niter<maxiter) { // Main iteration loop 1385 prev_step = b-a; // Distance from the last but one 1386 // to the last approximation 1387 1388 if ( Math.abs(fc) < Math.abs(fb) ) { 1389 // Swap data for b to be the 1390 a = b; b = c; c = a; // best approximation 1391 fa=fb; fb=fc; fc=fa; 1392 } 1393 tol_act = 2*eps*Math.abs(b) + tol*0.5; 1394 new_step = (c-b)*0.5; 1395 1396 if ( Math.abs(new_step) <= tol_act && Math.abs(fb) <= eps ) { 1397 //JXG.debug("nfev="+nfev); 1398 return b; // Acceptable approx. is found 1399 } 1400 1401 // Decide if the interpolation can be tried 1402 if ( Math.abs(prev_step) >= tol_act // If prev_step was large enough 1403 && Math.abs(fa) > Math.abs(fb) ) { // and was in true direction, 1404 // Interpolatiom may be tried 1405 cb = c-b; 1406 if ( a==c ) { // If we have only two distinct 1407 // points linear interpolation 1408 t1 = fb/fa; // can only be applied 1409 p = cb*t1; 1410 q = 1.0 - t1; 1411 } else { // Quadric inverse interpolation 1412 q = fa/fc; t1 = fb/fc; t2 = fb/fa; 1413 p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) ); 1414 q = (q-1.0) * (t1-1.0) * (t2-1.0); 1415 } 1416 if ( p>0 ) { // p was calculated with the op- 1417 q = -q; // posite sign; make p positive 1418 } else { // and assign possible minus to 1419 p = -p; // q 1420 } 1421 1422 if( p < (0.75*cb*q-Math.abs(tol_act*q)*0.5) // If b+p/q falls in [b,c] 1423 && p < Math.abs(prev_step*q*0.5) ) { // and isn't too large 1424 new_step = p/q; // it is accepted 1425 } // If p/q is too large then the 1426 // bissection procedure can 1427 // reduce [b,c] range to more 1428 // extent 1429 } 1430 1431 if ( Math.abs(new_step) < tol_act ) { // Adjust the step to be not less 1432 if ( new_step > 0 ) { // than tolerance 1433 new_step = tol_act; 1434 } else { 1435 new_step = -tol_act; 1436 } 1437 } 1438 1439 a = b; fa = fb; // Save the previous approx. 1440 b += new_step; 1441 fb = f.apply(object,[b]); nfev++; // Do step to a new approxim. 1442 if ( (fb>0 && fc>0) || (fb<0 && fc<0) ) { 1443 // Adjust c for it to have a sign 1444 c = a; fc = fa; // opposite to that of b 1445 } 1446 niter++; 1447 } // End while 1448 1449 //JXG.debug("fzero: maxiter="+maxiter+" reached."); 1450 return b; 1451 }, 1452 1453 1454 /** 1455 * 1456 * Find minimum of an univariate function f. 1457 * @param {function} f Function, whose minimum is to be found 1458 * @param {array} x0 Start interval enclosing the minimum 1459 * @param {object} object Parent object in case f is method of it 1460 * @return {number} the approximation of the minimum 1461 * Algorithm: 1462 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 1463 * computations. M., Mir, 1980, p.180 of the Russian edition 1464 * x0 1465 **/ 1466 fminbr: function(f, x0, object) { // An estimate to the min location 1467 var a, b, x, v, w, 1468 fx, fv, fw, 1469 r = (3.-Math.sqrt(5.0))*0.5, // Golden section ratio 1470 tol = JXG.Math.eps, 1471 sqrteps = Math.sqrt(JXG.Math.eps), 1472 maxiter = this.maxIterationsMinimize, 1473 niter = 0, 1474 range, middle_range, tol_act, new_step, 1475 p, q, t, ft, 1476 nfev = 0; 1477 1478 if (!JXG.isArray(x0) || x0.length<2) { 1479 throw new Error("JXG.Math.Numerics.fminbr: length of array x0 has to be at least two."); 1480 } 1481 a = x0[0]; 1482 b = x0[1]; 1483 v = a + r*(b-a); 1484 fv = f.apply(object,[v]); nfev++; // First step - always gold section 1485 x = v; w = v; 1486 fx=fv; fw=fv; 1487 1488 while (niter<maxiter) { 1489 range = b-a; // Range over which the minimum 1490 // is seeked for 1491 middle_range = (a+b)*0.5; 1492 tol_act = sqrteps*Math.abs(x) + tol/3; // Actual tolerance 1493 if( Math.abs(x-middle_range) + range*0.5 <= 2*tol_act ) { 1494 //JXG.debug(nfev); 1495 return x; // Acceptable approx. is found 1496 } 1497 // Obtain the golden section step 1498 new_step = r * ( x<middle_range ? b-x : a-x ); 1499 // Decide if the interpolation can be tried 1500 if ( Math.abs(x-w) >= tol_act ) { // If x and w are distinct 1501 // interpolatiom may be tried 1502 // Interpolation step is calculated as p/q; 1503 // division operation is delayed until last moment 1504 t = (x-w) * (fx-fv); 1505 q = (x-v) * (fx-fw); 1506 p = (x-v)*q - (x-w)*t; 1507 q = 2*(q-t); 1508 1509 if ( q>0 ) { // q was calculated with the op- 1510 p = -p; // posite sign; make q positive 1511 } else { // and assign possible minus to 1512 q = -q; // p 1513 } 1514 if ( Math.abs(p) < Math.abs(new_step*q) && // If x+p/q falls in [a,b] 1515 p > q*(a-x+2*tol_act) && // not too close to a and 1516 p < q*(b-x-2*tol_act) ) { // b, and isn't too large */ 1517 new_step = p/q; // it is accepted 1518 } 1519 // If p/q is too large then the 1520 // golden section procedure can 1521 // reduce [a,b] range to more 1522 // extent 1523 } 1524 1525 if ( Math.abs(new_step) < tol_act ) { // Adjust the step to be not less 1526 if( new_step > 0 ) { // than tolerance 1527 new_step = tol_act; 1528 } else { 1529 new_step = -tol_act; 1530 } 1531 } 1532 1533 // Obtain the next approximation to min 1534 // and reduce the enveloping range 1535 t = x + new_step; // Tentative point for the min 1536 ft = f.apply(object,[t]); nfev++; 1537 if ( ft <= fx ) { // t is a better approximation 1538 if ( t < x ) { // Reduce the range so that 1539 b = x; // t would fall within it 1540 } else { 1541 a = x; 1542 } 1543 v = w; w = x; x = t; // Assign the best approx to x 1544 fv=fw; fw=fx; fx=ft; 1545 } else { // x remains the better approx 1546 if ( t < x ) { // Reduce the range enclosing x 1547 a = t; 1548 } else { 1549 b = t; 1550 } 1551 if ( ft <= fw || w==x ) { 1552 v = w; w = t; 1553 fv=fw; fw=ft; 1554 } else if ( ft<=fv || v==x || v==w ) { 1555 v = t; 1556 fv=ft; 1557 } 1558 } 1559 niter++; 1560 } 1561 //JXG.debug("fminbr: maxiter="+maxiter+" reached."); 1562 return x; 1563 }, 1564 1565 /** 1566 * Helper function to create curve which displays Reuleaux polygons. 1567 * @param {array} points Array of points which should be the vertices of the Reuleaux polygon. Typically, 1568 * these point list is the array vrtices of a regular polygon. 1569 * @param {number} nr Number of vertices 1570 * @returns {array} An array containing the two functions defining the Reuleaux polygon and the two values 1571 * for the start and the end of the paramtric curve. 1572 * array may be used as parent array of a {@link JXG.Curve}. 1573 * 1574 * @example 1575 * var A = brd.create('point',[-2,-2]); 1576 * var B = brd.create('point',[0,1]); 1577 * var pol = brd.create('regularpolygon',[A,B,3], {withLines:false, fillColor:'none', highlightFillColor:'none', fillOpacity:0.0}); 1578 * var reuleauxTriangle = brd.create('curve', JXG.Math.Numerics.reuleauxPolygon(pol.vertices, 3), 1579 * {strokeWidth:6, strokeColor:'#d66d55', fillColor:'#ad5544', highlightFillColor:'#ad5544'}); 1580 * 1581 * </pre><div id="2543a843-46a9-4372-abc1-94d9ad2db7ac" style="width: 300px; height: 300px;"></div> 1582 * <script type="text/javascript"> 1583 * var brd = JXG.JSXGraph.initBoard('2543a843-46a9-4372-abc1-94d9ad2db7ac', {boundingbox: [-5, 5, 5, -5], axis: true, showcopyright:false, shownavigation: false}); 1584 * var A = brd.create('point',[-2,-2]); 1585 * var B = brd.create('point',[0,1]); 1586 * var pol = brd.create('regularpolygon',[A,B,3], {withLines:false, fillColor:'none', highlightFillColor:'none', fillOpacity:0.0}); 1587 * var reuleauxTriangle = brd.create('curve', JXG.Math.Numerics.reuleauxPolygon(pol.vertices, 3), 1588 * {strokeWidth:6, strokeColor:'#d66d55', fillColor:'#ad5544', highlightFillColor:'#ad5544'}); 1589 * </script><pre> 1590 */ 1591 reuleauxPolygon: function(points, nr) { 1592 var pi2 = Math.PI*2, 1593 pi2_n = pi2/nr, 1594 diag = (nr-1)/2, 1595 beta, d = 0, 1596 makeFct = function(which, trig) { 1597 return function(t, suspendUpdate) { 1598 if (!suspendUpdate) { 1599 d = points[0].Dist(points[diag]); 1600 beta = JXG.Math.Geometry.rad([points[0].X()+1,points[0].Y()],points[0],points[(diag)%nr]); 1601 } 1602 var t1 = (t%pi2 + pi2) % pi2; 1603 var j = Math.floor(t1 / pi2_n)%nr; 1604 if (isNaN(j)) return j; 1605 //t1 = (t1-j*pi2_n)*0.5 + beta+j*pi2_n; 1606 t1 = t1*0.5+j*pi2_n*0.5 + beta; 1607 return points[j][which]()+d*Math[trig](t1); 1608 }; 1609 }; 1610 return [ 1611 makeFct('X','cos'), 1612 makeFct('Y','sin'), 1613 0, 1614 Math.PI*2 1615 ]; 1616 }, 1617 1618 /** 1619 * Implements the Ramer-Douglas-Peuker algorithm. 1620 * It discards points which are not necessary from the polygonal line defined by the point array 1621 * pts. The computation is done in screen coordinates. 1622 * Average runtime is O(nlog(n)), worst case runtime is O(n^2), where n is the number of points. 1623 * @param {Array} pts Array of {@link JXG.Coords} 1624 * @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>. 1625 * @returns {Array} An array containing points which represent an apparently identical curve as the points of pts do, but contains fewer points. 1626 */ 1627 RamerDouglasPeuker: function(pts, eps) { 1628 var newPts = [], i, k, len, 1629 /** 1630 * RDP() is a private subroutine of {@link JXG.Math.Numerics#RamerDouglasPeuker}. 1631 * It runs recursively through the point set and searches the 1632 * point which has the largest distance from the line between the first point and 1633 * the last point. If the distance from the line is greater than eps, this point is 1634 * included in our new point set otherwise it is discarded. 1635 * If it is taken, we recursively apply the subroutine to the point set before 1636 * and after the chosen point. 1637 * @param {Array} pts Array of {@link JXG.Coords} 1638 * @param {Number} i Index of an element of pts 1639 * @param {Number} j Index of an element of pts 1640 * @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>. 1641 * @param {Array} newPts Array of {@link JXG.Coords} 1642 * @private 1643 */ 1644 RDP = function(pts, i, j, eps, newPts) { 1645 var result = findSplit(pts, i, j); 1646 1647 if (result[0] > eps) { 1648 RDP(pts, i, result[1], eps, newPts); 1649 RDP(pts, result[1], j, eps, newPts); 1650 } else { 1651 newPts.push(pts[j]); 1652 } 1653 }, 1654 /** 1655 * findSplit() is a subroutine of {@link JXG.Math.Numerics#RamerDouglasPeuker}. 1656 * It searches for the point between index i and j which 1657 * has the largest distance from the line between the points i and j. 1658 * @param {Array} pts Array of {@link JXG.Coords} 1659 * @param {Number} i Index of a point in pts 1660 * @param {Number} j Index of a point in pts 1661 **/ 1662 findSplit = function(pts, i, j) { 1663 var dist = 0, 1664 f = i, 1665 d, k, ci, cj, ck, 1666 x0, y0, x1, y1, 1667 den, lbda; 1668 1669 if (j - i < 2) return [-1.0,0]; 1670 1671 ci = pts[i].scrCoords; 1672 cj = pts[j].scrCoords; 1673 if (isNaN(ci[1] + ci[2] + cj[1] + cj[2])) return [NaN,j]; 1674 1675 for (k = i + 1; k < j; k++) { 1676 ck = pts[k].scrCoords; 1677 x0 = ck[1] - ci[1]; 1678 y0 = ck[2] - ci[2]; 1679 x1 = cj[1] - ci[1]; 1680 y1 = cj[2] - ci[2]; 1681 den = x1 * x1 + y1 * y1; 1682 /* 1683 if (den >= JXG.Math.eps) { 1684 lbda = (x0 * x1 + y0 * y1) / den; 1685 d = x0 * x0 + y0 * y0 - lbda * (x0 * x1 + y0 * y1); 1686 } else { 1687 lbda = 0.0; 1688 d = x0 * x0 + y0 * y0; 1689 } 1690 if (lbda < 0.0) { 1691 d = x0 * x0 + y0 * y0; 1692 } else if (lbda > 1.0) { 1693 x0 = ck[1] - cj[1]; 1694 y0 = ck[2] - cj[2]; 1695 d = x0 * x0 + y0 * y0; 1696 } 1697 */ 1698 if (den >= JXG.Math.eps) { 1699 lbda = (x0 * x1 + y0 * y1) / den; 1700 //d = x0 * x0 + y0 * y0 - lbda * (x0 * x1 + y0 * y1); 1701 if (lbda<0.0) { 1702 lbda = 0.0; 1703 } else if (lbda>1.0) { 1704 lbda = 1.0; 1705 } 1706 x0 = x0-lbda*x1; 1707 y0 = y0-lbda*y1; 1708 d = x0*x0+y0*y0; 1709 } else { 1710 lbda = 0.0; 1711 d = x0 * x0 + y0 * y0; 1712 } 1713 1714 if (d > dist) { 1715 dist = d; 1716 f = k; 1717 } 1718 } 1719 return [Math.sqrt(dist),f]; 1720 }; 1721 1722 len = pts.length; 1723 1724 // Search for the left most point woithout NaN coordinates 1725 i = 0; 1726 while (i < len && isNaN(pts[i].scrCoords[1] + pts[i].scrCoords[2])) { 1727 i++; 1728 } 1729 // Search for the right most point woithout NaN coordinates 1730 k = len - 1; 1731 while (k > i && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])) { 1732 k--; 1733 } 1734 1735 // Only proceed if something is left 1736 if (!(i > k || i == len)) { 1737 newPts[0] = pts[i]; 1738 RDP(pts, i, k, eps, newPts); 1739 } 1740 1741 return newPts; 1742 } 1743 1744 } 1745 })(JXG, Math); 1746 1747