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