1 /* 2 Copyright 2008-2025 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 dual licensed under the GNU LGPL or MIT License. 13 14 You can redistribute it and/or modify it under the terms of the 15 16 * GNU Lesser General Public License as published by 17 the Free Software Foundation, either version 3 of the License, or 18 (at your option) any later version 19 OR 20 * MIT License: https://github.com/jsxgraph/jsxgraph/blob/master/LICENSE.MIT 21 22 JSXGraph is distributed in the hope that it will be useful, 23 but WITHOUT ANY WARRANTY; without even the implied warranty of 24 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 25 GNU Lesser General Public License for more details. 26 27 You should have received a copy of the GNU Lesser General Public License and 28 the MIT License along with JSXGraph. If not, see <https://www.gnu.org/licenses/> 29 and <https://opensource.org/licenses/MIT/>. 30 */ 31 32 /*global JXG: true, define: true*/ 33 /*jslint nomen: true, plusplus: true*/ 34 /*eslint no-loss-of-precision: off */ 35 36 /** 37 * @fileoverview In this file the namespace Math.Numerics is defined, which holds numerical 38 * algorithms for solving linear equations etc. 39 */ 40 41 import JXG from "../jxg.js"; 42 import Type from "../utils/type.js"; 43 import Env from "../utils/env.js"; 44 import Mat from "./math.js"; 45 46 // Predefined butcher tableaus for the common Runge-Kutta method (fourth order), Heun method (second order), and Euler method (first order). 47 var predefinedButcher = { 48 rk4: { 49 s: 4, 50 A: [ 51 [0, 0, 0, 0], 52 [0.5, 0, 0, 0], 53 [0, 0.5, 0, 0], 54 [0, 0, 1, 0] 55 ], 56 b: [1.0 / 6.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 6.0], 57 c: [0, 0.5, 0.5, 1] 58 }, 59 heun: { 60 s: 2, 61 A: [ 62 [0, 0], 63 [1, 0] 64 ], 65 b: [0.5, 0.5], 66 c: [0, 1] 67 }, 68 euler: { 69 s: 1, 70 A: [[0]], 71 b: [1], 72 c: [0] 73 } 74 }; 75 76 /** 77 * The JXG.Math.Numerics namespace holds numerical algorithms, constants, and variables. 78 * @name JXG.Math.Numerics 79 * @exports Mat.Numerics as JXG.Math.Numerics 80 * @namespace 81 */ 82 Mat.Numerics = { 83 //JXG.extend(Mat.Numerics, /** @lends JXG.Math.Numerics */ { 84 /** 85 * Solves a system of linear equations given by A and b using the Gauss-Jordan-elimination. 86 * The algorithm runs in-place. I.e. the entries of A and b are changed. 87 * @param {Array} A Square matrix represented by an array of rows, containing the coefficients of the lineare equation system. 88 * @param {Array} b A vector containing the linear equation system's right hand side. 89 * @throws {Error} If a non-square-matrix is given or if b has not the right length or A's rank is not full. 90 * @returns {Array} A vector that solves the linear equation system. 91 * @memberof JXG.Math.Numerics 92 */ 93 Gauss: function (A, b) { 94 var i, 95 j, 96 k, 97 // copy the matrix to prevent changes in the original 98 Acopy, 99 // solution vector, to prevent changing b 100 x, 101 eps = Mat.eps, 102 // number of columns of A 103 n = A.length > 0 ? A[0].length : 0; 104 105 if (n !== b.length || n !== A.length) { 106 throw new Error( 107 "JXG.Math.Numerics.Gauss: Dimensions don't match. A must be a square matrix and b must be of the same length as A." 108 ); 109 } 110 111 // initialize solution vector 112 Acopy = []; 113 x = b.slice(0, n); 114 115 for (i = 0; i < n; i++) { 116 Acopy[i] = A[i].slice(0, n); 117 } 118 119 // Gauss-Jordan-elimination 120 for (j = 0; j < n; j++) { 121 for (i = n - 1; i > j; i--) { 122 // Is the element which is to eliminate greater than zero? 123 if (Math.abs(Acopy[i][j]) > eps) { 124 // Equals pivot element zero? 125 if (Math.abs(Acopy[j][j]) < eps) { 126 // At least numerically, so we have to exchange the rows 127 Type.swap(Acopy, i, j); 128 Type.swap(x, i, j); 129 } else { 130 // Saves the L matrix of the LR-decomposition. unnecessary. 131 Acopy[i][j] /= Acopy[j][j]; 132 // Transform right-hand-side b 133 x[i] -= Acopy[i][j] * x[j]; 134 135 // subtract the multiple of A[i][j] / A[j][j] of the j-th row from the i-th. 136 for (k = j + 1; k < n; k++) { 137 Acopy[i][k] -= Acopy[i][j] * Acopy[j][k]; 138 } 139 } 140 } 141 } 142 143 // The absolute values of all coefficients below the j-th row in the j-th column are smaller than JXG.Math.eps. 144 if (Math.abs(Acopy[j][j]) < eps) { 145 throw new Error( 146 "JXG.Math.Numerics.Gauss(): The given matrix seems to be singular." 147 ); 148 } 149 } 150 151 this.backwardSolve(Acopy, x, true); 152 153 return x; 154 }, 155 156 /** 157 * Solves a system of linear equations given by the right triangular matrix R and vector b. 158 * @param {Array} R Right triangular matrix represented by an array of rows. All entries a_(i,j) with i < j are ignored. 159 * @param {Array} b Right hand side of the linear equation system. 160 * @param {Boolean} [canModify=false] If true, the right hand side vector is allowed to be changed by this method. 161 * @returns {Array} An array representing a vector that solves the system of linear equations. 162 * @memberof JXG.Math.Numerics 163 */ 164 backwardSolve: function (R, b, canModify) { 165 var x, m, n, i, j; 166 167 if (canModify) { 168 x = b; 169 } else { 170 x = b.slice(0, b.length); 171 } 172 173 // m: number of rows of R 174 // n: number of columns of R 175 m = R.length; 176 n = R.length > 0 ? R[0].length : 0; 177 178 for (i = m - 1; i >= 0; i--) { 179 for (j = n - 1; j > i; j--) { 180 x[i] -= R[i][j] * x[j]; 181 } 182 x[i] /= R[i][i]; 183 } 184 185 return x; 186 }, 187 188 /** 189 * Gauss-Bareiss algorithm to compute the 190 * determinant of matrix without fractions. 191 * See Henri Cohen, "A Course in Computational 192 * Algebraic Number Theory (Graduate texts 193 * in mathematics; 138)", Springer-Verlag, 194 * ISBN 3-540-55640-0 / 0-387-55640-0 195 * Third, Corrected Printing 1996 196 * "Algorithm 2.2.6", pg. 52-53 197 * 198 * @param {Array} mat Matrix 199 * @returns Number 200 * @private 201 * @memberof JXG.Math.Numerics 202 */ 203 gaussBareiss: function (mat) { 204 var k, c, s, 205 i, j, p, 206 n, M, t, 207 eps = Mat.eps; 208 209 n = mat.length; 210 211 if (n <= 0) { 212 return 0; 213 } 214 215 if (mat[0].length < n) { 216 n = mat[0].length; 217 } 218 219 // Copy the input matrix to M 220 M = []; 221 222 for (i = 0; i < n; i++) { 223 M[i] = mat[i].slice(0, n); 224 } 225 226 c = 1; 227 s = 1; 228 229 for (k = 0; k < n - 1; k++) { 230 p = M[k][k]; 231 232 // Pivot step 233 if (Math.abs(p) < eps) { 234 for (i = k + 1; i < n; i++) { 235 if (Math.abs(M[i][k]) >= eps) { 236 break; 237 } 238 } 239 240 // No nonzero entry found in column k -> det(M) = 0 241 if (i === n) { 242 return 0.0; 243 } 244 245 // swap row i and k partially 246 for (j = k; j < n; j++) { 247 t = M[i][j]; 248 M[i][j] = M[k][j]; 249 M[k][j] = t; 250 } 251 s = -s; 252 p = M[k][k]; 253 } 254 255 // Main step 256 for (i = k + 1; i < n; i++) { 257 for (j = k + 1; j < n; j++) { 258 t = p * M[i][j] - M[i][k] * M[k][j]; 259 M[i][j] = t / c; 260 } 261 } 262 263 c = p; 264 } 265 266 return s * M[n - 1][n - 1]; 267 }, 268 269 /** 270 * Computes the determinant of a square nxn matrix with the 271 * Gauss-Bareiss algorithm. 272 * @param {Array} mat Matrix. 273 * @returns {Number} The determinant pf the matrix mat. 274 * The empty matrix returns 0. 275 * @memberof JXG.Math.Numerics 276 */ 277 det: function (mat) { 278 var n = mat.length; 279 280 if (n === 2 && mat[0].length === 2) { 281 return mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1]; 282 } 283 284 return this.gaussBareiss(mat); 285 }, 286 287 /** 288 * Compute the Eigenvalues and Eigenvectors of a symmetric 3x3 matrix with the Jacobi method 289 * Adaption of a FORTRAN program by Ed Wilson, Dec. 25, 1990 290 * @param {Array} Ain A symmetric 3x3 matrix. 291 * @returns {Array} [A,V] the matrices A and V. The diagonal of A contains the Eigenvalues, V contains the Eigenvectors. 292 * @memberof JXG.Math.Numerics 293 */ 294 Jacobi: function (Ain) { 295 var i, 296 j, 297 k, 298 aa, 299 si, 300 co, 301 tt, 302 ssum, 303 amax, 304 eps = Mat.eps * Mat.eps, 305 sum = 0.0, 306 n = Ain.length, 307 V = [ 308 [0, 0, 0], 309 [0, 0, 0], 310 [0, 0, 0] 311 ], 312 A = [ 313 [0, 0, 0], 314 [0, 0, 0], 315 [0, 0, 0] 316 ], 317 nloops = 0; 318 319 // Initialization. Set initial Eigenvectors. 320 for (i = 0; i < n; i++) { 321 for (j = 0; j < n; j++) { 322 V[i][j] = 0.0; 323 A[i][j] = Ain[i][j]; 324 sum += Math.abs(A[i][j]); 325 } 326 V[i][i] = 1.0; 327 } 328 329 // Trivial problems 330 if (n === 1) { 331 return [A, V]; 332 } 333 334 if (sum <= 0.0) { 335 return [A, V]; 336 } 337 338 sum /= n * n; 339 340 // Reduce matrix to diagonal 341 do { 342 ssum = 0.0; 343 amax = 0.0; 344 for (j = 1; j < n; j++) { 345 for (i = 0; i < j; i++) { 346 // Check if A[i][j] is to be reduced 347 aa = Math.abs(A[i][j]); 348 349 if (aa > amax) { 350 amax = aa; 351 } 352 353 ssum += aa; 354 355 if (aa >= eps) { 356 // calculate rotation angle 357 aa = Math.atan2(2.0 * A[i][j], A[i][i] - A[j][j]) * 0.5; 358 si = Math.sin(aa); 359 co = Math.cos(aa); 360 361 // Modify 'i' and 'j' columns 362 for (k = 0; k < n; k++) { 363 tt = A[k][i]; 364 A[k][i] = co * tt + si * A[k][j]; 365 A[k][j] = -si * tt + co * A[k][j]; 366 tt = V[k][i]; 367 V[k][i] = co * tt + si * V[k][j]; 368 V[k][j] = -si * tt + co * V[k][j]; 369 } 370 371 // Modify diagonal terms 372 A[i][i] = co * A[i][i] + si * A[j][i]; 373 A[j][j] = -si * A[i][j] + co * A[j][j]; 374 A[i][j] = 0.0; 375 376 // Make 'A' matrix symmetrical 377 for (k = 0; k < n; k++) { 378 A[i][k] = A[k][i]; 379 A[j][k] = A[k][j]; 380 } 381 // A[i][j] made zero by rotation 382 } 383 } 384 } 385 nloops += 1; 386 } while (Math.abs(ssum) / sum > eps && nloops < 2000); 387 388 return [A, V]; 389 }, 390 391 /** 392 * Calculates the integral of function f over interval using Newton-Cotes-algorithm. 393 * @param {Array} interval The integration interval, e.g. [0, 3]. 394 * @param {function} f A function which takes one argument of type number and returns a number. 395 * @param {Object} [config] The algorithm setup. Accepted properties are number_of_nodes of type number and integration_type 396 * with value being either 'trapez', 'simpson', or 'milne'. 397 * @param {Number} [config.number_of_nodes=28] 398 * @param {String} [config.integration_type='milne'] Possible values are 'milne', 'simpson', 'trapez' 399 * @returns {Number} Integral value of f over interval 400 * @throws {Error} If config.number_of_nodes doesn't match config.integration_type an exception is thrown. If you want to use 401 * simpson rule respectively milne rule config.number_of_nodes must be dividable by 2 respectively 4. 402 * @example 403 * function f(x) { 404 * return x*x; 405 * } 406 * 407 * // calculates integral of <tt>f</tt> from 0 to 2. 408 * var area1 = JXG.Math.Numerics.NewtonCotes([0, 2], f); 409 * 410 * // the same with an anonymous function 411 * var area2 = JXG.Math.Numerics.NewtonCotes([0, 2], function (x) { return x*x; }); 412 * 413 * // use trapez rule with 16 nodes 414 * var area3 = JXG.Math.Numerics.NewtonCotes([0, 2], f, 415 * {number_of_nodes: 16, integration_type: 'trapez'}); 416 * @memberof JXG.Math.Numerics 417 */ 418 NewtonCotes: function (interval, f, config) { 419 var evaluation_point, 420 i, 421 number_of_intervals, 422 integral_value = 0.0, 423 number_of_nodes = 424 config && Type.isNumber(config.number_of_nodes) ? config.number_of_nodes : 28, 425 available_types = { trapez: true, simpson: true, milne: true }, 426 integration_type = 427 config && 428 config.integration_type && 429 available_types.hasOwnProperty(config.integration_type) && 430 available_types[config.integration_type] 431 ? config.integration_type 432 : "milne", 433 step_size = (interval[1] - interval[0]) / number_of_nodes; 434 435 switch (integration_type) { 436 case "trapez": 437 integral_value = (f(interval[0]) + f(interval[1])) * 0.5; 438 evaluation_point = interval[0]; 439 440 for (i = 0; i < number_of_nodes - 1; i++) { 441 evaluation_point += step_size; 442 integral_value += f(evaluation_point); 443 } 444 445 integral_value *= step_size; 446 break; 447 case "simpson": 448 if (number_of_nodes % 2 > 0) { 449 throw new Error( 450 "JSXGraph: INT_SIMPSON requires config.number_of_nodes dividable by 2." 451 ); 452 } 453 454 number_of_intervals = number_of_nodes / 2.0; 455 integral_value = f(interval[0]) + f(interval[1]); 456 evaluation_point = interval[0]; 457 458 for (i = 0; i < number_of_intervals - 1; i++) { 459 evaluation_point += 2.0 * step_size; 460 integral_value += 2.0 * f(evaluation_point); 461 } 462 463 evaluation_point = interval[0] - step_size; 464 465 for (i = 0; i < number_of_intervals; i++) { 466 evaluation_point += 2.0 * step_size; 467 integral_value += 4.0 * f(evaluation_point); 468 } 469 470 integral_value *= step_size / 3.0; 471 break; 472 default: 473 if (number_of_nodes % 4 > 0) { 474 throw new Error( 475 "JSXGraph: Error in INT_MILNE: config.number_of_nodes must be a multiple of 4" 476 ); 477 } 478 479 number_of_intervals = number_of_nodes * 0.25; 480 integral_value = 7.0 * (f(interval[0]) + f(interval[1])); 481 evaluation_point = interval[0]; 482 483 for (i = 0; i < number_of_intervals - 1; i++) { 484 evaluation_point += 4.0 * step_size; 485 integral_value += 14.0 * f(evaluation_point); 486 } 487 488 evaluation_point = interval[0] - 3.0 * step_size; 489 490 for (i = 0; i < number_of_intervals; i++) { 491 evaluation_point += 4.0 * step_size; 492 integral_value += 493 32.0 * (f(evaluation_point) + f(evaluation_point + 2 * step_size)); 494 } 495 496 evaluation_point = interval[0] - 2.0 * step_size; 497 498 for (i = 0; i < number_of_intervals; i++) { 499 evaluation_point += 4.0 * step_size; 500 integral_value += 12.0 * f(evaluation_point); 501 } 502 503 integral_value *= (2.0 * step_size) / 45.0; 504 } 505 return integral_value; 506 }, 507 508 /** 509 * Calculates the integral of function f over interval using Romberg iteration. 510 * @param {Array} interval The integration interval, e.g. [0, 3]. 511 * @param {function} f A function which takes one argument of type number and returns a number. 512 * @param {Object} [config] The algorithm setup. Accepted properties are max_iterations of type number and precision eps. 513 * @param {Number} [config.max_iterations=20] 514 * @param {Number} [config.eps=0.0000001] 515 * @returns {Number} Integral value of f over interval 516 * @example 517 * function f(x) { 518 * return x*x; 519 * } 520 * 521 * // calculates integral of <tt>f</tt> from 0 to 2. 522 * var area1 = JXG.Math.Numerics.Romberg([0, 2], f); 523 * 524 * // the same with an anonymous function 525 * var area2 = JXG.Math.Numerics.Romberg([0, 2], function (x) { return x*x; }); 526 * 527 * // use trapez rule with maximum of 16 iterations or stop if the precision 0.0001 has been reached. 528 * var area3 = JXG.Math.Numerics.Romberg([0, 2], f, 529 * {max_iterations: 16, eps: 0.0001}); 530 * @memberof JXG.Math.Numerics 531 */ 532 Romberg: function (interval, f, config) { 533 var a, 534 b, 535 h, 536 s, 537 n, 538 k, 539 i, 540 q, 541 p = [], 542 integral = 0.0, 543 last = Infinity, 544 m = config && Type.isNumber(config.max_iterations) ? config.max_iterations : 20, 545 eps = config && Type.isNumber(config.eps) ? config.eps : config.eps || 0.0000001; 546 547 a = interval[0]; 548 b = interval[1]; 549 h = b - a; 550 n = 1; 551 552 p[0] = 0.5 * h * (f(a) + f(b)); 553 554 for (k = 0; k < m; ++k) { 555 s = 0; 556 h *= 0.5; 557 n *= 2; 558 q = 1; 559 560 for (i = 1; i < n; i += 2) { 561 s += f(a + i * h); 562 } 563 564 p[k + 1] = 0.5 * p[k] + s * h; 565 566 integral = p[k + 1]; 567 for (i = k - 1; i >= 0; --i) { 568 q *= 4; 569 p[i] = p[i + 1] + (p[i + 1] - p[i]) / (q - 1.0); 570 integral = p[i]; 571 } 572 573 if (Math.abs(integral - last) < eps * Math.abs(integral)) { 574 break; 575 } 576 last = integral; 577 } 578 579 return integral; 580 }, 581 582 /** 583 * Calculates the integral of function f over interval using Gauss-Legendre quadrature. 584 * @param {Array} interval The integration interval, e.g. [0, 3]. 585 * @param {function} f A function which takes one argument of type number and returns a number. 586 * @param {Object} [config] The algorithm setup. Accepted property is the order n of type number. n is allowed to take 587 * values between 2 and 18, default value is 12. 588 * @param {Number} [config.n=16] 589 * @returns {Number} Integral value of f over interval 590 * @example 591 * function f(x) { 592 * return x*x; 593 * } 594 * 595 * // calculates integral of <tt>f</tt> from 0 to 2. 596 * var area1 = JXG.Math.Numerics.GaussLegendre([0, 2], f); 597 * 598 * // the same with an anonymous function 599 * var area2 = JXG.Math.Numerics.GaussLegendre([0, 2], function (x) { return x*x; }); 600 * 601 * // use 16 point Gauss-Legendre rule. 602 * var area3 = JXG.Math.Numerics.GaussLegendre([0, 2], f, 603 * {n: 16}); 604 * @memberof JXG.Math.Numerics 605 */ 606 GaussLegendre: function (interval, f, config) { 607 var a, 608 b, 609 i, 610 m, 611 xp, 612 xm, 613 result = 0.0, 614 table_xi = [], 615 table_w = [], 616 xi, 617 w, 618 n = config && Type.isNumber(config.n) ? config.n : 12; 619 620 if (n > 18) { 621 n = 18; 622 } 623 624 /* n = 2 */ 625 table_xi[2] = [0.5773502691896257645091488]; 626 table_w[2] = [1.0]; 627 628 /* n = 4 */ 629 table_xi[4] = [0.3399810435848562648026658, 0.8611363115940525752239465]; 630 table_w[4] = [0.6521451548625461426269361, 0.3478548451374538573730639]; 631 632 /* n = 6 */ 633 table_xi[6] = [ 634 0.2386191860831969086305017, 0.6612093864662645136613996, 635 0.9324695142031520278123016 636 ]; 637 table_w[6] = [ 638 0.4679139345726910473898703, 0.3607615730481386075698335, 639 0.1713244923791703450402961 640 ]; 641 642 /* n = 8 */ 643 table_xi[8] = [ 644 0.1834346424956498049394761, 0.525532409916328985817739, 645 0.7966664774136267395915539, 0.9602898564975362316835609 646 ]; 647 table_w[8] = [ 648 0.3626837833783619829651504, 0.3137066458778872873379622, 649 0.222381034453374470544356, 0.1012285362903762591525314 650 ]; 651 652 /* n = 10 */ 653 table_xi[10] = [ 654 0.148874338981631210884826, 0.4333953941292471907992659, 655 0.6794095682990244062343274, 0.8650633666889845107320967, 0.973906528517171720077964 656 ]; 657 table_w[10] = [ 658 0.295524224714752870173893, 0.2692667193099963550912269, 659 0.2190863625159820439955349, 0.1494513491505805931457763, 660 0.0666713443086881375935688 661 ]; 662 663 /* n = 12 */ 664 table_xi[12] = [ 665 0.1252334085114689154724414, 0.3678314989981801937526915, 666 0.5873179542866174472967024, 0.7699026741943046870368938, 667 0.9041172563704748566784659, 0.9815606342467192506905491 668 ]; 669 table_w[12] = [ 670 0.2491470458134027850005624, 0.2334925365383548087608499, 671 0.2031674267230659217490645, 0.1600783285433462263346525, 672 0.1069393259953184309602547, 0.047175336386511827194616 673 ]; 674 675 /* n = 14 */ 676 table_xi[14] = [ 677 0.1080549487073436620662447, 0.3191123689278897604356718, 678 0.5152486363581540919652907, 0.6872929048116854701480198, 679 0.8272013150697649931897947, 0.9284348836635735173363911, 680 0.9862838086968123388415973 681 ]; 682 table_w[14] = [ 683 0.2152638534631577901958764, 0.2051984637212956039659241, 684 0.1855383974779378137417166, 0.1572031671581935345696019, 685 0.1215185706879031846894148, 0.0801580871597602098056333, 686 0.0351194603317518630318329 687 ]; 688 689 /* n = 16 */ 690 table_xi[16] = [ 691 0.0950125098376374401853193, 0.2816035507792589132304605, 692 0.4580167776572273863424194, 0.6178762444026437484466718, 693 0.7554044083550030338951012, 0.8656312023878317438804679, 694 0.9445750230732325760779884, 0.9894009349916499325961542 695 ]; 696 table_w[16] = [ 697 0.1894506104550684962853967, 0.1826034150449235888667637, 698 0.1691565193950025381893121, 0.1495959888165767320815017, 699 0.1246289712555338720524763, 0.0951585116824927848099251, 700 0.0622535239386478928628438, 0.0271524594117540948517806 701 ]; 702 703 /* n = 18 */ 704 table_xi[18] = [ 705 0.0847750130417353012422619, 0.2518862256915055095889729, 706 0.4117511614628426460359318, 0.5597708310739475346078715, 707 0.6916870430603532078748911, 0.8037049589725231156824175, 708 0.8926024664975557392060606, 0.9558239495713977551811959, 0.991565168420930946730016 709 ]; 710 table_w[18] = [ 711 0.1691423829631435918406565, 0.1642764837458327229860538, 712 0.154684675126265244925418, 0.1406429146706506512047313, 713 0.1225552067114784601845191, 0.100942044106287165562814, 714 0.0764257302548890565291297, 0.0497145488949697964533349, 715 0.0216160135264833103133427 716 ]; 717 718 /* n = 3 */ 719 table_xi[3] = [0.0, 0.7745966692414833770358531]; 720 table_w[3] = [0.8888888888888888888888889, 0.5555555555555555555555556]; 721 722 /* n = 5 */ 723 table_xi[5] = [0.0, 0.5384693101056830910363144, 0.9061798459386639927976269]; 724 table_w[5] = [ 725 0.5688888888888888888888889, 0.4786286704993664680412915, 0.236926885056189087514264 726 ]; 727 728 /* n = 7 */ 729 table_xi[7] = [ 730 0.0, 0.4058451513773971669066064, 0.7415311855993944398638648, 731 0.9491079123427585245261897 732 ]; 733 table_w[7] = [ 734 0.417959183673469387755102, 0.3818300505051189449503698, 735 0.2797053914892766679014678, 0.1294849661688696932706114 736 ]; 737 738 /* n = 9 */ 739 table_xi[9] = [ 740 0.0, 0.324253423403808929038538, 0.613371432700590397308702, 741 0.8360311073266357942994298, 0.9681602395076260898355762 742 ]; 743 table_w[9] = [ 744 0.3302393550012597631645251, 0.3123470770400028400686304, 745 0.2606106964029354623187429, 0.180648160694857404058472, 0.0812743883615744119718922 746 ]; 747 748 /* n = 11 */ 749 table_xi[11] = [ 750 0.0, 0.269543155952344972331532, 0.5190961292068118159257257, 751 0.7301520055740493240934163, 0.8870625997680952990751578, 0.978228658146056992803938 752 ]; 753 table_w[11] = [ 754 0.2729250867779006307144835, 0.2628045445102466621806889, 755 0.2331937645919904799185237, 0.1862902109277342514260976, 756 0.1255803694649046246346943, 0.0556685671161736664827537 757 ]; 758 759 /* n = 13 */ 760 table_xi[13] = [ 761 0.0, 0.2304583159551347940655281, 0.4484927510364468528779129, 762 0.6423493394403402206439846, 0.8015780907333099127942065, 763 0.9175983992229779652065478, 0.9841830547185881494728294 764 ]; 765 table_w[13] = [ 766 0.2325515532308739101945895, 0.2262831802628972384120902, 767 0.2078160475368885023125232, 0.1781459807619457382800467, 768 0.1388735102197872384636018, 0.0921214998377284479144218, 769 0.0404840047653158795200216 770 ]; 771 772 /* n = 15 */ 773 table_xi[15] = [ 774 0.0, 0.2011940939974345223006283, 0.3941513470775633698972074, 775 0.5709721726085388475372267, 0.7244177313601700474161861, 776 0.8482065834104272162006483, 0.9372733924007059043077589, 777 0.9879925180204854284895657 778 ]; 779 table_w[15] = [ 780 0.2025782419255612728806202, 0.1984314853271115764561183, 781 0.1861610000155622110268006, 0.1662692058169939335532009, 782 0.1395706779261543144478048, 0.1071592204671719350118695, 783 0.0703660474881081247092674, 0.0307532419961172683546284 784 ]; 785 786 /* n = 17 */ 787 table_xi[17] = [ 788 0.0, 0.1784841814958478558506775, 0.3512317634538763152971855, 789 0.5126905370864769678862466, 0.6576711592166907658503022, 790 0.7815140038968014069252301, 0.8802391537269859021229557, 791 0.950675521768767761222717, 0.990575475314417335675434 792 ]; 793 table_w[17] = [ 794 0.1794464703562065254582656, 0.176562705366992646325271, 795 0.1680041021564500445099707, 0.1540457610768102880814316, 0.13513636846852547328632, 796 0.1118838471934039710947884, 0.0850361483171791808835354, 797 0.0554595293739872011294402, 0.02414830286854793196011 798 ]; 799 800 a = interval[0]; 801 b = interval[1]; 802 803 //m = Math.ceil(n * 0.5); 804 m = (n + 1) >> 1; 805 806 xi = table_xi[n]; 807 w = table_w[n]; 808 809 xm = 0.5 * (b - a); 810 xp = 0.5 * (b + a); 811 812 if (n & (1 === 1)) { 813 // n odd 814 result = w[0] * f(xp); 815 for (i = 1; i < m; ++i) { 816 result += w[i] * (f(xp + xm * xi[i]) + f(xp - xm * xi[i])); 817 } 818 } else { 819 // n even 820 result = 0.0; 821 for (i = 0; i < m; ++i) { 822 result += w[i] * (f(xp + xm * xi[i]) + f(xp - xm * xi[i])); 823 } 824 } 825 826 return xm * result; 827 }, 828 829 /** 830 * Scale error in Gauss Kronrod quadrature. 831 * Internal method used in {@link JXG.Math.Numerics._gaussKronrod}. 832 * @private 833 */ 834 _rescale_error: function (err, result_abs, result_asc) { 835 var scale, 836 min_err, 837 DBL_MIN = 2.2250738585072014e-308, 838 DBL_EPS = 2.2204460492503131e-16; 839 840 err = Math.abs(err); 841 if (result_asc !== 0 && err !== 0) { 842 scale = Math.pow((200 * err) / result_asc, 1.5); 843 844 if (scale < 1.0) { 845 err = result_asc * scale; 846 } else { 847 err = result_asc; 848 } 849 } 850 if (result_abs > DBL_MIN / (50 * DBL_EPS)) { 851 min_err = 50 * DBL_EPS * result_abs; 852 853 if (min_err > err) { 854 err = min_err; 855 } 856 } 857 858 return err; 859 }, 860 861 /** 862 * Generic Gauss-Kronrod quadrature algorithm. 863 * Internal method used in {@link JXG.Math.Numerics.GaussKronrod15}, 864 * {@link JXG.Math.Numerics.GaussKronrod21}, 865 * {@link JXG.Math.Numerics.GaussKronrod31}. 866 * Taken from QUADPACK. 867 * 868 * @param {Array} interval The integration interval, e.g. [0, 3]. 869 * @param {function} f A function which takes one argument of type number and returns a number. 870 * @param {Number} n order of approximation. Actually, n is the length of the array xgk. For example, for the 15-point Kronrod rule, n is equal to 8. 871 * @param {Array} xgk Kronrod quadrature abscissae 872 * @param {Array} wg Weights of the Gauss rule 873 * @param {Array} wgk Weights of the Kronrod rule 874 * @param {Object} resultObj Object returning resultObj.abserr, resultObj.resabs, resultObj.resasc. 875 * See the library QUADPACK for an explanation. 876 * 877 * @returns {Number} Integral value of f over interval 878 * 879 * @private 880 */ 881 _gaussKronrod: function (interval, f, n, xgk, wg, wgk, resultObj) { 882 var a = interval[0], 883 b = interval[1], 884 up, 885 result, 886 center = 0.5 * (a + b), 887 half_length = 0.5 * (b - a), 888 abs_half_length = Math.abs(half_length), 889 f_center = f(center), 890 result_gauss = 0.0, 891 result_kronrod = f_center * wgk[n - 1], 892 result_abs = Math.abs(result_kronrod), 893 result_asc = 0.0, 894 mean = 0.0, 895 err = 0.0, 896 j, jtw, jtwm1, 897 abscissa, 898 fval1, fval2, fsum, 899 fv1 = [], 900 fv2 = []; 901 902 if (n % 2 === 0) { 903 result_gauss = f_center * wg[n / 2 - 1]; 904 } 905 906 up = Math.floor((n - 1) / 2); 907 for (j = 0; j < up; j++) { 908 jtw = j * 2 + 1; // in original fortran j=1,2,3 jtw=2,4,6 909 abscissa = half_length * xgk[jtw]; 910 fval1 = f(center - abscissa); 911 fval2 = f(center + abscissa); 912 fsum = fval1 + fval2; 913 fv1[jtw] = fval1; 914 fv2[jtw] = fval2; 915 result_gauss += wg[j] * fsum; 916 result_kronrod += wgk[jtw] * fsum; 917 result_abs += wgk[jtw] * (Math.abs(fval1) + Math.abs(fval2)); 918 } 919 920 up = Math.floor(n / 2); 921 for (j = 0; j < up; j++) { 922 jtwm1 = j * 2; 923 abscissa = half_length * xgk[jtwm1]; 924 fval1 = f(center - abscissa); 925 fval2 = f(center + abscissa); 926 fv1[jtwm1] = fval1; 927 fv2[jtwm1] = fval2; 928 result_kronrod += wgk[jtwm1] * (fval1 + fval2); 929 result_abs += wgk[jtwm1] * (Math.abs(fval1) + Math.abs(fval2)); 930 } 931 932 mean = result_kronrod * 0.5; 933 result_asc = wgk[n - 1] * Math.abs(f_center - mean); 934 935 for (j = 0; j < n - 1; j++) { 936 result_asc += wgk[j] * (Math.abs(fv1[j] - mean) + Math.abs(fv2[j] - mean)); 937 } 938 939 // scale by the width of the integration region 940 err = (result_kronrod - result_gauss) * half_length; 941 942 result_kronrod *= half_length; 943 result_abs *= abs_half_length; 944 result_asc *= abs_half_length; 945 result = result_kronrod; 946 947 resultObj.abserr = this._rescale_error(err, result_abs, result_asc); 948 resultObj.resabs = result_abs; 949 resultObj.resasc = result_asc; 950 951 return result; 952 }, 953 954 /** 955 * 15-point Gauss-Kronrod quadrature algorithm, see the library QUADPACK 956 * @param {Array} interval The integration interval, e.g. [0, 3]. 957 * @param {function} f A function which takes one argument of type number and returns a number. 958 * @param {Object} resultObj Object returning resultObj.abserr, resultObj.resabs, resultObj.resasc. See the library 959 * QUADPACK for an explanation. 960 * 961 * @returns {Number} Integral value of f over interval 962 * 963 * @memberof JXG.Math.Numerics 964 */ 965 GaussKronrod15: function (interval, f, resultObj) { 966 /* Gauss quadrature weights and kronrod quadrature abscissae and 967 weights as evaluated with 80 decimal digit arithmetic by 968 L. W. Fullerton, Bell Labs, Nov. 1981. */ 969 970 var xgk = 971 /* abscissae of the 15-point kronrod rule */ 972 [ 973 0.991455371120812639206854697526329, 0.949107912342758524526189684047851, 974 0.864864423359769072789712788640926, 0.741531185599394439863864773280788, 975 0.58608723546769113029414483825873, 0.405845151377397166906606412076961, 976 0.207784955007898467600689403773245, 0.0 977 ], 978 /* xgk[1], xgk[3], ... abscissae of the 7-point gauss rule. 979 xgk[0], xgk[2], ... abscissae to optimally extend the 7-point gauss rule */ 980 981 wg = 982 /* weights of the 7-point gauss rule */ 983 [ 984 0.129484966168869693270611432679082, 0.27970539148927666790146777142378, 985 0.381830050505118944950369775488975, 0.417959183673469387755102040816327 986 ], 987 wgk = 988 /* weights of the 15-point kronrod rule */ 989 [ 990 0.02293532201052922496373200805897, 0.063092092629978553290700663189204, 991 0.104790010322250183839876322541518, 0.140653259715525918745189590510238, 992 0.16900472663926790282658342659855, 0.190350578064785409913256402421014, 993 0.204432940075298892414161999234649, 0.209482141084727828012999174891714 994 ]; 995 996 return this._gaussKronrod(interval, f, 8, xgk, wg, wgk, resultObj); 997 }, 998 999 /** 1000 * 21 point Gauss-Kronrod quadrature algorithm, see the library QUADPACK 1001 * @param {Array} interval The integration interval, e.g. [0, 3]. 1002 * @param {function} f A function which takes one argument of type number and returns a number. 1003 * @param {Object} resultObj Object returning resultObj.abserr, resultObj.resabs, resultObj.resasc. See the library 1004 * QUADPACK for an explanation. 1005 * 1006 * @returns {Number} Integral value of f over interval 1007 * 1008 * @memberof JXG.Math.Numerics 1009 */ 1010 GaussKronrod21: function (interval, f, resultObj) { 1011 /* Gauss quadrature weights and kronrod quadrature abscissae and 1012 weights as evaluated with 80 decimal digit arithmetic by 1013 L. W. Fullerton, Bell Labs, Nov. 1981. */ 1014 1015 var xgk = 1016 /* abscissae of the 21-point kronrod rule */ 1017 [ 1018 0.995657163025808080735527280689003, 0.973906528517171720077964012084452, 1019 0.930157491355708226001207180059508, 0.865063366688984510732096688423493, 1020 0.780817726586416897063717578345042, 0.679409568299024406234327365114874, 1021 0.562757134668604683339000099272694, 0.433395394129247190799265943165784, 1022 0.294392862701460198131126603103866, 0.14887433898163121088482600112972, 0.0 1023 ], 1024 /* xgk[1], xgk[3], ... abscissae of the 10-point gauss rule. 1025 xgk[0], xgk[2], ... abscissae to optimally extend the 10-point gauss rule */ 1026 wg = 1027 /* weights of the 10-point gauss rule */ 1028 [ 1029 0.066671344308688137593568809893332, 0.149451349150580593145776339657697, 1030 0.219086362515982043995534934228163, 0.269266719309996355091226921569469, 1031 0.295524224714752870173892994651338 1032 ], 1033 wgk = 1034 /* weights of the 21-point kronrod rule */ 1035 [ 1036 0.011694638867371874278064396062192, 0.03255816230796472747881897245939, 1037 0.05475589657435199603138130024458, 0.07503967481091995276704314091619, 1038 0.093125454583697605535065465083366, 0.109387158802297641899210590325805, 1039 0.123491976262065851077958109831074, 0.134709217311473325928054001771707, 1040 0.142775938577060080797094273138717, 0.147739104901338491374841515972068, 1041 0.149445554002916905664936468389821 1042 ]; 1043 1044 return this._gaussKronrod(interval, f, 11, xgk, wg, wgk, resultObj); 1045 }, 1046 1047 /** 1048 * 31 point Gauss-Kronrod quadrature algorithm, see the library QUADPACK 1049 * @param {Array} interval The integration interval, e.g. [0, 3]. 1050 * @param {function} f A function which takes one argument of type number and returns a number. 1051 * @param {Object} resultObj Object returning resultObj.abserr, resultObj.resabs, resultObj.resasc. See the library 1052 * QUADPACK for an explanation. 1053 * 1054 * @returns {Number} Integral value of f over interval 1055 * 1056 * @memberof JXG.Math.Numerics 1057 */ 1058 GaussKronrod31: function (interval, f, resultObj) { 1059 /* Gauss quadrature weights and kronrod quadrature abscissae and 1060 weights as evaluated with 80 decimal digit arithmetic by 1061 L. W. Fullerton, Bell Labs, Nov. 1981. */ 1062 1063 var xgk = 1064 /* abscissae of the 21-point kronrod rule */ 1065 [ 1066 0.998002298693397060285172840152271, 0.987992518020485428489565718586613, 1067 0.967739075679139134257347978784337, 0.937273392400705904307758947710209, 1068 0.897264532344081900882509656454496, 0.848206583410427216200648320774217, 1069 0.790418501442465932967649294817947, 0.724417731360170047416186054613938, 1070 0.650996741297416970533735895313275, 0.570972172608538847537226737253911, 1071 0.485081863640239680693655740232351, 0.394151347077563369897207370981045, 1072 0.299180007153168812166780024266389, 0.201194093997434522300628303394596, 1073 0.101142066918717499027074231447392, 0.0 1074 ], 1075 /* xgk[1], xgk[3], ... abscissae of the 10-point gauss rule. 1076 xgk[0], xgk[2], ... abscissae to optimally extend the 10-point gauss rule */ 1077 wg = 1078 /* weights of the 10-point gauss rule */ 1079 [ 1080 0.030753241996117268354628393577204, 0.070366047488108124709267416450667, 1081 0.107159220467171935011869546685869, 0.139570677926154314447804794511028, 1082 0.166269205816993933553200860481209, 0.186161000015562211026800561866423, 1083 0.198431485327111576456118326443839, 0.202578241925561272880620199967519 1084 ], 1085 wgk = 1086 /* weights of the 21-point kronrod rule */ 1087 [ 1088 0.005377479872923348987792051430128, 0.015007947329316122538374763075807, 1089 0.025460847326715320186874001019653, 0.03534636079137584622203794847836, 1090 0.04458975132476487660822729937328, 0.05348152469092808726534314723943, 1091 0.062009567800670640285139230960803, 0.069854121318728258709520077099147, 1092 0.076849680757720378894432777482659, 0.083080502823133021038289247286104, 1093 0.088564443056211770647275443693774, 0.093126598170825321225486872747346, 1094 0.096642726983623678505179907627589, 0.099173598721791959332393173484603, 1095 0.10076984552387559504494666261757, 0.101330007014791549017374792767493 1096 ]; 1097 1098 return this._gaussKronrod(interval, f, 16, xgk, wg, wgk, resultObj); 1099 }, 1100 1101 /** 1102 * Generate workspace object for {@link JXG.Math.Numerics.Qag}. 1103 * @param {Array} interval The integration interval, e.g. [0, 3]. 1104 * @param {Number} n Max. limit 1105 * @returns {Object} Workspace object 1106 * 1107 * @private 1108 * @memberof JXG.Math.Numerics 1109 */ 1110 _workspace: function (interval, n) { 1111 return { 1112 limit: n, 1113 size: 0, 1114 nrmax: 0, 1115 i: 0, 1116 alist: [interval[0]], 1117 blist: [interval[1]], 1118 rlist: [0.0], 1119 elist: [0.0], 1120 order: [0], 1121 level: [0], 1122 1123 qpsrt: function () { 1124 var last = this.size - 1, 1125 limit = this.limit, 1126 errmax, 1127 errmin, 1128 i, 1129 k, 1130 top, 1131 i_nrmax = this.nrmax, 1132 i_maxerr = this.order[i_nrmax]; 1133 1134 /* Check whether the list contains more than two error estimates */ 1135 if (last < 2) { 1136 this.order[0] = 0; 1137 this.order[1] = 1; 1138 this.i = i_maxerr; 1139 return; 1140 } 1141 1142 errmax = this.elist[i_maxerr]; 1143 1144 /* This part of the routine is only executed if, due to a difficult 1145 integrand, subdivision increased the error estimate. In the normal 1146 case the insert procedure should start after the nrmax-th largest 1147 error estimate. */ 1148 while (i_nrmax > 0 && errmax > this.elist[this.order[i_nrmax - 1]]) { 1149 this.order[i_nrmax] = this.order[i_nrmax - 1]; 1150 i_nrmax--; 1151 } 1152 1153 /* Compute the number of elements in the list to be maintained in 1154 descending order. This number depends on the number of 1155 subdivisions still allowed. */ 1156 if (last < limit / 2 + 2) { 1157 top = last; 1158 } else { 1159 top = limit - last + 1; 1160 } 1161 1162 /* Insert errmax by traversing the list top-down, starting 1163 comparison from the element elist(order(i_nrmax+1)). */ 1164 i = i_nrmax + 1; 1165 1166 /* The order of the tests in the following line is important to 1167 prevent a segmentation fault */ 1168 while (i < top && errmax < this.elist[this.order[i]]) { 1169 this.order[i - 1] = this.order[i]; 1170 i++; 1171 } 1172 1173 this.order[i - 1] = i_maxerr; 1174 1175 /* Insert errmin by traversing the list bottom-up */ 1176 errmin = this.elist[last]; 1177 k = top - 1; 1178 1179 while (k > i - 2 && errmin >= this.elist[this.order[k]]) { 1180 this.order[k + 1] = this.order[k]; 1181 k--; 1182 } 1183 1184 this.order[k + 1] = last; 1185 1186 /* Set i_max and e_max */ 1187 i_maxerr = this.order[i_nrmax]; 1188 this.i = i_maxerr; 1189 this.nrmax = i_nrmax; 1190 }, 1191 1192 set_initial_result: function (result, error) { 1193 this.size = 1; 1194 this.rlist[0] = result; 1195 this.elist[0] = error; 1196 }, 1197 1198 update: function (a1, b1, area1, error1, a2, b2, area2, error2) { 1199 var i_max = this.i, 1200 i_new = this.size, 1201 new_level = this.level[this.i] + 1; 1202 1203 /* append the newly-created intervals to the list */ 1204 1205 if (error2 > error1) { 1206 this.alist[i_max] = a2; /* blist[maxerr] is already == b2 */ 1207 this.rlist[i_max] = area2; 1208 this.elist[i_max] = error2; 1209 this.level[i_max] = new_level; 1210 1211 this.alist[i_new] = a1; 1212 this.blist[i_new] = b1; 1213 this.rlist[i_new] = area1; 1214 this.elist[i_new] = error1; 1215 this.level[i_new] = new_level; 1216 } else { 1217 this.blist[i_max] = b1; /* alist[maxerr] is already == a1 */ 1218 this.rlist[i_max] = area1; 1219 this.elist[i_max] = error1; 1220 this.level[i_max] = new_level; 1221 1222 this.alist[i_new] = a2; 1223 this.blist[i_new] = b2; 1224 this.rlist[i_new] = area2; 1225 this.elist[i_new] = error2; 1226 this.level[i_new] = new_level; 1227 } 1228 1229 this.size++; 1230 1231 if (new_level > this.maximum_level) { 1232 this.maximum_level = new_level; 1233 } 1234 1235 this.qpsrt(); 1236 }, 1237 1238 retrieve: function () { 1239 var i = this.i; 1240 return { 1241 a: this.alist[i], 1242 b: this.blist[i], 1243 r: this.rlist[i], 1244 e: this.elist[i] 1245 }; 1246 }, 1247 1248 sum_results: function () { 1249 var nn = this.size, 1250 k, 1251 result_sum = 0.0; 1252 1253 for (k = 0; k < nn; k++) { 1254 result_sum += this.rlist[k]; 1255 } 1256 1257 return result_sum; 1258 }, 1259 1260 subinterval_too_small: function (a1, a2, b2) { 1261 var e = 2.2204460492503131e-16, 1262 u = 2.2250738585072014e-308, 1263 tmp = (1 + 100 * e) * (Math.abs(a2) + 1000 * u); 1264 1265 return Math.abs(a1) <= tmp && Math.abs(b2) <= tmp; 1266 } 1267 }; 1268 }, 1269 1270 /** 1271 * Quadrature algorithm qag from QUADPACK. 1272 * Internal method used in {@link JXG.Math.Numerics.GaussKronrod15}, 1273 * {@link JXG.Math.Numerics.GaussKronrod21}, 1274 * {@link JXG.Math.Numerics.GaussKronrod31}. 1275 * 1276 * @param {Array} interval The integration interval, e.g. [0, 3]. 1277 * @param {function} f A function which takes one argument of type number and returns a number. 1278 * @param {Object} [config] The algorithm setup. Accepted propert are max. recursion limit of type number, 1279 * and epsrel and epsabs, the relative and absolute required precision of type number. Further, 1280 * q the internal quadrature sub-algorithm of type function. 1281 * @param {Number} [config.limit=15] 1282 * @param {Number} [config.epsrel=0.0000001] 1283 * @param {Number} [config.epsabs=0.0000001] 1284 * @param {Number} [config.q=JXG.Math.Numerics.GaussKronrod15] 1285 * @returns {Number} Integral value of f over interval 1286 * 1287 * @example 1288 * function f(x) { 1289 * return x*x; 1290 * } 1291 * 1292 * // calculates integral of <tt>f</tt> from 0 to 2. 1293 * var area1 = JXG.Math.Numerics.Qag([0, 2], f); 1294 * 1295 * // the same with an anonymous function 1296 * var area2 = JXG.Math.Numerics.Qag([0, 2], function (x) { return x*x; }); 1297 * 1298 * // use JXG.Math.Numerics.GaussKronrod31 rule as sub-algorithm. 1299 * var area3 = JXG.Math.Numerics.Quag([0, 2], f, 1300 * {q: JXG.Math.Numerics.GaussKronrod31}); 1301 * @memberof JXG.Math.Numerics 1302 */ 1303 Qag: function (interval, f, config) { 1304 var DBL_EPS = 2.2204460492503131e-16, 1305 ws = this._workspace(interval, 1000), 1306 limit = config && Type.isNumber(config.limit) ? config.limit : 15, 1307 epsrel = config && Type.isNumber(config.epsrel) ? config.epsrel : 0.0000001, 1308 epsabs = config && Type.isNumber(config.epsabs) ? config.epsabs : 0.0000001, 1309 q = config && Type.isFunction(config.q) ? config.q : this.GaussKronrod15, 1310 resultObj = {}, 1311 area, 1312 errsum, 1313 result0, 1314 abserr0, 1315 resabs0, 1316 resasc0, 1317 result, 1318 tolerance, 1319 iteration = 0, 1320 roundoff_type1 = 0, 1321 roundoff_type2 = 0, 1322 error_type = 0, 1323 round_off, 1324 a1, 1325 b1, 1326 a2, 1327 b2, 1328 a_i, 1329 b_i, 1330 r_i, 1331 e_i, 1332 area1 = 0, 1333 area2 = 0, 1334 area12 = 0, 1335 error1 = 0, 1336 error2 = 0, 1337 error12 = 0, 1338 resasc1, 1339 resasc2, 1340 // resabs1, resabs2, 1341 wsObj, 1342 delta; 1343 1344 if (limit > ws.limit) { 1345 JXG.warn("iteration limit exceeds available workspace"); 1346 } 1347 if (epsabs <= 0 && (epsrel < 50 * Mat.eps || epsrel < 0.5e-28)) { 1348 JXG.warn("tolerance cannot be acheived with given epsabs and epsrel"); 1349 } 1350 1351 result0 = q.apply(this, [interval, f, resultObj]); 1352 abserr0 = resultObj.abserr; 1353 resabs0 = resultObj.resabs; 1354 resasc0 = resultObj.resasc; 1355 1356 ws.set_initial_result(result0, abserr0); 1357 tolerance = Math.max(epsabs, epsrel * Math.abs(result0)); 1358 round_off = 50 * DBL_EPS * resabs0; 1359 1360 if (abserr0 <= round_off && abserr0 > tolerance) { 1361 result = result0; 1362 // abserr = abserr0; 1363 1364 JXG.warn("cannot reach tolerance because of roundoff error on first attempt"); 1365 return -Infinity; 1366 } 1367 1368 if ((abserr0 <= tolerance && abserr0 !== resasc0) || abserr0 === 0.0) { 1369 result = result0; 1370 // abserr = abserr0; 1371 1372 return result; 1373 } 1374 1375 if (limit === 1) { 1376 result = result0; 1377 // abserr = abserr0; 1378 1379 JXG.warn("a maximum of one iteration was insufficient"); 1380 return -Infinity; 1381 } 1382 1383 area = result0; 1384 errsum = abserr0; 1385 iteration = 1; 1386 1387 do { 1388 area1 = 0; 1389 area2 = 0; 1390 area12 = 0; 1391 error1 = 0; 1392 error2 = 0; 1393 error12 = 0; 1394 1395 /* Bisect the subinterval with the largest error estimate */ 1396 wsObj = ws.retrieve(); 1397 a_i = wsObj.a; 1398 b_i = wsObj.b; 1399 r_i = wsObj.r; 1400 e_i = wsObj.e; 1401 1402 a1 = a_i; 1403 b1 = 0.5 * (a_i + b_i); 1404 a2 = b1; 1405 b2 = b_i; 1406 1407 area1 = q.apply(this, [[a1, b1], f, resultObj]); 1408 error1 = resultObj.abserr; 1409 // resabs1 = resultObj.resabs; 1410 resasc1 = resultObj.resasc; 1411 1412 area2 = q.apply(this, [[a2, b2], f, resultObj]); 1413 error2 = resultObj.abserr; 1414 // resabs2 = resultObj.resabs; 1415 resasc2 = resultObj.resasc; 1416 1417 area12 = area1 + area2; 1418 error12 = error1 + error2; 1419 1420 errsum += error12 - e_i; 1421 area += area12 - r_i; 1422 1423 if (resasc1 !== error1 && resasc2 !== error2) { 1424 delta = r_i - area12; 1425 if (Math.abs(delta) <= 1.0e-5 * Math.abs(area12) && error12 >= 0.99 * e_i) { 1426 roundoff_type1++; 1427 } 1428 if (iteration >= 10 && error12 > e_i) { 1429 roundoff_type2++; 1430 } 1431 } 1432 1433 tolerance = Math.max(epsabs, epsrel * Math.abs(area)); 1434 1435 if (errsum > tolerance) { 1436 if (roundoff_type1 >= 6 || roundoff_type2 >= 20) { 1437 error_type = 2; /* round off error */ 1438 } 1439 1440 /* set error flag in the case of bad integrand behaviour at 1441 a point of the integration range */ 1442 1443 if (ws.subinterval_too_small(a1, a2, b2)) { 1444 error_type = 3; 1445 } 1446 } 1447 1448 ws.update(a1, b1, area1, error1, a2, b2, area2, error2); 1449 wsObj = ws.retrieve(); 1450 a_i = wsObj.a_i; 1451 b_i = wsObj.b_i; 1452 r_i = wsObj.r_i; 1453 e_i = wsObj.e_i; 1454 1455 iteration++; 1456 } while (iteration < limit && !error_type && errsum > tolerance); 1457 1458 result = ws.sum_results(); 1459 // abserr = errsum; 1460 /* 1461 if (errsum <= tolerance) 1462 { 1463 return GSL_SUCCESS; 1464 } 1465 else if (error_type == 2) 1466 { 1467 GSL_ERROR ("roundoff error prevents tolerance from being achieved", 1468 GSL_EROUND); 1469 } 1470 else if (error_type == 3) 1471 { 1472 GSL_ERROR ("bad integrand behavior found in the integration interval", 1473 GSL_ESING); 1474 } 1475 else if (iteration == limit) 1476 { 1477 GSL_ERROR ("maximum number of subdivisions reached", GSL_EMAXITER); 1478 } 1479 else 1480 { 1481 GSL_ERROR ("could not integrate function", GSL_EFAILED); 1482 } 1483 */ 1484 1485 return result; 1486 }, 1487 1488 /** 1489 * Integral of function f over interval. 1490 * @param {Array} interval The integration interval, e.g. [0, 3]. 1491 * @param {function} f A function which takes one argument of type number and returns a number. 1492 * @returns {Number} The value of the integral of f over interval 1493 * @see JXG.Math.Numerics.NewtonCotes 1494 * @see JXG.Math.Numerics.Romberg 1495 * @see JXG.Math.Numerics.Qag 1496 * @memberof JXG.Math.Numerics 1497 */ 1498 I: function (interval, f) { 1499 // return this.NewtonCotes(interval, f, {number_of_nodes: 16, integration_type: 'milne'}); 1500 // return this.Romberg(interval, f, {max_iterations: 20, eps: 0.0000001}); 1501 return this.Qag(interval, f, { 1502 q: this.GaussKronrod15, 1503 limit: 15, 1504 epsrel: 0.0000001, 1505 epsabs: 0.0000001 1506 }); 1507 }, 1508 1509 /** 1510 * Newton's method to find roots of a funtion in one variable. 1511 * @param {function} f We search for a solution of f(x)=0. 1512 * @param {Number} x initial guess for the root, i.e. start value. 1513 * @param {Object} context optional object that is treated as "this" in the function body. This is useful if 1514 * the function is a method of an object and contains a reference to its parent object via "this". 1515 * @returns {Number} A root of the function f. 1516 * @memberof JXG.Math.Numerics 1517 */ 1518 Newton: function (f, x, context) { 1519 var df, 1520 i = 0, 1521 h = Mat.eps, 1522 newf = f.apply(context, [x]); 1523 // nfev = 1; 1524 1525 // For compatibility 1526 if (Type.isArray(x)) { 1527 x = x[0]; 1528 } 1529 1530 while (i < 50 && Math.abs(newf) > h) { 1531 df = this.D(f, context)(x); 1532 // nfev += 2; 1533 1534 if (Math.abs(df) > h) { 1535 x -= newf / df; 1536 } else { 1537 x += Math.random() * 0.2 - 1.0; 1538 } 1539 1540 newf = f.apply(context, [x]); 1541 // nfev += 1; 1542 i += 1; 1543 } 1544 1545 return x; 1546 }, 1547 1548 /** 1549 * Abstract method to find roots of univariate functions, which - for the time being - 1550 * is an alias for {@link JXG.Math.Numerics.chandrupatla}. 1551 * @param {function} f We search for a solution of f(x)=0. 1552 * @param {Number|Array} x initial guess for the root, i.e. starting value, or start interval enclosing the root. 1553 * If x is an interval [a,b], it is required that f(a)f(b) <= 0, otherwise the minimum of f in [a, b] will be returned. 1554 * If x is a number, the algorithms tries to enclose the root by an interval [a, b] containing x and the root and 1555 * f(a)f(b) <= 0. If this fails, the algorithm falls back to Newton's method. 1556 * 1557 * @param {Object} context optional object that is treated as "this" in the function body. This is useful if 1558 * the function is a method of an object and contains a reference to its parent object via "this". 1559 * @returns {Number} A root of the function f. 1560 * 1561 * @see JXG.Math.Numerics.chandrupatla 1562 * @see JXG.Math.Numerics.fzero 1563 * @see JXG.Math.Numerics.polzeros 1564 * @see JXG.Math.Numerics.Newton 1565 * @memberof JXG.Math.Numerics 1566 */ 1567 root: function (f, x, context) { 1568 //return this.fzero(f, x, context); 1569 return this.chandrupatla(f, x, context); 1570 }, 1571 1572 /** 1573 * Compute an intersection of the curves c1 and c2 1574 * with a generalized Newton method. 1575 * We want to find values t1, t2 such that 1576 * c1(t1) = c2(t2), i.e. 1577 * <br> 1578 * (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2)) = (0,0). 1579 * <p> 1580 * We set 1581 * (e,f) := (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2)) 1582 * <p> 1583 * The Jacobian J is defined by 1584 * <pre> 1585 * J = (a, b) 1586 * (c, d) 1587 * </pre> 1588 * where 1589 * <ul> 1590 * <li> a = c1_x'(t1) 1591 * <li> b = -c2_x'(t2) 1592 * <li> c = c1_y'(t1) 1593 * <li> d = -c2_y'(t2) 1594 * </ul> 1595 * The inverse J^(-1) of J is equal to 1596 * <pre> 1597 * (d, -b)/ (ad-bc) 1598 * (-c, a) / (ad-bc) 1599 * </pre> 1600 * 1601 * Then, (t1new, t2new) := (t1,t2) - J^(-1)*(e,f). 1602 * <p> 1603 * If the function meetCurveCurve has the properties 1604 * t1memo and t2memo then these are taken as start values 1605 * for the Newton algorithm. 1606 * After stopping of the Newton algorithm the values of t1 and t2 are stored in 1607 * t1memo and t2memo. 1608 * 1609 * @param {JXG.Curve} c1 Curve, Line or Circle 1610 * @param {JXG.Curve} c2 Curve, Line or Circle 1611 * @param {Number} t1ini start value for t1 1612 * @param {Number} t2ini start value for t2 1613 * @returns {JXG.Coords} intersection point 1614 * @memberof JXG.Math.Numerics 1615 */ 1616 generalizedNewton: function (c1, c2, t1ini, t2ini) { 1617 var t1, t2, 1618 a, b, c, d, e, f, 1619 disc, 1620 F, 1621 D00, D01, D10, D11, 1622 count = 0; 1623 1624 if (this.generalizedNewton.t1memo) { 1625 t1 = this.generalizedNewton.t1memo; 1626 t2 = this.generalizedNewton.t2memo; 1627 } else { 1628 t1 = t1ini; 1629 t2 = t2ini; 1630 } 1631 1632 e = c1.X(t1) - c2.X(t2); 1633 f = c1.Y(t1) - c2.Y(t2); 1634 F = e * e + f * f; 1635 1636 D00 = this.D(c1.X, c1); 1637 D01 = this.D(c2.X, c2); 1638 D10 = this.D(c1.Y, c1); 1639 D11 = this.D(c2.Y, c2); 1640 1641 while (F > Mat.eps && count < 10) { 1642 a = D00(t1); 1643 b = -D01(t2); 1644 c = D10(t1); 1645 d = -D11(t2); 1646 disc = a * d - b * c; 1647 t1 -= (d * e - b * f) / disc; 1648 t2 -= (a * f - c * e) / disc; 1649 e = c1.X(t1) - c2.X(t2); 1650 f = c1.Y(t1) - c2.Y(t2); 1651 F = e * e + f * f; 1652 count += 1; 1653 } 1654 1655 this.generalizedNewton.t1memo = t1; 1656 this.generalizedNewton.t2memo = t2; 1657 1658 if (Math.abs(t1) < Math.abs(t2)) { 1659 return [c1.X(t1), c1.Y(t1)]; 1660 } 1661 1662 return [c2.X(t2), c2.Y(t2)]; 1663 }, 1664 1665 /** 1666 * Returns the Lagrange polynomials for curves with equidistant nodes, see 1667 * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation, 1668 * SIAM Review, Vol 46, No 3, (2004) 501-517. 1669 * The graph of the parametric curve [x(t),y(t)] runs through the given points. 1670 * @param {Array} p Array of JXG.Points 1671 * @returns {Array} An array consisting of two functions x(t), y(t) which define a parametric curve 1672 * f(t) = (x(t), y(t)), a number x1 (which equals 0) and a function x2 defining the curve's domain. 1673 * That means the curve is defined between x1 and x2(). x2 returns the (length of array p minus one). 1674 * @memberof JXG.Math.Numerics 1675 * 1676 * @example 1677 * var p = []; 1678 * 1679 * p[0] = board.create('point', [0, -2], {size:2, name: 'C(a)'}); 1680 * p[1] = board.create('point', [-1.5, 5], {size:2, name: ''}); 1681 * p[2] = board.create('point', [1, 4], {size:2, name: ''}); 1682 * p[3] = board.create('point', [3, 3], {size:2, name: 'C(b)'}); 1683 * 1684 * // Curve 1685 * var fg = JXG.Math.Numerics.Neville(p); 1686 * var graph = board.create('curve', fg, {strokeWidth:3, strokeOpacity:0.5}); 1687 * 1688 * </pre><div id="JXG88a8b3a8-6561-44f5-a678-76bca13fd484" class="jxgbox" style="width: 300px; height: 300px;"></div> 1689 * <script type="text/javascript"> 1690 * (function() { 1691 * var board = JXG.JSXGraph.initBoard('JXG88a8b3a8-6561-44f5-a678-76bca13fd484', 1692 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 1693 * var p = []; 1694 * 1695 * p[0] = board.create('point', [0, -2], {size:2, name: 'C(a)'}); 1696 * p[1] = board.create('point', [-1.5, 5], {size:2, name: ''}); 1697 * p[2] = board.create('point', [1, 4], {size:2, name: ''}); 1698 * p[3] = board.create('point', [3, 3], {size:2, name: 'C(b)'}); 1699 * 1700 * // Curve 1701 * var fg = JXG.Math.Numerics.Neville(p); 1702 * var graph = board.create('curve', fg, {strokeWidth:3, strokeOpacity:0.5}); 1703 * 1704 * })(); 1705 * 1706 * </script><pre> 1707 * 1708 */ 1709 Neville: function (p) { 1710 var w = [], 1711 /** @ignore */ 1712 makeFct = function (fun) { 1713 return function (t, suspendedUpdate) { 1714 var i, 1715 d, 1716 s, 1717 bin = Mat.binomial, 1718 len = p.length, 1719 len1 = len - 1, 1720 num = 0.0, 1721 denom = 0.0; 1722 1723 if (!suspendedUpdate) { 1724 s = 1; 1725 for (i = 0; i < len; i++) { 1726 w[i] = bin(len1, i) * s; 1727 s *= -1; 1728 } 1729 } 1730 1731 d = t; 1732 1733 for (i = 0; i < len; i++) { 1734 if (d === 0) { 1735 return p[i][fun](); 1736 } 1737 s = w[i] / d; 1738 d -= 1; 1739 num += p[i][fun]() * s; 1740 denom += s; 1741 } 1742 return num / denom; 1743 }; 1744 }, 1745 xfct = makeFct("X"), 1746 yfct = makeFct("Y"); 1747 1748 return [ 1749 xfct, 1750 yfct, 1751 0, 1752 function () { 1753 return p.length - 1; 1754 } 1755 ]; 1756 }, 1757 1758 /** 1759 * Calculates second derivatives at the given knots. 1760 * @param {Array} x x values of knots 1761 * @param {Array} y y values of knots 1762 * @returns {Array} Second derivatives of the interpolated function at the knots. 1763 * @see JXG.Math.Numerics.splineEval 1764 * @memberof JXG.Math.Numerics 1765 */ 1766 splineDef: function (x, y) { 1767 var pair, 1768 i, 1769 l, 1770 n = Math.min(x.length, y.length), 1771 diag = [], 1772 z = [], 1773 data = [], 1774 dx = [], 1775 delta = [], 1776 F = []; 1777 1778 if (n === 2) { 1779 return [0, 0]; 1780 } 1781 1782 for (i = 0; i < n; i++) { 1783 pair = { X: x[i], Y: y[i] }; 1784 data.push(pair); 1785 } 1786 data.sort(function (a, b) { 1787 return a.X - b.X; 1788 }); 1789 for (i = 0; i < n; i++) { 1790 x[i] = data[i].X; 1791 y[i] = data[i].Y; 1792 } 1793 1794 for (i = 0; i < n - 1; i++) { 1795 dx.push(x[i + 1] - x[i]); 1796 } 1797 for (i = 0; i < n - 2; i++) { 1798 delta.push( 1799 (6 * (y[i + 2] - y[i + 1])) / dx[i + 1] - (6 * (y[i + 1] - y[i])) / dx[i] 1800 ); 1801 } 1802 1803 // ForwardSolve 1804 diag.push(2 * (dx[0] + dx[1])); 1805 z.push(delta[0]); 1806 1807 for (i = 0; i < n - 3; i++) { 1808 l = dx[i + 1] / diag[i]; 1809 diag.push(2 * (dx[i + 1] + dx[i + 2]) - l * dx[i + 1]); 1810 z.push(delta[i + 1] - l * z[i]); 1811 } 1812 1813 // BackwardSolve 1814 F[n - 3] = z[n - 3] / diag[n - 3]; 1815 for (i = n - 4; i >= 0; i--) { 1816 F[i] = (z[i] - dx[i + 1] * F[i + 1]) / diag[i]; 1817 } 1818 1819 // Generate f''-Vector 1820 for (i = n - 3; i >= 0; i--) { 1821 F[i + 1] = F[i]; 1822 } 1823 1824 // natural cubic spline 1825 F[0] = 0; 1826 F[n - 1] = 0; 1827 1828 return F; 1829 }, 1830 1831 /** 1832 * Evaluate points on spline. 1833 * @param {Number|Array} x0 A single float value or an array of values to evaluate 1834 * @param {Array} x x values of knots 1835 * @param {Array} y y values of knots 1836 * @param {Array} F Second derivatives at knots, calculated by {@link JXG.Math.Numerics.splineDef} 1837 * @see JXG.Math.Numerics.splineDef 1838 * @returns {Number|Array} A single value or an array, depending on what is given as x0. 1839 * @memberof JXG.Math.Numerics 1840 */ 1841 splineEval: function (x0, x, y, F) { 1842 var i, 1843 j, 1844 a, 1845 b, 1846 c, 1847 d, 1848 x_, 1849 n = Math.min(x.length, y.length), 1850 l = 1, 1851 asArray = false, 1852 y0 = []; 1853 1854 // number of points to be evaluated 1855 if (Type.isArray(x0)) { 1856 l = x0.length; 1857 asArray = true; 1858 } else { 1859 x0 = [x0]; 1860 } 1861 1862 for (i = 0; i < l; i++) { 1863 // is x0 in defining interval? 1864 if (x0[i] < x[0] || x[i] > x[n - 1]) { 1865 return NaN; 1866 } 1867 1868 // determine part of spline in which x0 lies 1869 for (j = 1; j < n; j++) { 1870 if (x0[i] <= x[j]) { 1871 break; 1872 } 1873 } 1874 1875 j -= 1; 1876 1877 // we're now in the j-th partial interval, i.e. x[j] < x0[i] <= x[j+1]; 1878 // determine the coefficients of the polynomial in this interval 1879 a = y[j]; 1880 b = 1881 (y[j + 1] - y[j]) / (x[j + 1] - x[j]) - 1882 ((x[j + 1] - x[j]) / 6) * (F[j + 1] + 2 * F[j]); 1883 c = F[j] / 2; 1884 d = (F[j + 1] - F[j]) / (6 * (x[j + 1] - x[j])); 1885 // evaluate x0[i] 1886 x_ = x0[i] - x[j]; 1887 //y0.push(a + b*x_ + c*x_*x_ + d*x_*x_*x_); 1888 y0.push(a + (b + (c + d * x_) * x_) * x_); 1889 } 1890 1891 if (asArray) { 1892 return y0; 1893 } 1894 1895 return y0[0]; 1896 }, 1897 1898 /** 1899 * Generate a string containing the function term of a polynomial. 1900 * @param {Array} coeffs Coefficients of the polynomial. The position i belongs to x^i. 1901 * @param {Number} deg Degree of the polynomial 1902 * @param {String} varname Name of the variable (usually 'x') 1903 * @param {Number} prec Precision 1904 * @returns {String} A string containing the function term of the polynomial. 1905 * @memberof JXG.Math.Numerics 1906 */ 1907 generatePolynomialTerm: function (coeffs, deg, varname, prec) { 1908 var i, 1909 t = []; 1910 1911 for (i = deg; i >= 0; i--) { 1912 Type.concat(t, ["(", coeffs[i].toPrecision(prec), ")"]); 1913 1914 if (i > 1) { 1915 Type.concat(t, ["*", varname, "<sup>", i, "<", "/sup> + "]); 1916 } else if (i === 1) { 1917 Type.concat(t, ["*", varname, " + "]); 1918 } 1919 } 1920 1921 return t.join(""); 1922 }, 1923 1924 /** 1925 * Computes the polynomial through a given set of coordinates in Lagrange form. 1926 * Returns the Lagrange polynomials, see 1927 * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation, 1928 * SIAM Review, Vol 46, No 3, (2004) 501-517. 1929 * <p> 1930 * It possesses the method getTerm() which returns the string containing the function term of the polynomial and 1931 * the method getCoefficients() which returns an array containing the coefficients of the polynomial. 1932 * @param {Array} p Array of JXG.Points 1933 * @returns {function} A function of one parameter which returns the value of the polynomial, whose graph runs through the given points. 1934 * @memberof JXG.Math.Numerics 1935 * 1936 * @example 1937 * var p = []; 1938 * p[0] = board.create('point', [-1,2], {size:4}); 1939 * p[1] = board.create('point', [0,3], {size:4}); 1940 * p[2] = board.create('point', [1,1], {size:4}); 1941 * p[3] = board.create('point', [3,-1], {size:4}); 1942 * var f = JXG.Math.Numerics.lagrangePolynomial(p); 1943 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 1944 * 1945 * </pre><div id="JXGc058aa6b-74d4-41e1-af94-df06169a2d89" class="jxgbox" style="width: 300px; height: 300px;"></div> 1946 * <script type="text/javascript"> 1947 * (function() { 1948 * var board = JXG.JSXGraph.initBoard('JXGc058aa6b-74d4-41e1-af94-df06169a2d89', 1949 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 1950 * var p = []; 1951 * p[0] = board.create('point', [-1,2], {size:4}); 1952 * p[1] = board.create('point', [0,3], {size:4}); 1953 * p[2] = board.create('point', [1,1], {size:4}); 1954 * p[3] = board.create('point', [3,-1], {size:4}); 1955 * var f = JXG.Math.Numerics.lagrangePolynomial(p); 1956 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 1957 * 1958 * })(); 1959 * 1960 * </script><pre> 1961 * 1962 * @example 1963 * var points = []; 1964 * points[0] = board.create('point', [-1,2], {size:4}); 1965 * points[1] = board.create('point', [0, 0], {size:4}); 1966 * points[2] = board.create('point', [2, 1], {size:4}); 1967 * 1968 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 1969 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 1970 * var txt = board.create('text', [-3, -4, () => f.getTerm(2, 't', ' * ')], {fontSize: 16}); 1971 * var txt2 = board.create('text', [-3, -6, () => f.getCoefficients()], {fontSize: 12}); 1972 * 1973 * </pre><div id="JXG73fdaf12-e257-4374-b488-ae063e4eecbb" class="jxgbox" style="width: 300px; height: 300px;"></div> 1974 * <script type="text/javascript"> 1975 * (function() { 1976 * var board = JXG.JSXGraph.initBoard('JXG73fdaf12-e257-4374-b488-ae063e4eecbb', 1977 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 1978 * var points = []; 1979 * points[0] = board.create('point', [-1,2], {size:4}); 1980 * points[1] = board.create('point', [0, 0], {size:4}); 1981 * points[2] = board.create('point', [2, 1], {size:4}); 1982 * 1983 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 1984 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 1985 * var txt = board.create('text', [-3, -4, () => f.getTerm(2, 't', ' * ')], {fontSize: 16}); 1986 * var txt2 = board.create('text', [-3, -6, () => f.getCoefficients()], {fontSize: 12}); 1987 * 1988 * })(); 1989 * 1990 * </script><pre> 1991 * 1992 */ 1993 lagrangePolynomial: function (p) { 1994 var w = [], 1995 that = this, 1996 /** @ignore */ 1997 fct = function (x, suspendedUpdate) { 1998 var i, // j, 1999 k, 2000 xi, 2001 s, //M, 2002 len = p.length, 2003 num = 0, 2004 denom = 0; 2005 2006 if (!suspendedUpdate) { 2007 for (i = 0; i < len; i++) { 2008 w[i] = 1.0; 2009 xi = p[i].X(); 2010 2011 for (k = 0; k < len; k++) { 2012 if (k !== i) { 2013 w[i] *= xi - p[k].X(); 2014 } 2015 } 2016 2017 w[i] = 1 / w[i]; 2018 } 2019 2020 // M = []; 2021 // for (k = 0; k < len; k++) { 2022 // M.push([1]); 2023 // } 2024 } 2025 2026 for (i = 0; i < len; i++) { 2027 xi = p[i].X(); 2028 2029 if (x === xi) { 2030 return p[i].Y(); 2031 } 2032 2033 s = w[i] / (x - xi); 2034 denom += s; 2035 num += s * p[i].Y(); 2036 } 2037 2038 return num / denom; 2039 }; 2040 2041 /** 2042 * Get the term of the Lagrange polynomial as string. 2043 * Calls {@link JXG.Math.Numerics#lagrangePolynomialTerm}. 2044 * 2045 * @name JXG.Math.Numerics.lagrangePolynomial#getTerm 2046 * @param {Number} digits Number of digits of the coefficients 2047 * @param {String} param Variable name 2048 * @param {String} dot Dot symbol 2049 * @returns {String} containing the term of Lagrange polynomial as string. 2050 * @see JXG.Math.Numerics.lagrangePolynomialTerm 2051 * @example 2052 * var points = []; 2053 * points[0] = board.create('point', [-1,2], {size:4}); 2054 * points[1] = board.create('point', [0, 0], {size:4}); 2055 * points[2] = board.create('point', [2, 1], {size:4}); 2056 * 2057 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2058 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2059 * var txt = board.create('text', [-3, -4, () => f.getTerm(2, 't', ' * ')], {fontSize: 16}); 2060 * 2061 * </pre><div id="JXG73fdaf12-e257-4374-b488-ae063e4eeccf" class="jxgbox" style="width: 300px; height: 300px;"></div> 2062 * <script type="text/javascript"> 2063 * (function() { 2064 * var board = JXG.JSXGraph.initBoard('JXG73fdaf12-e257-4374-b488-ae063e4eeccf', 2065 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 2066 * var points = []; 2067 * points[0] = board.create('point', [-1,2], {size:4}); 2068 * points[1] = board.create('point', [0, 0], {size:4}); 2069 * points[2] = board.create('point', [2, 1], {size:4}); 2070 * 2071 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2072 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2073 * var txt = board.create('text', [-3, -4, () => f.getTerm(2, 't', ' * ')], {fontSize: 16}); 2074 * 2075 * })(); 2076 * 2077 * </script><pre> 2078 * 2079 */ 2080 fct.getTerm = function (digits, param, dot) { 2081 return that.lagrangePolynomialTerm(p, digits, param, dot)(); 2082 }; 2083 2084 /** 2085 * Get the coefficients of the Lagrange polynomial as array. The leading 2086 * coefficient is at position 0. 2087 * Calls {@link JXG.Math.Numerics#lagrangePolynomialCoefficients}. 2088 * 2089 * @name JXG.Math.Numerics.lagrangePolynomial#getCoefficients 2090 * @returns {Array} containing the coefficients of the Lagrange polynomial. 2091 * @see JXG.Math.Numerics.lagrangePolynomial.getTerm 2092 * @see JXG.Math.Numerics.lagrangePolynomialTerm 2093 * @see JXG.Math.Numerics.lagrangePolynomialCoefficients 2094 * @example 2095 * var points = []; 2096 * points[0] = board.create('point', [-1,2], {size:4}); 2097 * points[1] = board.create('point', [0, 0], {size:4}); 2098 * points[2] = board.create('point', [2, 1], {size:4}); 2099 * 2100 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2101 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2102 * var txt = board.create('text', [1, -4, () => f.getCoefficients()], {fontSize: 10}); 2103 * 2104 * </pre><div id="JXG52a883a5-2e0c-4caf-8f84-8650c173c365" class="jxgbox" style="width: 300px; height: 300px;"></div> 2105 * <script type="text/javascript"> 2106 * (function() { 2107 * var board = JXG.JSXGraph.initBoard('JXG52a883a5-2e0c-4caf-8f84-8650c173c365', 2108 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 2109 * var points = []; 2110 * points[0] = board.create('point', [-1,2], {size:4}); 2111 * points[1] = board.create('point', [0, 0], {size:4}); 2112 * points[2] = board.create('point', [2, 1], {size:4}); 2113 * 2114 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2115 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2116 * var txt = board.create('text', [1, -4, () => f.getCoefficients()], {fontSize: 10}); 2117 * 2118 * })(); 2119 * 2120 * </script><pre> 2121 * 2122 */ 2123 fct.getCoefficients = function () { 2124 return that.lagrangePolynomialCoefficients(p)(); 2125 }; 2126 2127 return fct; 2128 }, 2129 2130 /** 2131 * Determine the Lagrange polynomial through an array of points and 2132 * return the term of the polynomial as string. 2133 * 2134 * @param {Array} points Array of JXG.Points 2135 * @param {Number} digits Number of decimal digits of the coefficients 2136 * @param {String} param Name of the parameter. Default: 'x'. 2137 * @param {String} dot Multiplication symbol. Default: ' * '. 2138 * @returns {Function} returning the Lagrange polynomial term through 2139 * the supplied points as string 2140 * @memberof JXG.Math.Numerics 2141 * 2142 * @example 2143 * var points = []; 2144 * points[0] = board.create('point', [-1,2], {size:4}); 2145 * points[1] = board.create('point', [0, 0], {size:4}); 2146 * points[2] = board.create('point', [2, 1], {size:4}); 2147 * 2148 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2149 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2150 * 2151 * var f_txt = JXG.Math.Numerics.lagrangePolynomialTerm(points, 2, 't', ' * '); 2152 * var txt = board.create('text', [-3, -4, f_txt], {fontSize: 16}); 2153 * 2154 * </pre><div id="JXGd45e9e96-7526-486d-aa43-e1178d5f2baa" class="jxgbox" style="width: 300px; height: 300px;"></div> 2155 * <script type="text/javascript"> 2156 * (function() { 2157 * var board = JXG.JSXGraph.initBoard('JXGd45e9e96-7526-486d-aa43-e1178d5f2baa', 2158 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 2159 * var points = []; 2160 * points[0] = board.create('point', [-1,2], {size:4}); 2161 * points[1] = board.create('point', [0, 0], {size:4}); 2162 * points[2] = board.create('point', [2, 1], {size:4}); 2163 * 2164 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2165 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2166 * 2167 * var f_txt = JXG.Math.Numerics.lagrangePolynomialTerm(points, 2, 't', ' * '); 2168 * var txt = board.create('text', [-3, -4, f_txt], {fontSize: 16}); 2169 * 2170 * })(); 2171 * 2172 * </script><pre> 2173 * 2174 */ 2175 lagrangePolynomialTerm: function (points, digits, param, dot) { 2176 var that = this; 2177 2178 return function () { 2179 var len = points.length, 2180 coeffs = [], 2181 isLeading = true, 2182 n, t, j, c; 2183 2184 param = param || "x"; 2185 if (dot === undefined) { 2186 dot = " * "; 2187 } 2188 2189 n = len - 1; // (Max) degree of the polynomial 2190 coeffs = that.lagrangePolynomialCoefficients(points)(); 2191 2192 t = ""; 2193 for (j = 0; j < coeffs.length; j++) { 2194 c = coeffs[j]; 2195 if (Math.abs(c) < Mat.eps) { 2196 continue; 2197 } 2198 if (JXG.exists(digits)) { 2199 c = Env._round10(c, -digits); 2200 } 2201 if (isLeading) { 2202 t += c > 0 ? c : "-" + -c; 2203 isLeading = false; 2204 } else { 2205 t += c > 0 ? " + " + c : " - " + -c; 2206 } 2207 2208 if (n - j > 1) { 2209 t += dot + param + "^" + (n - j); 2210 } else if (n - j === 1) { 2211 t += dot + param; 2212 } 2213 } 2214 return t; // board.jc.manipulate('f = map(x) -> ' + t + ';'); 2215 }; 2216 }, 2217 2218 /** 2219 * Determine the Lagrange polynomial through an array of points and 2220 * return the coefficients of the polynomial as array. 2221 * The leading coefficient is at position 0. 2222 * 2223 * @param {Array} points Array of JXG.Points 2224 * @returns {Function} returning the coefficients of the Lagrange polynomial through 2225 * the supplied points. 2226 * @memberof JXG.Math.Numerics 2227 * 2228 * @example 2229 * var points = []; 2230 * points[0] = board.create('point', [-1,2], {size:4}); 2231 * points[1] = board.create('point', [0, 0], {size:4}); 2232 * points[2] = board.create('point', [2, 1], {size:4}); 2233 * 2234 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2235 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2236 * 2237 * var f_arr = JXG.Math.Numerics.lagrangePolynomialCoefficients(points); 2238 * var txt = board.create('text', [1, -4, f_arr], {fontSize: 10}); 2239 * 2240 * </pre><div id="JXG1778f0d1-a420-473f-99e8-1755ef4be97e" class="jxgbox" style="width: 300px; height: 300px;"></div> 2241 * <script type="text/javascript"> 2242 * (function() { 2243 * var board = JXG.JSXGraph.initBoard('JXG1778f0d1-a420-473f-99e8-1755ef4be97e', 2244 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 2245 * var points = []; 2246 * points[0] = board.create('point', [-1,2], {size:4}); 2247 * points[1] = board.create('point', [0, 0], {size:4}); 2248 * points[2] = board.create('point', [2, 1], {size:4}); 2249 * 2250 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2251 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2252 * 2253 * var f_arr = JXG.Math.Numerics.lagrangePolynomialCoefficients(points); 2254 * var txt = board.create('text', [1, -4, f_arr], {fontSize: 10}); 2255 * 2256 * })(); 2257 * 2258 * </script><pre> 2259 * 2260 */ 2261 lagrangePolynomialCoefficients: function (points) { 2262 return function () { 2263 var len = points.length, 2264 zeroes = [], 2265 coeffs = [], 2266 coeffs_sum = [], 2267 i, j, c, p; 2268 2269 // n = len - 1; // (Max) degree of the polynomial 2270 for (j = 0; j < len; j++) { 2271 coeffs_sum[j] = 0; 2272 } 2273 2274 for (i = 0; i < len; i++) { 2275 c = points[i].Y(); 2276 p = points[i].X(); 2277 zeroes = []; 2278 for (j = 0; j < len; j++) { 2279 if (j !== i) { 2280 c /= p - points[j].X(); 2281 zeroes.push(points[j].X()); 2282 } 2283 } 2284 coeffs = [1].concat(Mat.Vieta(zeroes)); 2285 for (j = 0; j < coeffs.length; j++) { 2286 coeffs_sum[j] += (j % 2 === 1 ? -1 : 1) * coeffs[j] * c; 2287 } 2288 } 2289 2290 return coeffs_sum; 2291 }; 2292 }, 2293 2294 /** 2295 * Determine the coefficients of a cardinal spline polynom, See 2296 * https://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections 2297 * @param {Number} x1 point 1 2298 * @param {Number} x2 point 2 2299 * @param {Number} t1 tangent slope 1 2300 * @param {Number} t2 tangent slope 2 2301 * @return {Array} coefficents array c for the polynomial t maps to 2302 * c[0] + c[1]*t + c[2]*t*t + c[3]*t*t*t 2303 */ 2304 _initCubicPoly: function (x1, x2, t1, t2) { 2305 return [x1, t1, -3 * x1 + 3 * x2 - 2 * t1 - t2, 2 * x1 - 2 * x2 + t1 + t2]; 2306 }, 2307 2308 /** 2309 * Computes the cubic cardinal spline curve through a given set of points. The curve 2310 * is uniformly parametrized. 2311 * Two artificial control points at the beginning and the end are added. 2312 * 2313 * The implementation (especially the centripetal parametrization) is from 2314 * https://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections . 2315 * @param {Array} points Array consisting of JXG.Points. 2316 * @param {Number|Function} tau The tension parameter, either a constant number or a function returning a number. This number is between 0 and 1. 2317 * tau=1/2 give Catmull-Rom splines. 2318 * @param {String} type (Optional) parameter which allows to choose between "uniform" (default) and 2319 * "centripetal" parameterization. Thus the two possible values are "uniform" or "centripetal". 2320 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 2321 * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, 2322 * and a function simply returning the length of the points array 2323 * minus three. 2324 * @memberof JXG.Math.Numerics 2325 */ 2326 CardinalSpline: function (points, tau_param, type) { 2327 var p, 2328 coeffs = [], 2329 makeFct, 2330 tau, _tau, 2331 that = this; 2332 2333 if (Type.isFunction(tau_param)) { 2334 _tau = tau_param; 2335 } else { 2336 _tau = function () { 2337 return tau_param; 2338 }; 2339 } 2340 2341 if (type === undefined) { 2342 type = "uniform"; 2343 } 2344 2345 /** @ignore */ 2346 makeFct = function (which) { 2347 return function (t, suspendedUpdate) { 2348 var s, 2349 c, 2350 // control point at the beginning and at the end 2351 first, 2352 last, 2353 t1, 2354 t2, 2355 dt0, 2356 dt1, 2357 dt2, 2358 // dx, dy, 2359 len; 2360 2361 if (points.length < 2) { 2362 return NaN; 2363 } 2364 2365 if (!suspendedUpdate) { 2366 tau = _tau(); 2367 2368 // New point list p: [first, points ..., last] 2369 first = { 2370 X: function () { 2371 return 2 * points[0].X() - points[1].X(); 2372 }, 2373 Y: function () { 2374 return 2 * points[0].Y() - points[1].Y(); 2375 }, 2376 Dist: function (p) { 2377 var dx = this.X() - p.X(), 2378 dy = this.Y() - p.Y(); 2379 return Mat.hypot(dx, dy); 2380 } 2381 }; 2382 2383 last = { 2384 X: function () { 2385 return ( 2386 2 * points[points.length - 1].X() - 2387 points[points.length - 2].X() 2388 ); 2389 }, 2390 Y: function () { 2391 return ( 2392 2 * points[points.length - 1].Y() - 2393 points[points.length - 2].Y() 2394 ); 2395 }, 2396 Dist: function (p) { 2397 var dx = this.X() - p.X(), 2398 dy = this.Y() - p.Y(); 2399 return Mat.hypot(dx, dy); 2400 } 2401 }; 2402 2403 p = [first].concat(points, [last]); 2404 len = p.length; 2405 2406 coeffs[which] = []; 2407 2408 for (s = 0; s < len - 3; s++) { 2409 if (type === "centripetal") { 2410 // The order is important, since p[0].coords === undefined 2411 dt0 = p[s].Dist(p[s + 1]); 2412 dt1 = p[s + 2].Dist(p[s + 1]); 2413 dt2 = p[s + 3].Dist(p[s + 2]); 2414 2415 dt0 = Math.sqrt(dt0); 2416 dt1 = Math.sqrt(dt1); 2417 dt2 = Math.sqrt(dt2); 2418 2419 if (dt1 < Mat.eps) { 2420 dt1 = 1.0; 2421 } 2422 if (dt0 < Mat.eps) { 2423 dt0 = dt1; 2424 } 2425 if (dt2 < Mat.eps) { 2426 dt2 = dt1; 2427 } 2428 2429 t1 = 2430 (p[s + 1][which]() - p[s][which]()) / dt0 - 2431 (p[s + 2][which]() - p[s][which]()) / (dt1 + dt0) + 2432 (p[s + 2][which]() - p[s + 1][which]()) / dt1; 2433 2434 t2 = 2435 (p[s + 2][which]() - p[s + 1][which]()) / dt1 - 2436 (p[s + 3][which]() - p[s + 1][which]()) / (dt2 + dt1) + 2437 (p[s + 3][which]() - p[s + 2][which]()) / dt2; 2438 2439 t1 *= dt1; 2440 t2 *= dt1; 2441 2442 coeffs[which][s] = that._initCubicPoly( 2443 p[s + 1][which](), 2444 p[s + 2][which](), 2445 tau * t1, 2446 tau * t2 2447 ); 2448 } else { 2449 coeffs[which][s] = that._initCubicPoly( 2450 p[s + 1][which](), 2451 p[s + 2][which](), 2452 tau * (p[s + 2][which]() - p[s][which]()), 2453 tau * (p[s + 3][which]() - p[s + 1][which]()) 2454 ); 2455 } 2456 } 2457 } 2458 2459 if (isNaN(t)) { 2460 return NaN; 2461 } 2462 2463 len = points.length; 2464 // This is necessary for our advanced plotting algorithm: 2465 if (t <= 0.0) { 2466 return points[0][which](); 2467 } 2468 if (t >= len) { 2469 return points[len - 1][which](); 2470 } 2471 2472 s = Math.floor(t); 2473 if (s === t) { 2474 return points[s][which](); 2475 } 2476 2477 t -= s; 2478 c = coeffs[which][s]; 2479 if (c === undefined) { 2480 return NaN; 2481 } 2482 2483 return ((c[3] * t + c[2]) * t + c[1]) * t + c[0]; 2484 }; 2485 }; 2486 2487 return [ 2488 makeFct("X"), 2489 makeFct("Y"), 2490 0, 2491 function () { 2492 return points.length - 1; 2493 } 2494 ]; 2495 }, 2496 2497 /** 2498 * Computes the cubic Catmull-Rom spline curve through a given set of points. The curve 2499 * is uniformly parametrized. The curve is the cardinal spline curve for tau=0.5. 2500 * Two artificial control points at the beginning and the end are added. 2501 * @param {Array} points Array consisting of JXG.Points. 2502 * @param {String} type (Optional) parameter which allows to choose between "uniform" (default) and 2503 * "centripetal" parameterization. Thus the two possible values are "uniform" or "centripetal". 2504 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 2505 * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, and a function simply 2506 * returning the length of the points array minus three. 2507 * @memberof JXG.Math.Numerics 2508 */ 2509 CatmullRomSpline: function (points, type) { 2510 return this.CardinalSpline(points, 0.5, type); 2511 }, 2512 2513 /** 2514 * Computes the regression polynomial of a given degree through a given set of coordinates. 2515 * Returns the regression polynomial function. 2516 * @param {Number|function|Slider} degree number, function or slider. 2517 * Either 2518 * @param {Array} dataX Array containing either the x-coordinates of the data set or both coordinates in 2519 * an array of {@link JXG.Point}s or {@link JXG.Coords}. 2520 * In the latter case, the <tt>dataY</tt> parameter will be ignored. 2521 * @param {Array} dataY Array containing the y-coordinates of the data set, 2522 * @returns {function} A function of one parameter which returns the value of the regression polynomial of the given degree. 2523 * It possesses the method getTerm() which returns the string containing the function term of the polynomial. 2524 * The function returned will throw an exception, if the data set is malformed. 2525 * @memberof JXG.Math.Numerics 2526 */ 2527 regressionPolynomial: function (degree, dataX, dataY) { 2528 var coeffs, deg, dX, dY, inputType, fct, 2529 term = ""; 2530 2531 // Slider 2532 if (Type.isPoint(degree) && Type.isFunction(degree.Value)) { 2533 /** @ignore */ 2534 deg = function () { 2535 return degree.Value(); 2536 }; 2537 // function 2538 } else if (Type.isFunction(degree)) { 2539 deg = degree; 2540 // number 2541 } else if (Type.isNumber(degree)) { 2542 /** @ignore */ 2543 deg = function () { 2544 return degree; 2545 }; 2546 } else { 2547 throw new Error( 2548 "JSXGraph: Can't create regressionPolynomial from degree of type'" + 2549 typeof degree + 2550 "'." 2551 ); 2552 } 2553 2554 // Parameters degree, dataX, dataY 2555 if (arguments.length === 3 && Type.isArray(dataX) && Type.isArray(dataY)) { 2556 inputType = 0; 2557 // Parameters degree, point array 2558 } else if ( 2559 arguments.length === 2 && 2560 Type.isArray(dataX) && 2561 dataX.length > 0 && 2562 Type.isPoint(dataX[0]) 2563 ) { 2564 inputType = 1; 2565 } else if ( 2566 arguments.length === 2 && 2567 Type.isArray(dataX) && 2568 dataX.length > 0 && 2569 dataX[0].usrCoords && 2570 dataX[0].scrCoords 2571 ) { 2572 inputType = 2; 2573 } else { 2574 throw new Error("JSXGraph: Can't create regressionPolynomial. Wrong parameters."); 2575 } 2576 2577 /** @ignore */ 2578 fct = function (x, suspendedUpdate) { 2579 var i, j, 2580 M, MT, y, B, c, s, d, 2581 // input data 2582 len = dataX.length; 2583 2584 d = Math.floor(deg()); 2585 2586 if (!suspendedUpdate) { 2587 // point list as input 2588 if (inputType === 1) { 2589 dX = []; 2590 dY = []; 2591 2592 for (i = 0; i < len; i++) { 2593 dX[i] = dataX[i].X(); 2594 dY[i] = dataX[i].Y(); 2595 } 2596 } 2597 2598 if (inputType === 2) { 2599 dX = []; 2600 dY = []; 2601 2602 for (i = 0; i < len; i++) { 2603 dX[i] = dataX[i].usrCoords[1]; 2604 dY[i] = dataX[i].usrCoords[2]; 2605 } 2606 } 2607 2608 // check for functions 2609 if (inputType === 0) { 2610 dX = []; 2611 dY = []; 2612 2613 for (i = 0; i < len; i++) { 2614 if (Type.isFunction(dataX[i])) { 2615 dX.push(dataX[i]()); 2616 } else { 2617 dX.push(dataX[i]); 2618 } 2619 2620 if (Type.isFunction(dataY[i])) { 2621 dY.push(dataY[i]()); 2622 } else { 2623 dY.push(dataY[i]); 2624 } 2625 } 2626 } 2627 2628 M = []; 2629 for (j = 0; j < len; j++) { 2630 M.push([1]); 2631 } 2632 for (i = 1; i <= d; i++) { 2633 for (j = 0; j < len; j++) { 2634 M[j][i] = M[j][i - 1] * dX[j]; 2635 } 2636 } 2637 2638 y = dY; 2639 MT = Mat.transpose(M); 2640 B = Mat.matMatMult(MT, M); 2641 c = Mat.matVecMult(MT, y); 2642 coeffs = Mat.Numerics.Gauss(B, c); 2643 term = Mat.Numerics.generatePolynomialTerm(coeffs, d, "x", 3); 2644 } 2645 2646 // Horner's scheme to evaluate polynomial 2647 s = coeffs[d]; 2648 for (i = d - 1; i >= 0; i--) { 2649 s = s * x + coeffs[i]; 2650 } 2651 2652 return s; 2653 }; 2654 2655 /** @ignore */ 2656 fct.getTerm = function () { 2657 return term; 2658 }; 2659 2660 return fct; 2661 }, 2662 2663 /** 2664 * Computes the cubic Bezier curve through a given set of points. 2665 * @param {Array} points Array consisting of 3*k+1 {@link JXG.Points}. 2666 * The points at position k with k mod 3 = 0 are the data points, 2667 * points at position k with k mod 3 = 1 or 2 are the control points. 2668 * @returns {Array} An array consisting of two functions of one parameter t which return the 2669 * x resp. y coordinates of the Bezier curve in t, one zero value, and a third function accepting 2670 * no parameters and returning one third of the length of the points. 2671 * @memberof JXG.Math.Numerics 2672 */ 2673 bezier: function (points) { 2674 var len, 2675 flen, 2676 /** @ignore */ 2677 makeFct = function (which) { 2678 return function (t, suspendedUpdate) { 2679 var z = Math.floor(t) * 3, 2680 t0 = t % 1, 2681 t1 = 1 - t0; 2682 2683 if (!suspendedUpdate) { 2684 flen = 3 * Math.floor((points.length - 1) / 3); 2685 len = Math.floor(flen / 3); 2686 } 2687 2688 if (t < 0) { 2689 return points[0][which](); 2690 } 2691 2692 if (t >= len) { 2693 return points[flen][which](); 2694 } 2695 2696 if (isNaN(t)) { 2697 return NaN; 2698 } 2699 2700 return ( 2701 t1 * t1 * (t1 * points[z][which]() + 3 * t0 * points[z + 1][which]()) + 2702 (3 * t1 * points[z + 2][which]() + t0 * points[z + 3][which]()) * 2703 t0 * 2704 t0 2705 ); 2706 }; 2707 }; 2708 2709 return [ 2710 makeFct("X"), 2711 makeFct("Y"), 2712 0, 2713 function () { 2714 return Math.floor(points.length / 3); 2715 } 2716 ]; 2717 }, 2718 2719 /** 2720 * Computes the B-spline curve of order k (order = degree+1) through a given set of points. 2721 * @param {Array} points Array consisting of JXG.Points. 2722 * @param {Number} order Order of the B-spline curve. 2723 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 2724 * which return the x resp. y coordinates of the B-spline curve in t, a zero value, and a function simply 2725 * returning the length of the points array minus one. 2726 * @memberof JXG.Math.Numerics 2727 */ 2728 bspline: function (points, order) { 2729 var knots, 2730 _knotVector = function (n, k) { 2731 var j, 2732 kn = []; 2733 2734 for (j = 0; j < n + k + 1; j++) { 2735 if (j < k) { 2736 kn[j] = 0.0; 2737 } else if (j <= n) { 2738 kn[j] = j - k + 1; 2739 } else { 2740 kn[j] = n - k + 2; 2741 } 2742 } 2743 2744 return kn; 2745 }, 2746 _evalBasisFuncs = function (t, kn, k, s) { 2747 var i, 2748 j, 2749 a, 2750 b, 2751 den, 2752 N = []; 2753 2754 if (kn[s] <= t && t < kn[s + 1]) { 2755 N[s] = 1; 2756 } else { 2757 N[s] = 0; 2758 } 2759 2760 for (i = 2; i <= k; i++) { 2761 for (j = s - i + 1; j <= s; j++) { 2762 if (j <= s - i + 1 || j < 0) { 2763 a = 0.0; 2764 } else { 2765 a = N[j]; 2766 } 2767 2768 if (j >= s) { 2769 b = 0.0; 2770 } else { 2771 b = N[j + 1]; 2772 } 2773 2774 den = kn[j + i - 1] - kn[j]; 2775 2776 if (den === 0) { 2777 N[j] = 0; 2778 } else { 2779 N[j] = ((t - kn[j]) / den) * a; 2780 } 2781 2782 den = kn[j + i] - kn[j + 1]; 2783 2784 if (den !== 0) { 2785 N[j] += ((kn[j + i] - t) / den) * b; 2786 } 2787 } 2788 } 2789 return N; 2790 }, 2791 /** @ignore */ 2792 makeFct = function (which) { 2793 return function (t, suspendedUpdate) { 2794 var y, 2795 j, 2796 s, 2797 N = [], 2798 len = points.length, 2799 n = len - 1, 2800 k = order; 2801 2802 if (n <= 0) { 2803 return NaN; 2804 } 2805 2806 if (n + 2 <= k) { 2807 k = n + 1; 2808 } 2809 2810 if (t <= 0) { 2811 return points[0][which](); 2812 } 2813 2814 if (t >= n - k + 2) { 2815 return points[n][which](); 2816 } 2817 2818 s = Math.floor(t) + k - 1; 2819 knots = _knotVector(n, k); 2820 N = _evalBasisFuncs(t, knots, k, s); 2821 2822 y = 0.0; 2823 for (j = s - k + 1; j <= s; j++) { 2824 if (j < len && j >= 0) { 2825 y += points[j][which]() * N[j]; 2826 } 2827 } 2828 2829 return y; 2830 }; 2831 }; 2832 2833 return [ 2834 makeFct("X"), 2835 makeFct("Y"), 2836 0, 2837 function () { 2838 return points.length - 1; 2839 } 2840 ]; 2841 }, 2842 2843 /** 2844 * Numerical (symmetric) approximation of derivative. suspendUpdate is piped through, 2845 * see {@link JXG.Curve#updateCurve} 2846 * and {@link JXG.Curve#hasPoint}. 2847 * @param {function} f Function in one variable to be differentiated. 2848 * @param {object} [obj] Optional object that is treated as "this" in the function body. This is useful, if the function is a 2849 * method of an object and contains a reference to its parent object via "this". 2850 * @returns {function} Derivative function of a given function f. 2851 * @memberof JXG.Math.Numerics 2852 */ 2853 D: function (f, obj) { 2854 if (!Type.exists(obj)) { 2855 return function (x, suspendedUpdate) { 2856 var h = 0.00001, 2857 h2 = h * 2.0; 2858 2859 // Experiments with Richardsons rule 2860 /* 2861 var phi = (f(x + h, suspendedUpdate) - f(x - h, suspendedUpdate)) / h2; 2862 var phi2; 2863 h *= 0.5; 2864 h2 *= 0.5; 2865 phi2 = (f(x + h, suspendedUpdate) - f(x - h, suspendedUpdate)) / h2; 2866 2867 return phi2 + (phi2 - phi) / 3.0; 2868 */ 2869 return (f(x + h, suspendedUpdate) - f(x - h, suspendedUpdate)) / h2; 2870 }; 2871 } 2872 2873 return function (x, suspendedUpdate) { 2874 var h = 0.00001, 2875 h2 = h * 2.0; 2876 2877 return ( 2878 (f.apply(obj, [x + h, suspendedUpdate]) - 2879 f.apply(obj, [x - h, suspendedUpdate])) / 2880 h2 2881 ); 2882 }; 2883 }, 2884 2885 /** 2886 * Evaluate the function term for {@link JXG.Math.Numerics.riemann}. 2887 * @private 2888 * @param {Number} x function argument 2889 * @param {function} f JavaScript function returning a number 2890 * @param {String} type Name of the Riemann sum type, e.g. 'lower'. 2891 * @param {Number} delta Width of the bars in user coordinates 2892 * @returns {Number} Upper (delta > 0) or lower (delta < 0) value of the bar containing x of the Riemann sum. 2893 * @see JXG.Math.Numerics.riemann 2894 * @private 2895 * @memberof JXG.Math.Numerics 2896 */ 2897 _riemannValue: function (x, f, type, delta) { 2898 var y, y1, x1, delta1; 2899 2900 if (delta < 0) { 2901 // delta is negative if the lower function term is evaluated 2902 if (type !== "trapezoidal") { 2903 x = x + delta; 2904 } 2905 delta *= -1; 2906 if (type === "lower") { 2907 type = "upper"; 2908 } else if (type === "upper") { 2909 type = "lower"; 2910 } 2911 } 2912 2913 delta1 = delta * 0.01; // for 'lower' and 'upper' 2914 2915 if (type === "right") { 2916 y = f(x + delta); 2917 } else if (type === "middle") { 2918 y = f(x + delta * 0.5); 2919 } else if (type === "left" || type === "trapezoidal") { 2920 y = f(x); 2921 } else if (type === "lower") { 2922 y = f(x); 2923 2924 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 2925 y1 = f(x1); 2926 2927 if (y1 < y) { 2928 y = y1; 2929 } 2930 } 2931 2932 y1 = f(x + delta); 2933 if (y1 < y) { 2934 y = y1; 2935 } 2936 } else if (type === "upper") { 2937 y = f(x); 2938 2939 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 2940 y1 = f(x1); 2941 if (y1 > y) { 2942 y = y1; 2943 } 2944 } 2945 2946 y1 = f(x + delta); 2947 if (y1 > y) { 2948 y = y1; 2949 } 2950 } else if (type === "random") { 2951 y = f(x + delta * Math.random()); 2952 } else if (type === "simpson") { 2953 y = (f(x) + 4 * f(x + delta * 0.5) + f(x + delta)) / 6.0; 2954 } else { 2955 y = f(x); // default is lower 2956 } 2957 2958 return y; 2959 }, 2960 2961 /** 2962 * Helper function to create curve which displays Riemann sums. 2963 * Compute coordinates for the rectangles showing the Riemann sum. 2964 * <p> 2965 * In case of type "simpson" and "trapezoidal", the horizontal line approximating the function value 2966 * is replaced by a parabola or a secant. IN case of "simpson", 2967 * the parabola is approximated visually by a polygonal chain of fixed step width. 2968 * 2969 * @param {Function|Array} f Function or array of two functions. 2970 * If f is a function the integral of this function is approximated by the Riemann sum. 2971 * If f is an array consisting of two functions the area between the two functions is filled 2972 * by the Riemann sum bars. 2973 * @param {Number} n number of rectangles. 2974 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson', or 'trapezoidal'. 2975 * "simpson" is Simpson's 1/3 rule. 2976 * @param {Number} start Left border of the approximation interval 2977 * @param {Number} end Right border of the approximation interval 2978 * @returns {Array} An array of two arrays containing the x and y coordinates for the rectangles showing the Riemann sum. This 2979 * 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 2980 * rectangles. 2981 * @memberof JXG.Math.Numerics 2982 */ 2983 riemann: function (gf, n, type, start, end) { 2984 var i, delta, 2985 k, a, b, c, f0, f1, f2, xx, h, 2986 steps = 30, // Fixed step width for Simpson's rule 2987 xarr = [], 2988 yarr = [], 2989 x = start, 2990 sum = 0, 2991 y, f, g; 2992 2993 if (Type.isArray(gf)) { 2994 g = gf[0]; 2995 f = gf[1]; 2996 } else { 2997 f = gf; 2998 } 2999 3000 n = Math.floor(n); 3001 3002 if (n <= 0) { 3003 return [xarr, yarr, sum]; 3004 } 3005 3006 delta = (end - start) / n; 3007 3008 // "Upper" horizontal line defined by function 3009 for (i = 0; i < n; i++) { 3010 if (type === "simpson") { 3011 sum += this._riemannValue(x, f, type, delta) * delta; 3012 3013 h = delta * 0.5; 3014 f0 = f(x); 3015 f1 = f(x + h); 3016 f2 = f(x + 2 * h); 3017 3018 a = (f2 + f0 - 2 * f1) / (h * h) * 0.5; 3019 b = (f2 - f0) / (2 * h); 3020 c = f1; 3021 for (k = 0; k < steps; k++) { 3022 xx = k * delta / steps - h; 3023 xarr.push(x + xx + h); 3024 yarr.push(a * xx * xx + b * xx + c); 3025 } 3026 x += delta; 3027 y = f2; 3028 } else { 3029 y = this._riemannValue(x, f, type, delta); 3030 xarr.push(x); 3031 yarr.push(y); 3032 3033 x += delta; 3034 if (type === "trapezoidal") { 3035 f2 = f(x); 3036 sum += (y + f2) * 0.5 * delta; 3037 y = f2; 3038 } else { 3039 sum += y * delta; 3040 } 3041 3042 xarr.push(x); 3043 yarr.push(y); 3044 } 3045 xarr.push(x); 3046 yarr.push(y); 3047 } 3048 3049 // "Lower" horizontal line 3050 // Go backwards 3051 for (i = 0; i < n; i++) { 3052 if (type === "simpson" && g) { 3053 sum -= this._riemannValue(x, g, type, -delta) * delta; 3054 3055 h = delta * 0.5; 3056 f0 = g(x); 3057 f1 = g(x - h); 3058 f2 = g(x - 2 * h); 3059 3060 a = (f2 + f0 - 2 * f1) / (h * h) * 0.5; 3061 b = (f2 - f0) / (2 * h); 3062 c = f1; 3063 for (k = 0; k < steps; k++) { 3064 xx = k * delta / steps - h; 3065 xarr.push(x - xx - h); 3066 yarr.push(a * xx * xx + b * xx + c); 3067 } 3068 x -= delta; 3069 y = f2; 3070 } else { 3071 if (g) { 3072 y = this._riemannValue(x, g, type, -delta); 3073 } else { 3074 y = 0.0; 3075 } 3076 xarr.push(x); 3077 yarr.push(y); 3078 3079 x -= delta; 3080 if (g) { 3081 if (type === "trapezoidal") { 3082 f2 = g(x); 3083 sum -= (y + f2) * 0.5 * delta; 3084 y = f2; 3085 } else { 3086 sum -= y * delta; 3087 } 3088 } 3089 } 3090 xarr.push(x); 3091 yarr.push(y); 3092 3093 // Draw the vertical lines 3094 xarr.push(x); 3095 yarr.push(f(x)); 3096 } 3097 3098 return [xarr, yarr, sum]; 3099 }, 3100 3101 /** 3102 * Approximate the integral by Riemann sums. 3103 * Compute the area described by the riemann sum rectangles. 3104 * 3105 * If there is an element of type {@link Riemannsum}, then it is more efficient 3106 * to use the method JXG.Curve.Value() of this element instead. 3107 * 3108 * @param {Function_Array} f Function or array of two functions. 3109 * If f is a function the integral of this function is approximated by the Riemann sum. 3110 * If f is an array consisting of two functions the area between the two functions is approximated 3111 * by the Riemann sum. 3112 * @param {Number} n number of rectangles. 3113 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson' or 'trapezoidal'. 3114 * 3115 * @param {Number} start Left border of the approximation interval 3116 * @param {Number} end Right border of the approximation interval 3117 * @returns {Number} The sum of the areas of the rectangles. 3118 * @memberof JXG.Math.Numerics 3119 */ 3120 riemannsum: function (f, n, type, start, end) { 3121 JXG.deprecated("Numerics.riemannsum()", "Numerics.riemann()[2]"); 3122 return this.riemann(f, n, type, start, end)[2]; 3123 }, 3124 3125 /** 3126 * Solve initial value problems numerically using <i>explicit</i> Runge-Kutta methods. 3127 * See {@link https://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm. 3128 * @param {object|String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing 3129 * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure 3130 * <pre> 3131 * { 3132 * s: <Number>, 3133 * A: <matrix>, 3134 * b: <Array>, 3135 * c: <Array> 3136 * } 3137 * </pre> 3138 * which corresponds to the Butcher tableau structure 3139 * shown here: https://en.wikipedia.org/w/index.php?title=List_of_Runge%E2%80%93Kutta_methods&oldid=357796696 . 3140 * <i>Default</i> is 'euler'. 3141 * @param {Array} x0 Initial value vector. Even if the problem is one-dimensional, the initial value has to be given in an array. 3142 * @param {Array} I Interval on which to integrate. 3143 * @param {Number} N Number of integration intervals, i.e. there are <i>N+1</i> evaluation points. 3144 * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode 3145 * 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 3146 * vector <tt>x</tt>, and has to return a vector of the same length as <tt>x</tt> has. 3147 * @returns {Array} An array of vectors describing the solution of the ode on the given interval I. 3148 * @example 3149 * // A very simple autonomous system dx(t)/dt = x(t); 3150 * var f = function(t, x) { 3151 * return [x[0]]; 3152 * } 3153 * 3154 * // Solve it with initial value x(0) = 1 on the interval [0, 2] 3155 * // with 20 evaluation points. 3156 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 3157 * 3158 * // Prepare data for plotting the solution of the ode using a curve. 3159 * var dataX = []; 3160 * var dataY = []; 3161 * var h = 0.1; // (I[1] - I[0])/N = (2-0)/20 3162 * var i; 3163 * for(i=0; i<data.length; i++) { 3164 * dataX[i] = i*h; 3165 * dataY[i] = data[i][0]; 3166 * } 3167 * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'}); 3168 * </pre><div class="jxgbox" id="JXGd2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div> 3169 * <script type="text/javascript"> 3170 * var board = JXG.JSXGraph.initBoard('JXGd2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false}); 3171 * var f = function(t, x) { 3172 * // we have to copy the value. 3173 * // return x; would just return the reference. 3174 * return [x[0]]; 3175 * } 3176 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 3177 * var dataX = []; 3178 * var dataY = []; 3179 * var h = 0.1; 3180 * for(var i=0; i<data.length; i++) { 3181 * dataX[i] = i*h; 3182 * dataY[i] = data[i][0]; 3183 * } 3184 * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'}); 3185 * </script><pre> 3186 * @memberof JXG.Math.Numerics 3187 */ 3188 rungeKutta: function (butcher, x0, I, N, f) { 3189 var e, 3190 i, j, k, l, s, 3191 x = [], 3192 y = [], 3193 h = (I[1] - I[0]) / N, 3194 t = I[0], 3195 dim = x0.length, 3196 result = [], 3197 r = 0; 3198 3199 if (Type.isString(butcher)) { 3200 butcher = predefinedButcher[butcher] || predefinedButcher.euler; 3201 } 3202 s = butcher.s; 3203 3204 // Don't change x0, so copy it 3205 x = x0.slice(); 3206 3207 for (i = 0; i <= N; i++) { 3208 result[r] = x.slice(); 3209 3210 r++; 3211 k = []; 3212 3213 for (j = 0; j < s; j++) { 3214 // Init y = 0 3215 for (e = 0; e < dim; e++) { 3216 y[e] = 0.0; 3217 } 3218 3219 // Calculate linear combination of former k's and save it in y 3220 for (l = 0; l < j; l++) { 3221 for (e = 0; e < dim; e++) { 3222 y[e] += butcher.A[j][l] * h * k[l][e]; 3223 } 3224 } 3225 3226 // Add x(t) to y 3227 for (e = 0; e < dim; e++) { 3228 y[e] += x[e]; 3229 } 3230 3231 // Calculate new k and add it to the k matrix 3232 k.push(f(t + butcher.c[j] * h, y)); 3233 } 3234 3235 // Init y = 0 3236 for (e = 0; e < dim; e++) { 3237 y[e] = 0.0; 3238 } 3239 3240 for (l = 0; l < s; l++) { 3241 for (e = 0; e < dim; e++) { 3242 y[e] += butcher.b[l] * k[l][e]; 3243 } 3244 } 3245 3246 for (e = 0; e < dim; e++) { 3247 x[e] = x[e] + h * y[e]; 3248 } 3249 3250 t += h; 3251 } 3252 3253 return result; 3254 }, 3255 3256 /** 3257 * Maximum number of iterations in {@link JXG.Math.Numerics.fzero} and 3258 * {@link JXG.Math.Numerics.chandrupatla} 3259 * @type Number 3260 * @default 80 3261 * @memberof JXG.Math.Numerics 3262 */ 3263 maxIterationsRoot: 80, 3264 3265 /** 3266 * Maximum number of iterations in {@link JXG.Math.Numerics.fminbr} 3267 * @type Number 3268 * @default 500 3269 * @memberof JXG.Math.Numerics 3270 */ 3271 maxIterationsMinimize: 500, 3272 3273 /** 3274 * Given a number x_0, this function tries to find a second number x_1 such that 3275 * the function f has opposite signs at x_0 and x_1. 3276 * The return values have to be tested if the method succeeded. 3277 * 3278 * @param {Function} f Function, whose root is to be found 3279 * @param {Number} x0 Start value 3280 * @param {Object} [context] Parent object in case f is method of it 3281 * @returns {Array} [x_0, f(x_0), x_1, f(x_1)] in case that x_0 <= x_1 3282 * or [x_1, f(x_1), x_0, f(x_0)] in case that x_1 < x_0. 3283 * 3284 * @see JXG.Math.Numerics.fzero 3285 * @see JXG.Math.Numerics.chandrupatla 3286 * 3287 * @memberof JXG.Math.Numerics 3288 */ 3289 findBracket: function (f, x0, context) { 3290 var a, aa, fa, blist, b, fb, u, fu, i, len; 3291 3292 if (Type.isArray(x0)) { 3293 return x0; 3294 } 3295 3296 a = x0; 3297 fa = f.call(context, a); 3298 // nfev += 1; 3299 3300 // Try to get b, by trying several values related to a 3301 aa = a === 0 ? 1 : a; 3302 blist = [ 3303 a - 0.1 * aa, 3304 a + 0.1 * aa, 3305 a - 1, 3306 a + 1, 3307 a - 0.5 * aa, 3308 a + 0.5 * aa, 3309 a - 0.6 * aa, 3310 a + 0.6 * aa, 3311 a - 1 * aa, 3312 a + 1 * aa, 3313 a - 2 * aa, 3314 a + 2 * aa, 3315 a - 5 * aa, 3316 a + 5 * aa, 3317 a - 10 * aa, 3318 a + 10 * aa, 3319 a - 50 * aa, 3320 a + 50 * aa, 3321 a - 100 * aa, 3322 a + 100 * aa 3323 ]; 3324 len = blist.length; 3325 3326 for (i = 0; i < len; i++) { 3327 b = blist[i]; 3328 fb = f.call(context, b); 3329 // nfev += 1; 3330 3331 if (fa * fb <= 0) { 3332 break; 3333 } 3334 } 3335 if (b < a) { 3336 u = a; 3337 a = b; 3338 b = u; 3339 3340 fu = fa; 3341 fa = fb; 3342 fb = fu; 3343 } 3344 return [a, fa, b, fb]; 3345 }, 3346 3347 /** 3348 * 3349 * Find zero of an univariate function f. 3350 * @param {function} f Function, whose root is to be found 3351 * @param {Array|Number} x0 Start value or start interval enclosing the root. 3352 * If x0 is an interval [a,b], it is required that f(a)f(b) <= 0, otherwise the minimum of f in [a, b] will be returned. 3353 * If x0 is a number, the algorithms tries to enclose the root by an interval [a, b] containing x0 and the root and 3354 * f(a)f(b) <= 0. If this fails, the algorithm falls back to Newton's method. 3355 * @param {Object} [context] Parent object in case f is method of it 3356 * @returns {Number} the approximation of the root 3357 * Algorithm: 3358 * Brent's root finder from 3359 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 3360 * computations. M., Mir, 1980, p.180 of the Russian edition 3361 * https://www.netlib.org/c/brent.shar 3362 * 3363 * If x0 is an array containing lower and upper bound for the zero 3364 * algorithm 748 is applied. Otherwise, if x0 is a number, 3365 * the algorithm tries to bracket a zero of f starting from x0. 3366 * If this fails, we fall back to Newton's method. 3367 * 3368 * @see JXG.Math.Numerics.chandrupatla 3369 * @see JXG.Math.Numerics.root 3370 * @see JXG.Math.Numerics.findBracket 3371 * @see JXG.Math.Numerics.Newton 3372 * @see JXG.Math.Numerics.fminbr 3373 * @memberof JXG.Math.Numerics 3374 */ 3375 fzero: function (f, x0, context) { 3376 var a, b, c, 3377 fa, fb, fc, 3378 res, x00, 3379 prev_step, 3380 t1, t2, 3381 cb, 3382 tol_act, // Actual tolerance 3383 p, // Interpolation step is calculated in the form p/q; division 3384 q, // operations is delayed until the last moment 3385 new_step, // Step at this iteration 3386 eps = Mat.eps, 3387 maxiter = this.maxIterationsRoot, 3388 niter = 0; 3389 // nfev = 0; 3390 3391 if (Type.isArray(x0)) { 3392 if (x0.length < 2) { 3393 throw new Error( 3394 "JXG.Math.Numerics.fzero: length of array x0 has to be at least two." 3395 ); 3396 } 3397 3398 x00 = this.findDomain(f, x0, context); 3399 a = x00[0]; 3400 b = x00[1]; 3401 // a = x0[0]; 3402 // b = x0[1]; 3403 3404 fa = f.call(context, a); 3405 // nfev += 1; 3406 fb = f.call(context, b); 3407 // nfev += 1; 3408 } else { 3409 res = this.findBracket(f, x0, context); 3410 a = res[0]; 3411 fa = res[1]; 3412 b = res[2]; 3413 fb = res[3]; 3414 } 3415 3416 if (Math.abs(fa) <= eps) { 3417 return a; 3418 } 3419 if (Math.abs(fb) <= eps) { 3420 return b; 3421 } 3422 3423 if (fa * fb > 0) { 3424 // Bracketing not successful, fall back to Newton's method or to fminbr 3425 if (Type.isArray(x0)) { 3426 return this.fminbr(f, [a, b], context); 3427 } 3428 3429 return this.Newton(f, a, context); 3430 } 3431 3432 // OK, we have enclosed a zero of f. 3433 // Now we can start Brent's method 3434 c = a; 3435 fc = fa; 3436 3437 // Main iteration loop 3438 while (niter < maxiter) { 3439 // Distance from the last but one to the last approximation 3440 prev_step = b - a; 3441 3442 // Swap data for b to be the best approximation 3443 if (Math.abs(fc) < Math.abs(fb)) { 3444 a = b; 3445 b = c; 3446 c = a; 3447 3448 fa = fb; 3449 fb = fc; 3450 fc = fa; 3451 } 3452 tol_act = 2 * eps * Math.abs(b) + eps * 0.5; 3453 new_step = (c - b) * 0.5; 3454 3455 if (Math.abs(new_step) <= tol_act || Math.abs(fb) <= eps) { 3456 // Acceptable approx. is found 3457 return b; 3458 } 3459 3460 // Decide if the interpolation can be tried 3461 // If prev_step was large enough and was in true direction Interpolatiom may be tried 3462 if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) { 3463 cb = c - b; 3464 3465 // If we have only two distinct points linear interpolation can only be applied 3466 if (a === c) { 3467 t1 = fb / fa; 3468 p = cb * t1; 3469 q = 1.0 - t1; 3470 // Quadric inverse interpolation 3471 } else { 3472 q = fa / fc; 3473 t1 = fb / fc; 3474 t2 = fb / fa; 3475 3476 p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0)); 3477 q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0); 3478 } 3479 3480 // p was calculated with the opposite sign; make p positive 3481 if (p > 0) { 3482 q = -q; 3483 // and assign possible minus to q 3484 } else { 3485 p = -p; 3486 } 3487 3488 // If b+p/q falls in [b,c] and isn't too large it is accepted 3489 // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent 3490 if ( 3491 p < 0.75 * cb * q - Math.abs(tol_act * q) * 0.5 && 3492 p < Math.abs(prev_step * q * 0.5) 3493 ) { 3494 new_step = p / q; 3495 } 3496 } 3497 3498 // Adjust the step to be not less than tolerance 3499 if (Math.abs(new_step) < tol_act) { 3500 new_step = new_step > 0 ? tol_act : -tol_act; 3501 } 3502 3503 // Save the previous approx. 3504 a = b; 3505 fa = fb; 3506 b += new_step; 3507 fb = f.call(context, b); 3508 // Do step to a new approxim. 3509 // nfev += 1; 3510 3511 // Adjust c for it to have a sign opposite to that of b 3512 if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) { 3513 c = a; 3514 fc = fa; 3515 } 3516 niter++; 3517 } // End while 3518 3519 return b; 3520 }, 3521 3522 /** 3523 * Find zero of an univariate function f. 3524 * @param {function} f Function, whose root is to be found 3525 * @param {Array|Number} x0 Start value or start interval enclosing the root. 3526 * If x0 is an interval [a,b], it is required that f(a)f(b) <= 0, otherwise the minimum of f in [a, b] will be returned. 3527 * If x0 is a number, the algorithms tries to enclose the root by an interval [a, b] containing x0 and the root and 3528 * f(a)f(b) <= 0. If this fails, the algorithm falls back to Newton's method. 3529 * @param {Object} [context] Parent object in case f is method of it 3530 * @returns {Number} the approximation of the root 3531 * Algorithm: 3532 * Chandrupatla's method, see 3533 * Tirupathi R. Chandrupatla, 3534 * "A new hybrid quadratic/bisection algorithm for finding the zero of a nonlinear function without using derivatives", 3535 * Advances in Engineering Software, Volume 28, Issue 3, April 1997, Pages 145-149. 3536 * 3537 * If x0 is an array containing lower and upper bound for the zero 3538 * algorithm 748 is applied. Otherwise, if x0 is a number, 3539 * the algorithm tries to bracket a zero of f starting from x0. 3540 * If this fails, we fall back to Newton's method. 3541 * 3542 * @see JXG.Math.Numerics.root 3543 * @see JXG.Math.Numerics.fzero 3544 * @see JXG.Math.Numerics.findBracket 3545 * @see JXG.Math.Numerics.Newton 3546 * @see JXG.Math.Numerics.fminbr 3547 * @memberof JXG.Math.Numerics 3548 */ 3549 chandrupatla: function (f, x0, context) { 3550 var a, b, fa, fb, 3551 res, 3552 niter = 0, 3553 maxiter = this.maxIterationsRoot, 3554 rand = 1 + Math.random() * 0.001, 3555 t = 0.5 * rand, 3556 eps = Mat.eps, // 1.e-10, 3557 dlt = 0.00001, 3558 x1, x2, x3, x, 3559 f1, f2, f3, y, 3560 xm, fm, 3561 tol, tl, 3562 xi, ph, fl, fh, 3563 AL, A, B, C, D; 3564 3565 if (Type.isArray(x0)) { 3566 if (x0.length < 2) { 3567 throw new Error( 3568 "JXG.Math.Numerics.fzero: length of array x0 has to be at least two." 3569 ); 3570 } 3571 3572 a = x0[0]; 3573 fa = f.call(context, a); 3574 // nfev += 1; 3575 b = x0[1]; 3576 fb = f.call(context, b); 3577 // nfev += 1; 3578 } else { 3579 res = this.findBracket(f, x0, context); 3580 a = res[0]; 3581 fa = res[1]; 3582 b = res[2]; 3583 fb = res[3]; 3584 } 3585 3586 if (fa * fb > 0) { 3587 // Bracketing not successful, fall back to Newton's method or to fminbr 3588 if (Type.isArray(x0)) { 3589 return this.fminbr(f, [a, b], context); 3590 } 3591 3592 return this.Newton(f, a, context); 3593 } 3594 3595 x1 = a; 3596 x2 = b; 3597 f1 = fa; 3598 f2 = fb; 3599 do { 3600 x = x1 + t * (x2 - x1); 3601 y = f.call(context, x); 3602 3603 // Arrange 2-1-3: 2-1 interval, 1 middle, 3 discarded point 3604 if (Math.sign(y) === Math.sign(f1)) { 3605 x3 = x1; 3606 x1 = x; 3607 f3 = f1; 3608 f1 = y; 3609 } else { 3610 x3 = x2; 3611 x2 = x1; 3612 f3 = f2; 3613 f2 = f1; 3614 } 3615 x1 = x; 3616 f1 = y; 3617 3618 xm = x1; 3619 fm = f1; 3620 if (Math.abs(f2) < Math.abs(f1)) { 3621 xm = x2; 3622 fm = f2; 3623 } 3624 tol = 2 * eps * Math.abs(xm) + 0.5 * dlt; 3625 tl = tol / Math.abs(x2 - x1); 3626 if (tl > 0.5 || fm === 0) { 3627 break; 3628 } 3629 // If inverse quadratic interpolation holds, use it 3630 xi = (x1 - x2) / (x3 - x2); 3631 ph = (f1 - f2) / (f3 - f2); 3632 fl = 1 - Math.sqrt(1 - xi); 3633 fh = Math.sqrt(xi); 3634 if (fl < ph && ph < fh) { 3635 AL = (x3 - x1) / (x2 - x1); 3636 A = f1 / (f2 - f1); 3637 B = f3 / (f2 - f3); 3638 C = f1 / (f3 - f1); 3639 D = f2 / (f3 - f2); 3640 t = A * B + C * D * AL; 3641 } else { 3642 t = 0.5 * rand; 3643 } 3644 // Adjust t away from the interval boundary 3645 if (t < tl) { 3646 t = tl; 3647 } 3648 if (t > 1 - tl) { 3649 t = 1 - tl; 3650 } 3651 niter++; 3652 } while (niter <= maxiter); 3653 // console.log(niter); 3654 3655 return xm; 3656 }, 3657 3658 /** 3659 * Find a small enclosing interval of the domain of a function by 3660 * tightening the input interval x0. 3661 * <p> 3662 * This is a helper function which is used in {@link JXG.Math.Numerics.fminbr} 3663 * and {@link JXG.Math.Numerics.fzero} 3664 * to avoid search in an interval where the function is mostly undefined. 3665 * 3666 * @param {function} f 3667 * @param {Array} x0 Start interval 3668 * @param {Object} context Parent object in case f is method of it 3669 * @returns Array 3670 * 3671 * @example 3672 * var f = (x) => Math.sqrt(x); 3673 * console.log(JXG.Math.Numerics.findDomain(f, [-5, 5])); 3674 * 3675 * // Output: [ -0.00020428174445492973, 5 ] 3676 */ 3677 findDomain: function (f, x0, context) { 3678 var a, b, c, fc, 3679 x, 3680 gr = 1 - 1 / 1.61803398875, 3681 eps = 0.001, 3682 cnt, 3683 max_cnt = 20; 3684 3685 if (!Type.isArray(x0) || x0.length < 2) { 3686 throw new Error( 3687 "JXG.Math.Numerics.findDomain: length of array x0 has to be at least two." 3688 ); 3689 } 3690 3691 x = x0.slice(); 3692 a = x[0]; 3693 b = x[1]; 3694 fc = f.call(context, a); 3695 if (isNaN(fc)) { 3696 // Divide the interval with the golden ratio 3697 // and keep a such that f(a) = NaN 3698 cnt = 0; 3699 while (b - a > eps && cnt < max_cnt) { 3700 c = (b - a) * gr + a; 3701 fc = f.call(context, c); 3702 if (isNaN(fc)) { 3703 a = c; 3704 } else { 3705 b = c; 3706 } 3707 cnt++; 3708 } 3709 x[0] = a; 3710 } 3711 3712 a = x[0]; 3713 b = x[1]; 3714 fc = f.call(context, b); 3715 if (isNaN(fc)) { 3716 // Divide the interval with the golden ratio 3717 // and keep b such that f(b) = NaN 3718 cnt = 0; 3719 while (b - a > eps && cnt < max_cnt) { 3720 c = b - (b - a) * gr; 3721 fc = f.call(context, c); 3722 if (isNaN(fc)) { 3723 b = c; 3724 } else { 3725 a = c; 3726 } 3727 cnt++; 3728 } 3729 x[0] = b; 3730 } 3731 return x; 3732 }, 3733 3734 /** 3735 * 3736 * Find minimum of an univariate function f. 3737 * <p> 3738 * Algorithm: 3739 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 3740 * computations. M., Mir, 1980, p.180 of the Russian edition 3741 * 3742 * @param {function} f Function, whose minimum is to be found 3743 * @param {Array} x0 Start interval enclosing the minimum 3744 * @param {Object} [context] Parent object in case f is method of it 3745 * @returns {Number} the approximation of the minimum value position 3746 * @memberof JXG.Math.Numerics 3747 **/ 3748 fminbr: function (f, x0, context) { 3749 var a, b, x, v, w, 3750 fx, fv, fw, 3751 x00, 3752 range, middle_range, tol_act, new_step, 3753 p, q, t, ft, 3754 r = (3.0 - Math.sqrt(5.0)) * 0.5, // Golden section ratio 3755 tol = Mat.eps, 3756 sqrteps = Mat.eps, // Math.sqrt(Mat.eps), 3757 maxiter = this.maxIterationsMinimize, 3758 niter = 0; 3759 // nfev = 0; 3760 3761 if (!Type.isArray(x0) || x0.length < 2) { 3762 throw new Error( 3763 "JXG.Math.Numerics.fminbr: length of array x0 has to be at least two." 3764 ); 3765 } 3766 3767 x00 = this.findDomain(f, x0, context); 3768 a = x00[0]; 3769 b = x00[1]; 3770 v = a + r * (b - a); 3771 fv = f.call(context, v); 3772 3773 // First step - always gold section 3774 // nfev += 1; 3775 x = v; 3776 w = v; 3777 fx = fv; 3778 fw = fv; 3779 3780 while (niter < maxiter) { 3781 // Range over the interval in which we are looking for the minimum 3782 range = b - a; 3783 middle_range = (a + b) * 0.5; 3784 3785 // Actual tolerance 3786 tol_act = sqrteps * Math.abs(x) + tol / 3.0; 3787 3788 if (Math.abs(x - middle_range) + range * 0.5 <= 2.0 * tol_act) { 3789 // Acceptable approx. is found 3790 return x; 3791 } 3792 3793 // Obtain the golden section step 3794 new_step = r * (x < middle_range ? b - x : a - x); 3795 3796 // Decide if the interpolation can be tried. If x and w are distinct, interpolatiom may be tried 3797 if (Math.abs(x - w) >= tol_act) { 3798 // Interpolation step is calculated as p/q; 3799 // division operation is delayed until last moment 3800 t = (x - w) * (fx - fv); 3801 q = (x - v) * (fx - fw); 3802 p = (x - v) * q - (x - w) * t; 3803 q = 2 * (q - t); 3804 3805 if (q > 0) { 3806 p = -p; // q was calculated with the opposite sign; make q positive 3807 } else { 3808 q = -q; // // and assign possible minus to p 3809 } 3810 if ( 3811 Math.abs(p) < Math.abs(new_step * q) && // If x+p/q falls in [a,b] 3812 p > q * (a - x + 2 * tol_act) && // not too close to a and 3813 p < q * (b - x - 2 * tol_act) 3814 ) { 3815 // b, and isn't too large 3816 new_step = p / q; // it is accepted 3817 } 3818 // If p / q is too large then the 3819 // golden section procedure can 3820 // reduce [a,b] range to more 3821 // extent 3822 } 3823 3824 // Adjust the step to be not less than tolerance 3825 if (Math.abs(new_step) < tol_act) { 3826 if (new_step > 0) { 3827 new_step = tol_act; 3828 } else { 3829 new_step = -tol_act; 3830 } 3831 } 3832 3833 // Obtain the next approximation to min 3834 // and reduce the enveloping range 3835 3836 // Tentative point for the min 3837 t = x + new_step; 3838 ft = f.call(context, t); 3839 // nfev += 1; 3840 3841 // t is a better approximation 3842 if (ft <= fx) { 3843 // Reduce the range so that t would fall within it 3844 if (t < x) { 3845 b = x; 3846 } else { 3847 a = x; 3848 } 3849 3850 // Assign the best approx to x 3851 v = w; 3852 w = x; 3853 x = t; 3854 3855 fv = fw; 3856 fw = fx; 3857 fx = ft; 3858 // x remains the better approx 3859 } else { 3860 // Reduce the range enclosing x 3861 if (t < x) { 3862 a = t; 3863 } else { 3864 b = t; 3865 } 3866 3867 if (ft <= fw || w === x) { 3868 v = w; 3869 w = t; 3870 fv = fw; 3871 fw = ft; 3872 } else if (ft <= fv || v === x || v === w) { 3873 v = t; 3874 fv = ft; 3875 } 3876 } 3877 niter += 1; 3878 } 3879 3880 return x; 3881 }, 3882 3883 /** 3884 * GLOMIN seeks a global minimum of a function F(X) in an interval [A,B] 3885 * and is the adaption of the algorithm GLOMIN by Richard Brent. 3886 * 3887 * Here is the original documentation: 3888 * <pre> 3889 * 3890 * Discussion: 3891 * 3892 * This function assumes that F(X) is twice continuously differentiable over [A,B] 3893 * and that |F''(X)| <= M for all X in [A,B]. 3894 * 3895 * Licensing: 3896 * This code is distributed under the GNU LGPL license. 3897 * 3898 * Modified: 3899 * 3900 * 17 April 2008 3901 * 3902 * Author: 3903 * 3904 * Original FORTRAN77 version by Richard Brent. 3905 * C version by John Burkardt. 3906 * https://people.math.sc.edu/Burkardt/c_src/brent/brent.c 3907 * 3908 * Reference: 3909 * 3910 * Richard Brent, 3911 * Algorithms for Minimization Without Derivatives, 3912 * Dover, 2002, 3913 * ISBN: 0-486-41998-3, 3914 * LC: QA402.5.B74. 3915 * 3916 * Parameters: 3917 * 3918 * Input, double A, B, the endpoints of the interval. 3919 * It must be the case that A < B. 3920 * 3921 * Input, double C, an initial guess for the global 3922 * minimizer. If no good guess is known, C = A or B is acceptable. 3923 * 3924 * Input, double M, the bound on the second derivative. 3925 * 3926 * Input, double MACHEP, an estimate for the relative machine 3927 * precision. 3928 * 3929 * Input, double E, a positive tolerance, a bound for the 3930 * absolute error in the evaluation of F(X) for any X in [A,B]. 3931 * 3932 * Input, double T, a positive error tolerance. 3933 * 3934 * Input, double F (double x ), a user-supplied 3935 * function whose global minimum is being sought. 3936 * 3937 * Output, double *X, the estimated value of the abscissa 3938 * for which F attains its global minimum value in [A,B]. 3939 * 3940 * Output, double GLOMIN, the value F(X). 3941 * </pre> 3942 * 3943 * In JSXGraph, some parameters of the original algorithm are set to fixed values: 3944 * <ul> 3945 * <li> M = 10000000.0 3946 * <li> C = A or B, depending if f(A) <= f(B) 3947 * <li> T = JXG.Math.eps 3948 * <li> E = JXG.Math.eps * JXG.Math.eps 3949 * <li> MACHEP = JXG.Math.eps * JXG.Math.eps * JXG.Math.eps 3950 * </ul> 3951 * @param {function} f Function, whose global minimum is to be found 3952 * @param {Array} x0 Array of length 2 determining the interval [A, B] for which the global minimum is to be found 3953 * @returns {Array} [x, y] x is the position of the global minimum and y = f(x). 3954 */ 3955 glomin: function (f, x0) { 3956 var a0, a2, a3, d0, d1, d2, h, 3957 k, m2, 3958 p, q, qs, 3959 r, s, sc, 3960 y, y0, y1, y2, y3, yb, 3961 z0, z1, z2, 3962 a, b, c, x, 3963 m = 10000000.0, 3964 t = Mat.eps, // * Mat.eps, 3965 e = Mat.eps * Mat.eps, 3966 machep = Mat.eps * Mat.eps * Mat.eps; 3967 3968 a = x0[0]; 3969 b = x0[1]; 3970 c = (f(a) < f(b)) ? a : b; 3971 3972 a0 = b; 3973 x = a0; 3974 a2 = a; 3975 y0 = f(b); 3976 yb = y0; 3977 y2 = f(a); 3978 y = y2; 3979 3980 if (y0 < y) { 3981 y = y0; 3982 } else { 3983 x = a; 3984 } 3985 3986 if (m <= 0.0 || b <= a) { 3987 return y; 3988 } 3989 3990 m2 = 0.5 * (1.0 + 16.0 * machep) * m; 3991 3992 if (c <= a || b <= c) { 3993 sc = 0.5 * (a + b); 3994 } else { 3995 sc = c; 3996 } 3997 3998 y1 = f(sc); 3999 k = 3; 4000 d0 = a2 - sc; 4001 h = 9.0 / 11.0; 4002 4003 if (y1 < y) { 4004 x = sc; 4005 y = y1; 4006 } 4007 4008 for (; ;) { 4009 d1 = a2 - a0; 4010 d2 = sc - a0; 4011 z2 = b - a2; 4012 z0 = y2 - y1; 4013 z1 = y2 - y0; 4014 r = d1 * d1 * z0 - d0 * d0 * z1; 4015 p = r; 4016 qs = 2.0 * (d0 * z1 - d1 * z0); 4017 q = qs; 4018 4019 if (k < 1000000 || y2 <= y) { 4020 for (; ;) { 4021 if (q * (r * (yb - y2) + z2 * q * ((y2 - y) + t)) < 4022 z2 * m2 * r * (z2 * q - r)) { 4023 4024 a3 = a2 + r / q; 4025 y3 = f(a3); 4026 4027 if (y3 < y) { 4028 x = a3; 4029 y = y3; 4030 } 4031 } 4032 k = ((1611 * k) % 1048576); 4033 q = 1.0; 4034 r = (b - a) * 0.00001 * k; 4035 4036 if (z2 <= r) { 4037 break; 4038 } 4039 } 4040 } else { 4041 k = ((1611 * k) % 1048576); 4042 q = 1.0; 4043 r = (b - a) * 0.00001 * k; 4044 4045 while (r < z2) { 4046 if (q * (r * (yb - y2) + z2 * q * ((y2 - y) + t)) < 4047 z2 * m2 * r * (z2 * q - r)) { 4048 4049 a3 = a2 + r / q; 4050 y3 = f(a3); 4051 4052 if (y3 < y) { 4053 x = a3; 4054 y = y3; 4055 } 4056 } 4057 k = ((1611 * k) % 1048576); 4058 q = 1.0; 4059 r = (b - a) * 0.00001 * k; 4060 } 4061 } 4062 4063 r = m2 * d0 * d1 * d2; 4064 s = Math.sqrt(((y2 - y) + t) / m2); 4065 h = 0.5 * (1.0 + h); 4066 p = h * (p + 2.0 * r * s); 4067 q = q + 0.5 * qs; 4068 r = - 0.5 * (d0 + (z0 + 2.01 * e) / (d0 * m2)); 4069 4070 if (r < s || d0 < 0.0) { 4071 r = a2 + s; 4072 } else { 4073 r = a2 + r; 4074 } 4075 4076 if (0.0 < p * q) { 4077 a3 = a2 + p / q; 4078 } else { 4079 a3 = r; 4080 } 4081 4082 for (; ;) { 4083 a3 = Math.max(a3, r); 4084 4085 if (b <= a3) { 4086 a3 = b; 4087 y3 = yb; 4088 } else { 4089 y3 = f(a3); 4090 } 4091 4092 if (y3 < y) { 4093 x = a3; 4094 y = y3; 4095 } 4096 4097 d0 = a3 - a2; 4098 4099 if (a3 <= r) { 4100 break; 4101 } 4102 4103 p = 2.0 * (y2 - y3) / (m * d0); 4104 4105 if ((1.0 + 9.0 * machep) * d0 <= Math.abs(p)) { 4106 break; 4107 } 4108 4109 if (0.5 * m2 * (d0 * d0 + p * p) <= (y2 - y) + (y3 - y) + 2.0 * t) { 4110 break; 4111 } 4112 a3 = 0.5 * (a2 + a3); 4113 h = 0.9 * h; 4114 } 4115 4116 if (b <= a3) { 4117 break; 4118 } 4119 4120 a0 = sc; 4121 sc = a2; 4122 a2 = a3; 4123 y0 = y1; 4124 y1 = y2; 4125 y2 = y3; 4126 } 4127 4128 return [x, y]; 4129 }, 4130 4131 /** 4132 * Determine all roots of a polynomial with real or complex coefficients by using the 4133 * iterative method attributed to Weierstrass, Durand, Kerner, Aberth, and Ehrlich. In particular, 4134 * the iteration method with cubic convergence is used that is usually attributed to Ehrlich-Aberth. 4135 * <p> 4136 * The returned roots are sorted with respect to their real values. 4137 * <p> This method makes use of the JSXGraph classes {@link JXG.Complex} and {@link JXG.C} to handle 4138 * complex numbers. 4139 * 4140 * @param {Array} a Array of coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2... 4141 * The coefficients are of type Number or JXG.Complex. 4142 * @param {Number} [deg] Optional degree of the polynomial. Otherwise all entries are taken, with 4143 * leading zeros removed. 4144 * @param {Number} [tol=Number.EPSILON] Approximation tolerance 4145 * @param {Number} [max_it=30] Maximum number of iterations 4146 * @param {Array} [initial_values=null] Array of initial values for the roots. If not given, 4147 * starting values are determined by the method of Ozawa. 4148 * @returns {Array} Array of complex numbers (of JXG.Complex) approximating the roots of the polynomial. 4149 * @memberof JXG.Math.Numerics 4150 * @see JXG.Complex 4151 * @see JXG.C 4152 * 4153 * @example 4154 * // Polynomial p(z) = -1 + 1z^2 4155 * var i, roots, 4156 * p = [-1, 0, 1]; 4157 * 4158 * roots = JXG.Math.Numerics.polzeros(p); 4159 * for (i = 0; i < roots.length; i++) { 4160 * console.log(i, roots[i].toString()); 4161 * } 4162 * // Output: 4163 * 0 -1 + -3.308722450212111e-24i 4164 * 1 1 + 0i 4165 * 4166 * @example 4167 * // Polynomial p(z) = -1 + 3z - 9z^2 + z^3 - 8z^6 + 9z^7 - 9z^8 + z^9 4168 * var i, roots, 4169 * p = [-1, 3, -9, 1, 0, 0, -8, 9, -9, 1]; 4170 * 4171 * roots = JXG.Math.Numerics.polzeros(p); 4172 * for (i = 0; i < roots.length; i++) { 4173 * console.log(i, roots[i].toString()); 4174 * } 4175 * // Output: 4176 * 0 -0.7424155888401961 + 0.4950476539211721i 4177 * 1 -0.7424155888401961 + -0.4950476539211721i 4178 * 2 0.16674869833354108 + 0.2980502714610669i 4179 * 3 0.16674869833354108 + -0.29805027146106694i 4180 * 4 0.21429002063640837 + 1.0682775088132996i 4181 * 5 0.21429002063640842 + -1.0682775088132999i 4182 * 6 0.861375497926218 + -0.6259177003583295i 4183 * 7 0.8613754979262181 + 0.6259177003583295i 4184 * 8 8.000002743888055 + -1.8367099231598242e-40i 4185 * 4186 */ 4187 polzeros: function (coeffs, deg, tol, max_it, initial_values) { 4188 var i, le, off, it, 4189 debug = false, 4190 cc = [], 4191 obvious = [], 4192 roots = [], 4193 4194 /** 4195 * Horner method to evaluate polynomial or the derivative thereof for complex numbers, 4196 * i.e. coefficients and variable are complex. 4197 * @function 4198 * @param {Array} a Array of complex coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2... 4199 * @param {JXG.Complex} z Value for which the polynomial will be evaluated. 4200 * @param {Boolean} [derivative=false] If true the derivative will be evaluated. 4201 * @ignore 4202 */ 4203 hornerComplex = function (a, z, derivative) { 4204 var i, s, 4205 n = a.length - 1; 4206 4207 derivative = derivative || false; 4208 if (derivative) { 4209 // s = n * a_n 4210 s = JXG.C.mult(n, a[n]); 4211 for (i = n - 1; i > 0; i--) { 4212 // s = s * z + i * a_i 4213 s.mult(z); 4214 s.add(JXG.C.mult(a[i], i)); 4215 } 4216 } else { 4217 // s = a_n 4218 s = JXG.C.copy(a[n]); 4219 for (i = n - 1; i >= 0; i--) { 4220 // s = s * z + a_i 4221 s.mult(z); 4222 s.add(a[i]); 4223 } 4224 } 4225 return s; 4226 }, 4227 4228 /** 4229 * Horner method to evaluate reciprocal polynomial or the derivative thereof for complex numbers, 4230 * i.e. coefficients and variable are complex. 4231 * @function 4232 * @param {Array} a Array of complex coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2... 4233 * @param {JXG.Complex} z Value for which the reciprocal polynomial will be evaluated. 4234 * @param {Boolean} [derivative=false] If true the derivative will be evaluated. 4235 * @ignore 4236 */ 4237 hornerRec = function (a, x, derivative) { 4238 var i, s, 4239 n = a.length - 1; 4240 4241 derivative = derivative || false; 4242 if (derivative) { 4243 // s = n * a_0 4244 s = JXG.C.mult(n, a[0]); 4245 for (i = n - 1; i > 0; i--) { 4246 // s = s * x + i * a_{n-i} 4247 s.mult(x); 4248 s.add(JXG.C.mult(a[n - i], i)); 4249 } 4250 } else { 4251 // s = a_0 4252 s = JXG.C.copy(a[0]); 4253 for (i = n - 1; i >= 0; i--) { 4254 // s = s * x + a_{n-i} 4255 s.mult(x); 4256 s.add(a[n - i]); 4257 } 4258 } 4259 return s; 4260 }, 4261 4262 /** 4263 * Horner method to evaluate real polynomial at a real value. 4264 * @function 4265 * @param {Array} a Array of real coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2... 4266 * @param {Number} z Value for which the polynomial will be evaluated. 4267 * @ignore 4268 */ 4269 horner = function (a, x) { 4270 var i, s, 4271 n = a.length - 1; 4272 4273 s = a[n]; 4274 for (i = n - 1; i >= 0; i--) { 4275 s = s * x + a[i]; 4276 } 4277 return s; 4278 }, 4279 4280 /** 4281 * Determine start values for the Aberth iteration, see 4282 * Ozawa, "An experimental study of the starting values 4283 * of the Durand-Kerner-Aberth iteration" (1995). 4284 * 4285 * @function 4286 * @param {Array} a Array of complex coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2... 4287 * @returns {Array} Array Initial values for the roots. 4288 * @ignore 4289 */ 4290 initial_guess = function (a) { 4291 var i, r, 4292 n = a.length - 1, // degree 4293 alpha1 = Math.PI * 2 / n, 4294 alpha0 = Math.PI / n * 0.5, 4295 b, z, 4296 init = []; 4297 4298 4299 // From Ozawa, "An experimental study of the starting values 4300 // of the Durand-Kerner-Aberth iteration" (1995) 4301 4302 // b is the arithmetic mean of the roots. 4303 // With is Vieta's formula <https://en.wikipedia.org/wiki/Vieta%27s_formulas> 4304 // b = -a_{n-1} / (n * a_n) 4305 b = JXG.C.mult(-1, a[n - 1]); 4306 b.div(JXG.C.mult(n, a[n])); 4307 4308 // r is the geometric mean of the deviations |b - root_i|. 4309 // Using 4310 // p(z) = a_n prod(z - root_i) 4311 // and therefore 4312 // |p(b)| = |a_n| prod(|b - root_i|) 4313 // we arrive at: 4314 // r = |p(b)/a_n|^(1/n) 4315 z = JXG.C.div(hornerComplex(a, b), a[n]); 4316 r = Math.pow(JXG.C.abs(z), 1 / n); 4317 if (r === 0) { r = 1; } 4318 4319 for (i = 0; i < n; i++) { 4320 a = new JXG.Complex(r * Math.cos(alpha1 * i + alpha0), r * Math.sin(alpha1 * i + alpha0)); 4321 init[i] = JXG.C.add(b, a); 4322 } 4323 4324 return init; 4325 }, 4326 4327 /** 4328 * Ehrlich-Aberth iteration. The stopping criterion is from 4329 * D.A. Bini, "Numerical computation of polynomial zeros 4330 * by means of Aberths's method", Numerical Algorithms (1996). 4331 * 4332 * @function 4333 * @param {Array} a Array of complex coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2... 4334 * @param {Number} mu Machine precision 4335 * @param {Number} max_it Maximum number of iterations 4336 * @param {Array} z Initial guess for the roots. Will be changed in place. 4337 * @returns {Number} Number of iterations 4338 * @ignore 4339 */ 4340 aberthIteration = function (cc, mu, max_it, z) { 4341 var k, i, j, 4342 done = [], 4343 cr = [], 4344 gamma, x, 4345 done_sum = 0, 4346 num, denom, s, pp, 4347 n = z.length; 4348 4349 for (i = 0; i < n; i++) { 4350 done.push(false); 4351 } 4352 for (i = 0; i < cc.length; i++) { 4353 cr.push(JXG.C.abs(cc[i]) * (4 * i + 1)); 4354 } 4355 for (k = 0; k < max_it && done_sum < n; k++) { 4356 for (i = 0; i < n; i++) { 4357 if (done[i]) { 4358 continue; 4359 } 4360 num = hornerComplex(cc, z[i]); 4361 x = JXG.C.abs(z[i]); 4362 4363 // Stopping criterion by D.A. Bini 4364 // "Numerical computation of polynomial zeros 4365 // by means of Aberths's method", Numerical Algorithms (1996). 4366 // 4367 if (JXG.C.abs(num) < mu * horner(cr, x)) { 4368 done[i] = true; 4369 done_sum++; 4370 if (done_sum === n) { 4371 break; 4372 } 4373 continue; 4374 } 4375 4376 // num = P(z_i) / P'(z_i) 4377 if (x > 1) { 4378 gamma = JXG.C.div(1, z[i]); 4379 pp = hornerRec(cc, gamma, true); 4380 pp.div(hornerRec(cc, gamma)); 4381 pp.mult(gamma); 4382 num = JXG.C.sub(n, pp); 4383 num = JXG.C.div(z[i], num); 4384 } else { 4385 num.div(hornerComplex(cc, z[i], true)); 4386 } 4387 4388 // denom = sum_{i\neq j} 1 / (z_i - z_j) 4389 denom = new JXG.Complex(0); 4390 for (j = 0; j < n; j++) { 4391 if (j === i) { 4392 continue; 4393 } 4394 s = JXG.C.sub(z[i], z[j]); 4395 s = JXG.C.div(1, s); 4396 denom.add(s); 4397 } 4398 4399 // num = num / 1 - num * sum_{i\neq j} 1 / (z_i - z_j) 4400 denom.mult(num); 4401 denom = JXG.C.sub(1, denom); 4402 num.div(denom); 4403 // z_i = z_i - num 4404 z[i].sub(num); 4405 } 4406 } 4407 4408 return k; 4409 }; 4410 4411 4412 tol = tol || Number.EPSILON; 4413 max_it = max_it || 30; 4414 4415 le = coeffs.length; 4416 if (JXG.isNumber(deg) && deg >= 0 && deg < le - 1) { 4417 le = deg + 1; 4418 } 4419 4420 // Convert coefficient array to complex numbers 4421 for (i = 0; i < le; i++) { 4422 cc.push(new JXG.Complex(coeffs[i])); 4423 } 4424 4425 // Search for (multiple) roots at x=0 4426 for (i = 0; i < le; i++) { 4427 if (cc[i].real !== 0 || cc[i].imaginary !== 0) { 4428 off = i; 4429 break; 4430 } 4431 } 4432 4433 // Deflate root x=0, store roots at x=0 in obvious 4434 for (i = 0; i < off; i++) { 4435 obvious.push(new JXG.Complex(0)); 4436 } 4437 cc = cc.slice(off); 4438 le = cc.length; 4439 4440 // Remove leading zeros from the coefficient array 4441 for (i = le - 1; i >= 0; i--) { 4442 if (cc[i].real !== 0 || cc[i].imaginary !== 0) { 4443 break; 4444 } 4445 cc.pop(); 4446 } 4447 le = cc.length; 4448 if (le === 0) { 4449 return []; 4450 } 4451 4452 // From now on we can assume that the 4453 // constant coefficient and the leading coefficient 4454 // are not zero. 4455 if (initial_values) { 4456 for (i = 0; i < le - 1; i++) { 4457 roots.push(new JXG.Complex(initial_values[i])); 4458 } 4459 } else { 4460 roots = initial_guess(cc); 4461 } 4462 it = aberthIteration(cc, tol, max_it, roots); 4463 4464 // Append the roots at x=0 4465 roots = obvious.concat(roots); 4466 4467 if (debug) { 4468 console.log("Iterations:", it); 4469 console.log('Roots:'); 4470 for (i = 0; i < roots.length; i++) { 4471 console.log(i, roots[i].toString(), JXG.C.abs(hornerComplex(cc, roots[i]))); 4472 } 4473 } 4474 4475 // Sort roots according to their real part 4476 roots.sort(function (a, b) { 4477 if (a.real < b.real) { 4478 return -1; 4479 } 4480 if (a.real > b.real) { 4481 return 1; 4482 } 4483 return 0; 4484 }); 4485 4486 return roots; 4487 }, 4488 4489 /** 4490 * Implements the Ramer-Douglas-Peucker algorithm. 4491 * It discards points which are not necessary from the polygonal line defined by the point array 4492 * pts. The computation is done in screen coordinates. 4493 * Average runtime is O(nlog(n)), worst case runtime is O(n^2), where n is the number of points. 4494 * @param {Array} pts Array of {@link JXG.Coords} 4495 * @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>. 4496 * @returns {Array} An array containing points which represent an apparently identical curve as the points of pts do, but contains fewer points. 4497 * @memberof JXG.Math.Numerics 4498 */ 4499 RamerDouglasPeucker: function (pts, eps) { 4500 var allPts = [], 4501 newPts = [], 4502 i, k, len, 4503 endless = true, 4504 4505 /** 4506 * findSplit() is a subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}. 4507 * It searches for the point between index i and j which 4508 * has the largest distance from the line between the points i and j. 4509 * @param {Array} pts Array of {@link JXG.Coords} 4510 * @param {Number} i Index of a point in pts 4511 * @param {Number} j Index of a point in pts 4512 * @ignore 4513 * @private 4514 */ 4515 findSplit = function (pts, i, j) { 4516 var d, k, ci, cj, ck, 4517 x0, y0, x1, y1, 4518 den, lbda, 4519 eps = Mat.eps * Mat.eps, 4520 huge = 10000, 4521 dist = 0, 4522 f = i; 4523 4524 if (j - i < 2) { 4525 return [-1.0, 0]; 4526 } 4527 4528 ci = pts[i].scrCoords; 4529 cj = pts[j].scrCoords; 4530 4531 if (isNaN(ci[1]) || isNaN(ci[2])) { 4532 return [NaN, i]; 4533 } 4534 if (isNaN(cj[1]) || isNaN(cj[2])) { 4535 return [NaN, j]; 4536 } 4537 4538 for (k = i + 1; k < j; k++) { 4539 ck = pts[k].scrCoords; 4540 if (isNaN(ck[1]) || isNaN(ck[2])) { 4541 return [NaN, k]; 4542 } 4543 4544 x0 = ck[1] - ci[1]; 4545 y0 = ck[2] - ci[2]; 4546 x1 = cj[1] - ci[1]; 4547 y1 = cj[2] - ci[2]; 4548 x0 = x0 === Infinity ? huge : x0; 4549 y0 = y0 === Infinity ? huge : y0; 4550 x1 = x1 === Infinity ? huge : x1; 4551 y1 = y1 === Infinity ? huge : y1; 4552 x0 = x0 === -Infinity ? -huge : x0; 4553 y0 = y0 === -Infinity ? -huge : y0; 4554 x1 = x1 === -Infinity ? -huge : x1; 4555 y1 = y1 === -Infinity ? -huge : y1; 4556 den = x1 * x1 + y1 * y1; 4557 4558 if (den > eps) { 4559 lbda = (x0 * x1 + y0 * y1) / den; 4560 4561 if (lbda < 0.0) { 4562 lbda = 0.0; 4563 } else if (lbda > 1.0) { 4564 lbda = 1.0; 4565 } 4566 4567 x0 = x0 - lbda * x1; 4568 y0 = y0 - lbda * y1; 4569 d = x0 * x0 + y0 * y0; 4570 } else { 4571 lbda = 0.0; 4572 d = x0 * x0 + y0 * y0; 4573 } 4574 4575 if (d > dist) { 4576 dist = d; 4577 f = k; 4578 } 4579 } 4580 return [Math.sqrt(dist), f]; 4581 }, 4582 /** 4583 * RDP() is a private subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}. 4584 * It runs recursively through the point set and searches the 4585 * point which has the largest distance from the line between the first point and 4586 * the last point. If the distance from the line is greater than eps, this point is 4587 * included in our new point set otherwise it is discarded. 4588 * If it is taken, we recursively apply the subroutine to the point set before 4589 * and after the chosen point. 4590 * @param {Array} pts Array of {@link JXG.Coords} 4591 * @param {Number} i Index of an element of pts 4592 * @param {Number} j Index of an element of pts 4593 * @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>. 4594 * @param {Array} newPts Array of {@link JXG.Coords} 4595 * @ignore 4596 * @private 4597 */ 4598 RDP = function (pts, i, j, eps, newPts) { 4599 var result = findSplit(pts, i, j), 4600 k = result[1]; 4601 4602 if (isNaN(result[0])) { 4603 RDP(pts, i, k - 1, eps, newPts); 4604 newPts.push(pts[k]); 4605 do { 4606 ++k; 4607 } while (k <= j && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])); 4608 if (k <= j) { 4609 newPts.push(pts[k]); 4610 } 4611 RDP(pts, k + 1, j, eps, newPts); 4612 } else if (result[0] > eps) { 4613 RDP(pts, i, k, eps, newPts); 4614 RDP(pts, k, j, eps, newPts); 4615 } else { 4616 newPts.push(pts[j]); 4617 } 4618 }; 4619 4620 len = pts.length; 4621 4622 i = 0; 4623 while (endless) { 4624 // Search for the next point without NaN coordinates 4625 while (i < len && isNaN(pts[i].scrCoords[1] + pts[i].scrCoords[2])) { 4626 i += 1; 4627 } 4628 // Search for the next position of a NaN point 4629 k = i + 1; 4630 while (k < len && !isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])) { 4631 k += 1; 4632 } 4633 k--; 4634 4635 // Only proceed if something is left 4636 if (i < len && k > i) { 4637 newPts = []; 4638 newPts[0] = pts[i]; 4639 RDP(pts, i, k, eps, newPts); 4640 allPts = allPts.concat(newPts); 4641 } 4642 if (i >= len) { 4643 break; 4644 } 4645 // Push the NaN point 4646 if (k < len - 1) { 4647 allPts.push(pts[k + 1]); 4648 } 4649 i = k + 1; 4650 } 4651 4652 return allPts; 4653 }, 4654 4655 /** 4656 * Old name for the implementation of the Ramer-Douglas-Peucker algorithm. 4657 * @deprecated Use {@link JXG.Math.Numerics.RamerDouglasPeucker} 4658 * @memberof JXG.Math.Numerics 4659 */ 4660 RamerDouglasPeuker: function (pts, eps) { 4661 JXG.deprecated("Numerics.RamerDouglasPeuker()", "Numerics.RamerDouglasPeucker()"); 4662 return this.RamerDouglasPeucker(pts, eps); 4663 }, 4664 4665 /** 4666 * Implements the Visvalingam-Whyatt algorithm. 4667 * See M. Visvalingam, J. D. Whyatt: 4668 * "Line generalisation by repeated elimination of the smallest area", C.I.S.R.G Discussion paper 10, July 1992 4669 * 4670 * The algorithm discards points which are not necessary from the polygonal line defined by the point array 4671 * pts (consisting of type JXG.Coords). 4672 * @param {Array} pts Array of {@link JXG.Coords} 4673 * @param {Number} numPoints Number of remaining intermediate points. The first and the last point of the original points will 4674 * be taken in any case. 4675 * @returns {Array} An array containing points which approximates the curve defined by pts. 4676 * @memberof JXG.Math.Numerics 4677 * 4678 * @example 4679 * var i, p = []; 4680 * for (i = 0; i < 5; ++i) { 4681 * p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6])); 4682 * } 4683 * 4684 * // Plot a cardinal spline curve 4685 * var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5); 4686 * var cu1 = board.create('curve', splineArr, {strokeColor: 'green'}); 4687 * 4688 * var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'}); 4689 * c.updateDataArray = function() { 4690 * var i, len, points; 4691 * 4692 * // Reduce number of intermediate points with Visvakingam-Whyatt to 6 4693 * points = JXG.Math.Numerics.Visvalingam(cu1.points, 6); 4694 * // Plot the remaining points 4695 * len = points.length; 4696 * this.dataX = []; 4697 * this.dataY = []; 4698 * for (i = 0; i < len; i++) { 4699 * this.dataX.push(points[i].usrCoords[1]); 4700 * this.dataY.push(points[i].usrCoords[2]); 4701 * } 4702 * }; 4703 * board.update(); 4704 * 4705 * </pre><div id="JXGce0cc55c-b592-11e6-8270-104a7d3be7eb" class="jxgbox" style="width: 300px; height: 300px;"></div> 4706 * <script type="text/javascript"> 4707 * (function() { 4708 * var board = JXG.JSXGraph.initBoard('JXGce0cc55c-b592-11e6-8270-104a7d3be7eb', 4709 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 4710 * 4711 * var i, p = []; 4712 * for (i = 0; i < 5; ++i) { 4713 * p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6])); 4714 * } 4715 * 4716 * // Plot a cardinal spline curve 4717 * var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5); 4718 * var cu1 = board.create('curve', splineArr, {strokeColor: 'green'}); 4719 * 4720 * var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'}); 4721 * c.updateDataArray = function() { 4722 * var i, len, points; 4723 * 4724 * // Reduce number of intermediate points with Visvakingam-Whyatt to 6 4725 * points = JXG.Math.Numerics.Visvalingam(cu1.points, 6); 4726 * // Plot the remaining points 4727 * len = points.length; 4728 * this.dataX = []; 4729 * this.dataY = []; 4730 * for (i = 0; i < len; i++) { 4731 * this.dataX.push(points[i].usrCoords[1]); 4732 * this.dataY.push(points[i].usrCoords[2]); 4733 * } 4734 * }; 4735 * board.update(); 4736 * 4737 * })(); 4738 * 4739 * </script><pre> 4740 * 4741 */ 4742 Visvalingam: function (pts, numPoints) { 4743 var i, 4744 len, 4745 vol, 4746 lastVol, 4747 linkedList = [], 4748 heap = [], 4749 points = [], 4750 lft, 4751 rt, 4752 lft2, 4753 rt2, 4754 obj; 4755 4756 len = pts.length; 4757 4758 // At least one intermediate point is needed 4759 if (len <= 2) { 4760 return pts; 4761 } 4762 4763 // Fill the linked list 4764 // Add first point to the linked list 4765 linkedList[0] = { 4766 used: true, 4767 lft: null, 4768 node: null 4769 }; 4770 4771 // Add all intermediate points to the linked list, 4772 // whose triangle area is nonzero. 4773 lft = 0; 4774 for (i = 1; i < len - 1; i++) { 4775 vol = Math.abs( 4776 JXG.Math.Numerics.det([ 4777 pts[i - 1].usrCoords, 4778 pts[i].usrCoords, 4779 pts[i + 1].usrCoords 4780 ]) 4781 ); 4782 if (!isNaN(vol)) { 4783 obj = { 4784 v: vol, 4785 idx: i 4786 }; 4787 heap.push(obj); 4788 linkedList[i] = { 4789 used: true, 4790 lft: lft, 4791 node: obj 4792 }; 4793 linkedList[lft].rt = i; 4794 lft = i; 4795 } 4796 } 4797 4798 // Add last point to the linked list 4799 linkedList[len - 1] = { 4800 used: true, 4801 rt: null, 4802 lft: lft, 4803 node: null 4804 }; 4805 linkedList[lft].rt = len - 1; 4806 4807 // Remove points until only numPoints intermediate points remain 4808 lastVol = -Infinity; 4809 while (heap.length > numPoints) { 4810 // Sort the heap with the updated volume values 4811 heap.sort(function (a, b) { 4812 // descending sort 4813 return b.v - a.v; 4814 }); 4815 4816 // Remove the point with the smallest triangle 4817 i = heap.pop().idx; 4818 linkedList[i].used = false; 4819 lastVol = linkedList[i].node.v; 4820 4821 // Update the pointers of the linked list 4822 lft = linkedList[i].lft; 4823 rt = linkedList[i].rt; 4824 linkedList[lft].rt = rt; 4825 linkedList[rt].lft = lft; 4826 4827 // Update the values for the volumes in the linked list 4828 lft2 = linkedList[lft].lft; 4829 if (lft2 !== null) { 4830 vol = Math.abs( 4831 JXG.Math.Numerics.det([ 4832 pts[lft2].usrCoords, 4833 pts[lft].usrCoords, 4834 pts[rt].usrCoords 4835 ]) 4836 ); 4837 4838 linkedList[lft].node.v = vol >= lastVol ? vol : lastVol; 4839 } 4840 rt2 = linkedList[rt].rt; 4841 if (rt2 !== null) { 4842 vol = Math.abs( 4843 JXG.Math.Numerics.det([ 4844 pts[lft].usrCoords, 4845 pts[rt].usrCoords, 4846 pts[rt2].usrCoords 4847 ]) 4848 ); 4849 4850 linkedList[rt].node.v = vol >= lastVol ? vol : lastVol; 4851 } 4852 } 4853 4854 // Return an array with the remaining points 4855 i = 0; 4856 points = [pts[i]]; 4857 do { 4858 i = linkedList[i].rt; 4859 points.push(pts[i]); 4860 } while (linkedList[i].rt !== null); 4861 4862 return points; 4863 } 4864 }; 4865 4866 export default Mat.Numerics; 4867