1 /* 2 Copyright 2008-2023 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"; 42 import Type from "../utils/type"; 43 import Env from "../utils/env"; 44 import Mat from "./math"; 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 * @private 190 * Gauss-Bareiss algorithm to compute the 191 * determinant of matrix without fractions. 192 * See Henri Cohen, "A Course in Computational 193 * Algebraic Number Theory (Graduate texts 194 * in mathematics; 138)", Springer-Verlag, 195 * ISBN 3-540-55640-0 / 0-387-55640-0 196 * Third, Corrected Printing 1996 197 * "Algorithm 2.2.6", pg. 52-53 198 * @memberof JXG.Math.Numerics 199 */ 200 gaussBareiss: function (mat) { 201 var k, 202 c, 203 s, 204 i, 205 j, 206 p, 207 n, 208 M, 209 t, 210 eps = Mat.eps; 211 212 n = mat.length; 213 214 if (n <= 0) { 215 return 0; 216 } 217 218 if (mat[0].length < n) { 219 n = mat[0].length; 220 } 221 222 // Copy the input matrix to M 223 M = []; 224 225 for (i = 0; i < n; i++) { 226 M[i] = mat[i].slice(0, n); 227 } 228 229 c = 1; 230 s = 1; 231 232 for (k = 0; k < n - 1; k++) { 233 p = M[k][k]; 234 235 // Pivot step 236 if (Math.abs(p) < eps) { 237 for (i = k + 1; i < n; i++) { 238 if (Math.abs(M[i][k]) >= eps) { 239 break; 240 } 241 } 242 243 // No nonzero entry found in column k -> det(M) = 0 244 if (i === n) { 245 return 0.0; 246 } 247 248 // swap row i and k partially 249 for (j = k; j < n; j++) { 250 t = M[i][j]; 251 M[i][j] = M[k][j]; 252 M[k][j] = t; 253 } 254 s = -s; 255 p = M[k][k]; 256 } 257 258 // Main step 259 for (i = k + 1; i < n; i++) { 260 for (j = k + 1; j < n; j++) { 261 t = p * M[i][j] - M[i][k] * M[k][j]; 262 M[i][j] = t / c; 263 } 264 } 265 266 c = p; 267 } 268 269 return s * M[n - 1][n - 1]; 270 }, 271 272 /** 273 * Computes the determinant of a square nxn matrix with the 274 * Gauss-Bareiss algorithm. 275 * @param {Array} mat Matrix. 276 * @returns {Number} The determinant pf the matrix mat. 277 * The empty matrix returns 0. 278 * @memberof JXG.Math.Numerics 279 */ 280 det: function (mat) { 281 var n = mat.length; 282 283 if (n === 2 && mat[0].length === 2) { 284 return mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1]; 285 } 286 287 return this.gaussBareiss(mat); 288 }, 289 290 /** 291 * Compute the Eigenvalues and Eigenvectors of a symmetric 3x3 matrix with the Jacobi method 292 * Adaption of a FORTRAN program by Ed Wilson, Dec. 25, 1990 293 * @param {Array} Ain A symmetric 3x3 matrix. 294 * @returns {Array} [A,V] the matrices A and V. The diagonal of A contains the Eigenvalues, V contains the Eigenvectors. 295 * @memberof JXG.Math.Numerics 296 */ 297 Jacobi: function (Ain) { 298 var i, 299 j, 300 k, 301 aa, 302 si, 303 co, 304 tt, 305 ssum, 306 amax, 307 eps = Mat.eps * Mat.eps, 308 sum = 0.0, 309 n = Ain.length, 310 V = [ 311 [0, 0, 0], 312 [0, 0, 0], 313 [0, 0, 0] 314 ], 315 A = [ 316 [0, 0, 0], 317 [0, 0, 0], 318 [0, 0, 0] 319 ], 320 nloops = 0; 321 322 // Initialization. Set initial Eigenvectors. 323 for (i = 0; i < n; i++) { 324 for (j = 0; j < n; j++) { 325 V[i][j] = 0.0; 326 A[i][j] = Ain[i][j]; 327 sum += Math.abs(A[i][j]); 328 } 329 V[i][i] = 1.0; 330 } 331 332 // Trivial problems 333 if (n === 1) { 334 return [A, V]; 335 } 336 337 if (sum <= 0.0) { 338 return [A, V]; 339 } 340 341 sum /= n * n; 342 343 // Reduce matrix to diagonal 344 do { 345 ssum = 0.0; 346 amax = 0.0; 347 for (j = 1; j < n; j++) { 348 for (i = 0; i < j; i++) { 349 // Check if A[i][j] is to be reduced 350 aa = Math.abs(A[i][j]); 351 352 if (aa > amax) { 353 amax = aa; 354 } 355 356 ssum += aa; 357 358 if (aa >= eps) { 359 // calculate rotation angle 360 aa = Math.atan2(2.0 * A[i][j], A[i][i] - A[j][j]) * 0.5; 361 si = Math.sin(aa); 362 co = Math.cos(aa); 363 364 // Modify 'i' and 'j' columns 365 for (k = 0; k < n; k++) { 366 tt = A[k][i]; 367 A[k][i] = co * tt + si * A[k][j]; 368 A[k][j] = -si * tt + co * A[k][j]; 369 tt = V[k][i]; 370 V[k][i] = co * tt + si * V[k][j]; 371 V[k][j] = -si * tt + co * V[k][j]; 372 } 373 374 // Modify diagonal terms 375 A[i][i] = co * A[i][i] + si * A[j][i]; 376 A[j][j] = -si * A[i][j] + co * A[j][j]; 377 A[i][j] = 0.0; 378 379 // Make 'A' matrix symmetrical 380 for (k = 0; k < n; k++) { 381 A[i][k] = A[k][i]; 382 A[j][k] = A[k][j]; 383 } 384 // A[i][j] made zero by rotation 385 } 386 } 387 } 388 nloops += 1; 389 } while (Math.abs(ssum) / sum > eps && nloops < 2000); 390 391 return [A, V]; 392 }, 393 394 /** 395 * Calculates the integral of function f over interval using Newton-Cotes-algorithm. 396 * @param {Array} interval The integration interval, e.g. [0, 3]. 397 * @param {function} f A function which takes one argument of type number and returns a number. 398 * @param {Object} [config] The algorithm setup. Accepted properties are number_of_nodes of type number and integration_type 399 * with value being either 'trapez', 'simpson', or 'milne'. 400 * @param {Number} [config.number_of_nodes=28] 401 * @param {String} [config.integration_type='milne'] Possible values are 'milne', 'simpson', 'trapez' 402 * @returns {Number} Integral value of f over interval 403 * @throws {Error} If config.number_of_nodes doesn't match config.integration_type an exception is thrown. If you want to use 404 * simpson rule respectively milne rule config.number_of_nodes must be dividable by 2 respectively 4. 405 * @example 406 * function f(x) { 407 * return x*x; 408 * } 409 * 410 * // calculates integral of <tt>f</tt> from 0 to 2. 411 * var area1 = JXG.Math.Numerics.NewtonCotes([0, 2], f); 412 * 413 * // the same with an anonymous function 414 * var area2 = JXG.Math.Numerics.NewtonCotes([0, 2], function (x) { return x*x; }); 415 * 416 * // use trapez rule with 16 nodes 417 * var area3 = JXG.Math.Numerics.NewtonCotes([0, 2], f, 418 * {number_of_nodes: 16, integration_type: 'trapez'}); 419 * @memberof JXG.Math.Numerics 420 */ 421 NewtonCotes: function (interval, f, config) { 422 var evaluation_point, 423 i, 424 number_of_intervals, 425 integral_value = 0.0, 426 number_of_nodes = 427 config && Type.isNumber(config.number_of_nodes) ? config.number_of_nodes : 28, 428 available_types = { trapez: true, simpson: true, milne: true }, 429 integration_type = 430 config && 431 config.integration_type && 432 available_types.hasOwnProperty(config.integration_type) && 433 available_types[config.integration_type] 434 ? config.integration_type 435 : "milne", 436 step_size = (interval[1] - interval[0]) / number_of_nodes; 437 438 switch (integration_type) { 439 case "trapez": 440 integral_value = (f(interval[0]) + f(interval[1])) * 0.5; 441 evaluation_point = interval[0]; 442 443 for (i = 0; i < number_of_nodes - 1; i++) { 444 evaluation_point += step_size; 445 integral_value += f(evaluation_point); 446 } 447 448 integral_value *= step_size; 449 break; 450 case "simpson": 451 if (number_of_nodes % 2 > 0) { 452 throw new Error( 453 "JSXGraph: INT_SIMPSON requires config.number_of_nodes dividable by 2." 454 ); 455 } 456 457 number_of_intervals = number_of_nodes / 2.0; 458 integral_value = f(interval[0]) + f(interval[1]); 459 evaluation_point = interval[0]; 460 461 for (i = 0; i < number_of_intervals - 1; i++) { 462 evaluation_point += 2.0 * step_size; 463 integral_value += 2.0 * f(evaluation_point); 464 } 465 466 evaluation_point = interval[0] - step_size; 467 468 for (i = 0; i < number_of_intervals; i++) { 469 evaluation_point += 2.0 * step_size; 470 integral_value += 4.0 * f(evaluation_point); 471 } 472 473 integral_value *= step_size / 3.0; 474 break; 475 default: 476 if (number_of_nodes % 4 > 0) { 477 throw new Error( 478 "JSXGraph: Error in INT_MILNE: config.number_of_nodes must be a multiple of 4" 479 ); 480 } 481 482 number_of_intervals = number_of_nodes * 0.25; 483 integral_value = 7.0 * (f(interval[0]) + f(interval[1])); 484 evaluation_point = interval[0]; 485 486 for (i = 0; i < number_of_intervals - 1; i++) { 487 evaluation_point += 4.0 * step_size; 488 integral_value += 14.0 * f(evaluation_point); 489 } 490 491 evaluation_point = interval[0] - 3.0 * step_size; 492 493 for (i = 0; i < number_of_intervals; i++) { 494 evaluation_point += 4.0 * step_size; 495 integral_value += 496 32.0 * (f(evaluation_point) + f(evaluation_point + 2 * step_size)); 497 } 498 499 evaluation_point = interval[0] - 2.0 * step_size; 500 501 for (i = 0; i < number_of_intervals; i++) { 502 evaluation_point += 4.0 * step_size; 503 integral_value += 12.0 * f(evaluation_point); 504 } 505 506 integral_value *= (2.0 * step_size) / 45.0; 507 } 508 return integral_value; 509 }, 510 511 /** 512 * Calculates the integral of function f over interval using Romberg iteration. 513 * @param {Array} interval The integration interval, e.g. [0, 3]. 514 * @param {function} f A function which takes one argument of type number and returns a number. 515 * @param {Object} [config] The algorithm setup. Accepted properties are max_iterations of type number and precision eps. 516 * @param {Number} [config.max_iterations=20] 517 * @param {Number} [config.eps=0.0000001] 518 * @returns {Number} Integral value of f over interval 519 * @example 520 * function f(x) { 521 * return x*x; 522 * } 523 * 524 * // calculates integral of <tt>f</tt> from 0 to 2. 525 * var area1 = JXG.Math.Numerics.Romberg([0, 2], f); 526 * 527 * // the same with an anonymous function 528 * var area2 = JXG.Math.Numerics.Romberg([0, 2], function (x) { return x*x; }); 529 * 530 * // use trapez rule with maximum of 16 iterations or stop if the precision 0.0001 has been reached. 531 * var area3 = JXG.Math.Numerics.Romberg([0, 2], f, 532 * {max_iterations: 16, eps: 0.0001}); 533 * @memberof JXG.Math.Numerics 534 */ 535 Romberg: function (interval, f, config) { 536 var a, 537 b, 538 h, 539 s, 540 n, 541 k, 542 i, 543 q, 544 p = [], 545 integral = 0.0, 546 last = Infinity, 547 m = config && Type.isNumber(config.max_iterations) ? config.max_iterations : 20, 548 eps = config && Type.isNumber(config.eps) ? config.eps : config.eps || 0.0000001; 549 550 a = interval[0]; 551 b = interval[1]; 552 h = b - a; 553 n = 1; 554 555 p[0] = 0.5 * h * (f(a) + f(b)); 556 557 for (k = 0; k < m; ++k) { 558 s = 0; 559 h *= 0.5; 560 n *= 2; 561 q = 1; 562 563 for (i = 1; i < n; i += 2) { 564 s += f(a + i * h); 565 } 566 567 p[k + 1] = 0.5 * p[k] + s * h; 568 569 integral = p[k + 1]; 570 for (i = k - 1; i >= 0; --i) { 571 q *= 4; 572 p[i] = p[i + 1] + (p[i + 1] - p[i]) / (q - 1.0); 573 integral = p[i]; 574 } 575 576 if (Math.abs(integral - last) < eps * Math.abs(integral)) { 577 break; 578 } 579 last = integral; 580 } 581 582 return integral; 583 }, 584 585 /** 586 * Calculates the integral of function f over interval using Gauss-Legendre quadrature. 587 * @param {Array} interval The integration interval, e.g. [0, 3]. 588 * @param {function} f A function which takes one argument of type number and returns a number. 589 * @param {Object} [config] The algorithm setup. Accepted property is the order n of type number. n is allowed to take 590 * values between 2 and 18, default value is 12. 591 * @param {Number} [config.n=16] 592 * @returns {Number} Integral value of f over interval 593 * @example 594 * function f(x) { 595 * return x*x; 596 * } 597 * 598 * // calculates integral of <tt>f</tt> from 0 to 2. 599 * var area1 = JXG.Math.Numerics.GaussLegendre([0, 2], f); 600 * 601 * // the same with an anonymous function 602 * var area2 = JXG.Math.Numerics.GaussLegendre([0, 2], function (x) { return x*x; }); 603 * 604 * // use 16 point Gauss-Legendre rule. 605 * var area3 = JXG.Math.Numerics.GaussLegendre([0, 2], f, 606 * {n: 16}); 607 * @memberof JXG.Math.Numerics 608 */ 609 GaussLegendre: function (interval, f, config) { 610 var a, 611 b, 612 i, 613 m, 614 xp, 615 xm, 616 result = 0.0, 617 table_xi = [], 618 table_w = [], 619 xi, 620 w, 621 n = config && Type.isNumber(config.n) ? config.n : 12; 622 623 if (n > 18) { 624 n = 18; 625 } 626 627 /* n = 2 */ 628 table_xi[2] = [0.5773502691896257645091488]; 629 table_w[2] = [1.0]; 630 631 /* n = 4 */ 632 table_xi[4] = [0.3399810435848562648026658, 0.8611363115940525752239465]; 633 table_w[4] = [0.6521451548625461426269361, 0.3478548451374538573730639]; 634 635 /* n = 6 */ 636 table_xi[6] = [ 637 0.2386191860831969086305017, 0.6612093864662645136613996, 638 0.9324695142031520278123016 639 ]; 640 table_w[6] = [ 641 0.4679139345726910473898703, 0.3607615730481386075698335, 642 0.1713244923791703450402961 643 ]; 644 645 /* n = 8 */ 646 table_xi[8] = [ 647 0.1834346424956498049394761, 0.525532409916328985817739, 648 0.7966664774136267395915539, 0.9602898564975362316835609 649 ]; 650 table_w[8] = [ 651 0.3626837833783619829651504, 0.3137066458778872873379622, 652 0.222381034453374470544356, 0.1012285362903762591525314 653 ]; 654 655 /* n = 10 */ 656 table_xi[10] = [ 657 0.148874338981631210884826, 0.4333953941292471907992659, 658 0.6794095682990244062343274, 0.8650633666889845107320967, 0.973906528517171720077964 659 ]; 660 table_w[10] = [ 661 0.295524224714752870173893, 0.2692667193099963550912269, 662 0.2190863625159820439955349, 0.1494513491505805931457763, 663 0.0666713443086881375935688 664 ]; 665 666 /* n = 12 */ 667 table_xi[12] = [ 668 0.1252334085114689154724414, 0.3678314989981801937526915, 669 0.5873179542866174472967024, 0.7699026741943046870368938, 670 0.9041172563704748566784659, 0.9815606342467192506905491 671 ]; 672 table_w[12] = [ 673 0.2491470458134027850005624, 0.2334925365383548087608499, 674 0.2031674267230659217490645, 0.1600783285433462263346525, 675 0.1069393259953184309602547, 0.047175336386511827194616 676 ]; 677 678 /* n = 14 */ 679 table_xi[14] = [ 680 0.1080549487073436620662447, 0.3191123689278897604356718, 681 0.5152486363581540919652907, 0.6872929048116854701480198, 682 0.8272013150697649931897947, 0.9284348836635735173363911, 683 0.9862838086968123388415973 684 ]; 685 table_w[14] = [ 686 0.2152638534631577901958764, 0.2051984637212956039659241, 687 0.1855383974779378137417166, 0.1572031671581935345696019, 688 0.1215185706879031846894148, 0.0801580871597602098056333, 689 0.0351194603317518630318329 690 ]; 691 692 /* n = 16 */ 693 table_xi[16] = [ 694 0.0950125098376374401853193, 0.2816035507792589132304605, 695 0.4580167776572273863424194, 0.6178762444026437484466718, 696 0.7554044083550030338951012, 0.8656312023878317438804679, 697 0.9445750230732325760779884, 0.9894009349916499325961542 698 ]; 699 table_w[16] = [ 700 0.1894506104550684962853967, 0.1826034150449235888667637, 701 0.1691565193950025381893121, 0.1495959888165767320815017, 702 0.1246289712555338720524763, 0.0951585116824927848099251, 703 0.0622535239386478928628438, 0.0271524594117540948517806 704 ]; 705 706 /* n = 18 */ 707 table_xi[18] = [ 708 0.0847750130417353012422619, 0.2518862256915055095889729, 709 0.4117511614628426460359318, 0.5597708310739475346078715, 710 0.6916870430603532078748911, 0.8037049589725231156824175, 711 0.8926024664975557392060606, 0.9558239495713977551811959, 0.991565168420930946730016 712 ]; 713 table_w[18] = [ 714 0.1691423829631435918406565, 0.1642764837458327229860538, 715 0.154684675126265244925418, 0.1406429146706506512047313, 716 0.1225552067114784601845191, 0.100942044106287165562814, 717 0.0764257302548890565291297, 0.0497145488949697964533349, 718 0.0216160135264833103133427 719 ]; 720 721 /* n = 3 */ 722 table_xi[3] = [0.0, 0.7745966692414833770358531]; 723 table_w[3] = [0.8888888888888888888888889, 0.5555555555555555555555556]; 724 725 /* n = 5 */ 726 table_xi[5] = [0.0, 0.5384693101056830910363144, 0.9061798459386639927976269]; 727 table_w[5] = [ 728 0.5688888888888888888888889, 0.4786286704993664680412915, 0.236926885056189087514264 729 ]; 730 731 /* n = 7 */ 732 table_xi[7] = [ 733 0.0, 0.4058451513773971669066064, 0.7415311855993944398638648, 734 0.9491079123427585245261897 735 ]; 736 table_w[7] = [ 737 0.417959183673469387755102, 0.3818300505051189449503698, 738 0.2797053914892766679014678, 0.1294849661688696932706114 739 ]; 740 741 /* n = 9 */ 742 table_xi[9] = [ 743 0.0, 0.324253423403808929038538, 0.613371432700590397308702, 744 0.8360311073266357942994298, 0.9681602395076260898355762 745 ]; 746 table_w[9] = [ 747 0.3302393550012597631645251, 0.3123470770400028400686304, 748 0.2606106964029354623187429, 0.180648160694857404058472, 0.0812743883615744119718922 749 ]; 750 751 /* n = 11 */ 752 table_xi[11] = [ 753 0.0, 0.269543155952344972331532, 0.5190961292068118159257257, 754 0.7301520055740493240934163, 0.8870625997680952990751578, 0.978228658146056992803938 755 ]; 756 table_w[11] = [ 757 0.2729250867779006307144835, 0.2628045445102466621806889, 758 0.2331937645919904799185237, 0.1862902109277342514260976, 759 0.1255803694649046246346943, 0.0556685671161736664827537 760 ]; 761 762 /* n = 13 */ 763 table_xi[13] = [ 764 0.0, 0.2304583159551347940655281, 0.4484927510364468528779129, 765 0.6423493394403402206439846, 0.8015780907333099127942065, 766 0.9175983992229779652065478, 0.9841830547185881494728294 767 ]; 768 table_w[13] = [ 769 0.2325515532308739101945895, 0.2262831802628972384120902, 770 0.2078160475368885023125232, 0.1781459807619457382800467, 771 0.1388735102197872384636018, 0.0921214998377284479144218, 772 0.0404840047653158795200216 773 ]; 774 775 /* n = 15 */ 776 table_xi[15] = [ 777 0.0, 0.2011940939974345223006283, 0.3941513470775633698972074, 778 0.5709721726085388475372267, 0.7244177313601700474161861, 779 0.8482065834104272162006483, 0.9372733924007059043077589, 780 0.9879925180204854284895657 781 ]; 782 table_w[15] = [ 783 0.2025782419255612728806202, 0.1984314853271115764561183, 784 0.1861610000155622110268006, 0.1662692058169939335532009, 785 0.1395706779261543144478048, 0.1071592204671719350118695, 786 0.0703660474881081247092674, 0.0307532419961172683546284 787 ]; 788 789 /* n = 17 */ 790 table_xi[17] = [ 791 0.0, 0.1784841814958478558506775, 0.3512317634538763152971855, 792 0.5126905370864769678862466, 0.6576711592166907658503022, 793 0.7815140038968014069252301, 0.8802391537269859021229557, 794 0.950675521768767761222717, 0.990575475314417335675434 795 ]; 796 table_w[17] = [ 797 0.1794464703562065254582656, 0.176562705366992646325271, 798 0.1680041021564500445099707, 0.1540457610768102880814316, 0.13513636846852547328632, 799 0.1118838471934039710947884, 0.0850361483171791808835354, 800 0.0554595293739872011294402, 0.02414830286854793196011 801 ]; 802 803 a = interval[0]; 804 b = interval[1]; 805 806 //m = Math.ceil(n * 0.5); 807 m = (n + 1) >> 1; 808 809 xi = table_xi[n]; 810 w = table_w[n]; 811 812 xm = 0.5 * (b - a); 813 xp = 0.5 * (b + a); 814 815 if (n & (1 === 1)) { 816 // n odd 817 result = w[0] * f(xp); 818 for (i = 1; i < m; ++i) { 819 result += w[i] * (f(xp + xm * xi[i]) + f(xp - xm * xi[i])); 820 } 821 } else { 822 // n even 823 result = 0.0; 824 for (i = 0; i < m; ++i) { 825 result += w[i] * (f(xp + xm * xi[i]) + f(xp - xm * xi[i])); 826 } 827 } 828 829 return xm * result; 830 }, 831 832 /** 833 * Scale error in Gauss Kronrod quadrature. 834 * Internal method used in {@link JXG.Math.Numerics._gaussKronrod}. 835 * @private 836 */ 837 _rescale_error: function (err, result_abs, result_asc) { 838 var scale, 839 min_err, 840 DBL_MIN = 2.2250738585072014e-308, 841 DBL_EPS = 2.2204460492503131e-16; 842 843 err = Math.abs(err); 844 if (result_asc !== 0 && err !== 0) { 845 scale = Math.pow((200 * err) / result_asc, 1.5); 846 847 if (scale < 1.0) { 848 err = result_asc * scale; 849 } else { 850 err = result_asc; 851 } 852 } 853 if (result_abs > DBL_MIN / (50 * DBL_EPS)) { 854 min_err = 50 * DBL_EPS * result_abs; 855 856 if (min_err > err) { 857 err = min_err; 858 } 859 } 860 861 return err; 862 }, 863 864 /** 865 * Generic Gauss-Kronrod quadrature algorithm. 866 * Internal method used in {@link JXG.Math.Numerics.GaussKronrod15}, 867 * {@link JXG.Math.Numerics.GaussKronrod21}, 868 * {@link JXG.Math.Numerics.GaussKronrod31}. 869 * Taken from QUADPACK. 870 * 871 * @param {Array} interval The integration interval, e.g. [0, 3]. 872 * @param {function} f A function which takes one argument of type number and returns a number. 873 * @param {Number} n order 874 * @param {Array} xgk Kronrod quadrature abscissae 875 * @param {Array} wg Weights of the Gauss rule 876 * @param {Array} wgk Weights of the Kronrod rule 877 * @param {Object} resultObj Object returning resultObj.abserr, resultObj.resabs, resultObj.resasc. 878 * See the library QUADPACK for an explanation. 879 * 880 * @returns {Number} Integral value of f over interval 881 * 882 * @private 883 */ 884 _gaussKronrod: function (interval, f, n, xgk, wg, wgk, resultObj) { 885 var a = interval[0], 886 b = interval[1], 887 up, 888 result, 889 center = 0.5 * (a + b), 890 half_length = 0.5 * (b - a), 891 abs_half_length = Math.abs(half_length), 892 f_center = f(center), 893 result_gauss = 0.0, 894 result_kronrod = f_center * wgk[n - 1], 895 result_abs = Math.abs(result_kronrod), 896 result_asc = 0.0, 897 mean = 0.0, 898 err = 0.0, 899 j, 900 jtw, 901 abscissa, 902 fval1, 903 fval2, 904 fsum, 905 jtwm1, 906 fv1 = [], 907 fv2 = []; 908 909 if (n % 2 === 0) { 910 result_gauss = f_center * wg[n / 2 - 1]; 911 } 912 913 up = Math.floor((n - 1) / 2); 914 for (j = 0; j < up; j++) { 915 jtw = j * 2 + 1; // in original fortran j=1,2,3 jtw=2,4,6 916 abscissa = half_length * xgk[jtw]; 917 fval1 = f(center - abscissa); 918 fval2 = f(center + abscissa); 919 fsum = fval1 + fval2; 920 fv1[jtw] = fval1; 921 fv2[jtw] = fval2; 922 result_gauss += wg[j] * fsum; 923 result_kronrod += wgk[jtw] * fsum; 924 result_abs += wgk[jtw] * (Math.abs(fval1) + Math.abs(fval2)); 925 } 926 927 up = Math.floor(n / 2); 928 for (j = 0; j < up; j++) { 929 jtwm1 = j * 2; 930 abscissa = half_length * xgk[jtwm1]; 931 fval1 = f(center - abscissa); 932 fval2 = f(center + abscissa); 933 fv1[jtwm1] = fval1; 934 fv2[jtwm1] = fval2; 935 result_kronrod += wgk[jtwm1] * (fval1 + fval2); 936 result_abs += wgk[jtwm1] * (Math.abs(fval1) + Math.abs(fval2)); 937 } 938 939 mean = result_kronrod * 0.5; 940 result_asc = wgk[n - 1] * Math.abs(f_center - mean); 941 942 for (j = 0; j < n - 1; j++) { 943 result_asc += wgk[j] * (Math.abs(fv1[j] - mean) + Math.abs(fv2[j] - mean)); 944 } 945 946 // scale by the width of the integration region 947 err = (result_kronrod - result_gauss) * half_length; 948 949 result_kronrod *= half_length; 950 result_abs *= abs_half_length; 951 result_asc *= abs_half_length; 952 result = result_kronrod; 953 954 resultObj.abserr = this._rescale_error(err, result_abs, result_asc); 955 resultObj.resabs = result_abs; 956 resultObj.resasc = result_asc; 957 958 return result; 959 }, 960 961 /** 962 * 15 point Gauss-Kronrod quadrature algorithm, see the library QUADPACK 963 * @param {Array} interval The integration interval, e.g. [0, 3]. 964 * @param {function} f A function which takes one argument of type number and returns a number. 965 * @param {Object} resultObj Object returning resultObj.abserr, resultObj.resabs, resultObj.resasc. See the library 966 * QUADPACK for an explanation. 967 * 968 * @returns {Number} Integral value of f over interval 969 * 970 * @memberof JXG.Math.Numerics 971 */ 972 GaussKronrod15: function (interval, f, resultObj) { 973 /* Gauss quadrature weights and kronrod quadrature abscissae and 974 weights as evaluated with 80 decimal digit arithmetic by 975 L. W. Fullerton, Bell Labs, Nov. 1981. */ 976 977 var xgk = 978 /* abscissae of the 15-point kronrod rule */ 979 [ 980 0.991455371120812639206854697526329, 0.949107912342758524526189684047851, 981 0.864864423359769072789712788640926, 0.741531185599394439863864773280788, 982 0.58608723546769113029414483825873, 0.405845151377397166906606412076961, 983 0.207784955007898467600689403773245, 0.0 984 ], 985 /* xgk[1], xgk[3], ... abscissae of the 7-point gauss rule. 986 xgk[0], xgk[2], ... abscissae to optimally extend the 7-point gauss rule */ 987 988 wg = 989 /* weights of the 7-point gauss rule */ 990 [ 991 0.129484966168869693270611432679082, 0.27970539148927666790146777142378, 992 0.381830050505118944950369775488975, 0.417959183673469387755102040816327 993 ], 994 wgk = 995 /* weights of the 15-point kronrod rule */ 996 [ 997 0.02293532201052922496373200805897, 0.063092092629978553290700663189204, 998 0.104790010322250183839876322541518, 0.140653259715525918745189590510238, 999 0.16900472663926790282658342659855, 0.190350578064785409913256402421014, 1000 0.204432940075298892414161999234649, 0.209482141084727828012999174891714 1001 ]; 1002 1003 return this._gaussKronrod(interval, f, 8, xgk, wg, wgk, resultObj); 1004 }, 1005 1006 /** 1007 * 21 point Gauss-Kronrod quadrature algorithm, see the library QUADPACK 1008 * @param {Array} interval The integration interval, e.g. [0, 3]. 1009 * @param {function} f A function which takes one argument of type number and returns a number. 1010 * @param {Object} resultObj Object returning resultObj.abserr, resultObj.resabs, resultObj.resasc. See the library 1011 * QUADPACK for an explanation. 1012 * 1013 * @returns {Number} Integral value of f over interval 1014 * 1015 * @memberof JXG.Math.Numerics 1016 */ 1017 GaussKronrod21: function (interval, f, resultObj) { 1018 /* Gauss quadrature weights and kronrod quadrature abscissae and 1019 weights as evaluated with 80 decimal digit arithmetic by 1020 L. W. Fullerton, Bell Labs, Nov. 1981. */ 1021 1022 var xgk = 1023 /* abscissae of the 21-point kronrod rule */ 1024 [ 1025 0.995657163025808080735527280689003, 0.973906528517171720077964012084452, 1026 0.930157491355708226001207180059508, 0.865063366688984510732096688423493, 1027 0.780817726586416897063717578345042, 0.679409568299024406234327365114874, 1028 0.562757134668604683339000099272694, 0.433395394129247190799265943165784, 1029 0.294392862701460198131126603103866, 0.14887433898163121088482600112972, 0.0 1030 ], 1031 /* xgk[1], xgk[3], ... abscissae of the 10-point gauss rule. 1032 xgk[0], xgk[2], ... abscissae to optimally extend the 10-point gauss rule */ 1033 wg = 1034 /* weights of the 10-point gauss rule */ 1035 [ 1036 0.066671344308688137593568809893332, 0.149451349150580593145776339657697, 1037 0.219086362515982043995534934228163, 0.269266719309996355091226921569469, 1038 0.295524224714752870173892994651338 1039 ], 1040 wgk = 1041 /* weights of the 21-point kronrod rule */ 1042 [ 1043 0.011694638867371874278064396062192, 0.03255816230796472747881897245939, 1044 0.05475589657435199603138130024458, 0.07503967481091995276704314091619, 1045 0.093125454583697605535065465083366, 0.109387158802297641899210590325805, 1046 0.123491976262065851077958109831074, 0.134709217311473325928054001771707, 1047 0.142775938577060080797094273138717, 0.147739104901338491374841515972068, 1048 0.149445554002916905664936468389821 1049 ]; 1050 1051 return this._gaussKronrod(interval, f, 11, xgk, wg, wgk, resultObj); 1052 }, 1053 1054 /** 1055 * 31 point Gauss-Kronrod quadrature algorithm, see the library QUADPACK 1056 * @param {Array} interval The integration interval, e.g. [0, 3]. 1057 * @param {function} f A function which takes one argument of type number and returns a number. 1058 * @param {Object} resultObj Object returning resultObj.abserr, resultObj.resabs, resultObj.resasc. See the library 1059 * QUADPACK for an explanation. 1060 * 1061 * @returns {Number} Integral value of f over interval 1062 * 1063 * @memberof JXG.Math.Numerics 1064 */ 1065 GaussKronrod31: function (interval, f, resultObj) { 1066 /* Gauss quadrature weights and kronrod quadrature abscissae and 1067 weights as evaluated with 80 decimal digit arithmetic by 1068 L. W. Fullerton, Bell Labs, Nov. 1981. */ 1069 1070 var xgk = 1071 /* abscissae of the 21-point kronrod rule */ 1072 [ 1073 0.998002298693397060285172840152271, 0.987992518020485428489565718586613, 1074 0.967739075679139134257347978784337, 0.937273392400705904307758947710209, 1075 0.897264532344081900882509656454496, 0.848206583410427216200648320774217, 1076 0.790418501442465932967649294817947, 0.724417731360170047416186054613938, 1077 0.650996741297416970533735895313275, 0.570972172608538847537226737253911, 1078 0.485081863640239680693655740232351, 0.394151347077563369897207370981045, 1079 0.299180007153168812166780024266389, 0.201194093997434522300628303394596, 1080 0.101142066918717499027074231447392, 0.0 1081 ], 1082 /* xgk[1], xgk[3], ... abscissae of the 10-point gauss rule. 1083 xgk[0], xgk[2], ... abscissae to optimally extend the 10-point gauss rule */ 1084 wg = 1085 /* weights of the 10-point gauss rule */ 1086 [ 1087 0.030753241996117268354628393577204, 0.070366047488108124709267416450667, 1088 0.107159220467171935011869546685869, 0.139570677926154314447804794511028, 1089 0.166269205816993933553200860481209, 0.186161000015562211026800561866423, 1090 0.198431485327111576456118326443839, 0.202578241925561272880620199967519 1091 ], 1092 wgk = 1093 /* weights of the 21-point kronrod rule */ 1094 [ 1095 0.005377479872923348987792051430128, 0.015007947329316122538374763075807, 1096 0.025460847326715320186874001019653, 0.03534636079137584622203794847836, 1097 0.04458975132476487660822729937328, 0.05348152469092808726534314723943, 1098 0.062009567800670640285139230960803, 0.069854121318728258709520077099147, 1099 0.076849680757720378894432777482659, 0.083080502823133021038289247286104, 1100 0.088564443056211770647275443693774, 0.093126598170825321225486872747346, 1101 0.096642726983623678505179907627589, 0.099173598721791959332393173484603, 1102 0.10076984552387559504494666261757, 0.101330007014791549017374792767493 1103 ]; 1104 1105 return this._gaussKronrod(interval, f, 16, xgk, wg, wgk, resultObj); 1106 }, 1107 1108 /** 1109 * Generate workspace object for {@link JXG.Math.Numerics.Qag}. 1110 * @param {Array} interval The integration interval, e.g. [0, 3]. 1111 * @param {Number} n Max. limit 1112 * @returns {Object} Workspace object 1113 * 1114 * @private 1115 * @memberof JXG.Math.Numerics 1116 */ 1117 _workspace: function (interval, n) { 1118 return { 1119 limit: n, 1120 size: 0, 1121 nrmax: 0, 1122 i: 0, 1123 alist: [interval[0]], 1124 blist: [interval[1]], 1125 rlist: [0.0], 1126 elist: [0.0], 1127 order: [0], 1128 level: [0], 1129 1130 qpsrt: function () { 1131 var last = this.size - 1, 1132 limit = this.limit, 1133 errmax, 1134 errmin, 1135 i, 1136 k, 1137 top, 1138 i_nrmax = this.nrmax, 1139 i_maxerr = this.order[i_nrmax]; 1140 1141 /* Check whether the list contains more than two error estimates */ 1142 if (last < 2) { 1143 this.order[0] = 0; 1144 this.order[1] = 1; 1145 this.i = i_maxerr; 1146 return; 1147 } 1148 1149 errmax = this.elist[i_maxerr]; 1150 1151 /* This part of the routine is only executed if, due to a difficult 1152 integrand, subdivision increased the error estimate. In the normal 1153 case the insert procedure should start after the nrmax-th largest 1154 error estimate. */ 1155 while (i_nrmax > 0 && errmax > this.elist[this.order[i_nrmax - 1]]) { 1156 this.order[i_nrmax] = this.order[i_nrmax - 1]; 1157 i_nrmax--; 1158 } 1159 1160 /* Compute the number of elements in the list to be maintained in 1161 descending order. This number depends on the number of 1162 subdivisions still allowed. */ 1163 if (last < limit / 2 + 2) { 1164 top = last; 1165 } else { 1166 top = limit - last + 1; 1167 } 1168 1169 /* Insert errmax by traversing the list top-down, starting 1170 comparison from the element elist(order(i_nrmax+1)). */ 1171 i = i_nrmax + 1; 1172 1173 /* The order of the tests in the following line is important to 1174 prevent a segmentation fault */ 1175 while (i < top && errmax < this.elist[this.order[i]]) { 1176 this.order[i - 1] = this.order[i]; 1177 i++; 1178 } 1179 1180 this.order[i - 1] = i_maxerr; 1181 1182 /* Insert errmin by traversing the list bottom-up */ 1183 errmin = this.elist[last]; 1184 k = top - 1; 1185 1186 while (k > i - 2 && errmin >= this.elist[this.order[k]]) { 1187 this.order[k + 1] = this.order[k]; 1188 k--; 1189 } 1190 1191 this.order[k + 1] = last; 1192 1193 /* Set i_max and e_max */ 1194 i_maxerr = this.order[i_nrmax]; 1195 this.i = i_maxerr; 1196 this.nrmax = i_nrmax; 1197 }, 1198 1199 set_initial_result: function (result, error) { 1200 this.size = 1; 1201 this.rlist[0] = result; 1202 this.elist[0] = error; 1203 }, 1204 1205 update: function (a1, b1, area1, error1, a2, b2, area2, error2) { 1206 var i_max = this.i, 1207 i_new = this.size, 1208 new_level = this.level[this.i] + 1; 1209 1210 /* append the newly-created intervals to the list */ 1211 1212 if (error2 > error1) { 1213 this.alist[i_max] = a2; /* blist[maxerr] is already == b2 */ 1214 this.rlist[i_max] = area2; 1215 this.elist[i_max] = error2; 1216 this.level[i_max] = new_level; 1217 1218 this.alist[i_new] = a1; 1219 this.blist[i_new] = b1; 1220 this.rlist[i_new] = area1; 1221 this.elist[i_new] = error1; 1222 this.level[i_new] = new_level; 1223 } else { 1224 this.blist[i_max] = b1; /* alist[maxerr] is already == a1 */ 1225 this.rlist[i_max] = area1; 1226 this.elist[i_max] = error1; 1227 this.level[i_max] = new_level; 1228 1229 this.alist[i_new] = a2; 1230 this.blist[i_new] = b2; 1231 this.rlist[i_new] = area2; 1232 this.elist[i_new] = error2; 1233 this.level[i_new] = new_level; 1234 } 1235 1236 this.size++; 1237 1238 if (new_level > this.maximum_level) { 1239 this.maximum_level = new_level; 1240 } 1241 1242 this.qpsrt(); 1243 }, 1244 1245 retrieve: function () { 1246 var i = this.i; 1247 return { 1248 a: this.alist[i], 1249 b: this.blist[i], 1250 r: this.rlist[i], 1251 e: this.elist[i] 1252 }; 1253 }, 1254 1255 sum_results: function () { 1256 var nn = this.size, 1257 k, 1258 result_sum = 0.0; 1259 1260 for (k = 0; k < nn; k++) { 1261 result_sum += this.rlist[k]; 1262 } 1263 1264 return result_sum; 1265 }, 1266 1267 subinterval_too_small: function (a1, a2, b2) { 1268 var e = 2.2204460492503131e-16, 1269 u = 2.2250738585072014e-308, 1270 tmp = (1 + 100 * e) * (Math.abs(a2) + 1000 * u); 1271 1272 return Math.abs(a1) <= tmp && Math.abs(b2) <= tmp; 1273 } 1274 }; 1275 }, 1276 1277 /** 1278 * Quadrature algorithm qag from QUADPACK. 1279 * Internal method used in {@link JXG.Math.Numerics.GaussKronrod15}, 1280 * {@link JXG.Math.Numerics.GaussKronrod21}, 1281 * {@link JXG.Math.Numerics.GaussKronrod31}. 1282 * 1283 * @param {Array} interval The integration interval, e.g. [0, 3]. 1284 * @param {function} f A function which takes one argument of type number and returns a number. 1285 * @param {Object} [config] The algorithm setup. Accepted propert are max. recursion limit of type number, 1286 * and epsrel and epsabs, the relative and absolute required precision of type number. Further, 1287 * q the internal quadrature sub-algorithm of type function. 1288 * @param {Number} [config.limit=15] 1289 * @param {Number} [config.epsrel=0.0000001] 1290 * @param {Number} [config.epsabs=0.0000001] 1291 * @param {Number} [config.q=JXG.Math.Numerics.GaussKronrod15] 1292 * @returns {Number} Integral value of f over interval 1293 * 1294 * @example 1295 * function f(x) { 1296 * return x*x; 1297 * } 1298 * 1299 * // calculates integral of <tt>f</tt> from 0 to 2. 1300 * var area1 = JXG.Math.Numerics.Qag([0, 2], f); 1301 * 1302 * // the same with an anonymous function 1303 * var area2 = JXG.Math.Numerics.Qag([0, 2], function (x) { return x*x; }); 1304 * 1305 * // use JXG.Math.Numerics.GaussKronrod31 rule as sub-algorithm. 1306 * var area3 = JXG.Math.Numerics.Quag([0, 2], f, 1307 * {q: JXG.Math.Numerics.GaussKronrod31}); 1308 * @memberof JXG.Math.Numerics 1309 */ 1310 Qag: function (interval, f, config) { 1311 var DBL_EPS = 2.2204460492503131e-16, 1312 ws = this._workspace(interval, 1000), 1313 limit = config && Type.isNumber(config.limit) ? config.limit : 15, 1314 epsrel = config && Type.isNumber(config.epsrel) ? config.epsrel : 0.0000001, 1315 epsabs = config && Type.isNumber(config.epsabs) ? config.epsabs : 0.0000001, 1316 q = config && Type.isFunction(config.q) ? config.q : this.GaussKronrod15, 1317 resultObj = {}, 1318 area, 1319 errsum, 1320 result0, 1321 abserr0, 1322 resabs0, 1323 resasc0, 1324 result, 1325 tolerance, 1326 iteration = 0, 1327 roundoff_type1 = 0, 1328 roundoff_type2 = 0, 1329 error_type = 0, 1330 round_off, 1331 a1, 1332 b1, 1333 a2, 1334 b2, 1335 a_i, 1336 b_i, 1337 r_i, 1338 e_i, 1339 area1 = 0, 1340 area2 = 0, 1341 area12 = 0, 1342 error1 = 0, 1343 error2 = 0, 1344 error12 = 0, 1345 resasc1, 1346 resasc2, 1347 // resabs1, resabs2, 1348 wsObj, 1349 delta; 1350 1351 if (limit > ws.limit) { 1352 JXG.warn("iteration limit exceeds available workspace"); 1353 } 1354 if (epsabs <= 0 && (epsrel < 50 * Mat.eps || epsrel < 0.5e-28)) { 1355 JXG.warn("tolerance cannot be acheived with given epsabs and epsrel"); 1356 } 1357 1358 result0 = q.apply(this, [interval, f, resultObj]); 1359 abserr0 = resultObj.abserr; 1360 resabs0 = resultObj.resabs; 1361 resasc0 = resultObj.resasc; 1362 1363 ws.set_initial_result(result0, abserr0); 1364 tolerance = Math.max(epsabs, epsrel * Math.abs(result0)); 1365 round_off = 50 * DBL_EPS * resabs0; 1366 1367 if (abserr0 <= round_off && abserr0 > tolerance) { 1368 result = result0; 1369 // abserr = abserr0; 1370 1371 JXG.warn("cannot reach tolerance because of roundoff error on first attempt"); 1372 return -Infinity; 1373 } 1374 1375 if ((abserr0 <= tolerance && abserr0 !== resasc0) || abserr0 === 0.0) { 1376 result = result0; 1377 // abserr = abserr0; 1378 1379 return result; 1380 } 1381 1382 if (limit === 1) { 1383 result = result0; 1384 // abserr = abserr0; 1385 1386 JXG.warn("a maximum of one iteration was insufficient"); 1387 return -Infinity; 1388 } 1389 1390 area = result0; 1391 errsum = abserr0; 1392 iteration = 1; 1393 1394 do { 1395 area1 = 0; 1396 area2 = 0; 1397 area12 = 0; 1398 error1 = 0; 1399 error2 = 0; 1400 error12 = 0; 1401 1402 /* Bisect the subinterval with the largest error estimate */ 1403 wsObj = ws.retrieve(); 1404 a_i = wsObj.a; 1405 b_i = wsObj.b; 1406 r_i = wsObj.r; 1407 e_i = wsObj.e; 1408 1409 a1 = a_i; 1410 b1 = 0.5 * (a_i + b_i); 1411 a2 = b1; 1412 b2 = b_i; 1413 1414 area1 = q.apply(this, [[a1, b1], f, resultObj]); 1415 error1 = resultObj.abserr; 1416 // resabs1 = resultObj.resabs; 1417 resasc1 = resultObj.resasc; 1418 1419 area2 = q.apply(this, [[a2, b2], f, resultObj]); 1420 error2 = resultObj.abserr; 1421 // resabs2 = resultObj.resabs; 1422 resasc2 = resultObj.resasc; 1423 1424 area12 = area1 + area2; 1425 error12 = error1 + error2; 1426 1427 errsum += error12 - e_i; 1428 area += area12 - r_i; 1429 1430 if (resasc1 !== error1 && resasc2 !== error2) { 1431 delta = r_i - area12; 1432 if (Math.abs(delta) <= 1.0e-5 * Math.abs(area12) && error12 >= 0.99 * e_i) { 1433 roundoff_type1++; 1434 } 1435 if (iteration >= 10 && error12 > e_i) { 1436 roundoff_type2++; 1437 } 1438 } 1439 1440 tolerance = Math.max(epsabs, epsrel * Math.abs(area)); 1441 1442 if (errsum > tolerance) { 1443 if (roundoff_type1 >= 6 || roundoff_type2 >= 20) { 1444 error_type = 2; /* round off error */ 1445 } 1446 1447 /* set error flag in the case of bad integrand behaviour at 1448 a point of the integration range */ 1449 1450 if (ws.subinterval_too_small(a1, a2, b2)) { 1451 error_type = 3; 1452 } 1453 } 1454 1455 ws.update(a1, b1, area1, error1, a2, b2, area2, error2); 1456 wsObj = ws.retrieve(); 1457 a_i = wsObj.a_i; 1458 b_i = wsObj.b_i; 1459 r_i = wsObj.r_i; 1460 e_i = wsObj.e_i; 1461 1462 iteration++; 1463 } while (iteration < limit && !error_type && errsum > tolerance); 1464 1465 result = ws.sum_results(); 1466 // abserr = errsum; 1467 /* 1468 if (errsum <= tolerance) 1469 { 1470 return GSL_SUCCESS; 1471 } 1472 else if (error_type == 2) 1473 { 1474 GSL_ERROR ("roundoff error prevents tolerance from being achieved", 1475 GSL_EROUND); 1476 } 1477 else if (error_type == 3) 1478 { 1479 GSL_ERROR ("bad integrand behavior found in the integration interval", 1480 GSL_ESING); 1481 } 1482 else if (iteration == limit) 1483 { 1484 GSL_ERROR ("maximum number of subdivisions reached", GSL_EMAXITER); 1485 } 1486 else 1487 { 1488 GSL_ERROR ("could not integrate function", GSL_EFAILED); 1489 } 1490 */ 1491 1492 return result; 1493 }, 1494 1495 /** 1496 * Integral of function f over interval. 1497 * @param {Array} interval The integration interval, e.g. [0, 3]. 1498 * @param {function} f A function which takes one argument of type number and returns a number. 1499 * @returns {Number} The value of the integral of f over interval 1500 * @see JXG.Math.Numerics.NewtonCotes 1501 * @see JXG.Math.Numerics.Romberg 1502 * @see JXG.Math.Numerics.Qag 1503 * @memberof JXG.Math.Numerics 1504 */ 1505 I: function (interval, f) { 1506 // return this.NewtonCotes(interval, f, {number_of_nodes: 16, integration_type: 'milne'}); 1507 // return this.Romberg(interval, f, {max_iterations: 20, eps: 0.0000001}); 1508 return this.Qag(interval, f, { 1509 q: this.GaussKronrod15, 1510 limit: 15, 1511 epsrel: 0.0000001, 1512 epsabs: 0.0000001 1513 }); 1514 }, 1515 1516 /** 1517 * Newton's method to find roots of a funtion in one variable. 1518 * @param {function} f We search for a solution of f(x)=0. 1519 * @param {Number} x initial guess for the root, i.e. start value. 1520 * @param {Object} context optional object that is treated as "this" in the function body. This is useful if 1521 * the function is a method of an object and contains a reference to its parent object via "this". 1522 * @returns {Number} A root of the function f. 1523 * @memberof JXG.Math.Numerics 1524 */ 1525 Newton: function (f, x, context) { 1526 var df, 1527 i = 0, 1528 h = Mat.eps, 1529 newf = f.apply(context, [x]); 1530 // nfev = 1; 1531 1532 // For compatibility 1533 if (Type.isArray(x)) { 1534 x = x[0]; 1535 } 1536 1537 while (i < 50 && Math.abs(newf) > h) { 1538 df = this.D(f, context)(x); 1539 // nfev += 2; 1540 1541 if (Math.abs(df) > h) { 1542 x -= newf / df; 1543 } else { 1544 x += Math.random() * 0.2 - 1.0; 1545 } 1546 1547 newf = f.apply(context, [x]); 1548 // nfev += 1; 1549 i += 1; 1550 } 1551 1552 return x; 1553 }, 1554 1555 /** 1556 * Abstract method to find roots of univariate functions, which - for the time being - 1557 * is an alias for {@link JXG.Math.Numerics.chandrupatla}. 1558 * @param {function} f We search for a solution of f(x)=0. 1559 * @param {Number|Array} x initial guess for the root, i.e. starting value, or start interval enclosing the root. 1560 * @param {Object} context optional object that is treated as "this" in the function body. This is useful if 1561 * the function is a method of an object and contains a reference to its parent object via "this". 1562 * @returns {Number} A root of the function f. 1563 * 1564 * @see JXG.Math.Numerics.chandrupatla 1565 * @see JXG.Math.Numerics.fzero 1566 * @memberof JXG.Math.Numerics 1567 */ 1568 root: function (f, x, context) { 1569 //return this.fzero(f, x, context); 1570 return this.chandrupatla(f, x, context); 1571 }, 1572 1573 /** 1574 * Compute an intersection of the curves c1 and c2 1575 * with a generalized Newton method. 1576 * We want to find values t1, t2 such that 1577 * c1(t1) = c2(t2), i.e. 1578 * (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2)) = (0,0). 1579 * We set 1580 * (e,f) := (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2)) 1581 * 1582 * The Jacobian J is defined by 1583 * J = (a, b) 1584 * (c, d) 1585 * where 1586 * a = c1_x'(t1) 1587 * b = -c2_x'(t2) 1588 * c = c1_y'(t1) 1589 * d = -c2_y'(t2) 1590 * 1591 * The inverse J^(-1) of J is equal to 1592 * (d, -b)/ 1593 * (-c, a) / (ad-bc) 1594 * 1595 * Then, (t1new, t2new) := (t1,t2) - J^(-1)*(e,f). 1596 * If the function meetCurveCurve possesses the properties 1597 * t1memo and t2memo then these are taken as start values 1598 * for the Newton algorithm. 1599 * After stopping of the Newton algorithm the values of t1 and t2 are stored in 1600 * t1memo and t2memo. 1601 * 1602 * @param {JXG.Curve} c1 Curve, Line or Circle 1603 * @param {JXG.Curve} c2 Curve, Line or Circle 1604 * @param {Number} t1ini start value for t1 1605 * @param {Number} t2ini start value for t2 1606 * @returns {JXG.Coords} intersection point 1607 * @memberof JXG.Math.Numerics 1608 */ 1609 generalizedNewton: function (c1, c2, t1ini, t2ini) { 1610 var t1, 1611 t2, 1612 a, 1613 b, 1614 c, 1615 d, 1616 disc, 1617 e, 1618 f, 1619 F, 1620 D00, 1621 D01, 1622 D10, 1623 D11, 1624 count = 0; 1625 1626 if (this.generalizedNewton.t1memo) { 1627 t1 = this.generalizedNewton.t1memo; 1628 t2 = this.generalizedNewton.t2memo; 1629 } else { 1630 t1 = t1ini; 1631 t2 = t2ini; 1632 } 1633 1634 e = c1.X(t1) - c2.X(t2); 1635 f = c1.Y(t1) - c2.Y(t2); 1636 F = e * e + f * f; 1637 1638 D00 = this.D(c1.X, c1); 1639 D01 = this.D(c2.X, c2); 1640 D10 = this.D(c1.Y, c1); 1641 D11 = this.D(c2.Y, c2); 1642 1643 while (F > Mat.eps && count < 10) { 1644 a = D00(t1); 1645 b = -D01(t2); 1646 c = D10(t1); 1647 d = -D11(t2); 1648 disc = a * d - b * c; 1649 t1 -= (d * e - b * f) / disc; 1650 t2 -= (a * f - c * e) / disc; 1651 e = c1.X(t1) - c2.X(t2); 1652 f = c1.Y(t1) - c2.Y(t2); 1653 F = e * e + f * f; 1654 count += 1; 1655 } 1656 1657 this.generalizedNewton.t1memo = t1; 1658 this.generalizedNewton.t2memo = t2; 1659 1660 if (Math.abs(t1) < Math.abs(t2)) { 1661 return [c1.X(t1), c1.Y(t1)]; 1662 } 1663 1664 return [c2.X(t2), c2.Y(t2)]; 1665 }, 1666 1667 /** 1668 * Returns the Lagrange polynomials for curves with equidistant nodes, see 1669 * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation, 1670 * SIAM Review, Vol 46, No 3, (2004) 501-517. 1671 * The graph of the parametric curve [x(t),y(t)] runs through the given points. 1672 * @param {Array} p Array of JXG.Points 1673 * @returns {Array} An array consisting of two functions x(t), y(t) which define a parametric curve 1674 * f(t) = (x(t), y(t)), a number x1 (which equals 0) and a function x2 defining the curve's domain. 1675 * That means the curve is defined between x1 and x2(). x2 returns the (length of array p minus one). 1676 * @memberof JXG.Math.Numerics 1677 * 1678 * @example 1679 * var p = []; 1680 * 1681 * p[0] = board.create('point', [0, -2], {size:2, name: 'C(a)'}); 1682 * p[1] = board.create('point', [-1.5, 5], {size:2, name: ''}); 1683 * p[2] = board.create('point', [1, 4], {size:2, name: ''}); 1684 * p[3] = board.create('point', [3, 3], {size:2, name: 'C(b)'}); 1685 * 1686 * // Curve 1687 * var fg = JXG.Math.Numerics.Neville(p); 1688 * var graph = board.create('curve', fg, {strokeWidth:3, strokeOpacity:0.5}); 1689 * 1690 * </pre><div id="JXG88a8b3a8-6561-44f5-a678-76bca13fd484" class="jxgbox" style="width: 300px; height: 300px;"></div> 1691 * <script type="text/javascript"> 1692 * (function() { 1693 * var board = JXG.JSXGraph.initBoard('JXG88a8b3a8-6561-44f5-a678-76bca13fd484', 1694 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 1695 * var p = []; 1696 * 1697 * p[0] = board.create('point', [0, -2], {size:2, name: 'C(a)'}); 1698 * p[1] = board.create('point', [-1.5, 5], {size:2, name: ''}); 1699 * p[2] = board.create('point', [1, 4], {size:2, name: ''}); 1700 * p[3] = board.create('point', [3, 3], {size:2, name: 'C(b)'}); 1701 * 1702 * // Curve 1703 * var fg = JXG.Math.Numerics.Neville(p); 1704 * var graph = board.create('curve', fg, {strokeWidth:3, strokeOpacity:0.5}); 1705 * 1706 * })(); 1707 * 1708 * </script><pre> 1709 * 1710 */ 1711 Neville: function (p) { 1712 var w = [], 1713 /** @ignore */ 1714 makeFct = function (fun) { 1715 return function (t, suspendedUpdate) { 1716 var i, 1717 d, 1718 s, 1719 bin = Mat.binomial, 1720 len = p.length, 1721 len1 = len - 1, 1722 num = 0.0, 1723 denom = 0.0; 1724 1725 if (!suspendedUpdate) { 1726 s = 1; 1727 for (i = 0; i < len; i++) { 1728 w[i] = bin(len1, i) * s; 1729 s *= -1; 1730 } 1731 } 1732 1733 d = t; 1734 1735 for (i = 0; i < len; i++) { 1736 if (d === 0) { 1737 return p[i][fun](); 1738 } 1739 s = w[i] / d; 1740 d -= 1; 1741 num += p[i][fun]() * s; 1742 denom += s; 1743 } 1744 return num / denom; 1745 }; 1746 }, 1747 xfct = makeFct("X"), 1748 yfct = makeFct("Y"); 1749 1750 return [ 1751 xfct, 1752 yfct, 1753 0, 1754 function () { 1755 return p.length - 1; 1756 } 1757 ]; 1758 }, 1759 1760 /** 1761 * Calculates second derivatives at the given knots. 1762 * @param {Array} x x values of knots 1763 * @param {Array} y y values of knots 1764 * @returns {Array} Second derivatives of the interpolated function at the knots. 1765 * @see #splineEval 1766 * @memberof JXG.Math.Numerics 1767 */ 1768 splineDef: function (x, y) { 1769 var pair, 1770 i, 1771 l, 1772 n = Math.min(x.length, y.length), 1773 diag = [], 1774 z = [], 1775 data = [], 1776 dx = [], 1777 delta = [], 1778 F = []; 1779 1780 if (n === 2) { 1781 return [0, 0]; 1782 } 1783 1784 for (i = 0; i < n; i++) { 1785 pair = { X: x[i], Y: y[i] }; 1786 data.push(pair); 1787 } 1788 data.sort(function (a, b) { 1789 return a.X - b.X; 1790 }); 1791 for (i = 0; i < n; i++) { 1792 x[i] = data[i].X; 1793 y[i] = data[i].Y; 1794 } 1795 1796 for (i = 0; i < n - 1; i++) { 1797 dx.push(x[i + 1] - x[i]); 1798 } 1799 for (i = 0; i < n - 2; i++) { 1800 delta.push( 1801 (6 * (y[i + 2] - y[i + 1])) / dx[i + 1] - (6 * (y[i + 1] - y[i])) / dx[i] 1802 ); 1803 } 1804 1805 // ForwardSolve 1806 diag.push(2 * (dx[0] + dx[1])); 1807 z.push(delta[0]); 1808 1809 for (i = 0; i < n - 3; i++) { 1810 l = dx[i + 1] / diag[i]; 1811 diag.push(2 * (dx[i + 1] + dx[i + 2]) - l * dx[i + 1]); 1812 z.push(delta[i + 1] - l * z[i]); 1813 } 1814 1815 // BackwardSolve 1816 F[n - 3] = z[n - 3] / diag[n - 3]; 1817 for (i = n - 4; i >= 0; i--) { 1818 F[i] = (z[i] - dx[i + 1] * F[i + 1]) / diag[i]; 1819 } 1820 1821 // Generate f''-Vector 1822 for (i = n - 3; i >= 0; i--) { 1823 F[i + 1] = F[i]; 1824 } 1825 1826 // natural cubic spline 1827 F[0] = 0; 1828 F[n - 1] = 0; 1829 1830 return F; 1831 }, 1832 1833 /** 1834 * Evaluate points on spline. 1835 * @param {Number,Array} x0 A single float value or an array of values to evaluate 1836 * @param {Array} x x values of knots 1837 * @param {Array} y y values of knots 1838 * @param {Array} F Second derivatives at knots, calculated by {@link JXG.Math.Numerics.splineDef} 1839 * @see #splineDef 1840 * @returns {Number,Array} A single value or an array, depending on what is given as x0. 1841 * @memberof JXG.Math.Numerics 1842 */ 1843 splineEval: function (x0, x, y, F) { 1844 var i, 1845 j, 1846 a, 1847 b, 1848 c, 1849 d, 1850 x_, 1851 n = Math.min(x.length, y.length), 1852 l = 1, 1853 asArray = false, 1854 y0 = []; 1855 1856 // number of points to be evaluated 1857 if (Type.isArray(x0)) { 1858 l = x0.length; 1859 asArray = true; 1860 } else { 1861 x0 = [x0]; 1862 } 1863 1864 for (i = 0; i < l; i++) { 1865 // is x0 in defining interval? 1866 if (x0[i] < x[0] || x[i] > x[n - 1]) { 1867 return NaN; 1868 } 1869 1870 // determine part of spline in which x0 lies 1871 for (j = 1; j < n; j++) { 1872 if (x0[i] <= x[j]) { 1873 break; 1874 } 1875 } 1876 1877 j -= 1; 1878 1879 // we're now in the j-th partial interval, i.e. x[j] < x0[i] <= x[j+1]; 1880 // determine the coefficients of the polynomial in this interval 1881 a = y[j]; 1882 b = 1883 (y[j + 1] - y[j]) / (x[j + 1] - x[j]) - 1884 ((x[j + 1] - x[j]) / 6) * (F[j + 1] + 2 * F[j]); 1885 c = F[j] / 2; 1886 d = (F[j + 1] - F[j]) / (6 * (x[j + 1] - x[j])); 1887 // evaluate x0[i] 1888 x_ = x0[i] - x[j]; 1889 //y0.push(a + b*x_ + c*x_*x_ + d*x_*x_*x_); 1890 y0.push(a + (b + (c + d * x_) * x_) * x_); 1891 } 1892 1893 if (asArray) { 1894 return y0; 1895 } 1896 1897 return y0[0]; 1898 }, 1899 1900 /** 1901 * Generate a string containing the function term of a polynomial. 1902 * @param {Array} coeffs Coefficients of the polynomial. The position i belongs to x^i. 1903 * @param {Number} deg Degree of the polynomial 1904 * @param {String} varname Name of the variable (usually 'x') 1905 * @param {Number} prec Precision 1906 * @returns {String} A string containg the function term of the polynomial. 1907 * @memberof JXG.Math.Numerics 1908 */ 1909 generatePolynomialTerm: function (coeffs, deg, varname, prec) { 1910 var i, 1911 t = []; 1912 1913 for (i = deg; i >= 0; i--) { 1914 t = t.concat(["(", coeffs[i].toPrecision(prec), ")"]); 1915 1916 if (i > 1) { 1917 t = t.concat(["*", varname, "<sup>", i, "<", "/sup> + "]); 1918 } else if (i === 1) { 1919 t = t.concat(["*", varname, " + "]); 1920 } 1921 } 1922 1923 return t.join(""); 1924 }, 1925 1926 /** 1927 * Computes the polynomial through a given set of coordinates in Lagrange form. 1928 * Returns the Lagrange polynomials, see 1929 * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation, 1930 * SIAM Review, Vol 46, No 3, (2004) 501-517. 1931 * <p> 1932 * It possesses the method getTerm() which returns the string containing the function term of the polynomial. 1933 * @param {Array} p Array of JXG.Points 1934 * @returns {function} A function of one parameter which returns the value of the polynomial, whose graph runs through the given points. 1935 * @memberof JXG.Math.Numerics 1936 * 1937 * @example 1938 * var p = []; 1939 * p[0] = board.create('point', [-1,2], {size:4}); 1940 * p[1] = board.create('point', [0,3], {size:4}); 1941 * p[2] = board.create('point', [1,1], {size:4}); 1942 * p[3] = board.create('point', [3,-1], {size:4}); 1943 * var f = JXG.Math.Numerics.lagrangePolynomial(p); 1944 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 1945 * 1946 * </pre><div id="JXGc058aa6b-74d4-41e1-af94-df06169a2d89" class="jxgbox" style="width: 300px; height: 300px;"></div> 1947 * <script type="text/javascript"> 1948 * (function() { 1949 * var board = JXG.JSXGraph.initBoard('JXGc058aa6b-74d4-41e1-af94-df06169a2d89', 1950 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 1951 * var p = []; 1952 * p[0] = board.create('point', [-1,2], {size:4}); 1953 * p[1] = board.create('point', [0,3], {size:4}); 1954 * p[2] = board.create('point', [1,1], {size:4}); 1955 * p[3] = board.create('point', [3,-1], {size:4}); 1956 * var f = JXG.Math.Numerics.lagrangePolynomial(p); 1957 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 1958 * 1959 * })(); 1960 * 1961 * </script><pre> 1962 * 1963 * @example 1964 * var points = []; 1965 * points[0] = board.create('point', [-1,2], {size:4}); 1966 * points[1] = board.create('point', [0, 0], {size:4}); 1967 * points[2] = board.create('point', [2, 1], {size:4}); 1968 * 1969 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 1970 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 1971 * var txt = board.create('text', [-3, -4, () => f.getTerm(2, 't', ' * ')], {fontSize: 16}); 1972 * var txt2 = board.create('text', [-3, -6, () => f.getCoefficients()], {fontSize: 12}); 1973 * 1974 * </pre><div id="JXG73fdaf12-e257-4374-b488-ae063e4eecbb" class="jxgbox" style="width: 300px; height: 300px;"></div> 1975 * <script type="text/javascript"> 1976 * (function() { 1977 * var board = JXG.JSXGraph.initBoard('JXG73fdaf12-e257-4374-b488-ae063e4eecbb', 1978 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 1979 * var points = []; 1980 * points[0] = board.create('point', [-1,2], {size:4}); 1981 * points[1] = board.create('point', [0, 0], {size:4}); 1982 * points[2] = board.create('point', [2, 1], {size:4}); 1983 * 1984 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 1985 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 1986 * var txt = board.create('text', [-3, -4, () => f.getTerm(2, 't', ' * ')], {fontSize: 16}); 1987 * var txt2 = board.create('text', [-3, -6, () => f.getCoefficients()], {fontSize: 12}); 1988 * 1989 * })(); 1990 * 1991 * </script><pre> 1992 * 1993 */ 1994 lagrangePolynomial: function (p) { 1995 var w = [], 1996 that = this, 1997 /** @ignore */ 1998 fct = function (x, suspendedUpdate) { 1999 var i, // j, 2000 k, 2001 xi, 2002 s, //M, 2003 len = p.length, 2004 num = 0, 2005 denom = 0; 2006 2007 if (!suspendedUpdate) { 2008 for (i = 0; i < len; i++) { 2009 w[i] = 1.0; 2010 xi = p[i].X(); 2011 2012 for (k = 0; k < len; k++) { 2013 if (k !== i) { 2014 w[i] *= xi - p[k].X(); 2015 } 2016 } 2017 2018 w[i] = 1 / w[i]; 2019 } 2020 2021 // M = []; 2022 // for (k = 0; k < len; k++) { 2023 // M.push([1]); 2024 // } 2025 } 2026 2027 for (i = 0; i < len; i++) { 2028 xi = p[i].X(); 2029 2030 if (x === xi) { 2031 return p[i].Y(); 2032 } 2033 2034 s = w[i] / (x - xi); 2035 denom += s; 2036 num += s * p[i].Y(); 2037 } 2038 2039 return num / denom; 2040 }; 2041 2042 /** 2043 * Get the term of the Lagrange polynomial as string. 2044 * Calls {@link JXG.Math.Numerics#lagrangePolynomialTerm}. 2045 * 2046 * @name JXG.Math.Numerics#lagrangePolynomial.getTerm 2047 * @param {Number} digits Number of digits of the coefficients 2048 * @param {String} param Variable name 2049 * @param {String} dot Dot symbol 2050 * @returns {String} containing the term of Lagrange polynomial as string. 2051 * @see JXG.Math.Numerics#lagrangePolynomialTerm 2052 * @example 2053 * var points = []; 2054 * points[0] = board.create('point', [-1,2], {size:4}); 2055 * points[1] = board.create('point', [0, 0], {size:4}); 2056 * points[2] = board.create('point', [2, 1], {size:4}); 2057 * 2058 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2059 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2060 * var txt = board.create('text', [-3, -4, () => f.getTerm(2, 't', ' * ')], {fontSize: 16}); 2061 * 2062 * </pre><div id="JXG73fdaf12-e257-4374-b488-ae063e4eeccf" class="jxgbox" style="width: 300px; height: 300px;"></div> 2063 * <script type="text/javascript"> 2064 * (function() { 2065 * var board = JXG.JSXGraph.initBoard('JXG73fdaf12-e257-4374-b488-ae063e4eeccf', 2066 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 2067 * var points = []; 2068 * points[0] = board.create('point', [-1,2], {size:4}); 2069 * points[1] = board.create('point', [0, 0], {size:4}); 2070 * points[2] = board.create('point', [2, 1], {size:4}); 2071 * 2072 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2073 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2074 * var txt = board.create('text', [-3, -4, () => f.getTerm(2, 't', ' * ')], {fontSize: 16}); 2075 * 2076 * })(); 2077 * 2078 * </script><pre> 2079 * 2080 */ 2081 fct.getTerm = function (digits, param, dot) { 2082 return that.lagrangePolynomialTerm(p, digits, param, dot)(); 2083 }; 2084 2085 /** 2086 * Get the coefficients of the Lagrange polynomial as array. The leading 2087 * coefficient is at position 0. 2088 * Calls {@link JXG.Math.Numerics#lagrangePolynomialCoefficients}. 2089 * 2090 * @name JXG.Math.Numerics#lagrangePolynomial.getCoefficients 2091 * @returns {Array} containing the coefficients of the Lagrange polynomial. 2092 * @see JXG.Math.Numerics#lagrangePolynomial.getTerm 2093 * @see JXG.Math.Numerics#lagrangePolynomialTerm 2094 * @see JXG.Math.Numerics#lagrangePolynomialCoefficients 2095 * @example 2096 * var points = []; 2097 * points[0] = board.create('point', [-1,2], {size:4}); 2098 * points[1] = board.create('point', [0, 0], {size:4}); 2099 * points[2] = board.create('point', [2, 1], {size:4}); 2100 * 2101 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2102 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2103 * var txt = board.create('text', [1, -4, () => f.getCoefficients()], {fontSize: 10}); 2104 * 2105 * </pre><div id="JXG52a883a5-2e0c-4caf-8f84-8650c173c365" class="jxgbox" style="width: 300px; height: 300px;"></div> 2106 * <script type="text/javascript"> 2107 * (function() { 2108 * var board = JXG.JSXGraph.initBoard('JXG52a883a5-2e0c-4caf-8f84-8650c173c365', 2109 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 2110 * var points = []; 2111 * points[0] = board.create('point', [-1,2], {size:4}); 2112 * points[1] = board.create('point', [0, 0], {size:4}); 2113 * points[2] = board.create('point', [2, 1], {size:4}); 2114 * 2115 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2116 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2117 * var txt = board.create('text', [1, -4, () => f.getCoefficients()], {fontSize: 10}); 2118 * 2119 * })(); 2120 * 2121 * </script><pre> 2122 * 2123 */ 2124 fct.getCoefficients = function () { 2125 return that.lagrangePolynomialCoefficients(p)(); 2126 }; 2127 2128 return fct; 2129 }, 2130 2131 /** 2132 * Determine the Lagrange polynomial through an array of points and 2133 * return the term of the polynomial as string. 2134 * 2135 * @param {Array} points Array of JXG.Points 2136 * @param {Number} digits Number of decimal digits of the coefficients 2137 * @param {String} param Name of the parameter. Default: 'x'. 2138 * @param {String} dot Multiplication symbol. Default: ' * '. 2139 * @returns {Function} returning the Lagrange polynomial term through 2140 * the supplied points as string 2141 * @memberof JXG.Math.Numerics 2142 * 2143 * @example 2144 * var points = []; 2145 * points[0] = board.create('point', [-1,2], {size:4}); 2146 * points[1] = board.create('point', [0, 0], {size:4}); 2147 * points[2] = board.create('point', [2, 1], {size:4}); 2148 * 2149 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2150 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2151 * 2152 * var f_txt = JXG.Math.Numerics.lagrangePolynomialTerm(points, 2, 't', ' * '); 2153 * var txt = board.create('text', [-3, -4, f_txt], {fontSize: 16}); 2154 * 2155 * </pre><div id="JXGd45e9e96-7526-486d-aa43-e1178d5f2baa" class="jxgbox" style="width: 300px; height: 300px;"></div> 2156 * <script type="text/javascript"> 2157 * (function() { 2158 * var board = JXG.JSXGraph.initBoard('JXGd45e9e96-7526-486d-aa43-e1178d5f2baa', 2159 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 2160 * var points = []; 2161 * points[0] = board.create('point', [-1,2], {size:4}); 2162 * points[1] = board.create('point', [0, 0], {size:4}); 2163 * points[2] = board.create('point', [2, 1], {size:4}); 2164 * 2165 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2166 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2167 * 2168 * var f_txt = JXG.Math.Numerics.lagrangePolynomialTerm(points, 2, 't', ' * '); 2169 * var txt = board.create('text', [-3, -4, f_txt], {fontSize: 16}); 2170 * 2171 * })(); 2172 * 2173 * </script><pre> 2174 * 2175 */ 2176 lagrangePolynomialTerm: function (points, digits, param, dot) { 2177 var that = this; 2178 2179 return function () { 2180 var len = points.length, 2181 coeffs = [], 2182 isLeading = true, 2183 n, t, j, c; 2184 2185 param = param || "x"; 2186 if (dot === undefined) { 2187 dot = " * "; 2188 } 2189 2190 n = len - 1; // (Max) degree of the polynomial 2191 coeffs = that.lagrangePolynomialCoefficients(points)(); 2192 2193 t = ""; 2194 for (j = 0; j < coeffs.length; j++) { 2195 c = coeffs[j]; 2196 if (Math.abs(c) < Mat.eps) { 2197 continue; 2198 } 2199 if (JXG.exists(digits)) { 2200 c = Env._round10(c, -digits); 2201 } 2202 if (isLeading) { 2203 t += c > 0 ? c : "-" + -c; 2204 isLeading = false; 2205 } else { 2206 t += c > 0 ? " + " + c : " - " + -c; 2207 } 2208 2209 if (n - j > 1) { 2210 t += dot + param + "^" + (n - j); 2211 } else if (n - j === 1) { 2212 t += dot + param; 2213 } 2214 } 2215 return t; // board.jc.manipulate('f = map(x) -> ' + t + ';'); 2216 }; 2217 }, 2218 2219 /** 2220 * Determine the Lagrange polynomial through an array of points and 2221 * return the coefficients of the polynomial as array. 2222 * The leading coefficient is at position 0. 2223 * 2224 * @param {Array} points Array of JXG.Points 2225 * @returns {Function} returning the coefficients of the Lagrange polynomial through 2226 * the supplied points. 2227 * @memberof JXG.Math.Numerics 2228 * 2229 * @example 2230 * var points = []; 2231 * points[0] = board.create('point', [-1,2], {size:4}); 2232 * points[1] = board.create('point', [0, 0], {size:4}); 2233 * points[2] = board.create('point', [2, 1], {size:4}); 2234 * 2235 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2236 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2237 * 2238 * var f_arr = JXG.Math.Numerics.lagrangePolynomialCoefficients(points); 2239 * var txt = board.create('text', [1, -4, f_arr], {fontSize: 10}); 2240 * 2241 * </pre><div id="JXG1778f0d1-a420-473f-99e8-1755ef4be97e" class="jxgbox" style="width: 300px; height: 300px;"></div> 2242 * <script type="text/javascript"> 2243 * (function() { 2244 * var board = JXG.JSXGraph.initBoard('JXG1778f0d1-a420-473f-99e8-1755ef4be97e', 2245 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 2246 * var points = []; 2247 * points[0] = board.create('point', [-1,2], {size:4}); 2248 * points[1] = board.create('point', [0, 0], {size:4}); 2249 * points[2] = board.create('point', [2, 1], {size:4}); 2250 * 2251 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2252 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2253 * 2254 * var f_arr = JXG.Math.Numerics.lagrangePolynomialCoefficients(points); 2255 * var txt = board.create('text', [1, -4, f_arr], {fontSize: 10}); 2256 * 2257 * })(); 2258 * 2259 * </script><pre> 2260 * 2261 */ 2262 lagrangePolynomialCoefficients: function (points) { 2263 return function () { 2264 var len = points.length, 2265 zeroes = [], 2266 coeffs = [], 2267 coeffs_sum = [], 2268 n, i, j, c, p; 2269 2270 n = len - 1; // (Max) degree of the polynomial 2271 for (j = 0; j < len; j++) { 2272 coeffs_sum[j] = 0; 2273 } 2274 2275 for (i = 0; i < len; i++) { 2276 c = points[i].Y(); 2277 p = points[i].X(); 2278 zeroes = []; 2279 for (j = 0; j < len; j++) { 2280 if (j !== i) { 2281 c /= p - points[j].X(); 2282 zeroes.push(points[j].X()); 2283 } 2284 } 2285 coeffs = [1].concat(Mat.Vieta(zeroes)); 2286 for (j = 0; j < coeffs.length; j++) { 2287 coeffs_sum[j] += (j % 2 === 1 ? -1 : 1) * coeffs[j] * c; 2288 } 2289 } 2290 2291 return coeffs_sum; 2292 }; 2293 }, 2294 2295 /** 2296 * Determine the coefficients of a cardinal spline polynom, See 2297 * https://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections 2298 * @param {Number} x1 point 1 2299 * @param {Number} x2 point 2 2300 * @param {Number} t1 tangent slope 1 2301 * @param {Number} t2 tangent slope 2 2302 * @return {Array} coefficents array c for the polynomial t maps to 2303 * c[0] + c[1]*t + c[2]*t*t + c[3]*t*t*t 2304 */ 2305 _initCubicPoly: function (x1, x2, t1, t2) { 2306 return [x1, t1, -3 * x1 + 3 * x2 - 2 * t1 - t2, 2 * x1 - 2 * x2 + t1 + t2]; 2307 }, 2308 2309 /** 2310 * Computes the cubic cardinal spline curve through a given set of points. The curve 2311 * is uniformly parametrized. 2312 * Two artificial control points at the beginning and the end are added. 2313 * 2314 * The implementation (especially the centripetal parametrization) is from 2315 * https://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections . 2316 * @param {Array} points Array consisting of JXG.Points. 2317 * @param {Number|Function} tau The tension parameter, either a constant number or a function returning a number. This number is between 0 and 1. 2318 * tau=1/2 give Catmull-Rom splines. 2319 * @param {String} type (Optional) parameter which allows to choose between "uniform" (default) and 2320 * "centripetal" parameterization. Thus the two possible values are "uniform" or "centripetal". 2321 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 2322 * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, 2323 * and a function simply returning the length of the points array 2324 * minus three. 2325 * @memberof JXG.Math.Numerics 2326 */ 2327 CardinalSpline: function (points, tau_param, type) { 2328 var p, 2329 coeffs = [], 2330 makeFct, 2331 tau, _tau, 2332 that = this; 2333 2334 if (Type.isFunction(tau_param)) { 2335 _tau = tau_param; 2336 } else { 2337 _tau = function () { 2338 return tau_param; 2339 }; 2340 } 2341 2342 if (type === undefined) { 2343 type = "uniform"; 2344 } 2345 2346 /** @ignore */ 2347 makeFct = function (which) { 2348 return function (t, suspendedUpdate) { 2349 var s, 2350 c, 2351 // control point at the beginning and at the end 2352 first, 2353 last, 2354 t1, 2355 t2, 2356 dt0, 2357 dt1, 2358 dt2, 2359 // dx, dy, 2360 len; 2361 2362 if (points.length < 2) { 2363 return NaN; 2364 } 2365 2366 if (!suspendedUpdate) { 2367 tau = _tau(); 2368 2369 // New point list p: [first, points ..., last] 2370 first = { 2371 X: function () { 2372 return 2 * points[0].X() - points[1].X(); 2373 }, 2374 Y: function () { 2375 return 2 * points[0].Y() - points[1].Y(); 2376 }, 2377 Dist: function (p) { 2378 var dx = this.X() - p.X(), 2379 dy = this.Y() - p.Y(); 2380 return Math.sqrt(dx * dx + dy * dy); 2381 } 2382 }; 2383 2384 last = { 2385 X: function () { 2386 return ( 2387 2 * points[points.length - 1].X() - 2388 points[points.length - 2].X() 2389 ); 2390 }, 2391 Y: function () { 2392 return ( 2393 2 * points[points.length - 1].Y() - 2394 points[points.length - 2].Y() 2395 ); 2396 }, 2397 Dist: function (p) { 2398 var dx = this.X() - p.X(), 2399 dy = this.Y() - p.Y(); 2400 return Math.sqrt(dx * dx + dy * dy); 2401 } 2402 }; 2403 2404 p = [first].concat(points, [last]); 2405 len = p.length; 2406 2407 coeffs[which] = []; 2408 2409 for (s = 0; s < len - 3; s++) { 2410 if (type === "centripetal") { 2411 // The order is important, since p[0].coords === undefined 2412 dt0 = p[s].Dist(p[s + 1]); 2413 dt1 = p[s + 2].Dist(p[s + 1]); 2414 dt2 = p[s + 3].Dist(p[s + 2]); 2415 2416 dt0 = Math.sqrt(dt0); 2417 dt1 = Math.sqrt(dt1); 2418 dt2 = Math.sqrt(dt2); 2419 2420 if (dt1 < Mat.eps) { 2421 dt1 = 1.0; 2422 } 2423 if (dt0 < Mat.eps) { 2424 dt0 = dt1; 2425 } 2426 if (dt2 < Mat.eps) { 2427 dt2 = dt1; 2428 } 2429 2430 t1 = 2431 (p[s + 1][which]() - p[s][which]()) / dt0 - 2432 (p[s + 2][which]() - p[s][which]()) / (dt1 + dt0) + 2433 (p[s + 2][which]() - p[s + 1][which]()) / dt1; 2434 2435 t2 = 2436 (p[s + 2][which]() - p[s + 1][which]()) / dt1 - 2437 (p[s + 3][which]() - p[s + 1][which]()) / (dt2 + dt1) + 2438 (p[s + 3][which]() - p[s + 2][which]()) / dt2; 2439 2440 t1 *= dt1; 2441 t2 *= dt1; 2442 2443 coeffs[which][s] = that._initCubicPoly( 2444 p[s + 1][which](), 2445 p[s + 2][which](), 2446 tau * t1, 2447 tau * t2 2448 ); 2449 } else { 2450 coeffs[which][s] = that._initCubicPoly( 2451 p[s + 1][which](), 2452 p[s + 2][which](), 2453 tau * (p[s + 2][which]() - p[s][which]()), 2454 tau * (p[s + 3][which]() - p[s + 1][which]()) 2455 ); 2456 } 2457 } 2458 } 2459 2460 if (isNaN(t)) { 2461 return NaN; 2462 } 2463 2464 len = points.length; 2465 // This is necessary for our advanced plotting algorithm: 2466 if (t <= 0.0) { 2467 return points[0][which](); 2468 } 2469 if (t >= len) { 2470 return points[len - 1][which](); 2471 } 2472 2473 s = Math.floor(t); 2474 if (s === t) { 2475 return points[s][which](); 2476 } 2477 2478 t -= s; 2479 c = coeffs[which][s]; 2480 if (c === undefined) { 2481 return NaN; 2482 } 2483 2484 return ((c[3] * t + c[2]) * t + c[1]) * t + c[0]; 2485 }; 2486 }; 2487 2488 return [ 2489 makeFct("X"), 2490 makeFct("Y"), 2491 0, 2492 function () { 2493 return points.length - 1; 2494 } 2495 ]; 2496 }, 2497 2498 /** 2499 * Computes the cubic Catmull-Rom spline curve through a given set of points. The curve 2500 * is uniformly parametrized. The curve is the cardinal spline curve for tau=0.5. 2501 * Two artificial control points at the beginning and the end are added. 2502 * @param {Array} points Array consisting of JXG.Points. 2503 * @param {String} type (Optional) parameter which allows to choose between "uniform" (default) and 2504 * "centripetal" parameterization. Thus the two possible values are "uniform" or "centripetal". 2505 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 2506 * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, and a function simply 2507 * returning the length of the points array minus three. 2508 * @memberof JXG.Math.Numerics 2509 */ 2510 CatmullRomSpline: function (points, type) { 2511 return this.CardinalSpline(points, 0.5, type); 2512 }, 2513 2514 /** 2515 * Computes the regression polynomial of a given degree through a given set of coordinates. 2516 * Returns the regression polynomial function. 2517 * @param {Number,function,Slider} degree number, function or slider. 2518 * Either 2519 * @param {Array} dataX Array containing either the x-coordinates of the data set or both coordinates in 2520 * an array of {@link JXG.Point}s or {@link JXG.Coords}. 2521 * In the latter case, the <tt>dataY</tt> parameter will be ignored. 2522 * @param {Array} dataY Array containing the y-coordinates of the data set, 2523 * @returns {function} A function of one parameter which returns the value of the regression polynomial of the given degree. 2524 * It possesses the method getTerm() which returns the string containing the function term of the polynomial. 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 2649 for (i = d - 1; i >= 0; i--) { 2650 s = s * x + coeffs[i]; 2651 } 2652 2653 return s; 2654 }; 2655 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 {@see #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', see {@see #riemann}. 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 * 2894 * @memberof JXG.Math.Numerics 2895 */ 2896 _riemannValue: function (x, f, type, delta) { 2897 var y, y1, x1, delta1; 2898 2899 if (delta < 0) { 2900 // delta is negative if the lower function term is evaluated 2901 if (type !== "trapezoidal") { 2902 x = x + delta; 2903 } 2904 delta *= -1; 2905 if (type === "lower") { 2906 type = "upper"; 2907 } else if (type === "upper") { 2908 type = "lower"; 2909 } 2910 } 2911 2912 delta1 = delta * 0.01; // for 'lower' and 'upper' 2913 2914 if (type === "right") { 2915 y = f(x + delta); 2916 } else if (type === "middle") { 2917 y = f(x + delta * 0.5); 2918 } else if (type === "left" || type === "trapezoidal") { 2919 y = f(x); 2920 } else if (type === "lower") { 2921 y = f(x); 2922 2923 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 2924 y1 = f(x1); 2925 2926 if (y1 < y) { 2927 y = y1; 2928 } 2929 } 2930 2931 y1 = f(x + delta); 2932 if (y1 < y) { 2933 y = y1; 2934 } 2935 } else if (type === "upper") { 2936 y = f(x); 2937 2938 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 2939 y1 = f(x1); 2940 if (y1 > y) { 2941 y = y1; 2942 } 2943 } 2944 2945 y1 = f(x + delta); 2946 if (y1 > y) { 2947 y = y1; 2948 } 2949 } else if (type === "random") { 2950 y = f(x + delta * Math.random()); 2951 } else if (type === "simpson") { 2952 y = (f(x) + 4 * f(x + delta * 0.5) + f(x + delta)) / 6.0; 2953 } else { 2954 y = f(x); // default is lower 2955 } 2956 2957 return y; 2958 }, 2959 2960 /** 2961 * Helper function to create curve which displays Riemann sums. 2962 * Compute coordinates for the rectangles showing the Riemann sum. 2963 * @param {Function,Array} f Function or array of two functions. 2964 * If f is a function the integral of this function is approximated by the Riemann sum. 2965 * If f is an array consisting of two functions the area between the two functions is filled 2966 * by the Riemann sum bars. 2967 * @param {Number} n number of rectangles. 2968 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson', or 'trapezoidal'. 2969 * @param {Number} start Left border of the approximation interval 2970 * @param {Number} end Right border of the approximation interval 2971 * @returns {Array} An array of two arrays containing the x and y coordinates for the rectangles showing the Riemann sum. This 2972 * 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 2973 * rectangles. 2974 * @memberof JXG.Math.Numerics 2975 */ 2976 riemann: function (gf, n, type, start, end) { 2977 var i, 2978 delta, 2979 xarr = [], 2980 yarr = [], 2981 j = 0, 2982 x = start, 2983 y, 2984 sum = 0, 2985 f, 2986 g, 2987 ylow, 2988 yup; 2989 2990 if (Type.isArray(gf)) { 2991 g = gf[0]; 2992 f = gf[1]; 2993 } else { 2994 f = gf; 2995 } 2996 2997 n = Math.floor(n); 2998 2999 if (n <= 0) { 3000 return [xarr, yarr, sum]; 3001 } 3002 3003 delta = (end - start) / n; 3004 3005 // Upper bar ends 3006 for (i = 0; i < n; i++) { 3007 y = this._riemannValue(x, f, type, delta); 3008 xarr[j] = x; 3009 yarr[j] = y; 3010 3011 j += 1; 3012 x += delta; 3013 if (type === "trapezoidal") { 3014 y = f(x); 3015 } 3016 xarr[j] = x; 3017 yarr[j] = y; 3018 3019 j += 1; 3020 } 3021 3022 // Lower bar ends 3023 for (i = 0; i < n; i++) { 3024 if (g) { 3025 y = this._riemannValue(x, g, type, -delta); 3026 } else { 3027 y = 0.0; 3028 } 3029 xarr[j] = x; 3030 yarr[j] = y; 3031 3032 j += 1; 3033 x -= delta; 3034 if (type === "trapezoidal" && g) { 3035 y = g(x); 3036 } 3037 xarr[j] = x; 3038 yarr[j] = y; 3039 3040 // Add the area of the bar to 'sum' 3041 if (type !== "trapezoidal") { 3042 ylow = y; 3043 yup = yarr[2 * (n - 1) - 2 * i]; 3044 } else { 3045 yup = 0.5 * (f(x + delta) + f(x)); 3046 if (g) { 3047 ylow = 0.5 * (g(x + delta) + g(x)); 3048 } else { 3049 ylow = 0.0; 3050 } 3051 } 3052 sum += (yup - ylow) * delta; 3053 3054 // Draw the vertical lines 3055 j += 1; 3056 xarr[j] = x; 3057 yarr[j] = yarr[2 * (n - 1) - 2 * i]; 3058 3059 j += 1; 3060 } 3061 3062 return [xarr, yarr, sum]; 3063 }, 3064 3065 /** 3066 * Approximate the integral by Riemann sums. 3067 * Compute the area described by the riemann sum rectangles. 3068 * 3069 * If there is an element of type {@link Riemannsum}, then it is more efficient 3070 * to use the method JXG.Curve.Value() of this element instead. 3071 * 3072 * @param {Function_Array} f Function or array of two functions. 3073 * If f is a function the integral of this function is approximated by the Riemann sum. 3074 * If f is an array consisting of two functions the area between the two functions is approximated 3075 * by the Riemann sum. 3076 * @param {Number} n number of rectangles. 3077 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson' or 'trapezoidal'. 3078 * 3079 * @param {Number} start Left border of the approximation interval 3080 * @param {Number} end Right border of the approximation interval 3081 * @returns {Number} The sum of the areas of the rectangles. 3082 * @memberof JXG.Math.Numerics 3083 */ 3084 riemannsum: function (f, n, type, start, end) { 3085 JXG.deprecated("Numerics.riemannsum()", "Numerics.riemann()"); 3086 return this.riemann(f, n, type, start, end)[2]; 3087 }, 3088 3089 /** 3090 * Solve initial value problems numerically using Runge-Kutta-methods. 3091 * See {@link https://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm. 3092 * @param {object,String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing 3093 * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure 3094 * <pre> 3095 * { 3096 * s: <Number>, 3097 * A: <matrix>, 3098 * b: <Array>, 3099 * c: <Array> 3100 * } 3101 * </pre> 3102 * which corresponds to the Butcher tableau structure shown here: https://en.wikipedia.org/w/index.php?title=List_of_Runge%E2%80%93Kutta_methods&oldid=357796696 3103 * @param {Array} x0 Initial value vector. If the problem is of one-dimensional, the initial value also has to be given in an array. 3104 * @param {Array} I Interval on which to integrate. 3105 * @param {Number} N Number of evaluation points. 3106 * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode 3107 * 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 3108 * vector <tt>x</tt>, and has to return a vector of the same dimension as <tt>x</tt> has. 3109 * @returns {Array} An array of vectors describing the solution of the ode on the given interval I. 3110 * @example 3111 * // A very simple autonomous system dx(t)/dt = x(t); 3112 * function f(t, x) { 3113 * return x; 3114 * } 3115 * 3116 * // Solve it with initial value x(0) = 1 on the interval [0, 2] 3117 * // with 20 evaluation points. 3118 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 3119 * 3120 * // Prepare data for plotting the solution of the ode using a curve. 3121 * var dataX = []; 3122 * var dataY = []; 3123 * var h = 0.1; // (I[1] - I[0])/N = (2-0)/20 3124 * for(var i=0; i<data.length; i++) { 3125 * dataX[i] = i*h; 3126 * dataY[i] = data[i][0]; 3127 * } 3128 * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'}); 3129 * </pre><div class="jxgbox" id="JXGd2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div> 3130 * <script type="text/javascript"> 3131 * var board = JXG.JSXGraph.initBoard('JXGd2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false}); 3132 * function f(t, x) { 3133 * // we have to copy the value. 3134 * // return x; would just return the reference. 3135 * return [x[0]]; 3136 * } 3137 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 3138 * var dataX = []; 3139 * var dataY = []; 3140 * var h = 0.1; 3141 * for(var i=0; i<data.length; i++) { 3142 * dataX[i] = i*h; 3143 * dataY[i] = data[i][0]; 3144 * } 3145 * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'}); 3146 * </script><pre> 3147 * @memberof JXG.Math.Numerics 3148 */ 3149 rungeKutta: function (butcher, x0, I, N, f) { 3150 var e, 3151 i, 3152 j, 3153 k, 3154 l, 3155 s, 3156 x = [], 3157 y = [], 3158 h = (I[1] - I[0]) / N, 3159 t = I[0], 3160 dim = x0.length, 3161 result = [], 3162 r = 0; 3163 3164 if (Type.isString(butcher)) { 3165 butcher = predefinedButcher[butcher] || predefinedButcher.euler; 3166 } 3167 s = butcher.s; 3168 3169 // don't change x0, so copy it 3170 for (e = 0; e < dim; e++) { 3171 x[e] = x0[e]; 3172 } 3173 3174 for (i = 0; i < N; i++) { 3175 // Optimization doesn't work for ODEs plotted using time 3176 // if((i % quotient == 0) || (i == N-1)) { 3177 result[r] = []; 3178 for (e = 0; e < dim; e++) { 3179 result[r][e] = x[e]; 3180 } 3181 3182 r += 1; 3183 k = []; 3184 3185 for (j = 0; j < s; j++) { 3186 // init y = 0 3187 for (e = 0; e < dim; e++) { 3188 y[e] = 0.0; 3189 } 3190 3191 // Calculate linear combination of former k's and save it in y 3192 for (l = 0; l < j; l++) { 3193 for (e = 0; e < dim; e++) { 3194 y[e] += butcher.A[j][l] * h * k[l][e]; 3195 } 3196 } 3197 3198 // add x(t) to y 3199 for (e = 0; e < dim; e++) { 3200 y[e] += x[e]; 3201 } 3202 3203 // calculate new k and add it to the k matrix 3204 k.push(f(t + butcher.c[j] * h, y)); 3205 } 3206 3207 // init y = 0 3208 for (e = 0; e < dim; e++) { 3209 y[e] = 0.0; 3210 } 3211 3212 for (l = 0; l < s; l++) { 3213 for (e = 0; e < dim; e++) { 3214 y[e] += butcher.b[l] * k[l][e]; 3215 } 3216 } 3217 3218 for (e = 0; e < dim; e++) { 3219 x[e] = x[e] + h * y[e]; 3220 } 3221 3222 t += h; 3223 } 3224 3225 return result; 3226 }, 3227 3228 /** 3229 * Maximum number of iterations in {@link JXG.Math.Numerics.fzero} and 3230 * {@link JXG.Math.Numerics.chandrupatla} 3231 * @type Number 3232 * @default 80 3233 * @memberof JXG.Math.Numerics 3234 */ 3235 maxIterationsRoot: 80, 3236 3237 /** 3238 * Maximum number of iterations in {@link JXG.Math.Numerics.fminbr} 3239 * @type Number 3240 * @default 500 3241 * @memberof JXG.Math.Numerics 3242 */ 3243 maxIterationsMinimize: 500, 3244 3245 /** 3246 * Given a value x_0, this function tries to find a second value x_1 such that 3247 * the function f has opposite signs at x_0 and x_1. 3248 * The return values have to be tested if the method succeeded. 3249 * 3250 * @param {Function} f Function, whose root is to be found 3251 * @param {Number} x0 Start value 3252 * @param {Object} object Parent object in case f is method of it 3253 * @returns {Array} [x_0, f(x_0), x_1, f(x_1)] in case that x_0 <= x_1 3254 * or [x_1, f(x_1), x_0, f(x_0)] in case that x_1 < x_0. 3255 * 3256 * @see JXG.Math.Numerics.fzero 3257 * @see JXG.Math.Numerics.chandrupatla 3258 * 3259 * @memberof JXG.Math.Numerics 3260 */ 3261 findBracket: function (f, x0, object) { 3262 var a, aa, fa, blist, b, fb, u, fu, i, len; 3263 3264 if (Type.isArray(x0)) { 3265 return x0; 3266 } 3267 3268 a = x0; 3269 fa = f.call(object, a); 3270 // nfev += 1; 3271 3272 // Try to get b, by trying several values related to a 3273 aa = a === 0 ? 1 : a; 3274 blist = [ 3275 a - 0.1 * aa, 3276 a + 0.1 * aa, 3277 a - 1, 3278 a + 1, 3279 a - 0.5 * aa, 3280 a + 0.5 * aa, 3281 a - 0.6 * aa, 3282 a + 0.6 * aa, 3283 a - 1 * aa, 3284 a + 1 * aa, 3285 a - 2 * aa, 3286 a + 2 * aa, 3287 a - 5 * aa, 3288 a + 5 * aa, 3289 a - 10 * aa, 3290 a + 10 * aa, 3291 a - 50 * aa, 3292 a + 50 * aa, 3293 a - 100 * aa, 3294 a + 100 * aa 3295 ]; 3296 len = blist.length; 3297 3298 for (i = 0; i < len; i++) { 3299 b = blist[i]; 3300 fb = f.call(object, b); 3301 // nfev += 1; 3302 3303 if (fa * fb <= 0) { 3304 break; 3305 } 3306 } 3307 if (b < a) { 3308 u = a; 3309 a = b; 3310 b = u; 3311 3312 fu = fa; 3313 fa = fb; 3314 fb = fu; 3315 } 3316 return [a, fa, b, fb]; 3317 }, 3318 3319 /** 3320 * 3321 * Find zero of an univariate function f. 3322 * @param {function} f Function, whose root is to be found 3323 * @param {Array,Number} x0 Start value or start interval enclosing the root 3324 * @param {Object} object Parent object in case f is method of it 3325 * @returns {Number} the approximation of the root 3326 * Algorithm: 3327 * Brent's root finder from 3328 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 3329 * computations. M., Mir, 1980, p.180 of the Russian edition 3330 * https://www.netlib.org/c/brent.shar 3331 * 3332 * If x0 is an array containing lower and upper bound for the zero 3333 * algorithm 748 is applied. Otherwise, if x0 is a number, 3334 * the algorithm tries to bracket a zero of f starting from x0. 3335 * If this fails, we fall back to Newton's method. 3336 * 3337 * @see JXG.Math.Numerics.chandrupatla 3338 * @see JXG.Math.Numerics.root 3339 * @memberof JXG.Math.Numerics 3340 */ 3341 fzero: function (f, x0, object) { 3342 var a, 3343 b, 3344 c, 3345 d, 3346 e, 3347 fa, 3348 fb, 3349 fc, 3350 res, 3351 prev_step, 3352 t1, 3353 cb, 3354 t2, 3355 // Actual tolerance 3356 tol_act, 3357 // Interpolation step is calculated in the form p/q; division 3358 // operations is delayed until the last moment 3359 p, 3360 q, 3361 // Step at this iteration 3362 new_step, 3363 eps = Mat.eps, 3364 maxiter = this.maxIterationsRoot, 3365 niter = 0; 3366 // nfev = 0; 3367 3368 if (Type.isArray(x0)) { 3369 if (x0.length < 2) { 3370 throw new Error( 3371 "JXG.Math.Numerics.fzero: length of array x0 has to be at least two." 3372 ); 3373 } 3374 3375 a = x0[0]; 3376 fa = f.call(object, a); 3377 // nfev += 1; 3378 b = x0[1]; 3379 fb = f.call(object, b); 3380 // nfev += 1; 3381 } else { 3382 res = this.findBracket(f, x0, object); 3383 a = res[0]; 3384 fa = res[1]; 3385 b = res[2]; 3386 fb = res[3]; 3387 } 3388 3389 if (Math.abs(fa) <= eps) { 3390 return a; 3391 } 3392 if (Math.abs(fb) <= eps) { 3393 return b; 3394 } 3395 3396 if (fa * fb > 0) { 3397 // Bracketing not successful, fall back to Newton's method or to fminbr 3398 if (Type.isArray(x0)) { 3399 return this.fminbr(f, [a, b], object); 3400 } 3401 3402 return this.Newton(f, a, object); 3403 } 3404 3405 // OK, we have enclosed a zero of f. 3406 // Now we can start Brent's method 3407 c = a; 3408 fc = fa; 3409 3410 // Main iteration loop 3411 while (niter < maxiter) { 3412 // Distance from the last but one to the last approximation 3413 prev_step = b - a; 3414 3415 // Swap data for b to be the best approximation 3416 if (Math.abs(fc) < Math.abs(fb)) { 3417 a = b; 3418 b = c; 3419 c = a; 3420 3421 fa = fb; 3422 fb = fc; 3423 fc = fa; 3424 } 3425 tol_act = 2 * eps * Math.abs(b) + eps * 0.5; 3426 new_step = (c - b) * 0.5; 3427 3428 if (Math.abs(new_step) <= tol_act || Math.abs(fb) <= eps) { 3429 // Acceptable approx. is found 3430 return b; 3431 } 3432 3433 // Decide if the interpolation can be tried 3434 // If prev_step was large enough and was in true direction Interpolatiom may be tried 3435 if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) { 3436 cb = c - b; 3437 3438 // If we have only two distinct points linear interpolation can only be applied 3439 if (a === c) { 3440 t1 = fb / fa; 3441 p = cb * t1; 3442 q = 1.0 - t1; 3443 // Quadric inverse interpolation 3444 } else { 3445 q = fa / fc; 3446 t1 = fb / fc; 3447 t2 = fb / fa; 3448 3449 p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0)); 3450 q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0); 3451 } 3452 3453 // p was calculated with the opposite sign; make p positive 3454 if (p > 0) { 3455 q = -q; 3456 // and assign possible minus to q 3457 } else { 3458 p = -p; 3459 } 3460 3461 // If b+p/q falls in [b,c] and isn't too large it is accepted 3462 // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent 3463 if ( 3464 p < 0.75 * cb * q - Math.abs(tol_act * q) * 0.5 && 3465 p < Math.abs(prev_step * q * 0.5) 3466 ) { 3467 new_step = p / q; 3468 } 3469 } 3470 3471 // Adjust the step to be not less than tolerance 3472 if (Math.abs(new_step) < tol_act) { 3473 new_step = new_step > 0 ? tol_act : -tol_act; 3474 } 3475 3476 // Save the previous approx. 3477 a = b; 3478 fa = fb; 3479 b += new_step; 3480 fb = f.call(object, b); 3481 // Do step to a new approxim. 3482 // nfev += 1; 3483 3484 // Adjust c for it to have a sign opposite to that of b 3485 if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) { 3486 c = a; 3487 fc = fa; 3488 } 3489 niter++; 3490 } // End while 3491 3492 return b; 3493 }, 3494 3495 /** 3496 * Find zero of an univariate function f. 3497 * @param {function} f Function, whose root is to be found 3498 * @param {Array,Number} x0 Start value or start interval enclosing the root 3499 * @param {Object} object Parent object in case f is method of it 3500 * @returns {Number} the approximation of the root 3501 * Algorithm: 3502 * Chandrupatla's method, see 3503 * Tirupathi R. Chandrupatla, 3504 * "A new hybrid quadratic/bisection algorithm for finding the zero of a nonlinear function without using derivatives", 3505 * Advances in Engineering Software, Volume 28, Issue 3, April 1997, Pages 145-149. 3506 * 3507 * If x0 is an array containing lower and upper bound for the zero 3508 * algorithm 748 is applied. Otherwise, if x0 is a number, 3509 * the algorithm tries to bracket a zero of f starting from x0. 3510 * If this fails, we fall back to Newton's method. 3511 * 3512 * @see JXG.Math.Numerics.root 3513 * @see JXG.Math.Numerics.fzero 3514 * @memberof JXG.Math.Numerics 3515 */ 3516 chandrupatla: function (f, x0, object) { 3517 var a, 3518 fa, 3519 b, 3520 fb, 3521 res, 3522 niter = 0, 3523 maxiter = this.maxIterationsRoot, 3524 rand = 1 + Math.random() * 0.001, 3525 t = 0.5 * rand, 3526 eps = Mat.eps, // 1.e-10, 3527 dlt = 0.00001, 3528 x1, 3529 x2, 3530 x3, 3531 x, 3532 f1, 3533 f2, 3534 f3, 3535 y, 3536 xm, 3537 fm, 3538 tol, 3539 tl, 3540 xi, 3541 ph, 3542 fl, 3543 fh, 3544 AL, 3545 A, 3546 B, 3547 C, 3548 D; 3549 3550 if (Type.isArray(x0)) { 3551 if (x0.length < 2) { 3552 throw new Error( 3553 "JXG.Math.Numerics.fzero: length of array x0 has to be at least two." 3554 ); 3555 } 3556 3557 a = x0[0]; 3558 fa = f.call(object, a); 3559 // nfev += 1; 3560 b = x0[1]; 3561 fb = f.call(object, b); 3562 // nfev += 1; 3563 } else { 3564 res = this.findBracket(f, x0, object); 3565 a = res[0]; 3566 fa = res[1]; 3567 b = res[2]; 3568 fb = res[3]; 3569 } 3570 3571 if (fa * fb > 0) { 3572 // Bracketing not successful, fall back to Newton's method or to fminbr 3573 if (Type.isArray(x0)) { 3574 return this.fminbr(f, [a, b], object); 3575 } 3576 3577 return this.Newton(f, a, object); 3578 } 3579 3580 x1 = a; 3581 x2 = b; 3582 f1 = fa; 3583 f2 = fb; 3584 do { 3585 x = x1 + t * (x2 - x1); 3586 y = f.call(object, x); 3587 3588 // Arrange 2-1-3: 2-1 interval, 1 middle, 3 discarded point 3589 if (Math.sign(y) === Math.sign(f1)) { 3590 x3 = x1; 3591 x1 = x; 3592 f3 = f1; 3593 f1 = y; 3594 } else { 3595 x3 = x2; 3596 x2 = x1; 3597 f3 = f2; 3598 f2 = f1; 3599 } 3600 x1 = x; 3601 f1 = y; 3602 3603 xm = x1; 3604 fm = f1; 3605 if (Math.abs(f2) < Math.abs(f1)) { 3606 xm = x2; 3607 fm = f2; 3608 } 3609 tol = 2 * eps * Math.abs(xm) + 0.5 * dlt; 3610 tl = tol / Math.abs(x2 - x1); 3611 if (tl > 0.5 || fm === 0) { 3612 break; 3613 } 3614 // If inverse quadratic interpolation holds, use it 3615 xi = (x1 - x2) / (x3 - x2); 3616 ph = (f1 - f2) / (f3 - f2); 3617 fl = 1 - Math.sqrt(1 - xi); 3618 fh = Math.sqrt(xi); 3619 if (fl < ph && ph < fh) { 3620 AL = (x3 - x1) / (x2 - x1); 3621 A = f1 / (f2 - f1); 3622 B = f3 / (f2 - f3); 3623 C = f1 / (f3 - f1); 3624 D = f2 / (f3 - f2); 3625 t = A * B + C * D * AL; 3626 } else { 3627 t = 0.5 * rand; 3628 } 3629 // Adjust t away from the interval boundary 3630 if (t < tl) { 3631 t = tl; 3632 } 3633 if (t > 1 - tl) { 3634 t = 1 - tl; 3635 } 3636 niter++; 3637 } while (niter <= maxiter); 3638 // console.log(niter); 3639 3640 return xm; 3641 }, 3642 3643 /** 3644 * 3645 * Find minimum of an univariate function f. 3646 * <p> 3647 * Algorithm: 3648 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 3649 * computations. M., Mir, 1980, p.180 of the Russian edition 3650 * 3651 * @param {function} f Function, whose minimum is to be found 3652 * @param {Array} x0 Start interval enclosing the minimum 3653 * @param {Object} context Parent object in case f is method of it 3654 * @returns {Number} the approximation of the minimum value position 3655 * @memberof JXG.Math.Numerics 3656 **/ 3657 fminbr: function (f, x0, context) { 3658 var a, b, x, v, w, 3659 fx, fv, fw, 3660 range, middle_range, tol_act, new_step, 3661 p, q, t, ft, 3662 // Golden section ratio 3663 r = (3.0 - Math.sqrt(5.0)) * 0.5, 3664 tol = Mat.eps, 3665 sqrteps = Mat.eps, //Math.sqrt(Mat.eps), 3666 maxiter = this.maxIterationsMinimize, 3667 niter = 0; 3668 // nfev = 0; 3669 3670 if (!Type.isArray(x0) || x0.length < 2) { 3671 throw new Error( 3672 "JXG.Math.Numerics.fminbr: length of array x0 has to be at least two." 3673 ); 3674 } 3675 3676 a = x0[0]; 3677 b = x0[1]; 3678 v = a + r * (b - a); 3679 fv = f.call(context, v); 3680 3681 // First step - always gold section 3682 // nfev += 1; 3683 x = v; 3684 w = v; 3685 fx = fv; 3686 fw = fv; 3687 3688 while (niter < maxiter) { 3689 // Range over the interval in which we are looking for the minimum 3690 range = b - a; 3691 middle_range = (a + b) * 0.5; 3692 3693 // Actual tolerance 3694 tol_act = sqrteps * Math.abs(x) + tol / 3.0; 3695 3696 if (Math.abs(x - middle_range) + range * 0.5 <= 2.0 * tol_act) { 3697 // Acceptable approx. is found 3698 return x; 3699 } 3700 3701 // Obtain the golden section step 3702 new_step = r * (x < middle_range ? b - x : a - x); 3703 3704 // Decide if the interpolation can be tried. If x and w are distinct interpolatiom may be tried 3705 if (Math.abs(x - w) >= tol_act) { 3706 // Interpolation step is calculated as p/q; 3707 // division operation is delayed until last moment 3708 t = (x - w) * (fx - fv); 3709 q = (x - v) * (fx - fw); 3710 p = (x - v) * q - (x - w) * t; 3711 q = 2 * (q - t); 3712 3713 if (q > 0) { 3714 // q was calculated with the op- 3715 p = -p; // posite sign; make q positive 3716 } else { 3717 // and assign possible minus to 3718 q = -q; // p 3719 } 3720 if ( 3721 Math.abs(p) < Math.abs(new_step * q) && // If x+p/q falls in [a,b] 3722 p > q * (a - x + 2 * tol_act) && // not too close to a and 3723 p < q * (b - x - 2 * tol_act) 3724 ) { 3725 // b, and isn't too large 3726 new_step = p / q; // it is accepted 3727 } 3728 // If p/q is too large then the 3729 // golden section procedure can 3730 // reduce [a,b] range to more 3731 // extent 3732 } 3733 3734 // Adjust the step to be not less than tolerance 3735 if (Math.abs(new_step) < tol_act) { 3736 if (new_step > 0) { 3737 new_step = tol_act; 3738 } else { 3739 new_step = -tol_act; 3740 } 3741 } 3742 3743 // Obtain the next approximation to min 3744 // and reduce the enveloping range 3745 3746 // Tentative point for the min 3747 t = x + new_step; 3748 ft = f.call(context, t); 3749 // nfev += 1; 3750 3751 // t is a better approximation 3752 if (ft <= fx) { 3753 // Reduce the range so that t would fall within it 3754 if (t < x) { 3755 b = x; 3756 } else { 3757 a = x; 3758 } 3759 3760 // Assign the best approx to x 3761 v = w; 3762 w = x; 3763 x = t; 3764 3765 fv = fw; 3766 fw = fx; 3767 fx = ft; 3768 // x remains the better approx 3769 } else { 3770 // Reduce the range enclosing x 3771 if (t < x) { 3772 a = t; 3773 } else { 3774 b = t; 3775 } 3776 3777 if (ft <= fw || w === x) { 3778 v = w; 3779 w = t; 3780 fv = fw; 3781 fw = ft; 3782 } else if (ft <= fv || v === x || v === w) { 3783 v = t; 3784 fv = ft; 3785 } 3786 } 3787 niter += 1; 3788 } 3789 3790 return x; 3791 }, 3792 3793 /** 3794 * Implements the Ramer-Douglas-Peucker algorithm. 3795 * It discards points which are not necessary from the polygonal line defined by the point array 3796 * pts. The computation is done in screen coordinates. 3797 * Average runtime is O(nlog(n)), worst case runtime is O(n^2), where n is the number of points. 3798 * @param {Array} pts Array of {@link JXG.Coords} 3799 * @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>. 3800 * @returns {Array} An array containing points which represent an apparently identical curve as the points of pts do, but contains fewer points. 3801 * @memberof JXG.Math.Numerics 3802 */ 3803 RamerDouglasPeucker: function (pts, eps) { 3804 var allPts = [], 3805 newPts = [], 3806 i, k, len, 3807 endless = true, 3808 3809 /** 3810 * findSplit() is a subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}. 3811 * It searches for the point between index i and j which 3812 * has the largest distance from the line between the points i and j. 3813 * @param {Array} pts Array of {@link JXG.Coords} 3814 * @param {Number} i Index of a point in pts 3815 * @param {Number} j Index of a point in pts 3816 * @ignore 3817 * @private 3818 */ 3819 findSplit = function (pts, i, j) { 3820 var d, k, ci, cj, ck, 3821 x0, y0, x1, y1, 3822 den, lbda, 3823 huge = 10000, 3824 dist = 0, 3825 f = i; 3826 3827 if (j - i < 2) { 3828 return [-1.0, 0]; 3829 } 3830 3831 ci = pts[i].scrCoords; 3832 cj = pts[j].scrCoords; 3833 3834 if (isNaN(ci[1]) || isNaN(ci[2])) { 3835 return [NaN, i]; 3836 } 3837 if (isNaN(cj[1]) || isNaN(cj[2])) { 3838 return [NaN, j]; 3839 } 3840 3841 for (k = i + 1; k < j; k++) { 3842 ck = pts[k].scrCoords; 3843 if (isNaN(ck[1]) || isNaN(ck[2])) { 3844 return [NaN, k]; 3845 } 3846 3847 x0 = ck[1] - ci[1]; 3848 y0 = ck[2] - ci[2]; 3849 x1 = cj[1] - ci[1]; 3850 y1 = cj[2] - ci[2]; 3851 x0 = x0 === Infinity ? huge : x0; 3852 y0 = y0 === Infinity ? huge : y0; 3853 x1 = x1 === Infinity ? huge : x1; 3854 y1 = y1 === Infinity ? huge : y1; 3855 x0 = x0 === -Infinity ? -huge : x0; 3856 y0 = y0 === -Infinity ? -huge : y0; 3857 x1 = x1 === -Infinity ? -huge : x1; 3858 y1 = y1 === -Infinity ? -huge : y1; 3859 den = x1 * x1 + y1 * y1; 3860 3861 if (den >= Mat.eps) { 3862 lbda = (x0 * x1 + y0 * y1) / den; 3863 3864 if (lbda < 0.0) { 3865 lbda = 0.0; 3866 } else if (lbda > 1.0) { 3867 lbda = 1.0; 3868 } 3869 3870 x0 = x0 - lbda * x1; 3871 y0 = y0 - lbda * y1; 3872 d = x0 * x0 + y0 * y0; 3873 } else { 3874 lbda = 0.0; 3875 d = x0 * x0 + y0 * y0; 3876 } 3877 3878 if (d > dist) { 3879 dist = d; 3880 f = k; 3881 } 3882 } 3883 return [Math.sqrt(dist), f]; 3884 }, 3885 /** 3886 * RDP() is a private subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}. 3887 * It runs recursively through the point set and searches the 3888 * point which has the largest distance from the line between the first point and 3889 * the last point. If the distance from the line is greater than eps, this point is 3890 * included in our new point set otherwise it is discarded. 3891 * If it is taken, we recursively apply the subroutine to the point set before 3892 * and after the chosen point. 3893 * @param {Array} pts Array of {@link JXG.Coords} 3894 * @param {Number} i Index of an element of pts 3895 * @param {Number} j Index of an element of pts 3896 * @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>. 3897 * @param {Array} newPts Array of {@link JXG.Coords} 3898 * @ignore 3899 * @private 3900 */ 3901 RDP = function (pts, i, j, eps, newPts) { 3902 var result = findSplit(pts, i, j), 3903 k = result[1]; 3904 3905 if (isNaN(result[0])) { 3906 RDP(pts, i, k - 1, eps, newPts); 3907 newPts.push(pts[k]); 3908 do { 3909 ++k; 3910 } while (k <= j && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])); 3911 if (k <= j) { 3912 newPts.push(pts[k]); 3913 } 3914 RDP(pts, k + 1, j, eps, newPts); 3915 } else if (result[0] > eps) { 3916 RDP(pts, i, k, eps, newPts); 3917 RDP(pts, k, j, eps, newPts); 3918 } else { 3919 newPts.push(pts[j]); 3920 } 3921 }; 3922 3923 len = pts.length; 3924 3925 i = 0; 3926 while (endless) { 3927 // Search for the next point without NaN coordinates 3928 while (i < len && isNaN(pts[i].scrCoords[1] + pts[i].scrCoords[2])) { 3929 i += 1; 3930 } 3931 // Search for the next position of a NaN point 3932 k = i + 1; 3933 while (k < len && !isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])) { 3934 k += 1; 3935 } 3936 k--; 3937 3938 // Only proceed if something is left 3939 if (i < len && k > i) { 3940 newPts = []; 3941 newPts[0] = pts[i]; 3942 RDP(pts, i, k, eps, newPts); 3943 allPts = allPts.concat(newPts); 3944 } 3945 if (i >= len) { 3946 break; 3947 } 3948 // Push the NaN point 3949 if (k < len - 1) { 3950 allPts.push(pts[k + 1]); 3951 } 3952 i = k + 1; 3953 } 3954 3955 return allPts; 3956 }, 3957 3958 /** 3959 * Old name for the implementation of the Ramer-Douglas-Peucker algorithm. 3960 * @deprecated Use {@link JXG.Math.Numerics.RamerDouglasPeucker} 3961 * @memberof JXG.Math.Numerics 3962 */ 3963 RamerDouglasPeuker: function (pts, eps) { 3964 JXG.deprecated("Numerics.RamerDouglasPeuker()", "Numerics.RamerDouglasPeucker()"); 3965 return this.RamerDouglasPeucker(pts, eps); 3966 }, 3967 3968 /** 3969 * Implements the Visvalingam-Whyatt algorithm. 3970 * See M. Visvalingam, J. D. Whyatt: 3971 * "Line generalisation by repeated elimination of the smallest area", C.I.S.R.G Discussion paper 10, July 1992 3972 * 3973 * The algorithm discards points which are not necessary from the polygonal line defined by the point array 3974 * pts (consisting of type JXG.Coords). 3975 * @param {Array} pts Array of {@link JXG.Coords} 3976 * @param {Number} numPoints Number of remaining intermediate points. The first and the last point of the original points will 3977 * be taken in any case. 3978 * @returns {Array} An array containing points which approximates the curve defined by pts. 3979 * @memberof JXG.Math.Numerics 3980 * 3981 * @example 3982 * var i, p = []; 3983 * for (i = 0; i < 5; ++i) { 3984 * p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6])); 3985 * } 3986 * 3987 * // Plot a cardinal spline curve 3988 * var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5); 3989 * var cu1 = board.create('curve', splineArr, {strokeColor: 'green'}); 3990 * 3991 * var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'}); 3992 * c.updateDataArray = function() { 3993 * var i, len, points; 3994 * 3995 * // Reduce number of intermediate points with Visvakingam-Whyatt to 6 3996 * points = JXG.Math.Numerics.Visvalingam(cu1.points, 6); 3997 * // Plot the remaining points 3998 * len = points.length; 3999 * this.dataX = []; 4000 * this.dataY = []; 4001 * for (i = 0; i < len; i++) { 4002 * this.dataX.push(points[i].usrCoords[1]); 4003 * this.dataY.push(points[i].usrCoords[2]); 4004 * } 4005 * }; 4006 * board.update(); 4007 * 4008 * </pre><div id="JXGce0cc55c-b592-11e6-8270-104a7d3be7eb" class="jxgbox" style="width: 300px; height: 300px;"></div> 4009 * <script type="text/javascript"> 4010 * (function() { 4011 * var board = JXG.JSXGraph.initBoard('JXGce0cc55c-b592-11e6-8270-104a7d3be7eb', 4012 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 4013 * 4014 * var i, p = []; 4015 * for (i = 0; i < 5; ++i) { 4016 * p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6])); 4017 * } 4018 * 4019 * // Plot a cardinal spline curve 4020 * var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5); 4021 * var cu1 = board.create('curve', splineArr, {strokeColor: 'green'}); 4022 * 4023 * var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'}); 4024 * c.updateDataArray = function() { 4025 * var i, len, points; 4026 * 4027 * // Reduce number of intermediate points with Visvakingam-Whyatt to 6 4028 * points = JXG.Math.Numerics.Visvalingam(cu1.points, 6); 4029 * // Plot the remaining points 4030 * len = points.length; 4031 * this.dataX = []; 4032 * this.dataY = []; 4033 * for (i = 0; i < len; i++) { 4034 * this.dataX.push(points[i].usrCoords[1]); 4035 * this.dataY.push(points[i].usrCoords[2]); 4036 * } 4037 * }; 4038 * board.update(); 4039 * 4040 * })(); 4041 * 4042 * </script><pre> 4043 * 4044 */ 4045 Visvalingam: function (pts, numPoints) { 4046 var i, 4047 len, 4048 vol, 4049 lastVol, 4050 linkedList = [], 4051 heap = [], 4052 points = [], 4053 lft, 4054 rt, 4055 lft2, 4056 rt2, 4057 obj; 4058 4059 len = pts.length; 4060 4061 // At least one intermediate point is needed 4062 if (len <= 2) { 4063 return pts; 4064 } 4065 4066 // Fill the linked list 4067 // Add first point to the linked list 4068 linkedList[0] = { 4069 used: true, 4070 lft: null, 4071 node: null 4072 }; 4073 4074 // Add all intermediate points to the linked list, 4075 // whose triangle area is nonzero. 4076 lft = 0; 4077 for (i = 1; i < len - 1; i++) { 4078 vol = Math.abs( 4079 JXG.Math.Numerics.det([ 4080 pts[i - 1].usrCoords, 4081 pts[i].usrCoords, 4082 pts[i + 1].usrCoords 4083 ]) 4084 ); 4085 if (!isNaN(vol)) { 4086 obj = { 4087 v: vol, 4088 idx: i 4089 }; 4090 heap.push(obj); 4091 linkedList[i] = { 4092 used: true, 4093 lft: lft, 4094 node: obj 4095 }; 4096 linkedList[lft].rt = i; 4097 lft = i; 4098 } 4099 } 4100 4101 // Add last point to the linked list 4102 linkedList[len - 1] = { 4103 used: true, 4104 rt: null, 4105 lft: lft, 4106 node: null 4107 }; 4108 linkedList[lft].rt = len - 1; 4109 4110 // Remove points until only numPoints intermediate points remain 4111 lastVol = -Infinity; 4112 while (heap.length > numPoints) { 4113 // Sort the heap with the updated volume values 4114 heap.sort(function (a, b) { 4115 // descending sort 4116 return b.v - a.v; 4117 }); 4118 4119 // Remove the point with the smallest triangle 4120 i = heap.pop().idx; 4121 linkedList[i].used = false; 4122 lastVol = linkedList[i].node.v; 4123 4124 // Update the pointers of the linked list 4125 lft = linkedList[i].lft; 4126 rt = linkedList[i].rt; 4127 linkedList[lft].rt = rt; 4128 linkedList[rt].lft = lft; 4129 4130 // Update the values for the volumes in the linked list 4131 lft2 = linkedList[lft].lft; 4132 if (lft2 !== null) { 4133 vol = Math.abs( 4134 JXG.Math.Numerics.det([ 4135 pts[lft2].usrCoords, 4136 pts[lft].usrCoords, 4137 pts[rt].usrCoords 4138 ]) 4139 ); 4140 4141 linkedList[lft].node.v = vol >= lastVol ? vol : lastVol; 4142 } 4143 rt2 = linkedList[rt].rt; 4144 if (rt2 !== null) { 4145 vol = Math.abs( 4146 JXG.Math.Numerics.det([ 4147 pts[lft].usrCoords, 4148 pts[rt].usrCoords, 4149 pts[rt2].usrCoords 4150 ]) 4151 ); 4152 4153 linkedList[rt].node.v = vol >= lastVol ? vol : lastVol; 4154 } 4155 } 4156 4157 // Return an array with the remaining points 4158 i = 0; 4159 points = [pts[i]]; 4160 do { 4161 i = linkedList[i].rt; 4162 points.push(pts[i]); 4163 } while (linkedList[i].rt !== null); 4164 4165 return points; 4166 } 4167 }; 4168 4169 export default Mat.Numerics; 4170