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 * <p> 2964 * In case of type "simpson" and "trapezoidal", the horizontal line approximating the function value 2965 * is replaced by a parabola or a secant. IN case of "simpson", 2966 * the parabola is approximated visually by a polygonal chain of fixed step width. 2967 * 2968 * @param {Function,Array} f Function or array of two functions. 2969 * If f is a function the integral of this function is approximated by the Riemann sum. 2970 * If f is an array consisting of two functions the area between the two functions is filled 2971 * by the Riemann sum bars. 2972 * @param {Number} n number of rectangles. 2973 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson', or 'trapezoidal'. 2974 * "simpson" is Simpson's 1/3 rule. 2975 * @param {Number} start Left border of the approximation interval 2976 * @param {Number} end Right border of the approximation interval 2977 * @returns {Array} An array of two arrays containing the x and y coordinates for the rectangles showing the Riemann sum. This 2978 * 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 2979 * rectangles. 2980 * @memberof JXG.Math.Numerics 2981 */ 2982 riemann: function (gf, n, type, start, end) { 2983 var i, delta, 2984 k, a, b, c, f0, f1, f2, xx, h, 2985 steps = 30, // Fixed step width for Simpson's rule 2986 xarr = [], 2987 yarr = [], 2988 x = start, 2989 sum = 0, 2990 y, f, g; 2991 2992 if (Type.isArray(gf)) { 2993 g = gf[0]; 2994 f = gf[1]; 2995 } else { 2996 f = gf; 2997 } 2998 2999 n = Math.floor(n); 3000 3001 if (n <= 0) { 3002 return [xarr, yarr, sum]; 3003 } 3004 3005 delta = (end - start) / n; 3006 3007 // "Upper" horizontal line defined by function 3008 for (i = 0; i < n; i++) { 3009 if (type === "simpson") { 3010 sum += this._riemannValue(x, f, type, delta) * delta; 3011 3012 h = delta * 0.5; 3013 f0 = f(x); 3014 f1 = f(x + h); 3015 f2 = f(x + 2 * h); 3016 3017 a = (f2 + f0 - 2 * f1) / (h * h) * 0.5; 3018 b = (f2 - f0) / (2 * h); 3019 c = f1; 3020 for (k = 0; k < steps; k++) { 3021 xx = k * delta / steps - h; 3022 xarr.push(x + xx + h); 3023 yarr.push(a * xx * xx + b * xx + c); 3024 } 3025 x += delta; 3026 y = f2; 3027 } else { 3028 y = this._riemannValue(x, f, type, delta); 3029 xarr.push(x); 3030 yarr.push(y); 3031 3032 x += delta; 3033 if (type === "trapezoidal") { 3034 f2 = f(x); 3035 sum += (y + f2) * 0.5 * delta; 3036 y = f2; 3037 } else { 3038 sum += y * delta; 3039 } 3040 3041 xarr.push(x); 3042 yarr.push(y); 3043 } 3044 xarr.push(x); 3045 yarr.push(y); 3046 } 3047 3048 // "Lower" horizontal line 3049 // Go backwards 3050 for (i = 0; i < n; i++) { 3051 if (type === "simpson" && g) { 3052 sum -= this._riemannValue(x, g, type, -delta) * delta; 3053 3054 h = delta * 0.5; 3055 f0 = g(x); 3056 f1 = g(x - h); 3057 f2 = g(x - 2 * h); 3058 3059 a = (f2 + f0 - 2 * f1) / (h * h) * 0.5; 3060 b = (f2 - f0) / (2 * h); 3061 c = f1; 3062 for (k = 0; k < steps; k++) { 3063 xx = k * delta / steps - h; 3064 xarr.push(x - xx - h); 3065 yarr.push(a * xx * xx + b * xx + c); 3066 } 3067 x -= delta; 3068 y = f2; 3069 } else { 3070 if (g) { 3071 y = this._riemannValue(x, g, type, -delta); 3072 } else { 3073 y = 0.0; 3074 } 3075 xarr.push(x); 3076 yarr.push(y); 3077 3078 x -= delta; 3079 if (g) { 3080 if (type === "trapezoidal") { 3081 f2 = g(x); 3082 sum -= (y + f2) * 0.5 * delta; 3083 y = f2; 3084 } else { 3085 sum -= y * delta; 3086 } 3087 } 3088 } 3089 xarr.push(x); 3090 yarr.push(y); 3091 3092 // Draw the vertical lines 3093 xarr.push(x); 3094 yarr.push(f(x)); 3095 } 3096 3097 return [xarr, yarr, sum]; 3098 }, 3099 3100 /** 3101 * Approximate the integral by Riemann sums. 3102 * Compute the area described by the riemann sum rectangles. 3103 * 3104 * If there is an element of type {@link Riemannsum}, then it is more efficient 3105 * to use the method JXG.Curve.Value() of this element instead. 3106 * 3107 * @param {Function_Array} f Function or array of two functions. 3108 * If f is a function the integral of this function is approximated by the Riemann sum. 3109 * If f is an array consisting of two functions the area between the two functions is approximated 3110 * by the Riemann sum. 3111 * @param {Number} n number of rectangles. 3112 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson' or 'trapezoidal'. 3113 * 3114 * @param {Number} start Left border of the approximation interval 3115 * @param {Number} end Right border of the approximation interval 3116 * @returns {Number} The sum of the areas of the rectangles. 3117 * @memberof JXG.Math.Numerics 3118 */ 3119 riemannsum: function (f, n, type, start, end) { 3120 JXG.deprecated("Numerics.riemannsum()", "Numerics.riemann()[2]"); 3121 return this.riemann(f, n, type, start, end)[2]; 3122 }, 3123 3124 /** 3125 * Solve initial value problems numerically using <i>explicit</i> Runge-Kutta methods. 3126 * See {@link https://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm. 3127 * @param {object,String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing 3128 * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure 3129 * <pre> 3130 * { 3131 * s: <Number>, 3132 * A: <matrix>, 3133 * b: <Array>, 3134 * c: <Array> 3135 * } 3136 * </pre> 3137 * which corresponds to the Butcher tableau structure 3138 * shown here: https://en.wikipedia.org/w/index.php?title=List_of_Runge%E2%80%93Kutta_methods&oldid=357796696 . 3139 * <i>Default</i> is 'euler'. 3140 * @param {Array} x0 Initial value vector. Even if the problem is one-dimensional, the initial value has to be given in an array. 3141 * @param {Array} I Interval on which to integrate. 3142 * @param {Number} N Number of integration intervals, i.e. there are <i>N+1</i> evaluation points. 3143 * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode 3144 * 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 3145 * vector <tt>x</tt>, and has to return a vector of the same length as <tt>x</tt> has. 3146 * @returns {Array} An array of vectors describing the solution of the ode on the given interval I. 3147 * @example 3148 * // A very simple autonomous system dx(t)/dt = x(t); 3149 * var f = function(t, x) { 3150 * return [x[0]]; 3151 * } 3152 * 3153 * // Solve it with initial value x(0) = 1 on the interval [0, 2] 3154 * // with 20 evaluation points. 3155 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 3156 * 3157 * // Prepare data for plotting the solution of the ode using a curve. 3158 * var dataX = []; 3159 * var dataY = []; 3160 * var h = 0.1; // (I[1] - I[0])/N = (2-0)/20 3161 * var i; 3162 * for(i=0; i<data.length; i++) { 3163 * dataX[i] = i*h; 3164 * dataY[i] = data[i][0]; 3165 * } 3166 * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'}); 3167 * </pre><div class="jxgbox" id="JXGd2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div> 3168 * <script type="text/javascript"> 3169 * var board = JXG.JSXGraph.initBoard('JXGd2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false}); 3170 * var f = function(t, x) { 3171 * // we have to copy the value. 3172 * // return x; would just return the reference. 3173 * return [x[0]]; 3174 * } 3175 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 3176 * var dataX = []; 3177 * var dataY = []; 3178 * var h = 0.1; 3179 * for(var i=0; i<data.length; i++) { 3180 * dataX[i] = i*h; 3181 * dataY[i] = data[i][0]; 3182 * } 3183 * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'}); 3184 * </script><pre> 3185 * @memberof JXG.Math.Numerics 3186 */ 3187 rungeKutta: function (butcher, x0, I, N, f) { 3188 var e, 3189 i, j, k, l, s, 3190 x = [], 3191 y = [], 3192 h = (I[1] - I[0]) / N, 3193 t = I[0], 3194 dim = x0.length, 3195 result = [], 3196 r = 0; 3197 3198 if (Type.isString(butcher)) { 3199 butcher = predefinedButcher[butcher] || predefinedButcher.euler; 3200 } 3201 s = butcher.s; 3202 3203 // don't change x0, so copy it 3204 // for (e = 0; e < dim; e++) { 3205 // x[e] = x0[e]; 3206 // } 3207 x = x0.slice(); 3208 3209 for (i = 0; i <= N; i++) { 3210 // Optimization doesn't work for ODEs plotted using time 3211 // if((i % quotient == 0) || (i == N-1)) { 3212 // result[r] = []; 3213 // for (e = 0; e < dim; e++) { 3214 // result[r][e] = x[e]; 3215 // } 3216 result[r] = x.slice(); 3217 3218 r += 1; 3219 k = []; 3220 3221 for (j = 0; j < s; j++) { 3222 // init y = 0 3223 for (e = 0; e < dim; e++) { 3224 y[e] = 0.0; 3225 } 3226 3227 // Calculate linear combination of former k's and save it in y 3228 for (l = 0; l < j; l++) { 3229 for (e = 0; e < dim; e++) { 3230 y[e] += butcher.A[j][l] * h * k[l][e]; 3231 } 3232 } 3233 3234 // add x(t) to y 3235 for (e = 0; e < dim; e++) { 3236 y[e] += x[e]; 3237 } 3238 3239 // calculate new k and add it to the k matrix 3240 k.push(f(t + butcher.c[j] * h, y)); 3241 } 3242 3243 // init y = 0 3244 for (e = 0; e < dim; e++) { 3245 y[e] = 0.0; 3246 } 3247 3248 for (l = 0; l < s; l++) { 3249 for (e = 0; e < dim; e++) { 3250 y[e] += butcher.b[l] * k[l][e]; 3251 } 3252 } 3253 3254 for (e = 0; e < dim; e++) { 3255 x[e] = x[e] + h * y[e]; 3256 } 3257 3258 t += h; 3259 } 3260 3261 return result; 3262 }, 3263 3264 /** 3265 * Maximum number of iterations in {@link JXG.Math.Numerics.fzero} and 3266 * {@link JXG.Math.Numerics.chandrupatla} 3267 * @type Number 3268 * @default 80 3269 * @memberof JXG.Math.Numerics 3270 */ 3271 maxIterationsRoot: 80, 3272 3273 /** 3274 * Maximum number of iterations in {@link JXG.Math.Numerics.fminbr} 3275 * @type Number 3276 * @default 500 3277 * @memberof JXG.Math.Numerics 3278 */ 3279 maxIterationsMinimize: 500, 3280 3281 /** 3282 * Given a value x_0, this function tries to find a second value x_1 such that 3283 * the function f has opposite signs at x_0 and x_1. 3284 * The return values have to be tested if the method succeeded. 3285 * 3286 * @param {Function} f Function, whose root is to be found 3287 * @param {Number} x0 Start value 3288 * @param {Object} object Parent object in case f is method of it 3289 * @returns {Array} [x_0, f(x_0), x_1, f(x_1)] in case that x_0 <= x_1 3290 * or [x_1, f(x_1), x_0, f(x_0)] in case that x_1 < x_0. 3291 * 3292 * @see JXG.Math.Numerics.fzero 3293 * @see JXG.Math.Numerics.chandrupatla 3294 * 3295 * @memberof JXG.Math.Numerics 3296 */ 3297 findBracket: function (f, x0, object) { 3298 var a, aa, fa, blist, b, fb, u, fu, i, len; 3299 3300 if (Type.isArray(x0)) { 3301 return x0; 3302 } 3303 3304 a = x0; 3305 fa = f.call(object, a); 3306 // nfev += 1; 3307 3308 // Try to get b, by trying several values related to a 3309 aa = a === 0 ? 1 : a; 3310 blist = [ 3311 a - 0.1 * aa, 3312 a + 0.1 * aa, 3313 a - 1, 3314 a + 1, 3315 a - 0.5 * aa, 3316 a + 0.5 * aa, 3317 a - 0.6 * aa, 3318 a + 0.6 * aa, 3319 a - 1 * aa, 3320 a + 1 * aa, 3321 a - 2 * aa, 3322 a + 2 * aa, 3323 a - 5 * aa, 3324 a + 5 * aa, 3325 a - 10 * aa, 3326 a + 10 * aa, 3327 a - 50 * aa, 3328 a + 50 * aa, 3329 a - 100 * aa, 3330 a + 100 * aa 3331 ]; 3332 len = blist.length; 3333 3334 for (i = 0; i < len; i++) { 3335 b = blist[i]; 3336 fb = f.call(object, b); 3337 // nfev += 1; 3338 3339 if (fa * fb <= 0) { 3340 break; 3341 } 3342 } 3343 if (b < a) { 3344 u = a; 3345 a = b; 3346 b = u; 3347 3348 fu = fa; 3349 fa = fb; 3350 fb = fu; 3351 } 3352 return [a, fa, b, fb]; 3353 }, 3354 3355 /** 3356 * 3357 * Find zero of an univariate function f. 3358 * @param {function} f Function, whose root is to be found 3359 * @param {Array,Number} x0 Start value or start interval enclosing the root 3360 * @param {Object} object Parent object in case f is method of it 3361 * @returns {Number} the approximation of the root 3362 * Algorithm: 3363 * Brent's root finder from 3364 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 3365 * computations. M., Mir, 1980, p.180 of the Russian edition 3366 * https://www.netlib.org/c/brent.shar 3367 * 3368 * If x0 is an array containing lower and upper bound for the zero 3369 * algorithm 748 is applied. Otherwise, if x0 is a number, 3370 * the algorithm tries to bracket a zero of f starting from x0. 3371 * If this fails, we fall back to Newton's method. 3372 * 3373 * @see JXG.Math.Numerics.chandrupatla 3374 * @see JXG.Math.Numerics.root 3375 * @memberof JXG.Math.Numerics 3376 */ 3377 fzero: function (f, x0, object) { 3378 var a, 3379 b, 3380 c, 3381 d, 3382 e, 3383 fa, 3384 fb, 3385 fc, 3386 res, 3387 prev_step, 3388 t1, 3389 cb, 3390 t2, 3391 // Actual tolerance 3392 tol_act, 3393 // Interpolation step is calculated in the form p/q; division 3394 // operations is delayed until the last moment 3395 p, 3396 q, 3397 // Step at this iteration 3398 new_step, 3399 eps = Mat.eps, 3400 maxiter = this.maxIterationsRoot, 3401 niter = 0; 3402 // nfev = 0; 3403 3404 if (Type.isArray(x0)) { 3405 if (x0.length < 2) { 3406 throw new Error( 3407 "JXG.Math.Numerics.fzero: length of array x0 has to be at least two." 3408 ); 3409 } 3410 3411 a = x0[0]; 3412 fa = f.call(object, a); 3413 // nfev += 1; 3414 b = x0[1]; 3415 fb = f.call(object, b); 3416 // nfev += 1; 3417 } else { 3418 res = this.findBracket(f, x0, object); 3419 a = res[0]; 3420 fa = res[1]; 3421 b = res[2]; 3422 fb = res[3]; 3423 } 3424 3425 if (Math.abs(fa) <= eps) { 3426 return a; 3427 } 3428 if (Math.abs(fb) <= eps) { 3429 return b; 3430 } 3431 3432 if (fa * fb > 0) { 3433 // Bracketing not successful, fall back to Newton's method or to fminbr 3434 if (Type.isArray(x0)) { 3435 return this.fminbr(f, [a, b], object); 3436 } 3437 3438 return this.Newton(f, a, object); 3439 } 3440 3441 // OK, we have enclosed a zero of f. 3442 // Now we can start Brent's method 3443 c = a; 3444 fc = fa; 3445 3446 // Main iteration loop 3447 while (niter < maxiter) { 3448 // Distance from the last but one to the last approximation 3449 prev_step = b - a; 3450 3451 // Swap data for b to be the best approximation 3452 if (Math.abs(fc) < Math.abs(fb)) { 3453 a = b; 3454 b = c; 3455 c = a; 3456 3457 fa = fb; 3458 fb = fc; 3459 fc = fa; 3460 } 3461 tol_act = 2 * eps * Math.abs(b) + eps * 0.5; 3462 new_step = (c - b) * 0.5; 3463 3464 if (Math.abs(new_step) <= tol_act || Math.abs(fb) <= eps) { 3465 // Acceptable approx. is found 3466 return b; 3467 } 3468 3469 // Decide if the interpolation can be tried 3470 // If prev_step was large enough and was in true direction Interpolatiom may be tried 3471 if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) { 3472 cb = c - b; 3473 3474 // If we have only two distinct points linear interpolation can only be applied 3475 if (a === c) { 3476 t1 = fb / fa; 3477 p = cb * t1; 3478 q = 1.0 - t1; 3479 // Quadric inverse interpolation 3480 } else { 3481 q = fa / fc; 3482 t1 = fb / fc; 3483 t2 = fb / fa; 3484 3485 p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0)); 3486 q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0); 3487 } 3488 3489 // p was calculated with the opposite sign; make p positive 3490 if (p > 0) { 3491 q = -q; 3492 // and assign possible minus to q 3493 } else { 3494 p = -p; 3495 } 3496 3497 // If b+p/q falls in [b,c] and isn't too large it is accepted 3498 // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent 3499 if ( 3500 p < 0.75 * cb * q - Math.abs(tol_act * q) * 0.5 && 3501 p < Math.abs(prev_step * q * 0.5) 3502 ) { 3503 new_step = p / q; 3504 } 3505 } 3506 3507 // Adjust the step to be not less than tolerance 3508 if (Math.abs(new_step) < tol_act) { 3509 new_step = new_step > 0 ? tol_act : -tol_act; 3510 } 3511 3512 // Save the previous approx. 3513 a = b; 3514 fa = fb; 3515 b += new_step; 3516 fb = f.call(object, b); 3517 // Do step to a new approxim. 3518 // nfev += 1; 3519 3520 // Adjust c for it to have a sign opposite to that of b 3521 if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) { 3522 c = a; 3523 fc = fa; 3524 } 3525 niter++; 3526 } // End while 3527 3528 return b; 3529 }, 3530 3531 /** 3532 * Find zero of an univariate function f. 3533 * @param {function} f Function, whose root is to be found 3534 * @param {Array,Number} x0 Start value or start interval enclosing the root 3535 * @param {Object} object Parent object in case f is method of it 3536 * @returns {Number} the approximation of the root 3537 * Algorithm: 3538 * Chandrupatla's method, see 3539 * Tirupathi R. Chandrupatla, 3540 * "A new hybrid quadratic/bisection algorithm for finding the zero of a nonlinear function without using derivatives", 3541 * Advances in Engineering Software, Volume 28, Issue 3, April 1997, Pages 145-149. 3542 * 3543 * If x0 is an array containing lower and upper bound for the zero 3544 * algorithm 748 is applied. Otherwise, if x0 is a number, 3545 * the algorithm tries to bracket a zero of f starting from x0. 3546 * If this fails, we fall back to Newton's method. 3547 * 3548 * @see JXG.Math.Numerics.root 3549 * @see JXG.Math.Numerics.fzero 3550 * @memberof JXG.Math.Numerics 3551 */ 3552 chandrupatla: function (f, x0, object) { 3553 var a, 3554 fa, 3555 b, 3556 fb, 3557 res, 3558 niter = 0, 3559 maxiter = this.maxIterationsRoot, 3560 rand = 1 + Math.random() * 0.001, 3561 t = 0.5 * rand, 3562 eps = Mat.eps, // 1.e-10, 3563 dlt = 0.00001, 3564 x1, 3565 x2, 3566 x3, 3567 x, 3568 f1, 3569 f2, 3570 f3, 3571 y, 3572 xm, 3573 fm, 3574 tol, 3575 tl, 3576 xi, 3577 ph, 3578 fl, 3579 fh, 3580 AL, 3581 A, 3582 B, 3583 C, 3584 D; 3585 3586 if (Type.isArray(x0)) { 3587 if (x0.length < 2) { 3588 throw new Error( 3589 "JXG.Math.Numerics.fzero: length of array x0 has to be at least two." 3590 ); 3591 } 3592 3593 a = x0[0]; 3594 fa = f.call(object, a); 3595 // nfev += 1; 3596 b = x0[1]; 3597 fb = f.call(object, b); 3598 // nfev += 1; 3599 } else { 3600 res = this.findBracket(f, x0, object); 3601 a = res[0]; 3602 fa = res[1]; 3603 b = res[2]; 3604 fb = res[3]; 3605 } 3606 3607 if (fa * fb > 0) { 3608 // Bracketing not successful, fall back to Newton's method or to fminbr 3609 if (Type.isArray(x0)) { 3610 return this.fminbr(f, [a, b], object); 3611 } 3612 3613 return this.Newton(f, a, object); 3614 } 3615 3616 x1 = a; 3617 x2 = b; 3618 f1 = fa; 3619 f2 = fb; 3620 do { 3621 x = x1 + t * (x2 - x1); 3622 y = f.call(object, x); 3623 3624 // Arrange 2-1-3: 2-1 interval, 1 middle, 3 discarded point 3625 if (Math.sign(y) === Math.sign(f1)) { 3626 x3 = x1; 3627 x1 = x; 3628 f3 = f1; 3629 f1 = y; 3630 } else { 3631 x3 = x2; 3632 x2 = x1; 3633 f3 = f2; 3634 f2 = f1; 3635 } 3636 x1 = x; 3637 f1 = y; 3638 3639 xm = x1; 3640 fm = f1; 3641 if (Math.abs(f2) < Math.abs(f1)) { 3642 xm = x2; 3643 fm = f2; 3644 } 3645 tol = 2 * eps * Math.abs(xm) + 0.5 * dlt; 3646 tl = tol / Math.abs(x2 - x1); 3647 if (tl > 0.5 || fm === 0) { 3648 break; 3649 } 3650 // If inverse quadratic interpolation holds, use it 3651 xi = (x1 - x2) / (x3 - x2); 3652 ph = (f1 - f2) / (f3 - f2); 3653 fl = 1 - Math.sqrt(1 - xi); 3654 fh = Math.sqrt(xi); 3655 if (fl < ph && ph < fh) { 3656 AL = (x3 - x1) / (x2 - x1); 3657 A = f1 / (f2 - f1); 3658 B = f3 / (f2 - f3); 3659 C = f1 / (f3 - f1); 3660 D = f2 / (f3 - f2); 3661 t = A * B + C * D * AL; 3662 } else { 3663 t = 0.5 * rand; 3664 } 3665 // Adjust t away from the interval boundary 3666 if (t < tl) { 3667 t = tl; 3668 } 3669 if (t > 1 - tl) { 3670 t = 1 - tl; 3671 } 3672 niter++; 3673 } while (niter <= maxiter); 3674 // console.log(niter); 3675 3676 return xm; 3677 }, 3678 3679 /** 3680 * 3681 * Find minimum of an univariate function f. 3682 * <p> 3683 * Algorithm: 3684 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 3685 * computations. M., Mir, 1980, p.180 of the Russian edition 3686 * 3687 * @param {function} f Function, whose minimum is to be found 3688 * @param {Array} x0 Start interval enclosing the minimum 3689 * @param {Object} context Parent object in case f is method of it 3690 * @returns {Number} the approximation of the minimum value position 3691 * @memberof JXG.Math.Numerics 3692 **/ 3693 fminbr: function (f, x0, context) { 3694 var a, b, x, v, w, 3695 fx, fv, fw, 3696 range, middle_range, tol_act, new_step, 3697 p, q, t, ft, 3698 // Golden section ratio 3699 r = (3.0 - Math.sqrt(5.0)) * 0.5, 3700 tol = Mat.eps, 3701 sqrteps = Mat.eps, //Math.sqrt(Mat.eps), 3702 maxiter = this.maxIterationsMinimize, 3703 niter = 0; 3704 // nfev = 0; 3705 3706 if (!Type.isArray(x0) || x0.length < 2) { 3707 throw new Error( 3708 "JXG.Math.Numerics.fminbr: length of array x0 has to be at least two." 3709 ); 3710 } 3711 3712 a = x0[0]; 3713 b = x0[1]; 3714 v = a + r * (b - a); 3715 fv = f.call(context, v); 3716 3717 // First step - always gold section 3718 // nfev += 1; 3719 x = v; 3720 w = v; 3721 fx = fv; 3722 fw = fv; 3723 3724 while (niter < maxiter) { 3725 // Range over the interval in which we are looking for the minimum 3726 range = b - a; 3727 middle_range = (a + b) * 0.5; 3728 3729 // Actual tolerance 3730 tol_act = sqrteps * Math.abs(x) + tol / 3.0; 3731 3732 if (Math.abs(x - middle_range) + range * 0.5 <= 2.0 * tol_act) { 3733 // Acceptable approx. is found 3734 return x; 3735 } 3736 3737 // Obtain the golden section step 3738 new_step = r * (x < middle_range ? b - x : a - x); 3739 3740 // Decide if the interpolation can be tried. If x and w are distinct interpolatiom may be tried 3741 if (Math.abs(x - w) >= tol_act) { 3742 // Interpolation step is calculated as p/q; 3743 // division operation is delayed until last moment 3744 t = (x - w) * (fx - fv); 3745 q = (x - v) * (fx - fw); 3746 p = (x - v) * q - (x - w) * t; 3747 q = 2 * (q - t); 3748 3749 if (q > 0) { 3750 // q was calculated with the op- 3751 p = -p; // posite sign; make q positive 3752 } else { 3753 // and assign possible minus to 3754 q = -q; // p 3755 } 3756 if ( 3757 Math.abs(p) < Math.abs(new_step * q) && // If x+p/q falls in [a,b] 3758 p > q * (a - x + 2 * tol_act) && // not too close to a and 3759 p < q * (b - x - 2 * tol_act) 3760 ) { 3761 // b, and isn't too large 3762 new_step = p / q; // it is accepted 3763 } 3764 // If p/q is too large then the 3765 // golden section procedure can 3766 // reduce [a,b] range to more 3767 // extent 3768 } 3769 3770 // Adjust the step to be not less than tolerance 3771 if (Math.abs(new_step) < tol_act) { 3772 if (new_step > 0) { 3773 new_step = tol_act; 3774 } else { 3775 new_step = -tol_act; 3776 } 3777 } 3778 3779 // Obtain the next approximation to min 3780 // and reduce the enveloping range 3781 3782 // Tentative point for the min 3783 t = x + new_step; 3784 ft = f.call(context, t); 3785 // nfev += 1; 3786 3787 // t is a better approximation 3788 if (ft <= fx) { 3789 // Reduce the range so that t would fall within it 3790 if (t < x) { 3791 b = x; 3792 } else { 3793 a = x; 3794 } 3795 3796 // Assign the best approx to x 3797 v = w; 3798 w = x; 3799 x = t; 3800 3801 fv = fw; 3802 fw = fx; 3803 fx = ft; 3804 // x remains the better approx 3805 } else { 3806 // Reduce the range enclosing x 3807 if (t < x) { 3808 a = t; 3809 } else { 3810 b = t; 3811 } 3812 3813 if (ft <= fw || w === x) { 3814 v = w; 3815 w = t; 3816 fv = fw; 3817 fw = ft; 3818 } else if (ft <= fv || v === x || v === w) { 3819 v = t; 3820 fv = ft; 3821 } 3822 } 3823 niter += 1; 3824 } 3825 3826 return x; 3827 }, 3828 3829 /** 3830 * Implements the Ramer-Douglas-Peucker algorithm. 3831 * It discards points which are not necessary from the polygonal line defined by the point array 3832 * pts. The computation is done in screen coordinates. 3833 * Average runtime is O(nlog(n)), worst case runtime is O(n^2), where n is the number of points. 3834 * @param {Array} pts Array of {@link JXG.Coords} 3835 * @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>. 3836 * @returns {Array} An array containing points which represent an apparently identical curve as the points of pts do, but contains fewer points. 3837 * @memberof JXG.Math.Numerics 3838 */ 3839 RamerDouglasPeucker: function (pts, eps) { 3840 var allPts = [], 3841 newPts = [], 3842 i, k, len, 3843 endless = true, 3844 3845 /** 3846 * findSplit() is a subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}. 3847 * It searches for the point between index i and j which 3848 * has the largest distance from the line between the points i and j. 3849 * @param {Array} pts Array of {@link JXG.Coords} 3850 * @param {Number} i Index of a point in pts 3851 * @param {Number} j Index of a point in pts 3852 * @ignore 3853 * @private 3854 */ 3855 findSplit = function (pts, i, j) { 3856 var d, k, ci, cj, ck, 3857 x0, y0, x1, y1, 3858 den, lbda, 3859 huge = 10000, 3860 dist = 0, 3861 f = i; 3862 3863 if (j - i < 2) { 3864 return [-1.0, 0]; 3865 } 3866 3867 ci = pts[i].scrCoords; 3868 cj = pts[j].scrCoords; 3869 3870 if (isNaN(ci[1]) || isNaN(ci[2])) { 3871 return [NaN, i]; 3872 } 3873 if (isNaN(cj[1]) || isNaN(cj[2])) { 3874 return [NaN, j]; 3875 } 3876 3877 for (k = i + 1; k < j; k++) { 3878 ck = pts[k].scrCoords; 3879 if (isNaN(ck[1]) || isNaN(ck[2])) { 3880 return [NaN, k]; 3881 } 3882 3883 x0 = ck[1] - ci[1]; 3884 y0 = ck[2] - ci[2]; 3885 x1 = cj[1] - ci[1]; 3886 y1 = cj[2] - ci[2]; 3887 x0 = x0 === Infinity ? huge : x0; 3888 y0 = y0 === Infinity ? huge : y0; 3889 x1 = x1 === Infinity ? huge : x1; 3890 y1 = y1 === Infinity ? huge : y1; 3891 x0 = x0 === -Infinity ? -huge : x0; 3892 y0 = y0 === -Infinity ? -huge : y0; 3893 x1 = x1 === -Infinity ? -huge : x1; 3894 y1 = y1 === -Infinity ? -huge : y1; 3895 den = x1 * x1 + y1 * y1; 3896 3897 if (den >= Mat.eps) { 3898 lbda = (x0 * x1 + y0 * y1) / den; 3899 3900 if (lbda < 0.0) { 3901 lbda = 0.0; 3902 } else if (lbda > 1.0) { 3903 lbda = 1.0; 3904 } 3905 3906 x0 = x0 - lbda * x1; 3907 y0 = y0 - lbda * y1; 3908 d = x0 * x0 + y0 * y0; 3909 } else { 3910 lbda = 0.0; 3911 d = x0 * x0 + y0 * y0; 3912 } 3913 3914 if (d > dist) { 3915 dist = d; 3916 f = k; 3917 } 3918 } 3919 return [Math.sqrt(dist), f]; 3920 }, 3921 /** 3922 * RDP() is a private subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}. 3923 * It runs recursively through the point set and searches the 3924 * point which has the largest distance from the line between the first point and 3925 * the last point. If the distance from the line is greater than eps, this point is 3926 * included in our new point set otherwise it is discarded. 3927 * If it is taken, we recursively apply the subroutine to the point set before 3928 * and after the chosen point. 3929 * @param {Array} pts Array of {@link JXG.Coords} 3930 * @param {Number} i Index of an element of pts 3931 * @param {Number} j Index of an element of pts 3932 * @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>. 3933 * @param {Array} newPts Array of {@link JXG.Coords} 3934 * @ignore 3935 * @private 3936 */ 3937 RDP = function (pts, i, j, eps, newPts) { 3938 var result = findSplit(pts, i, j), 3939 k = result[1]; 3940 3941 if (isNaN(result[0])) { 3942 RDP(pts, i, k - 1, eps, newPts); 3943 newPts.push(pts[k]); 3944 do { 3945 ++k; 3946 } while (k <= j && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])); 3947 if (k <= j) { 3948 newPts.push(pts[k]); 3949 } 3950 RDP(pts, k + 1, j, eps, newPts); 3951 } else if (result[0] > eps) { 3952 RDP(pts, i, k, eps, newPts); 3953 RDP(pts, k, j, eps, newPts); 3954 } else { 3955 newPts.push(pts[j]); 3956 } 3957 }; 3958 3959 len = pts.length; 3960 3961 i = 0; 3962 while (endless) { 3963 // Search for the next point without NaN coordinates 3964 while (i < len && isNaN(pts[i].scrCoords[1] + pts[i].scrCoords[2])) { 3965 i += 1; 3966 } 3967 // Search for the next position of a NaN point 3968 k = i + 1; 3969 while (k < len && !isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])) { 3970 k += 1; 3971 } 3972 k--; 3973 3974 // Only proceed if something is left 3975 if (i < len && k > i) { 3976 newPts = []; 3977 newPts[0] = pts[i]; 3978 RDP(pts, i, k, eps, newPts); 3979 allPts = allPts.concat(newPts); 3980 } 3981 if (i >= len) { 3982 break; 3983 } 3984 // Push the NaN point 3985 if (k < len - 1) { 3986 allPts.push(pts[k + 1]); 3987 } 3988 i = k + 1; 3989 } 3990 3991 return allPts; 3992 }, 3993 3994 /** 3995 * Old name for the implementation of the Ramer-Douglas-Peucker algorithm. 3996 * @deprecated Use {@link JXG.Math.Numerics.RamerDouglasPeucker} 3997 * @memberof JXG.Math.Numerics 3998 */ 3999 RamerDouglasPeuker: function (pts, eps) { 4000 JXG.deprecated("Numerics.RamerDouglasPeuker()", "Numerics.RamerDouglasPeucker()"); 4001 return this.RamerDouglasPeucker(pts, eps); 4002 }, 4003 4004 /** 4005 * Implements the Visvalingam-Whyatt algorithm. 4006 * See M. Visvalingam, J. D. Whyatt: 4007 * "Line generalisation by repeated elimination of the smallest area", C.I.S.R.G Discussion paper 10, July 1992 4008 * 4009 * The algorithm discards points which are not necessary from the polygonal line defined by the point array 4010 * pts (consisting of type JXG.Coords). 4011 * @param {Array} pts Array of {@link JXG.Coords} 4012 * @param {Number} numPoints Number of remaining intermediate points. The first and the last point of the original points will 4013 * be taken in any case. 4014 * @returns {Array} An array containing points which approximates the curve defined by pts. 4015 * @memberof JXG.Math.Numerics 4016 * 4017 * @example 4018 * var i, p = []; 4019 * for (i = 0; i < 5; ++i) { 4020 * p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6])); 4021 * } 4022 * 4023 * // Plot a cardinal spline curve 4024 * var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5); 4025 * var cu1 = board.create('curve', splineArr, {strokeColor: 'green'}); 4026 * 4027 * var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'}); 4028 * c.updateDataArray = function() { 4029 * var i, len, points; 4030 * 4031 * // Reduce number of intermediate points with Visvakingam-Whyatt to 6 4032 * points = JXG.Math.Numerics.Visvalingam(cu1.points, 6); 4033 * // Plot the remaining points 4034 * len = points.length; 4035 * this.dataX = []; 4036 * this.dataY = []; 4037 * for (i = 0; i < len; i++) { 4038 * this.dataX.push(points[i].usrCoords[1]); 4039 * this.dataY.push(points[i].usrCoords[2]); 4040 * } 4041 * }; 4042 * board.update(); 4043 * 4044 * </pre><div id="JXGce0cc55c-b592-11e6-8270-104a7d3be7eb" class="jxgbox" style="width: 300px; height: 300px;"></div> 4045 * <script type="text/javascript"> 4046 * (function() { 4047 * var board = JXG.JSXGraph.initBoard('JXGce0cc55c-b592-11e6-8270-104a7d3be7eb', 4048 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 4049 * 4050 * var i, p = []; 4051 * for (i = 0; i < 5; ++i) { 4052 * p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6])); 4053 * } 4054 * 4055 * // Plot a cardinal spline curve 4056 * var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5); 4057 * var cu1 = board.create('curve', splineArr, {strokeColor: 'green'}); 4058 * 4059 * var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'}); 4060 * c.updateDataArray = function() { 4061 * var i, len, points; 4062 * 4063 * // Reduce number of intermediate points with Visvakingam-Whyatt to 6 4064 * points = JXG.Math.Numerics.Visvalingam(cu1.points, 6); 4065 * // Plot the remaining points 4066 * len = points.length; 4067 * this.dataX = []; 4068 * this.dataY = []; 4069 * for (i = 0; i < len; i++) { 4070 * this.dataX.push(points[i].usrCoords[1]); 4071 * this.dataY.push(points[i].usrCoords[2]); 4072 * } 4073 * }; 4074 * board.update(); 4075 * 4076 * })(); 4077 * 4078 * </script><pre> 4079 * 4080 */ 4081 Visvalingam: function (pts, numPoints) { 4082 var i, 4083 len, 4084 vol, 4085 lastVol, 4086 linkedList = [], 4087 heap = [], 4088 points = [], 4089 lft, 4090 rt, 4091 lft2, 4092 rt2, 4093 obj; 4094 4095 len = pts.length; 4096 4097 // At least one intermediate point is needed 4098 if (len <= 2) { 4099 return pts; 4100 } 4101 4102 // Fill the linked list 4103 // Add first point to the linked list 4104 linkedList[0] = { 4105 used: true, 4106 lft: null, 4107 node: null 4108 }; 4109 4110 // Add all intermediate points to the linked list, 4111 // whose triangle area is nonzero. 4112 lft = 0; 4113 for (i = 1; i < len - 1; i++) { 4114 vol = Math.abs( 4115 JXG.Math.Numerics.det([ 4116 pts[i - 1].usrCoords, 4117 pts[i].usrCoords, 4118 pts[i + 1].usrCoords 4119 ]) 4120 ); 4121 if (!isNaN(vol)) { 4122 obj = { 4123 v: vol, 4124 idx: i 4125 }; 4126 heap.push(obj); 4127 linkedList[i] = { 4128 used: true, 4129 lft: lft, 4130 node: obj 4131 }; 4132 linkedList[lft].rt = i; 4133 lft = i; 4134 } 4135 } 4136 4137 // Add last point to the linked list 4138 linkedList[len - 1] = { 4139 used: true, 4140 rt: null, 4141 lft: lft, 4142 node: null 4143 }; 4144 linkedList[lft].rt = len - 1; 4145 4146 // Remove points until only numPoints intermediate points remain 4147 lastVol = -Infinity; 4148 while (heap.length > numPoints) { 4149 // Sort the heap with the updated volume values 4150 heap.sort(function (a, b) { 4151 // descending sort 4152 return b.v - a.v; 4153 }); 4154 4155 // Remove the point with the smallest triangle 4156 i = heap.pop().idx; 4157 linkedList[i].used = false; 4158 lastVol = linkedList[i].node.v; 4159 4160 // Update the pointers of the linked list 4161 lft = linkedList[i].lft; 4162 rt = linkedList[i].rt; 4163 linkedList[lft].rt = rt; 4164 linkedList[rt].lft = lft; 4165 4166 // Update the values for the volumes in the linked list 4167 lft2 = linkedList[lft].lft; 4168 if (lft2 !== null) { 4169 vol = Math.abs( 4170 JXG.Math.Numerics.det([ 4171 pts[lft2].usrCoords, 4172 pts[lft].usrCoords, 4173 pts[rt].usrCoords 4174 ]) 4175 ); 4176 4177 linkedList[lft].node.v = vol >= lastVol ? vol : lastVol; 4178 } 4179 rt2 = linkedList[rt].rt; 4180 if (rt2 !== null) { 4181 vol = Math.abs( 4182 JXG.Math.Numerics.det([ 4183 pts[lft].usrCoords, 4184 pts[rt].usrCoords, 4185 pts[rt2].usrCoords 4186 ]) 4187 ); 4188 4189 linkedList[rt].node.v = vol >= lastVol ? vol : lastVol; 4190 } 4191 } 4192 4193 // Return an array with the remaining points 4194 i = 0; 4195 points = [pts[i]]; 4196 do { 4197 i = linkedList[i].rt; 4198 points.push(pts[i]); 4199 } while (linkedList[i].rt !== null); 4200 4201 return points; 4202 } 4203 }; 4204 4205 export default Mat.Numerics; 4206