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 * @param {Object} context optional object that is treated as "this" in the function body. This is useful if 1558 * the function is a method of an object and contains a reference to its parent object via "this". 1559 * @returns {Number} A root of the function f. 1560 * 1561 * @see JXG.Math.Numerics.chandrupatla 1562 * @see JXG.Math.Numerics.fzero 1563 * @memberof JXG.Math.Numerics 1564 */ 1565 root: function (f, x, context) { 1566 //return this.fzero(f, x, context); 1567 return this.chandrupatla(f, x, context); 1568 }, 1569 1570 /** 1571 * Compute an intersection of the curves c1 and c2 1572 * with a generalized Newton method. 1573 * We want to find values t1, t2 such that 1574 * c1(t1) = c2(t2), i.e. 1575 * <br> 1576 * (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2)) = (0,0). 1577 * <p> 1578 * We set 1579 * (e,f) := (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2)) 1580 * <p> 1581 * The Jacobian J is defined by 1582 * <pre> 1583 * J = (a, b) 1584 * (c, d) 1585 * </pre> 1586 * where 1587 * <ul> 1588 * <li> a = c1_x'(t1) 1589 * <li> b = -c2_x'(t2) 1590 * <li> c = c1_y'(t1) 1591 * <li> d = -c2_y'(t2) 1592 * </ul> 1593 * The inverse J^(-1) of J is equal to 1594 * <pre> 1595 * (d, -b)/ (ad-bc) 1596 * (-c, a) / (ad-bc) 1597 * </pre> 1598 * 1599 * Then, (t1new, t2new) := (t1,t2) - J^(-1)*(e,f). 1600 * <p> 1601 * If the function meetCurveCurve has the properties 1602 * t1memo and t2memo then these are taken as start values 1603 * for the Newton algorithm. 1604 * After stopping of the Newton algorithm the values of t1 and t2 are stored in 1605 * t1memo and t2memo. 1606 * 1607 * @param {JXG.Curve} c1 Curve, Line or Circle 1608 * @param {JXG.Curve} c2 Curve, Line or Circle 1609 * @param {Number} t1ini start value for t1 1610 * @param {Number} t2ini start value for t2 1611 * @returns {JXG.Coords} intersection point 1612 * @memberof JXG.Math.Numerics 1613 */ 1614 generalizedNewton: function (c1, c2, t1ini, t2ini) { 1615 var t1, t2, 1616 a, b, c, d, e, f, 1617 disc, 1618 F, 1619 D00, D01, D10, D11, 1620 count = 0; 1621 1622 if (this.generalizedNewton.t1memo) { 1623 t1 = this.generalizedNewton.t1memo; 1624 t2 = this.generalizedNewton.t2memo; 1625 } else { 1626 t1 = t1ini; 1627 t2 = t2ini; 1628 } 1629 1630 e = c1.X(t1) - c2.X(t2); 1631 f = c1.Y(t1) - c2.Y(t2); 1632 F = e * e + f * f; 1633 1634 D00 = this.D(c1.X, c1); 1635 D01 = this.D(c2.X, c2); 1636 D10 = this.D(c1.Y, c1); 1637 D11 = this.D(c2.Y, c2); 1638 1639 while (F > Mat.eps && count < 10) { 1640 a = D00(t1); 1641 b = -D01(t2); 1642 c = D10(t1); 1643 d = -D11(t2); 1644 disc = a * d - b * c; 1645 t1 -= (d * e - b * f) / disc; 1646 t2 -= (a * f - c * e) / disc; 1647 e = c1.X(t1) - c2.X(t2); 1648 f = c1.Y(t1) - c2.Y(t2); 1649 F = e * e + f * f; 1650 count += 1; 1651 } 1652 1653 this.generalizedNewton.t1memo = t1; 1654 this.generalizedNewton.t2memo = t2; 1655 1656 if (Math.abs(t1) < Math.abs(t2)) { 1657 return [c1.X(t1), c1.Y(t1)]; 1658 } 1659 1660 return [c2.X(t2), c2.Y(t2)]; 1661 }, 1662 1663 /** 1664 * Returns the Lagrange polynomials for curves with equidistant nodes, see 1665 * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation, 1666 * SIAM Review, Vol 46, No 3, (2004) 501-517. 1667 * The graph of the parametric curve [x(t),y(t)] runs through the given points. 1668 * @param {Array} p Array of JXG.Points 1669 * @returns {Array} An array consisting of two functions x(t), y(t) which define a parametric curve 1670 * f(t) = (x(t), y(t)), a number x1 (which equals 0) and a function x2 defining the curve's domain. 1671 * That means the curve is defined between x1 and x2(). x2 returns the (length of array p minus one). 1672 * @memberof JXG.Math.Numerics 1673 * 1674 * @example 1675 * var p = []; 1676 * 1677 * p[0] = board.create('point', [0, -2], {size:2, name: 'C(a)'}); 1678 * p[1] = board.create('point', [-1.5, 5], {size:2, name: ''}); 1679 * p[2] = board.create('point', [1, 4], {size:2, name: ''}); 1680 * p[3] = board.create('point', [3, 3], {size:2, name: 'C(b)'}); 1681 * 1682 * // Curve 1683 * var fg = JXG.Math.Numerics.Neville(p); 1684 * var graph = board.create('curve', fg, {strokeWidth:3, strokeOpacity:0.5}); 1685 * 1686 * </pre><div id="JXG88a8b3a8-6561-44f5-a678-76bca13fd484" class="jxgbox" style="width: 300px; height: 300px;"></div> 1687 * <script type="text/javascript"> 1688 * (function() { 1689 * var board = JXG.JSXGraph.initBoard('JXG88a8b3a8-6561-44f5-a678-76bca13fd484', 1690 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 1691 * var p = []; 1692 * 1693 * p[0] = board.create('point', [0, -2], {size:2, name: 'C(a)'}); 1694 * p[1] = board.create('point', [-1.5, 5], {size:2, name: ''}); 1695 * p[2] = board.create('point', [1, 4], {size:2, name: ''}); 1696 * p[3] = board.create('point', [3, 3], {size:2, name: 'C(b)'}); 1697 * 1698 * // Curve 1699 * var fg = JXG.Math.Numerics.Neville(p); 1700 * var graph = board.create('curve', fg, {strokeWidth:3, strokeOpacity:0.5}); 1701 * 1702 * })(); 1703 * 1704 * </script><pre> 1705 * 1706 */ 1707 Neville: function (p) { 1708 var w = [], 1709 /** @ignore */ 1710 makeFct = function (fun) { 1711 return function (t, suspendedUpdate) { 1712 var i, 1713 d, 1714 s, 1715 bin = Mat.binomial, 1716 len = p.length, 1717 len1 = len - 1, 1718 num = 0.0, 1719 denom = 0.0; 1720 1721 if (!suspendedUpdate) { 1722 s = 1; 1723 for (i = 0; i < len; i++) { 1724 w[i] = bin(len1, i) * s; 1725 s *= -1; 1726 } 1727 } 1728 1729 d = t; 1730 1731 for (i = 0; i < len; i++) { 1732 if (d === 0) { 1733 return p[i][fun](); 1734 } 1735 s = w[i] / d; 1736 d -= 1; 1737 num += p[i][fun]() * s; 1738 denom += s; 1739 } 1740 return num / denom; 1741 }; 1742 }, 1743 xfct = makeFct("X"), 1744 yfct = makeFct("Y"); 1745 1746 return [ 1747 xfct, 1748 yfct, 1749 0, 1750 function () { 1751 return p.length - 1; 1752 } 1753 ]; 1754 }, 1755 1756 /** 1757 * Calculates second derivatives at the given knots. 1758 * @param {Array} x x values of knots 1759 * @param {Array} y y values of knots 1760 * @returns {Array} Second derivatives of the interpolated function at the knots. 1761 * @see #splineEval 1762 * @memberof JXG.Math.Numerics 1763 */ 1764 splineDef: function (x, y) { 1765 var pair, 1766 i, 1767 l, 1768 n = Math.min(x.length, y.length), 1769 diag = [], 1770 z = [], 1771 data = [], 1772 dx = [], 1773 delta = [], 1774 F = []; 1775 1776 if (n === 2) { 1777 return [0, 0]; 1778 } 1779 1780 for (i = 0; i < n; i++) { 1781 pair = { X: x[i], Y: y[i] }; 1782 data.push(pair); 1783 } 1784 data.sort(function (a, b) { 1785 return a.X - b.X; 1786 }); 1787 for (i = 0; i < n; i++) { 1788 x[i] = data[i].X; 1789 y[i] = data[i].Y; 1790 } 1791 1792 for (i = 0; i < n - 1; i++) { 1793 dx.push(x[i + 1] - x[i]); 1794 } 1795 for (i = 0; i < n - 2; i++) { 1796 delta.push( 1797 (6 * (y[i + 2] - y[i + 1])) / dx[i + 1] - (6 * (y[i + 1] - y[i])) / dx[i] 1798 ); 1799 } 1800 1801 // ForwardSolve 1802 diag.push(2 * (dx[0] + dx[1])); 1803 z.push(delta[0]); 1804 1805 for (i = 0; i < n - 3; i++) { 1806 l = dx[i + 1] / diag[i]; 1807 diag.push(2 * (dx[i + 1] + dx[i + 2]) - l * dx[i + 1]); 1808 z.push(delta[i + 1] - l * z[i]); 1809 } 1810 1811 // BackwardSolve 1812 F[n - 3] = z[n - 3] / diag[n - 3]; 1813 for (i = n - 4; i >= 0; i--) { 1814 F[i] = (z[i] - dx[i + 1] * F[i + 1]) / diag[i]; 1815 } 1816 1817 // Generate f''-Vector 1818 for (i = n - 3; i >= 0; i--) { 1819 F[i + 1] = F[i]; 1820 } 1821 1822 // natural cubic spline 1823 F[0] = 0; 1824 F[n - 1] = 0; 1825 1826 return F; 1827 }, 1828 1829 /** 1830 * Evaluate points on spline. 1831 * @param {Number|Array} x0 A single float value or an array of values to evaluate 1832 * @param {Array} x x values of knots 1833 * @param {Array} y y values of knots 1834 * @param {Array} F Second derivatives at knots, calculated by {@link JXG.Math.Numerics.splineDef} 1835 * @see #splineDef 1836 * @returns {Number|Array} A single value or an array, depending on what is given as x0. 1837 * @memberof JXG.Math.Numerics 1838 */ 1839 splineEval: function (x0, x, y, F) { 1840 var i, 1841 j, 1842 a, 1843 b, 1844 c, 1845 d, 1846 x_, 1847 n = Math.min(x.length, y.length), 1848 l = 1, 1849 asArray = false, 1850 y0 = []; 1851 1852 // number of points to be evaluated 1853 if (Type.isArray(x0)) { 1854 l = x0.length; 1855 asArray = true; 1856 } else { 1857 x0 = [x0]; 1858 } 1859 1860 for (i = 0; i < l; i++) { 1861 // is x0 in defining interval? 1862 if (x0[i] < x[0] || x[i] > x[n - 1]) { 1863 return NaN; 1864 } 1865 1866 // determine part of spline in which x0 lies 1867 for (j = 1; j < n; j++) { 1868 if (x0[i] <= x[j]) { 1869 break; 1870 } 1871 } 1872 1873 j -= 1; 1874 1875 // we're now in the j-th partial interval, i.e. x[j] < x0[i] <= x[j+1]; 1876 // determine the coefficients of the polynomial in this interval 1877 a = y[j]; 1878 b = 1879 (y[j + 1] - y[j]) / (x[j + 1] - x[j]) - 1880 ((x[j + 1] - x[j]) / 6) * (F[j + 1] + 2 * F[j]); 1881 c = F[j] / 2; 1882 d = (F[j + 1] - F[j]) / (6 * (x[j + 1] - x[j])); 1883 // evaluate x0[i] 1884 x_ = x0[i] - x[j]; 1885 //y0.push(a + b*x_ + c*x_*x_ + d*x_*x_*x_); 1886 y0.push(a + (b + (c + d * x_) * x_) * x_); 1887 } 1888 1889 if (asArray) { 1890 return y0; 1891 } 1892 1893 return y0[0]; 1894 }, 1895 1896 /** 1897 * Generate a string containing the function term of a polynomial. 1898 * @param {Array} coeffs Coefficients of the polynomial. The position i belongs to x^i. 1899 * @param {Number} deg Degree of the polynomial 1900 * @param {String} varname Name of the variable (usually 'x') 1901 * @param {Number} prec Precision 1902 * @returns {String} A string containg the function term of the polynomial. 1903 * @memberof JXG.Math.Numerics 1904 */ 1905 generatePolynomialTerm: function (coeffs, deg, varname, prec) { 1906 var i, 1907 t = []; 1908 1909 for (i = deg; i >= 0; i--) { 1910 Type.concat(t, ["(", coeffs[i].toPrecision(prec), ")"]); 1911 1912 if (i > 1) { 1913 Type.concat(t, ["*", varname, "<sup>", i, "<", "/sup> + "]); 1914 } else if (i === 1) { 1915 Type.concat(t, ["*", varname, " + "]); 1916 } 1917 } 1918 1919 return t.join(""); 1920 }, 1921 1922 /** 1923 * Computes the polynomial through a given set of coordinates in Lagrange form. 1924 * Returns the Lagrange polynomials, see 1925 * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation, 1926 * SIAM Review, Vol 46, No 3, (2004) 501-517. 1927 * <p> 1928 * It possesses the method getTerm() which returns the string containing the function term of the polynomial. 1929 * @param {Array} p Array of JXG.Points 1930 * @returns {function} A function of one parameter which returns the value of the polynomial, whose graph runs through the given points. 1931 * @memberof JXG.Math.Numerics 1932 * 1933 * @example 1934 * var p = []; 1935 * p[0] = board.create('point', [-1,2], {size:4}); 1936 * p[1] = board.create('point', [0,3], {size:4}); 1937 * p[2] = board.create('point', [1,1], {size:4}); 1938 * p[3] = board.create('point', [3,-1], {size:4}); 1939 * var f = JXG.Math.Numerics.lagrangePolynomial(p); 1940 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 1941 * 1942 * </pre><div id="JXGc058aa6b-74d4-41e1-af94-df06169a2d89" class="jxgbox" style="width: 300px; height: 300px;"></div> 1943 * <script type="text/javascript"> 1944 * (function() { 1945 * var board = JXG.JSXGraph.initBoard('JXGc058aa6b-74d4-41e1-af94-df06169a2d89', 1946 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 1947 * var p = []; 1948 * p[0] = board.create('point', [-1,2], {size:4}); 1949 * p[1] = board.create('point', [0,3], {size:4}); 1950 * p[2] = board.create('point', [1,1], {size:4}); 1951 * p[3] = board.create('point', [3,-1], {size:4}); 1952 * var f = JXG.Math.Numerics.lagrangePolynomial(p); 1953 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 1954 * 1955 * })(); 1956 * 1957 * </script><pre> 1958 * 1959 * @example 1960 * var points = []; 1961 * points[0] = board.create('point', [-1,2], {size:4}); 1962 * points[1] = board.create('point', [0, 0], {size:4}); 1963 * points[2] = board.create('point', [2, 1], {size:4}); 1964 * 1965 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 1966 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 1967 * var txt = board.create('text', [-3, -4, () => f.getTerm(2, 't', ' * ')], {fontSize: 16}); 1968 * var txt2 = board.create('text', [-3, -6, () => f.getCoefficients()], {fontSize: 12}); 1969 * 1970 * </pre><div id="JXG73fdaf12-e257-4374-b488-ae063e4eecbb" class="jxgbox" style="width: 300px; height: 300px;"></div> 1971 * <script type="text/javascript"> 1972 * (function() { 1973 * var board = JXG.JSXGraph.initBoard('JXG73fdaf12-e257-4374-b488-ae063e4eecbb', 1974 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 1975 * var points = []; 1976 * points[0] = board.create('point', [-1,2], {size:4}); 1977 * points[1] = board.create('point', [0, 0], {size:4}); 1978 * points[2] = board.create('point', [2, 1], {size:4}); 1979 * 1980 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 1981 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 1982 * var txt = board.create('text', [-3, -4, () => f.getTerm(2, 't', ' * ')], {fontSize: 16}); 1983 * var txt2 = board.create('text', [-3, -6, () => f.getCoefficients()], {fontSize: 12}); 1984 * 1985 * })(); 1986 * 1987 * </script><pre> 1988 * 1989 */ 1990 lagrangePolynomial: function (p) { 1991 var w = [], 1992 that = this, 1993 /** @ignore */ 1994 fct = function (x, suspendedUpdate) { 1995 var i, // j, 1996 k, 1997 xi, 1998 s, //M, 1999 len = p.length, 2000 num = 0, 2001 denom = 0; 2002 2003 if (!suspendedUpdate) { 2004 for (i = 0; i < len; i++) { 2005 w[i] = 1.0; 2006 xi = p[i].X(); 2007 2008 for (k = 0; k < len; k++) { 2009 if (k !== i) { 2010 w[i] *= xi - p[k].X(); 2011 } 2012 } 2013 2014 w[i] = 1 / w[i]; 2015 } 2016 2017 // M = []; 2018 // for (k = 0; k < len; k++) { 2019 // M.push([1]); 2020 // } 2021 } 2022 2023 for (i = 0; i < len; i++) { 2024 xi = p[i].X(); 2025 2026 if (x === xi) { 2027 return p[i].Y(); 2028 } 2029 2030 s = w[i] / (x - xi); 2031 denom += s; 2032 num += s * p[i].Y(); 2033 } 2034 2035 return num / denom; 2036 }; 2037 2038 /** 2039 * Get the term of the Lagrange polynomial as string. 2040 * Calls {@link JXG.Math.Numerics#lagrangePolynomialTerm}. 2041 * 2042 * @name JXG.Math.Numerics.lagrangePolynomial#getTerm 2043 * @param {Number} digits Number of digits of the coefficients 2044 * @param {String} param Variable name 2045 * @param {String} dot Dot symbol 2046 * @returns {String} containing the term of Lagrange polynomial as string. 2047 * @see JXG.Math.Numerics#lagrangePolynomialTerm 2048 * @example 2049 * var points = []; 2050 * points[0] = board.create('point', [-1,2], {size:4}); 2051 * points[1] = board.create('point', [0, 0], {size:4}); 2052 * points[2] = board.create('point', [2, 1], {size:4}); 2053 * 2054 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2055 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2056 * var txt = board.create('text', [-3, -4, () => f.getTerm(2, 't', ' * ')], {fontSize: 16}); 2057 * 2058 * </pre><div id="JXG73fdaf12-e257-4374-b488-ae063e4eeccf" class="jxgbox" style="width: 300px; height: 300px;"></div> 2059 * <script type="text/javascript"> 2060 * (function() { 2061 * var board = JXG.JSXGraph.initBoard('JXG73fdaf12-e257-4374-b488-ae063e4eeccf', 2062 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 2063 * var points = []; 2064 * points[0] = board.create('point', [-1,2], {size:4}); 2065 * points[1] = board.create('point', [0, 0], {size:4}); 2066 * points[2] = board.create('point', [2, 1], {size:4}); 2067 * 2068 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2069 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2070 * var txt = board.create('text', [-3, -4, () => f.getTerm(2, 't', ' * ')], {fontSize: 16}); 2071 * 2072 * })(); 2073 * 2074 * </script><pre> 2075 * 2076 */ 2077 fct.getTerm = function (digits, param, dot) { 2078 return that.lagrangePolynomialTerm(p, digits, param, dot)(); 2079 }; 2080 2081 /** 2082 * Get the coefficients of the Lagrange polynomial as array. The leading 2083 * coefficient is at position 0. 2084 * Calls {@link JXG.Math.Numerics#lagrangePolynomialCoefficients}. 2085 * 2086 * @name JXG.Math.Numerics.lagrangePolynomial#getCoefficients 2087 * @returns {Array} containing the coefficients of the Lagrange polynomial. 2088 * @see JXG.Math.Numerics.lagrangePolynomial#getTerm 2089 * @see JXG.Math.Numerics#lagrangePolynomialTerm 2090 * @see JXG.Math.Numerics#lagrangePolynomialCoefficients 2091 * @example 2092 * var points = []; 2093 * points[0] = board.create('point', [-1,2], {size:4}); 2094 * points[1] = board.create('point', [0, 0], {size:4}); 2095 * points[2] = board.create('point', [2, 1], {size:4}); 2096 * 2097 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2098 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2099 * var txt = board.create('text', [1, -4, () => f.getCoefficients()], {fontSize: 10}); 2100 * 2101 * </pre><div id="JXG52a883a5-2e0c-4caf-8f84-8650c173c365" class="jxgbox" style="width: 300px; height: 300px;"></div> 2102 * <script type="text/javascript"> 2103 * (function() { 2104 * var board = JXG.JSXGraph.initBoard('JXG52a883a5-2e0c-4caf-8f84-8650c173c365', 2105 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 2106 * var points = []; 2107 * points[0] = board.create('point', [-1,2], {size:4}); 2108 * points[1] = board.create('point', [0, 0], {size:4}); 2109 * points[2] = board.create('point', [2, 1], {size:4}); 2110 * 2111 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2112 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2113 * var txt = board.create('text', [1, -4, () => f.getCoefficients()], {fontSize: 10}); 2114 * 2115 * })(); 2116 * 2117 * </script><pre> 2118 * 2119 */ 2120 fct.getCoefficients = function () { 2121 return that.lagrangePolynomialCoefficients(p)(); 2122 }; 2123 2124 return fct; 2125 }, 2126 2127 /** 2128 * Determine the Lagrange polynomial through an array of points and 2129 * return the term of the polynomial as string. 2130 * 2131 * @param {Array} points Array of JXG.Points 2132 * @param {Number} digits Number of decimal digits of the coefficients 2133 * @param {String} param Name of the parameter. Default: 'x'. 2134 * @param {String} dot Multiplication symbol. Default: ' * '. 2135 * @returns {Function} returning the Lagrange polynomial term through 2136 * the supplied points as string 2137 * @memberof JXG.Math.Numerics 2138 * 2139 * @example 2140 * var points = []; 2141 * points[0] = board.create('point', [-1,2], {size:4}); 2142 * points[1] = board.create('point', [0, 0], {size:4}); 2143 * points[2] = board.create('point', [2, 1], {size:4}); 2144 * 2145 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2146 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2147 * 2148 * var f_txt = JXG.Math.Numerics.lagrangePolynomialTerm(points, 2, 't', ' * '); 2149 * var txt = board.create('text', [-3, -4, f_txt], {fontSize: 16}); 2150 * 2151 * </pre><div id="JXGd45e9e96-7526-486d-aa43-e1178d5f2baa" class="jxgbox" style="width: 300px; height: 300px;"></div> 2152 * <script type="text/javascript"> 2153 * (function() { 2154 * var board = JXG.JSXGraph.initBoard('JXGd45e9e96-7526-486d-aa43-e1178d5f2baa', 2155 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 2156 * var points = []; 2157 * points[0] = board.create('point', [-1,2], {size:4}); 2158 * points[1] = board.create('point', [0, 0], {size:4}); 2159 * points[2] = board.create('point', [2, 1], {size:4}); 2160 * 2161 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2162 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2163 * 2164 * var f_txt = JXG.Math.Numerics.lagrangePolynomialTerm(points, 2, 't', ' * '); 2165 * var txt = board.create('text', [-3, -4, f_txt], {fontSize: 16}); 2166 * 2167 * })(); 2168 * 2169 * </script><pre> 2170 * 2171 */ 2172 lagrangePolynomialTerm: function (points, digits, param, dot) { 2173 var that = this; 2174 2175 return function () { 2176 var len = points.length, 2177 coeffs = [], 2178 isLeading = true, 2179 n, t, j, c; 2180 2181 param = param || "x"; 2182 if (dot === undefined) { 2183 dot = " * "; 2184 } 2185 2186 n = len - 1; // (Max) degree of the polynomial 2187 coeffs = that.lagrangePolynomialCoefficients(points)(); 2188 2189 t = ""; 2190 for (j = 0; j < coeffs.length; j++) { 2191 c = coeffs[j]; 2192 if (Math.abs(c) < Mat.eps) { 2193 continue; 2194 } 2195 if (JXG.exists(digits)) { 2196 c = Env._round10(c, -digits); 2197 } 2198 if (isLeading) { 2199 t += c > 0 ? c : "-" + -c; 2200 isLeading = false; 2201 } else { 2202 t += c > 0 ? " + " + c : " - " + -c; 2203 } 2204 2205 if (n - j > 1) { 2206 t += dot + param + "^" + (n - j); 2207 } else if (n - j === 1) { 2208 t += dot + param; 2209 } 2210 } 2211 return t; // board.jc.manipulate('f = map(x) -> ' + t + ';'); 2212 }; 2213 }, 2214 2215 /** 2216 * Determine the Lagrange polynomial through an array of points and 2217 * return the coefficients of the polynomial as array. 2218 * The leading coefficient is at position 0. 2219 * 2220 * @param {Array} points Array of JXG.Points 2221 * @returns {Function} returning the coefficients of the Lagrange polynomial through 2222 * the supplied points. 2223 * @memberof JXG.Math.Numerics 2224 * 2225 * @example 2226 * var points = []; 2227 * points[0] = board.create('point', [-1,2], {size:4}); 2228 * points[1] = board.create('point', [0, 0], {size:4}); 2229 * points[2] = board.create('point', [2, 1], {size:4}); 2230 * 2231 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2232 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2233 * 2234 * var f_arr = JXG.Math.Numerics.lagrangePolynomialCoefficients(points); 2235 * var txt = board.create('text', [1, -4, f_arr], {fontSize: 10}); 2236 * 2237 * </pre><div id="JXG1778f0d1-a420-473f-99e8-1755ef4be97e" class="jxgbox" style="width: 300px; height: 300px;"></div> 2238 * <script type="text/javascript"> 2239 * (function() { 2240 * var board = JXG.JSXGraph.initBoard('JXG1778f0d1-a420-473f-99e8-1755ef4be97e', 2241 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 2242 * var points = []; 2243 * points[0] = board.create('point', [-1,2], {size:4}); 2244 * points[1] = board.create('point', [0, 0], {size:4}); 2245 * points[2] = board.create('point', [2, 1], {size:4}); 2246 * 2247 * var f = JXG.Math.Numerics.lagrangePolynomial(points); 2248 * var graph = board.create('functiongraph', [f,-10, 10], {strokeWidth:3}); 2249 * 2250 * var f_arr = JXG.Math.Numerics.lagrangePolynomialCoefficients(points); 2251 * var txt = board.create('text', [1, -4, f_arr], {fontSize: 10}); 2252 * 2253 * })(); 2254 * 2255 * </script><pre> 2256 * 2257 */ 2258 lagrangePolynomialCoefficients: function (points) { 2259 return function () { 2260 var len = points.length, 2261 zeroes = [], 2262 coeffs = [], 2263 coeffs_sum = [], 2264 i, j, c, p; 2265 2266 // n = len - 1; // (Max) degree of the polynomial 2267 for (j = 0; j < len; j++) { 2268 coeffs_sum[j] = 0; 2269 } 2270 2271 for (i = 0; i < len; i++) { 2272 c = points[i].Y(); 2273 p = points[i].X(); 2274 zeroes = []; 2275 for (j = 0; j < len; j++) { 2276 if (j !== i) { 2277 c /= p - points[j].X(); 2278 zeroes.push(points[j].X()); 2279 } 2280 } 2281 coeffs = [1].concat(Mat.Vieta(zeroes)); 2282 for (j = 0; j < coeffs.length; j++) { 2283 coeffs_sum[j] += (j % 2 === 1 ? -1 : 1) * coeffs[j] * c; 2284 } 2285 } 2286 2287 return coeffs_sum; 2288 }; 2289 }, 2290 2291 /** 2292 * Determine the coefficients of a cardinal spline polynom, See 2293 * https://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections 2294 * @param {Number} x1 point 1 2295 * @param {Number} x2 point 2 2296 * @param {Number} t1 tangent slope 1 2297 * @param {Number} t2 tangent slope 2 2298 * @return {Array} coefficents array c for the polynomial t maps to 2299 * c[0] + c[1]*t + c[2]*t*t + c[3]*t*t*t 2300 */ 2301 _initCubicPoly: function (x1, x2, t1, t2) { 2302 return [x1, t1, -3 * x1 + 3 * x2 - 2 * t1 - t2, 2 * x1 - 2 * x2 + t1 + t2]; 2303 }, 2304 2305 /** 2306 * Computes the cubic cardinal spline curve through a given set of points. The curve 2307 * is uniformly parametrized. 2308 * Two artificial control points at the beginning and the end are added. 2309 * 2310 * The implementation (especially the centripetal parametrization) is from 2311 * https://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections . 2312 * @param {Array} points Array consisting of JXG.Points. 2313 * @param {Number|Function} tau The tension parameter, either a constant number or a function returning a number. This number is between 0 and 1. 2314 * tau=1/2 give Catmull-Rom splines. 2315 * @param {String} type (Optional) parameter which allows to choose between "uniform" (default) and 2316 * "centripetal" parameterization. Thus the two possible values are "uniform" or "centripetal". 2317 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 2318 * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, 2319 * and a function simply returning the length of the points array 2320 * minus three. 2321 * @memberof JXG.Math.Numerics 2322 */ 2323 CardinalSpline: function (points, tau_param, type) { 2324 var p, 2325 coeffs = [], 2326 makeFct, 2327 tau, _tau, 2328 that = this; 2329 2330 if (Type.isFunction(tau_param)) { 2331 _tau = tau_param; 2332 } else { 2333 _tau = function () { 2334 return tau_param; 2335 }; 2336 } 2337 2338 if (type === undefined) { 2339 type = "uniform"; 2340 } 2341 2342 /** @ignore */ 2343 makeFct = function (which) { 2344 return function (t, suspendedUpdate) { 2345 var s, 2346 c, 2347 // control point at the beginning and at the end 2348 first, 2349 last, 2350 t1, 2351 t2, 2352 dt0, 2353 dt1, 2354 dt2, 2355 // dx, dy, 2356 len; 2357 2358 if (points.length < 2) { 2359 return NaN; 2360 } 2361 2362 if (!suspendedUpdate) { 2363 tau = _tau(); 2364 2365 // New point list p: [first, points ..., last] 2366 first = { 2367 X: function () { 2368 return 2 * points[0].X() - points[1].X(); 2369 }, 2370 Y: function () { 2371 return 2 * points[0].Y() - points[1].Y(); 2372 }, 2373 Dist: function (p) { 2374 var dx = this.X() - p.X(), 2375 dy = this.Y() - p.Y(); 2376 return Mat.hypot(dx, dy); 2377 } 2378 }; 2379 2380 last = { 2381 X: function () { 2382 return ( 2383 2 * points[points.length - 1].X() - 2384 points[points.length - 2].X() 2385 ); 2386 }, 2387 Y: function () { 2388 return ( 2389 2 * points[points.length - 1].Y() - 2390 points[points.length - 2].Y() 2391 ); 2392 }, 2393 Dist: function (p) { 2394 var dx = this.X() - p.X(), 2395 dy = this.Y() - p.Y(); 2396 return Mat.hypot(dx, dy); 2397 } 2398 }; 2399 2400 p = [first].concat(points, [last]); 2401 len = p.length; 2402 2403 coeffs[which] = []; 2404 2405 for (s = 0; s < len - 3; s++) { 2406 if (type === "centripetal") { 2407 // The order is important, since p[0].coords === undefined 2408 dt0 = p[s].Dist(p[s + 1]); 2409 dt1 = p[s + 2].Dist(p[s + 1]); 2410 dt2 = p[s + 3].Dist(p[s + 2]); 2411 2412 dt0 = Math.sqrt(dt0); 2413 dt1 = Math.sqrt(dt1); 2414 dt2 = Math.sqrt(dt2); 2415 2416 if (dt1 < Mat.eps) { 2417 dt1 = 1.0; 2418 } 2419 if (dt0 < Mat.eps) { 2420 dt0 = dt1; 2421 } 2422 if (dt2 < Mat.eps) { 2423 dt2 = dt1; 2424 } 2425 2426 t1 = 2427 (p[s + 1][which]() - p[s][which]()) / dt0 - 2428 (p[s + 2][which]() - p[s][which]()) / (dt1 + dt0) + 2429 (p[s + 2][which]() - p[s + 1][which]()) / dt1; 2430 2431 t2 = 2432 (p[s + 2][which]() - p[s + 1][which]()) / dt1 - 2433 (p[s + 3][which]() - p[s + 1][which]()) / (dt2 + dt1) + 2434 (p[s + 3][which]() - p[s + 2][which]()) / dt2; 2435 2436 t1 *= dt1; 2437 t2 *= dt1; 2438 2439 coeffs[which][s] = that._initCubicPoly( 2440 p[s + 1][which](), 2441 p[s + 2][which](), 2442 tau * t1, 2443 tau * t2 2444 ); 2445 } else { 2446 coeffs[which][s] = that._initCubicPoly( 2447 p[s + 1][which](), 2448 p[s + 2][which](), 2449 tau * (p[s + 2][which]() - p[s][which]()), 2450 tau * (p[s + 3][which]() - p[s + 1][which]()) 2451 ); 2452 } 2453 } 2454 } 2455 2456 if (isNaN(t)) { 2457 return NaN; 2458 } 2459 2460 len = points.length; 2461 // This is necessary for our advanced plotting algorithm: 2462 if (t <= 0.0) { 2463 return points[0][which](); 2464 } 2465 if (t >= len) { 2466 return points[len - 1][which](); 2467 } 2468 2469 s = Math.floor(t); 2470 if (s === t) { 2471 return points[s][which](); 2472 } 2473 2474 t -= s; 2475 c = coeffs[which][s]; 2476 if (c === undefined) { 2477 return NaN; 2478 } 2479 2480 return ((c[3] * t + c[2]) * t + c[1]) * t + c[0]; 2481 }; 2482 }; 2483 2484 return [ 2485 makeFct("X"), 2486 makeFct("Y"), 2487 0, 2488 function () { 2489 return points.length - 1; 2490 } 2491 ]; 2492 }, 2493 2494 /** 2495 * Computes the cubic Catmull-Rom spline curve through a given set of points. The curve 2496 * is uniformly parametrized. The curve is the cardinal spline curve for tau=0.5. 2497 * Two artificial control points at the beginning and the end are added. 2498 * @param {Array} points Array consisting of JXG.Points. 2499 * @param {String} type (Optional) parameter which allows to choose between "uniform" (default) and 2500 * "centripetal" parameterization. Thus the two possible values are "uniform" or "centripetal". 2501 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 2502 * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, and a function simply 2503 * returning the length of the points array minus three. 2504 * @memberof JXG.Math.Numerics 2505 */ 2506 CatmullRomSpline: function (points, type) { 2507 return this.CardinalSpline(points, 0.5, type); 2508 }, 2509 2510 /** 2511 * Computes the regression polynomial of a given degree through a given set of coordinates. 2512 * Returns the regression polynomial function. 2513 * @param {Number|function|Slider} degree number, function or slider. 2514 * Either 2515 * @param {Array} dataX Array containing either the x-coordinates of the data set or both coordinates in 2516 * an array of {@link JXG.Point}s or {@link JXG.Coords}. 2517 * In the latter case, the <tt>dataY</tt> parameter will be ignored. 2518 * @param {Array} dataY Array containing the y-coordinates of the data set, 2519 * @returns {function} A function of one parameter which returns the value of the regression polynomial of the given degree. 2520 * It possesses the method getTerm() which returns the string containing the function term of the polynomial. 2521 * The function returned will throw an exception, if the data set is malformed. 2522 * @memberof JXG.Math.Numerics 2523 */ 2524 regressionPolynomial: function (degree, dataX, dataY) { 2525 var coeffs, deg, dX, dY, inputType, fct, 2526 term = ""; 2527 2528 // Slider 2529 if (Type.isPoint(degree) && Type.isFunction(degree.Value)) { 2530 /** @ignore */ 2531 deg = function () { 2532 return degree.Value(); 2533 }; 2534 // function 2535 } else if (Type.isFunction(degree)) { 2536 deg = degree; 2537 // number 2538 } else if (Type.isNumber(degree)) { 2539 /** @ignore */ 2540 deg = function () { 2541 return degree; 2542 }; 2543 } else { 2544 throw new Error( 2545 "JSXGraph: Can't create regressionPolynomial from degree of type'" + 2546 typeof degree + 2547 "'." 2548 ); 2549 } 2550 2551 // Parameters degree, dataX, dataY 2552 if (arguments.length === 3 && Type.isArray(dataX) && Type.isArray(dataY)) { 2553 inputType = 0; 2554 // Parameters degree, point array 2555 } else if ( 2556 arguments.length === 2 && 2557 Type.isArray(dataX) && 2558 dataX.length > 0 && 2559 Type.isPoint(dataX[0]) 2560 ) { 2561 inputType = 1; 2562 } else if ( 2563 arguments.length === 2 && 2564 Type.isArray(dataX) && 2565 dataX.length > 0 && 2566 dataX[0].usrCoords && 2567 dataX[0].scrCoords 2568 ) { 2569 inputType = 2; 2570 } else { 2571 throw new Error("JSXGraph: Can't create regressionPolynomial. Wrong parameters."); 2572 } 2573 2574 /** @ignore */ 2575 fct = function (x, suspendedUpdate) { 2576 var i, j, 2577 M, MT, y, B, c, s, d, 2578 // input data 2579 len = dataX.length; 2580 2581 d = Math.floor(deg()); 2582 2583 if (!suspendedUpdate) { 2584 // point list as input 2585 if (inputType === 1) { 2586 dX = []; 2587 dY = []; 2588 2589 for (i = 0; i < len; i++) { 2590 dX[i] = dataX[i].X(); 2591 dY[i] = dataX[i].Y(); 2592 } 2593 } 2594 2595 if (inputType === 2) { 2596 dX = []; 2597 dY = []; 2598 2599 for (i = 0; i < len; i++) { 2600 dX[i] = dataX[i].usrCoords[1]; 2601 dY[i] = dataX[i].usrCoords[2]; 2602 } 2603 } 2604 2605 // check for functions 2606 if (inputType === 0) { 2607 dX = []; 2608 dY = []; 2609 2610 for (i = 0; i < len; i++) { 2611 if (Type.isFunction(dataX[i])) { 2612 dX.push(dataX[i]()); 2613 } else { 2614 dX.push(dataX[i]); 2615 } 2616 2617 if (Type.isFunction(dataY[i])) { 2618 dY.push(dataY[i]()); 2619 } else { 2620 dY.push(dataY[i]); 2621 } 2622 } 2623 } 2624 2625 M = []; 2626 for (j = 0; j < len; j++) { 2627 M.push([1]); 2628 } 2629 for (i = 1; i <= d; i++) { 2630 for (j = 0; j < len; j++) { 2631 M[j][i] = M[j][i - 1] * dX[j]; 2632 } 2633 } 2634 2635 y = dY; 2636 MT = Mat.transpose(M); 2637 B = Mat.matMatMult(MT, M); 2638 c = Mat.matVecMult(MT, y); 2639 coeffs = Mat.Numerics.Gauss(B, c); 2640 term = Mat.Numerics.generatePolynomialTerm(coeffs, d, "x", 3); 2641 } 2642 2643 // Horner's scheme to evaluate polynomial 2644 s = coeffs[d]; 2645 for (i = d - 1; i >= 0; i--) { 2646 s = s * x + coeffs[i]; 2647 } 2648 2649 return s; 2650 }; 2651 2652 /** @ignore */ 2653 fct.getTerm = function () { 2654 return term; 2655 }; 2656 2657 return fct; 2658 }, 2659 2660 /** 2661 * Computes the cubic Bezier curve through a given set of points. 2662 * @param {Array} points Array consisting of 3*k+1 {@link JXG.Points}. 2663 * The points at position k with k mod 3 = 0 are the data points, 2664 * points at position k with k mod 3 = 1 or 2 are the control points. 2665 * @returns {Array} An array consisting of two functions of one parameter t which return the 2666 * x resp. y coordinates of the Bezier curve in t, one zero value, and a third function accepting 2667 * no parameters and returning one third of the length of the points. 2668 * @memberof JXG.Math.Numerics 2669 */ 2670 bezier: function (points) { 2671 var len, 2672 flen, 2673 /** @ignore */ 2674 makeFct = function (which) { 2675 return function (t, suspendedUpdate) { 2676 var z = Math.floor(t) * 3, 2677 t0 = t % 1, 2678 t1 = 1 - t0; 2679 2680 if (!suspendedUpdate) { 2681 flen = 3 * Math.floor((points.length - 1) / 3); 2682 len = Math.floor(flen / 3); 2683 } 2684 2685 if (t < 0) { 2686 return points[0][which](); 2687 } 2688 2689 if (t >= len) { 2690 return points[flen][which](); 2691 } 2692 2693 if (isNaN(t)) { 2694 return NaN; 2695 } 2696 2697 return ( 2698 t1 * t1 * (t1 * points[z][which]() + 3 * t0 * points[z + 1][which]()) + 2699 (3 * t1 * points[z + 2][which]() + t0 * points[z + 3][which]()) * 2700 t0 * 2701 t0 2702 ); 2703 }; 2704 }; 2705 2706 return [ 2707 makeFct("X"), 2708 makeFct("Y"), 2709 0, 2710 function () { 2711 return Math.floor(points.length / 3); 2712 } 2713 ]; 2714 }, 2715 2716 /** 2717 * Computes the B-spline curve of order k (order = degree+1) through a given set of points. 2718 * @param {Array} points Array consisting of JXG.Points. 2719 * @param {Number} order Order of the B-spline curve. 2720 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 2721 * which return the x resp. y coordinates of the B-spline curve in t, a zero value, and a function simply 2722 * returning the length of the points array minus one. 2723 * @memberof JXG.Math.Numerics 2724 */ 2725 bspline: function (points, order) { 2726 var knots, 2727 _knotVector = function (n, k) { 2728 var j, 2729 kn = []; 2730 2731 for (j = 0; j < n + k + 1; j++) { 2732 if (j < k) { 2733 kn[j] = 0.0; 2734 } else if (j <= n) { 2735 kn[j] = j - k + 1; 2736 } else { 2737 kn[j] = n - k + 2; 2738 } 2739 } 2740 2741 return kn; 2742 }, 2743 _evalBasisFuncs = function (t, kn, k, s) { 2744 var i, 2745 j, 2746 a, 2747 b, 2748 den, 2749 N = []; 2750 2751 if (kn[s] <= t && t < kn[s + 1]) { 2752 N[s] = 1; 2753 } else { 2754 N[s] = 0; 2755 } 2756 2757 for (i = 2; i <= k; i++) { 2758 for (j = s - i + 1; j <= s; j++) { 2759 if (j <= s - i + 1 || j < 0) { 2760 a = 0.0; 2761 } else { 2762 a = N[j]; 2763 } 2764 2765 if (j >= s) { 2766 b = 0.0; 2767 } else { 2768 b = N[j + 1]; 2769 } 2770 2771 den = kn[j + i - 1] - kn[j]; 2772 2773 if (den === 0) { 2774 N[j] = 0; 2775 } else { 2776 N[j] = ((t - kn[j]) / den) * a; 2777 } 2778 2779 den = kn[j + i] - kn[j + 1]; 2780 2781 if (den !== 0) { 2782 N[j] += ((kn[j + i] - t) / den) * b; 2783 } 2784 } 2785 } 2786 return N; 2787 }, 2788 /** @ignore */ 2789 makeFct = function (which) { 2790 return function (t, suspendedUpdate) { 2791 var y, 2792 j, 2793 s, 2794 N = [], 2795 len = points.length, 2796 n = len - 1, 2797 k = order; 2798 2799 if (n <= 0) { 2800 return NaN; 2801 } 2802 2803 if (n + 2 <= k) { 2804 k = n + 1; 2805 } 2806 2807 if (t <= 0) { 2808 return points[0][which](); 2809 } 2810 2811 if (t >= n - k + 2) { 2812 return points[n][which](); 2813 } 2814 2815 s = Math.floor(t) + k - 1; 2816 knots = _knotVector(n, k); 2817 N = _evalBasisFuncs(t, knots, k, s); 2818 2819 y = 0.0; 2820 for (j = s - k + 1; j <= s; j++) { 2821 if (j < len && j >= 0) { 2822 y += points[j][which]() * N[j]; 2823 } 2824 } 2825 2826 return y; 2827 }; 2828 }; 2829 2830 return [ 2831 makeFct("X"), 2832 makeFct("Y"), 2833 0, 2834 function () { 2835 return points.length - 1; 2836 } 2837 ]; 2838 }, 2839 2840 /** 2841 * Numerical (symmetric) approximation of derivative. suspendUpdate is piped through, 2842 * see {@link JXG.Curve#updateCurve} 2843 * and {@link JXG.Curve#hasPoint}. 2844 * @param {function} f Function in one variable to be differentiated. 2845 * @param {object} [obj] Optional object that is treated as "this" in the function body. This is useful, if the function is a 2846 * method of an object and contains a reference to its parent object via "this". 2847 * @returns {function} Derivative function of a given function f. 2848 * @memberof JXG.Math.Numerics 2849 */ 2850 D: function (f, obj) { 2851 if (!Type.exists(obj)) { 2852 return function (x, suspendedUpdate) { 2853 var h = 0.00001, 2854 h2 = h * 2.0; 2855 2856 // Experiments with Richardsons rule 2857 /* 2858 var phi = (f(x + h, suspendedUpdate) - f(x - h, suspendedUpdate)) / h2; 2859 var phi2; 2860 h *= 0.5; 2861 h2 *= 0.5; 2862 phi2 = (f(x + h, suspendedUpdate) - f(x - h, suspendedUpdate)) / h2; 2863 2864 return phi2 + (phi2 - phi) / 3.0; 2865 */ 2866 return (f(x + h, suspendedUpdate) - f(x - h, suspendedUpdate)) / h2; 2867 }; 2868 } 2869 2870 return function (x, suspendedUpdate) { 2871 var h = 0.00001, 2872 h2 = h * 2.0; 2873 2874 return ( 2875 (f.apply(obj, [x + h, suspendedUpdate]) - 2876 f.apply(obj, [x - h, suspendedUpdate])) / 2877 h2 2878 ); 2879 }; 2880 }, 2881 2882 /** 2883 * Evaluate the function term for {@see #riemann}. 2884 * @private 2885 * @param {Number} x function argument 2886 * @param {function} f JavaScript function returning a number 2887 * @param {String} type Name of the Riemann sum type, e.g. 'lower', see {@see #riemann}. 2888 * @param {Number} delta Width of the bars in user coordinates 2889 * @returns {Number} Upper (delta > 0) or lower (delta < 0) value of the bar containing x of the Riemann sum. 2890 * 2891 * @memberof JXG.Math.Numerics 2892 */ 2893 _riemannValue: function (x, f, type, delta) { 2894 var y, y1, x1, delta1; 2895 2896 if (delta < 0) { 2897 // delta is negative if the lower function term is evaluated 2898 if (type !== "trapezoidal") { 2899 x = x + delta; 2900 } 2901 delta *= -1; 2902 if (type === "lower") { 2903 type = "upper"; 2904 } else if (type === "upper") { 2905 type = "lower"; 2906 } 2907 } 2908 2909 delta1 = delta * 0.01; // for 'lower' and 'upper' 2910 2911 if (type === "right") { 2912 y = f(x + delta); 2913 } else if (type === "middle") { 2914 y = f(x + delta * 0.5); 2915 } else if (type === "left" || type === "trapezoidal") { 2916 y = f(x); 2917 } else if (type === "lower") { 2918 y = f(x); 2919 2920 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 2921 y1 = f(x1); 2922 2923 if (y1 < y) { 2924 y = y1; 2925 } 2926 } 2927 2928 y1 = f(x + delta); 2929 if (y1 < y) { 2930 y = y1; 2931 } 2932 } else if (type === "upper") { 2933 y = f(x); 2934 2935 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 2936 y1 = f(x1); 2937 if (y1 > y) { 2938 y = y1; 2939 } 2940 } 2941 2942 y1 = f(x + delta); 2943 if (y1 > y) { 2944 y = y1; 2945 } 2946 } else if (type === "random") { 2947 y = f(x + delta * Math.random()); 2948 } else if (type === "simpson") { 2949 y = (f(x) + 4 * f(x + delta * 0.5) + f(x + delta)) / 6.0; 2950 } else { 2951 y = f(x); // default is lower 2952 } 2953 2954 return y; 2955 }, 2956 2957 /** 2958 * Helper function to create curve which displays Riemann sums. 2959 * Compute coordinates for the rectangles showing the Riemann sum. 2960 * <p> 2961 * In case of type "simpson" and "trapezoidal", the horizontal line approximating the function value 2962 * is replaced by a parabola or a secant. IN case of "simpson", 2963 * the parabola is approximated visually by a polygonal chain of fixed step width. 2964 * 2965 * @param {Function|Array} f Function or array of two functions. 2966 * If f is a function the integral of this function is approximated by the Riemann sum. 2967 * If f is an array consisting of two functions the area between the two functions is filled 2968 * by the Riemann sum bars. 2969 * @param {Number} n number of rectangles. 2970 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson', or 'trapezoidal'. 2971 * "simpson" is Simpson's 1/3 rule. 2972 * @param {Number} start Left border of the approximation interval 2973 * @param {Number} end Right border of the approximation interval 2974 * @returns {Array} An array of two arrays containing the x and y coordinates for the rectangles showing the Riemann sum. This 2975 * 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 2976 * rectangles. 2977 * @memberof JXG.Math.Numerics 2978 */ 2979 riemann: function (gf, n, type, start, end) { 2980 var i, delta, 2981 k, a, b, c, f0, f1, f2, xx, h, 2982 steps = 30, // Fixed step width for Simpson's rule 2983 xarr = [], 2984 yarr = [], 2985 x = start, 2986 sum = 0, 2987 y, f, g; 2988 2989 if (Type.isArray(gf)) { 2990 g = gf[0]; 2991 f = gf[1]; 2992 } else { 2993 f = gf; 2994 } 2995 2996 n = Math.floor(n); 2997 2998 if (n <= 0) { 2999 return [xarr, yarr, sum]; 3000 } 3001 3002 delta = (end - start) / n; 3003 3004 // "Upper" horizontal line defined by function 3005 for (i = 0; i < n; i++) { 3006 if (type === "simpson") { 3007 sum += this._riemannValue(x, f, type, delta) * delta; 3008 3009 h = delta * 0.5; 3010 f0 = f(x); 3011 f1 = f(x + h); 3012 f2 = f(x + 2 * h); 3013 3014 a = (f2 + f0 - 2 * f1) / (h * h) * 0.5; 3015 b = (f2 - f0) / (2 * h); 3016 c = f1; 3017 for (k = 0; k < steps; k++) { 3018 xx = k * delta / steps - h; 3019 xarr.push(x + xx + h); 3020 yarr.push(a * xx * xx + b * xx + c); 3021 } 3022 x += delta; 3023 y = f2; 3024 } else { 3025 y = this._riemannValue(x, f, type, delta); 3026 xarr.push(x); 3027 yarr.push(y); 3028 3029 x += delta; 3030 if (type === "trapezoidal") { 3031 f2 = f(x); 3032 sum += (y + f2) * 0.5 * delta; 3033 y = f2; 3034 } else { 3035 sum += y * delta; 3036 } 3037 3038 xarr.push(x); 3039 yarr.push(y); 3040 } 3041 xarr.push(x); 3042 yarr.push(y); 3043 } 3044 3045 // "Lower" horizontal line 3046 // Go backwards 3047 for (i = 0; i < n; i++) { 3048 if (type === "simpson" && g) { 3049 sum -= this._riemannValue(x, g, type, -delta) * delta; 3050 3051 h = delta * 0.5; 3052 f0 = g(x); 3053 f1 = g(x - h); 3054 f2 = g(x - 2 * h); 3055 3056 a = (f2 + f0 - 2 * f1) / (h * h) * 0.5; 3057 b = (f2 - f0) / (2 * h); 3058 c = f1; 3059 for (k = 0; k < steps; k++) { 3060 xx = k * delta / steps - h; 3061 xarr.push(x - xx - h); 3062 yarr.push(a * xx * xx + b * xx + c); 3063 } 3064 x -= delta; 3065 y = f2; 3066 } else { 3067 if (g) { 3068 y = this._riemannValue(x, g, type, -delta); 3069 } else { 3070 y = 0.0; 3071 } 3072 xarr.push(x); 3073 yarr.push(y); 3074 3075 x -= delta; 3076 if (g) { 3077 if (type === "trapezoidal") { 3078 f2 = g(x); 3079 sum -= (y + f2) * 0.5 * delta; 3080 y = f2; 3081 } else { 3082 sum -= y * delta; 3083 } 3084 } 3085 } 3086 xarr.push(x); 3087 yarr.push(y); 3088 3089 // Draw the vertical lines 3090 xarr.push(x); 3091 yarr.push(f(x)); 3092 } 3093 3094 return [xarr, yarr, sum]; 3095 }, 3096 3097 /** 3098 * Approximate the integral by Riemann sums. 3099 * Compute the area described by the riemann sum rectangles. 3100 * 3101 * If there is an element of type {@link Riemannsum}, then it is more efficient 3102 * to use the method JXG.Curve.Value() of this element instead. 3103 * 3104 * @param {Function_Array} f Function or array of two functions. 3105 * If f is a function the integral of this function is approximated by the Riemann sum. 3106 * If f is an array consisting of two functions the area between the two functions is approximated 3107 * by the Riemann sum. 3108 * @param {Number} n number of rectangles. 3109 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson' or 'trapezoidal'. 3110 * 3111 * @param {Number} start Left border of the approximation interval 3112 * @param {Number} end Right border of the approximation interval 3113 * @returns {Number} The sum of the areas of the rectangles. 3114 * @memberof JXG.Math.Numerics 3115 */ 3116 riemannsum: function (f, n, type, start, end) { 3117 JXG.deprecated("Numerics.riemannsum()", "Numerics.riemann()[2]"); 3118 return this.riemann(f, n, type, start, end)[2]; 3119 }, 3120 3121 /** 3122 * Solve initial value problems numerically using <i>explicit</i> Runge-Kutta methods. 3123 * See {@link https://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm. 3124 * @param {object|String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing 3125 * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure 3126 * <pre> 3127 * { 3128 * s: <Number>, 3129 * A: <matrix>, 3130 * b: <Array>, 3131 * c: <Array> 3132 * } 3133 * </pre> 3134 * which corresponds to the Butcher tableau structure 3135 * shown here: https://en.wikipedia.org/w/index.php?title=List_of_Runge%E2%80%93Kutta_methods&oldid=357796696 . 3136 * <i>Default</i> is 'euler'. 3137 * @param {Array} x0 Initial value vector. Even if the problem is one-dimensional, the initial value has to be given in an array. 3138 * @param {Array} I Interval on which to integrate. 3139 * @param {Number} N Number of integration intervals, i.e. there are <i>N+1</i> evaluation points. 3140 * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode 3141 * 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 3142 * vector <tt>x</tt>, and has to return a vector of the same length as <tt>x</tt> has. 3143 * @returns {Array} An array of vectors describing the solution of the ode on the given interval I. 3144 * @example 3145 * // A very simple autonomous system dx(t)/dt = x(t); 3146 * var f = function(t, x) { 3147 * return [x[0]]; 3148 * } 3149 * 3150 * // Solve it with initial value x(0) = 1 on the interval [0, 2] 3151 * // with 20 evaluation points. 3152 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 3153 * 3154 * // Prepare data for plotting the solution of the ode using a curve. 3155 * var dataX = []; 3156 * var dataY = []; 3157 * var h = 0.1; // (I[1] - I[0])/N = (2-0)/20 3158 * var i; 3159 * for(i=0; i<data.length; i++) { 3160 * dataX[i] = i*h; 3161 * dataY[i] = data[i][0]; 3162 * } 3163 * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'}); 3164 * </pre><div class="jxgbox" id="JXGd2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div> 3165 * <script type="text/javascript"> 3166 * var board = JXG.JSXGraph.initBoard('JXGd2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false}); 3167 * var f = function(t, x) { 3168 * // we have to copy the value. 3169 * // return x; would just return the reference. 3170 * return [x[0]]; 3171 * } 3172 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 3173 * var dataX = []; 3174 * var dataY = []; 3175 * var h = 0.1; 3176 * for(var i=0; i<data.length; i++) { 3177 * dataX[i] = i*h; 3178 * dataY[i] = data[i][0]; 3179 * } 3180 * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'}); 3181 * </script><pre> 3182 * @memberof JXG.Math.Numerics 3183 */ 3184 rungeKutta: function (butcher, x0, I, N, f) { 3185 var e, 3186 i, j, k, l, s, 3187 x = [], 3188 y = [], 3189 h = (I[1] - I[0]) / N, 3190 t = I[0], 3191 dim = x0.length, 3192 result = [], 3193 r = 0; 3194 3195 if (Type.isString(butcher)) { 3196 butcher = predefinedButcher[butcher] || predefinedButcher.euler; 3197 } 3198 s = butcher.s; 3199 3200 // Don't change x0, so copy it 3201 x = x0.slice(); 3202 3203 for (i = 0; i <= N; i++) { 3204 result[r] = x.slice(); 3205 3206 r++; 3207 k = []; 3208 3209 for (j = 0; j < s; j++) { 3210 // Init y = 0 3211 for (e = 0; e < dim; e++) { 3212 y[e] = 0.0; 3213 } 3214 3215 // Calculate linear combination of former k's and save it in y 3216 for (l = 0; l < j; l++) { 3217 for (e = 0; e < dim; e++) { 3218 y[e] += butcher.A[j][l] * h * k[l][e]; 3219 } 3220 } 3221 3222 // Add x(t) to y 3223 for (e = 0; e < dim; e++) { 3224 y[e] += x[e]; 3225 } 3226 3227 // Calculate new k and add it to the k matrix 3228 k.push(f(t + butcher.c[j] * h, y)); 3229 } 3230 3231 // Init y = 0 3232 for (e = 0; e < dim; e++) { 3233 y[e] = 0.0; 3234 } 3235 3236 for (l = 0; l < s; l++) { 3237 for (e = 0; e < dim; e++) { 3238 y[e] += butcher.b[l] * k[l][e]; 3239 } 3240 } 3241 3242 for (e = 0; e < dim; e++) { 3243 x[e] = x[e] + h * y[e]; 3244 } 3245 3246 t += h; 3247 } 3248 3249 return result; 3250 }, 3251 3252 /** 3253 * Maximum number of iterations in {@link JXG.Math.Numerics.fzero} and 3254 * {@link JXG.Math.Numerics.chandrupatla} 3255 * @type Number 3256 * @default 80 3257 * @memberof JXG.Math.Numerics 3258 */ 3259 maxIterationsRoot: 80, 3260 3261 /** 3262 * Maximum number of iterations in {@link JXG.Math.Numerics.fminbr} 3263 * @type Number 3264 * @default 500 3265 * @memberof JXG.Math.Numerics 3266 */ 3267 maxIterationsMinimize: 500, 3268 3269 /** 3270 * Given a number x_0, this function tries to find a second number x_1 such that 3271 * the function f has opposite signs at x_0 and x_1. 3272 * The return values have to be tested if the method succeeded. 3273 * 3274 * @param {Function} f Function, whose root is to be found 3275 * @param {Number} x0 Start value 3276 * @param {Object} [context] Parent object in case f is method of it 3277 * @returns {Array} [x_0, f(x_0), x_1, f(x_1)] in case that x_0 <= x_1 3278 * or [x_1, f(x_1), x_0, f(x_0)] in case that x_1 < x_0. 3279 * 3280 * @see JXG.Math.Numerics.fzero 3281 * @see JXG.Math.Numerics.chandrupatla 3282 * 3283 * @memberof JXG.Math.Numerics 3284 */ 3285 findBracket: function (f, x0, context) { 3286 var a, aa, fa, blist, b, fb, u, fu, i, len; 3287 3288 if (Type.isArray(x0)) { 3289 return x0; 3290 } 3291 3292 a = x0; 3293 fa = f.call(context, a); 3294 // nfev += 1; 3295 3296 // Try to get b, by trying several values related to a 3297 aa = a === 0 ? 1 : a; 3298 blist = [ 3299 a - 0.1 * aa, 3300 a + 0.1 * aa, 3301 a - 1, 3302 a + 1, 3303 a - 0.5 * aa, 3304 a + 0.5 * aa, 3305 a - 0.6 * aa, 3306 a + 0.6 * aa, 3307 a - 1 * aa, 3308 a + 1 * aa, 3309 a - 2 * aa, 3310 a + 2 * aa, 3311 a - 5 * aa, 3312 a + 5 * aa, 3313 a - 10 * aa, 3314 a + 10 * aa, 3315 a - 50 * aa, 3316 a + 50 * aa, 3317 a - 100 * aa, 3318 a + 100 * aa 3319 ]; 3320 len = blist.length; 3321 3322 for (i = 0; i < len; i++) { 3323 b = blist[i]; 3324 fb = f.call(context, b); 3325 // nfev += 1; 3326 3327 if (fa * fb <= 0) { 3328 break; 3329 } 3330 } 3331 if (b < a) { 3332 u = a; 3333 a = b; 3334 b = u; 3335 3336 fu = fa; 3337 fa = fb; 3338 fb = fu; 3339 } 3340 return [a, fa, b, fb]; 3341 }, 3342 3343 /** 3344 * 3345 * Find zero of an univariate function f. 3346 * @param {function} f Function, whose root is to be found 3347 * @param {Array|Number} x0 Start value or start interval enclosing the root 3348 * @param {Object} [context] Parent object in case f is method of it 3349 * @returns {Number} the approximation of the root 3350 * Algorithm: 3351 * Brent's root finder from 3352 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 3353 * computations. M., Mir, 1980, p.180 of the Russian edition 3354 * https://www.netlib.org/c/brent.shar 3355 * 3356 * If x0 is an array containing lower and upper bound for the zero 3357 * algorithm 748 is applied. Otherwise, if x0 is a number, 3358 * the algorithm tries to bracket a zero of f starting from x0. 3359 * If this fails, we fall back to Newton's method. 3360 * 3361 * @see JXG.Math.Numerics.chandrupatla 3362 * @see JXG.Math.Numerics.root 3363 * @memberof JXG.Math.Numerics 3364 */ 3365 fzero: function (f, x0, context) { 3366 var a, b, c, 3367 fa, fb, fc, 3368 res, x00, 3369 prev_step, 3370 t1, t2, 3371 cb, 3372 tol_act, // Actual tolerance 3373 p, // Interpolation step is calculated in the form p/q; division 3374 q, // operations is delayed until the last moment 3375 new_step, // Step at this iteration 3376 eps = Mat.eps, 3377 maxiter = this.maxIterationsRoot, 3378 niter = 0; 3379 // nfev = 0; 3380 3381 if (Type.isArray(x0)) { 3382 if (x0.length < 2) { 3383 throw new Error( 3384 "JXG.Math.Numerics.fzero: length of array x0 has to be at least two." 3385 ); 3386 } 3387 3388 x00 = this.findDomain(f, x0, context); 3389 a = x00[0]; 3390 b = x00[1]; 3391 // a = x0[0]; 3392 // b = x0[1]; 3393 3394 fa = f.call(context, a); 3395 // nfev += 1; 3396 fb = f.call(context, b); 3397 // nfev += 1; 3398 } else { 3399 res = this.findBracket(f, x0, context); 3400 a = res[0]; 3401 fa = res[1]; 3402 b = res[2]; 3403 fb = res[3]; 3404 } 3405 3406 if (Math.abs(fa) <= eps) { 3407 return a; 3408 } 3409 if (Math.abs(fb) <= eps) { 3410 return b; 3411 } 3412 3413 if (fa * fb > 0) { 3414 // Bracketing not successful, fall back to Newton's method or to fminbr 3415 if (Type.isArray(x0)) { 3416 return this.fminbr(f, [a, b], context); 3417 } 3418 3419 return this.Newton(f, a, context); 3420 } 3421 3422 // OK, we have enclosed a zero of f. 3423 // Now we can start Brent's method 3424 c = a; 3425 fc = fa; 3426 3427 // Main iteration loop 3428 while (niter < maxiter) { 3429 // Distance from the last but one to the last approximation 3430 prev_step = b - a; 3431 3432 // Swap data for b to be the best approximation 3433 if (Math.abs(fc) < Math.abs(fb)) { 3434 a = b; 3435 b = c; 3436 c = a; 3437 3438 fa = fb; 3439 fb = fc; 3440 fc = fa; 3441 } 3442 tol_act = 2 * eps * Math.abs(b) + eps * 0.5; 3443 new_step = (c - b) * 0.5; 3444 3445 if (Math.abs(new_step) <= tol_act || Math.abs(fb) <= eps) { 3446 // Acceptable approx. is found 3447 return b; 3448 } 3449 3450 // Decide if the interpolation can be tried 3451 // If prev_step was large enough and was in true direction Interpolatiom may be tried 3452 if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) { 3453 cb = c - b; 3454 3455 // If we have only two distinct points linear interpolation can only be applied 3456 if (a === c) { 3457 t1 = fb / fa; 3458 p = cb * t1; 3459 q = 1.0 - t1; 3460 // Quadric inverse interpolation 3461 } else { 3462 q = fa / fc; 3463 t1 = fb / fc; 3464 t2 = fb / fa; 3465 3466 p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0)); 3467 q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0); 3468 } 3469 3470 // p was calculated with the opposite sign; make p positive 3471 if (p > 0) { 3472 q = -q; 3473 // and assign possible minus to q 3474 } else { 3475 p = -p; 3476 } 3477 3478 // If b+p/q falls in [b,c] and isn't too large it is accepted 3479 // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent 3480 if ( 3481 p < 0.75 * cb * q - Math.abs(tol_act * q) * 0.5 && 3482 p < Math.abs(prev_step * q * 0.5) 3483 ) { 3484 new_step = p / q; 3485 } 3486 } 3487 3488 // Adjust the step to be not less than tolerance 3489 if (Math.abs(new_step) < tol_act) { 3490 new_step = new_step > 0 ? tol_act : -tol_act; 3491 } 3492 3493 // Save the previous approx. 3494 a = b; 3495 fa = fb; 3496 b += new_step; 3497 fb = f.call(context, b); 3498 // Do step to a new approxim. 3499 // nfev += 1; 3500 3501 // Adjust c for it to have a sign opposite to that of b 3502 if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) { 3503 c = a; 3504 fc = fa; 3505 } 3506 niter++; 3507 } // End while 3508 3509 return b; 3510 }, 3511 3512 /** 3513 * Find zero of an univariate function f. 3514 * @param {function} f Function, whose root is to be found 3515 * @param {Array|Number} x0 Start value or start interval enclosing the root 3516 * @param {Object} [context] Parent object in case f is method of it 3517 * @returns {Number} the approximation of the root 3518 * Algorithm: 3519 * Chandrupatla's method, see 3520 * Tirupathi R. Chandrupatla, 3521 * "A new hybrid quadratic/bisection algorithm for finding the zero of a nonlinear function without using derivatives", 3522 * Advances in Engineering Software, Volume 28, Issue 3, April 1997, Pages 145-149. 3523 * 3524 * If x0 is an array containing lower and upper bound for the zero 3525 * algorithm 748 is applied. Otherwise, if x0 is a number, 3526 * the algorithm tries to bracket a zero of f starting from x0. 3527 * If this fails, we fall back to Newton's method. 3528 * 3529 * @see JXG.Math.Numerics.root 3530 * @see JXG.Math.Numerics.fzero 3531 * @memberof JXG.Math.Numerics 3532 */ 3533 chandrupatla: function (f, x0, context) { 3534 var a, b, fa, fb, 3535 res, 3536 niter = 0, 3537 maxiter = this.maxIterationsRoot, 3538 rand = 1 + Math.random() * 0.001, 3539 t = 0.5 * rand, 3540 eps = Mat.eps, // 1.e-10, 3541 dlt = 0.00001, 3542 x1, x2, x3, x, 3543 f1, f2, f3, y, 3544 xm, fm, 3545 tol, tl, 3546 xi, ph, fl, fh, 3547 AL, A, B, C, D; 3548 3549 if (Type.isArray(x0)) { 3550 if (x0.length < 2) { 3551 throw new Error( 3552 "JXG.Math.Numerics.fzero: length of array x0 has to be at least two." 3553 ); 3554 } 3555 3556 a = x0[0]; 3557 fa = f.call(context, a); 3558 // nfev += 1; 3559 b = x0[1]; 3560 fb = f.call(context, b); 3561 // nfev += 1; 3562 } else { 3563 res = this.findBracket(f, x0, context); 3564 a = res[0]; 3565 fa = res[1]; 3566 b = res[2]; 3567 fb = res[3]; 3568 } 3569 3570 if (fa * fb > 0) { 3571 // Bracketing not successful, fall back to Newton's method or to fminbr 3572 if (Type.isArray(x0)) { 3573 return this.fminbr(f, [a, b], context); 3574 } 3575 3576 return this.Newton(f, a, context); 3577 } 3578 3579 x1 = a; 3580 x2 = b; 3581 f1 = fa; 3582 f2 = fb; 3583 do { 3584 x = x1 + t * (x2 - x1); 3585 y = f.call(context, x); 3586 3587 // Arrange 2-1-3: 2-1 interval, 1 middle, 3 discarded point 3588 if (Math.sign(y) === Math.sign(f1)) { 3589 x3 = x1; 3590 x1 = x; 3591 f3 = f1; 3592 f1 = y; 3593 } else { 3594 x3 = x2; 3595 x2 = x1; 3596 f3 = f2; 3597 f2 = f1; 3598 } 3599 x1 = x; 3600 f1 = y; 3601 3602 xm = x1; 3603 fm = f1; 3604 if (Math.abs(f2) < Math.abs(f1)) { 3605 xm = x2; 3606 fm = f2; 3607 } 3608 tol = 2 * eps * Math.abs(xm) + 0.5 * dlt; 3609 tl = tol / Math.abs(x2 - x1); 3610 if (tl > 0.5 || fm === 0) { 3611 break; 3612 } 3613 // If inverse quadratic interpolation holds, use it 3614 xi = (x1 - x2) / (x3 - x2); 3615 ph = (f1 - f2) / (f3 - f2); 3616 fl = 1 - Math.sqrt(1 - xi); 3617 fh = Math.sqrt(xi); 3618 if (fl < ph && ph < fh) { 3619 AL = (x3 - x1) / (x2 - x1); 3620 A = f1 / (f2 - f1); 3621 B = f3 / (f2 - f3); 3622 C = f1 / (f3 - f1); 3623 D = f2 / (f3 - f2); 3624 t = A * B + C * D * AL; 3625 } else { 3626 t = 0.5 * rand; 3627 } 3628 // Adjust t away from the interval boundary 3629 if (t < tl) { 3630 t = tl; 3631 } 3632 if (t > 1 - tl) { 3633 t = 1 - tl; 3634 } 3635 niter++; 3636 } while (niter <= maxiter); 3637 // console.log(niter); 3638 3639 return xm; 3640 }, 3641 3642 /** 3643 * Find a small enclosing interval of the domain of a function by 3644 * tightening the input interval x0. 3645 * <p> 3646 * This is a helper function which is used in {@link JXG.Math.Numerics.fminbr} 3647 * and {@link JXG.Math.Numerics.fzero} 3648 * to avoid search in an interval where the function is mostly undefined. 3649 * 3650 * @param {function} f 3651 * @param {Array} x0 Start interval 3652 * @param {Object} context Parent object in case f is method of it 3653 * @returns Array 3654 * 3655 * @example 3656 * var f = (x) => Math.sqrt(x); 3657 * console.log(JXG.Math.Numerics.findDomain(f, [-5, 5])); 3658 * 3659 * // Output: [ -0.00020428174445492973, 5 ] 3660 */ 3661 findDomain: function (f, x0, context) { 3662 var a, b, c, fc, 3663 x, 3664 gr = 1 - 1 / 1.61803398875, 3665 eps = 0.001, 3666 cnt, 3667 max_cnt = 20; 3668 3669 if (!Type.isArray(x0) || x0.length < 2) { 3670 throw new Error( 3671 "JXG.Math.Numerics.findDomain: length of array x0 has to be at least two." 3672 ); 3673 } 3674 3675 x = x0.slice(); 3676 a = x[0]; 3677 b = x[1]; 3678 fc = f.call(context, a); 3679 if (isNaN(fc)) { 3680 // Divide the interval with the golden ratio 3681 // and keep a such that f(a) = NaN 3682 cnt = 0; 3683 while (b - a > eps && cnt < max_cnt) { 3684 c = (b - a) * gr + a; 3685 fc = f.call(context, c); 3686 if (isNaN(fc)) { 3687 a = c; 3688 } else { 3689 b = c; 3690 } 3691 cnt++; 3692 } 3693 x[0] = a; 3694 } 3695 3696 a = x[0]; 3697 b = x[1]; 3698 fc = f.call(context, b); 3699 if (isNaN(fc)) { 3700 // Divide the interval with the golden ratio 3701 // and keep b such that f(b) = NaN 3702 cnt = 0; 3703 while (b - a > eps && cnt < max_cnt) { 3704 c = b - (b - a) * gr; 3705 fc = f.call(context, c); 3706 if (isNaN(fc)) { 3707 b = c; 3708 } else { 3709 a = c; 3710 } 3711 cnt++; 3712 } 3713 x[0] = b; 3714 } 3715 return x; 3716 }, 3717 3718 /** 3719 * 3720 * Find minimum of an univariate function f. 3721 * <p> 3722 * Algorithm: 3723 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 3724 * computations. M., Mir, 1980, p.180 of the Russian edition 3725 * 3726 * @param {function} f Function, whose minimum is to be found 3727 * @param {Array} x0 Start interval enclosing the minimum 3728 * @param {Object} [context] Parent object in case f is method of it 3729 * @returns {Number} the approximation of the minimum value position 3730 * @memberof JXG.Math.Numerics 3731 **/ 3732 fminbr: function (f, x0, context) { 3733 var a, b, x, v, w, 3734 fx, fv, fw, 3735 x00, 3736 range, middle_range, tol_act, new_step, 3737 p, q, t, ft, 3738 r = (3.0 - Math.sqrt(5.0)) * 0.5, // Golden section ratio 3739 tol = Mat.eps, 3740 sqrteps = Mat.eps, // Math.sqrt(Mat.eps), 3741 maxiter = this.maxIterationsMinimize, 3742 niter = 0; 3743 // nfev = 0; 3744 3745 if (!Type.isArray(x0) || x0.length < 2) { 3746 throw new Error( 3747 "JXG.Math.Numerics.fminbr: length of array x0 has to be at least two." 3748 ); 3749 } 3750 3751 x00 = this.findDomain(f, x0, context); 3752 a = x00[0]; 3753 b = x00[1]; 3754 v = a + r * (b - a); 3755 fv = f.call(context, v); 3756 3757 // First step - always gold section 3758 // nfev += 1; 3759 x = v; 3760 w = v; 3761 fx = fv; 3762 fw = fv; 3763 3764 while (niter < maxiter) { 3765 // Range over the interval in which we are looking for the minimum 3766 range = b - a; 3767 middle_range = (a + b) * 0.5; 3768 3769 // Actual tolerance 3770 tol_act = sqrteps * Math.abs(x) + tol / 3.0; 3771 3772 if (Math.abs(x - middle_range) + range * 0.5 <= 2.0 * tol_act) { 3773 // Acceptable approx. is found 3774 return x; 3775 } 3776 3777 // Obtain the golden section step 3778 new_step = r * (x < middle_range ? b - x : a - x); 3779 3780 // Decide if the interpolation can be tried. If x and w are distinct, interpolatiom may be tried 3781 if (Math.abs(x - w) >= tol_act) { 3782 // Interpolation step is calculated as p/q; 3783 // division operation is delayed until last moment 3784 t = (x - w) * (fx - fv); 3785 q = (x - v) * (fx - fw); 3786 p = (x - v) * q - (x - w) * t; 3787 q = 2 * (q - t); 3788 3789 if (q > 0) { 3790 p = -p; // q was calculated with the opposite sign; make q positive 3791 } else { 3792 q = -q; // // and assign possible minus to p 3793 } 3794 if ( 3795 Math.abs(p) < Math.abs(new_step * q) && // If x+p/q falls in [a,b] 3796 p > q * (a - x + 2 * tol_act) && // not too close to a and 3797 p < q * (b - x - 2 * tol_act) 3798 ) { 3799 // b, and isn't too large 3800 new_step = p / q; // it is accepted 3801 } 3802 // If p / q is too large then the 3803 // golden section procedure can 3804 // reduce [a,b] range to more 3805 // extent 3806 } 3807 3808 // Adjust the step to be not less than tolerance 3809 if (Math.abs(new_step) < tol_act) { 3810 if (new_step > 0) { 3811 new_step = tol_act; 3812 } else { 3813 new_step = -tol_act; 3814 } 3815 } 3816 3817 // Obtain the next approximation to min 3818 // and reduce the enveloping range 3819 3820 // Tentative point for the min 3821 t = x + new_step; 3822 ft = f.call(context, t); 3823 // nfev += 1; 3824 3825 // t is a better approximation 3826 if (ft <= fx) { 3827 // Reduce the range so that t would fall within it 3828 if (t < x) { 3829 b = x; 3830 } else { 3831 a = x; 3832 } 3833 3834 // Assign the best approx to x 3835 v = w; 3836 w = x; 3837 x = t; 3838 3839 fv = fw; 3840 fw = fx; 3841 fx = ft; 3842 // x remains the better approx 3843 } else { 3844 // Reduce the range enclosing x 3845 if (t < x) { 3846 a = t; 3847 } else { 3848 b = t; 3849 } 3850 3851 if (ft <= fw || w === x) { 3852 v = w; 3853 w = t; 3854 fv = fw; 3855 fw = ft; 3856 } else if (ft <= fv || v === x || v === w) { 3857 v = t; 3858 fv = ft; 3859 } 3860 } 3861 niter += 1; 3862 } 3863 3864 return x; 3865 }, 3866 3867 /** 3868 * GLOMIN seeks a global minimum of a function F(X) in an interval [A,B]. 3869 * <pre> 3870 * 3871 * Discussion: 3872 * 3873 * This function assumes that F(X) is twice continuously differentiable over [A,B] 3874 * and that F''(X) <= M for all X in [A,B]. 3875 * 3876 * Licensing: 3877 * This code is distributed under the GNU LGPL license. 3878 * 3879 * Modified: 3880 * 3881 * 17 April 2008 3882 * 3883 * Author: 3884 * 3885 * Original FORTRAN77 version by Richard Brent. 3886 * C version by John Burkardt. 3887 * https://people.math.sc.edu/Burkardt/c_src/brent/brent.c 3888 * 3889 * Reference: 3890 * 3891 * Richard Brent, 3892 * Algorithms for Minimization Without Derivatives, 3893 * Dover, 2002, 3894 * ISBN: 0-486-41998-3, 3895 * LC: QA402.5.B74. 3896 * 3897 * Parameters: 3898 * 3899 * Input, double A, B, the endpoints of the interval. 3900 * It must be the case that A < B. 3901 * 3902 * Input, double C, an initial guess for the global 3903 * minimizer. If no good guess is known, C = A or B is acceptable. 3904 * 3905 * Input, double M, the bound on the second derivative. 3906 * 3907 * Input, double MACHEP, an estimate for the relative machine 3908 * precision. 3909 * 3910 * Input, double E, a positive tolerance, a bound for the 3911 * absolute error in the evaluation of F(X) for any X in [A,B]. 3912 * 3913 * Input, double T, a positive error tolerance. 3914 * 3915 * Input, double F (double x ), a user-supplied 3916 * function whose global minimum is being sought. 3917 * 3918 * Output, double *X, the estimated value of the abscissa 3919 * for which F attains its global minimum value in [A,B]. 3920 * 3921 * Output, double GLOMIN, the value F(X). 3922 * </pre> 3923 */ 3924 glomin: function (f, x0) { 3925 var a0, a2, a3, d0, d1, d2, h, 3926 k, m2, 3927 p, q, qs, 3928 r, s, sc, 3929 y, y0, y1, y2, y3, yb, 3930 z0, z1, z2, 3931 a, b, c, x, 3932 m = 10000000.0, 3933 t = Mat.eps, // * Mat.eps, 3934 e = Mat.eps * Mat.eps, 3935 machep = Mat.eps * Mat.eps * Mat.eps; 3936 3937 a = x0[0]; 3938 b = x0[1]; 3939 c = (f(a) < f(b)) ? a : b; 3940 3941 a0 = b; 3942 x = a0; 3943 a2 = a; 3944 y0 = f(b); 3945 yb = y0; 3946 y2 = f(a); 3947 y = y2; 3948 3949 if (y0 < y) { 3950 y = y0; 3951 } else { 3952 x = a; 3953 } 3954 3955 if (m <= 0.0 || b <= a) { 3956 return y; 3957 } 3958 3959 m2 = 0.5 * (1.0 + 16.0 * machep) * m; 3960 3961 if (c <= a || b <= c) { 3962 sc = 0.5 * (a + b); 3963 } else { 3964 sc = c; 3965 } 3966 3967 y1 = f(sc); 3968 k = 3; 3969 d0 = a2 - sc; 3970 h = 9.0 / 11.0; 3971 3972 if (y1 < y) { 3973 x = sc; 3974 y = y1; 3975 } 3976 3977 for (; ;) { 3978 d1 = a2 - a0; 3979 d2 = sc - a0; 3980 z2 = b - a2; 3981 z0 = y2 - y1; 3982 z1 = y2 - y0; 3983 r = d1 * d1 * z0 - d0 * d0 * z1; 3984 p = r; 3985 qs = 2.0 * (d0 * z1 - d1 * z0); 3986 q = qs; 3987 3988 if (k < 1000000 || y2 <= y) { 3989 for (; ;) { 3990 if (q * (r * (yb - y2) + z2 * q * ((y2 - y) + t)) < 3991 z2 * m2 * r * (z2 * q - r)) { 3992 3993 a3 = a2 + r / q; 3994 y3 = f(a3); 3995 3996 if (y3 < y) { 3997 x = a3; 3998 y = y3; 3999 } 4000 } 4001 k = ((1611 * k) % 1048576); 4002 q = 1.0; 4003 r = (b - a) * 0.00001 * k; 4004 4005 if (z2 <= r) { 4006 break; 4007 } 4008 } 4009 } else { 4010 k = ((1611 * k) % 1048576); 4011 q = 1.0; 4012 r = (b - a) * 0.00001 * k; 4013 4014 while (r < z2) { 4015 if (q * (r * (yb - y2) + z2 * q * ((y2 - y) + t)) < 4016 z2 * m2 * r * (z2 * q - r)) { 4017 4018 a3 = a2 + r / q; 4019 y3 = f(a3); 4020 4021 if (y3 < y) { 4022 x = a3; 4023 y = y3; 4024 } 4025 } 4026 k = ((1611 * k) % 1048576); 4027 q = 1.0; 4028 r = (b - a) * 0.00001 * k; 4029 } 4030 } 4031 4032 r = m2 * d0 * d1 * d2; 4033 s = Math.sqrt(((y2 - y) + t) / m2); 4034 h = 0.5 * (1.0 + h); 4035 p = h * (p + 2.0 * r * s); 4036 q = q + 0.5 * qs; 4037 r = - 0.5 * (d0 + (z0 + 2.01 * e) / (d0 * m2)); 4038 4039 if (r < s || d0 < 0.0) { 4040 r = a2 + s; 4041 } else { 4042 r = a2 + r; 4043 } 4044 4045 if (0.0 < p * q) { 4046 a3 = a2 + p / q; 4047 } else { 4048 a3 = r; 4049 } 4050 4051 for (; ;) { 4052 a3 = Math.max(a3, r); 4053 4054 if (b <= a3) { 4055 a3 = b; 4056 y3 = yb; 4057 } else { 4058 y3 = f(a3); 4059 } 4060 4061 if (y3 < y) { 4062 x = a3; 4063 y = y3; 4064 } 4065 4066 d0 = a3 - a2; 4067 4068 if (a3 <= r) { 4069 break; 4070 } 4071 4072 p = 2.0 * (y2 - y3) / (m * d0); 4073 4074 if ((1.0 + 9.0 * machep) * d0 <= Math.abs(p)) { 4075 break; 4076 } 4077 4078 if (0.5 * m2 * (d0 * d0 + p * p) <= (y2 - y) + (y3 - y) + 2.0 * t) { 4079 break; 4080 } 4081 a3 = 0.5 * (a2 + a3); 4082 h = 0.9 * h; 4083 } 4084 4085 if (b <= a3) { 4086 break; 4087 } 4088 4089 a0 = sc; 4090 sc = a2; 4091 a2 = a3; 4092 y0 = y1; 4093 y1 = y2; 4094 y2 = y3; 4095 } 4096 4097 return [x, y]; 4098 }, 4099 4100 /** 4101 * Determine all roots of a polynomial with real or complex coefficients by using the 4102 * iterative method attributed to Weierstrass, Durand, Kerner, Aberth, and Ehrlich. In particular, 4103 * the iteration method with cubic convergence is used that is usually attributed to Ehrlich-Aberth. 4104 * <p> 4105 * The returned roots are sorted with respect to their real values. 4106 * <p> This method makes use of the JSXGraph classes {@link JXG.Complex} and {@link JXG.C} to handle 4107 * complex numbers. 4108 * 4109 * @param {Array} a Array of coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2... 4110 * The coefficients are of type Number or JXG.Complex. 4111 * @param {Number} [deg] Optional degree of the polynomial. Otherwise all entries are taken, with 4112 * leading zeros removed. 4113 * @param {Number} [tol=Number.EPSILON] Approximation tolerance 4114 * @param {Number} [max_it=30] Maximum number of iterations 4115 * @param {Array} [initial_values=null] Array of initial values for the roots. If not given, 4116 * starting values are determined by the method of Ozawa. 4117 * @returns {Array} Array of complex numbers (of JXG.Complex) approximating the roots of the polynomial. 4118 * @memberof JXG.Math.Numerics 4119 * @see JXG.Complex 4120 * @see JXG.C 4121 * 4122 * @example 4123 * // Polynomial p(z) = -1 + 1z^2 4124 * var i, roots, 4125 * p = [-1, 0, 1]; 4126 * 4127 * roots = JXG.Math.Numerics.polzeros(p); 4128 * for (i = 0; i < roots.length; i++) { 4129 * console.log(i, roots[i].toString()); 4130 * } 4131 * // Output: 4132 * 0 -1 + -3.308722450212111e-24i 4133 * 1 1 + 0i 4134 * 4135 * @example 4136 * // Polynomial p(z) = -1 + 3z - 9z^2 + z^3 - 8z^6 + 9z^7 - 9z^8 + z^9 4137 * var i, roots, 4138 * p = [-1, 3, -9, 1, 0, 0, -8, 9, -9, 1]; 4139 * 4140 * roots = JXG.Math.Numerics.polzeros(p); 4141 * for (i = 0; i < roots.length; i++) { 4142 * console.log(i, roots[i].toString()); 4143 * } 4144 * // Output: 4145 * 0 -0.7424155888401961 + 0.4950476539211721i 4146 * 1 -0.7424155888401961 + -0.4950476539211721i 4147 * 2 0.16674869833354108 + 0.2980502714610669i 4148 * 3 0.16674869833354108 + -0.29805027146106694i 4149 * 4 0.21429002063640837 + 1.0682775088132996i 4150 * 5 0.21429002063640842 + -1.0682775088132999i 4151 * 6 0.861375497926218 + -0.6259177003583295i 4152 * 7 0.8613754979262181 + 0.6259177003583295i 4153 * 8 8.000002743888055 + -1.8367099231598242e-40i 4154 * 4155 */ 4156 polzeros: function (coeffs, deg, tol, max_it, initial_values) { 4157 var i, le, off, it, 4158 debug = false, 4159 cc = [], 4160 obvious = [], 4161 roots = [], 4162 4163 /** 4164 * Horner method to evaluate polynomial or the derivative thereof for complex numbers, 4165 * i.e. coefficients and variable are complex. 4166 * @function 4167 * @param {Array} a Array of complex coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2... 4168 * @param {JXG.Complex} z Value for which the polynomial will be evaluated. 4169 * @param {Boolean} [derivative=false] If true the derivative will be evaluated. 4170 * @ignore 4171 */ 4172 hornerComplex = function (a, z, derivative) { 4173 var i, s, 4174 n = a.length - 1; 4175 4176 derivative = derivative || false; 4177 if (derivative) { 4178 // s = n * a_n 4179 s = JXG.C.mult(n, a[n]); 4180 for (i = n - 1; i > 0; i--) { 4181 // s = s * z + i * a_i 4182 s.mult(z); 4183 s.add(JXG.C.mult(a[i], i)); 4184 } 4185 } else { 4186 // s = a_n 4187 s = JXG.C.copy(a[n]); 4188 for (i = n - 1; i >= 0; i--) { 4189 // s = s * z + a_i 4190 s.mult(z); 4191 s.add(a[i]); 4192 } 4193 } 4194 return s; 4195 }, 4196 4197 /** 4198 * Horner method to evaluate reciprocal polynomial or the derivative thereof for complex numbers, 4199 * i.e. coefficients and variable are complex. 4200 * @function 4201 * @param {Array} a Array of complex coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2... 4202 * @param {JXG.Complex} z Value for which the reciprocal polynomial will be evaluated. 4203 * @param {Boolean} [derivative=false] If true the derivative will be evaluated. 4204 * @ignore 4205 */ 4206 hornerRec = function (a, x, derivative) { 4207 var i, s, 4208 n = a.length - 1; 4209 4210 derivative = derivative || false; 4211 if (derivative) { 4212 // s = n * a_0 4213 s = JXG.C.mult(n, a[0]); 4214 for (i = n - 1; i > 0; i--) { 4215 // s = s * x + i * a_{n-i} 4216 s.mult(x); 4217 s.add(JXG.C.mult(a[n - i], i)); 4218 } 4219 } else { 4220 // s = a_0 4221 s = JXG.C.copy(a[0]); 4222 for (i = n - 1; i >= 0; i--) { 4223 // s = s * x + a_{n-i} 4224 s.mult(x); 4225 s.add(a[n - i]); 4226 } 4227 } 4228 return s; 4229 }, 4230 4231 /** 4232 * Horner method to evaluate real polynomial at a real value. 4233 * @function 4234 * @param {Array} a Array of real coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2... 4235 * @param {Number} z Value for which the polynomial will be evaluated. 4236 * @ignore 4237 */ 4238 horner = function (a, x) { 4239 var i, s, 4240 n = a.length - 1; 4241 4242 s = a[n]; 4243 for (i = n - 1; i >= 0; i--) { 4244 s = s * x + a[i]; 4245 } 4246 return s; 4247 }, 4248 4249 /** 4250 * Determine start values for the Aberth iteration, see 4251 * Ozawa, "An experimental study of the starting values 4252 * of the Durand-Kerner-Aberth iteration" (1995). 4253 * 4254 * @function 4255 * @param {Array} a Array of complex coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2... 4256 * @returns {Array} Array Initial values for the roots. 4257 * @ignore 4258 */ 4259 initial_guess = function (a) { 4260 var i, r, 4261 n = a.length - 1, // degree 4262 alpha1 = Math.PI * 2 / n, 4263 alpha0 = Math.PI / n * 0.5, 4264 b, z, 4265 init = []; 4266 4267 4268 // From Ozawa, "An experimental study of the starting values 4269 // of the Durand-Kerner-Aberth iteration" (1995) 4270 4271 // b is the arithmetic mean of the roots. 4272 // With is Vieta's formula <https://en.wikipedia.org/wiki/Vieta%27s_formulas> 4273 // b = -a_{n-1} / (n * a_n) 4274 b = JXG.C.mult(-1, a[n - 1]); 4275 b.div(JXG.C.mult(n, a[n])); 4276 4277 // r is the geometric mean of the deviations |b - root_i|. 4278 // Using 4279 // p(z) = a_n prod(z - root_i) 4280 // and therefore 4281 // |p(b)| = |a_n| prod(|b - root_i|) 4282 // we arrive at: 4283 // r = |p(b)/a_n|^(1/n) 4284 z = JXG.C.div(hornerComplex(a, b), a[n]); 4285 r = Math.pow(JXG.C.abs(z), 1 / n); 4286 if (r === 0) { r = 1; } 4287 4288 for (i = 0; i < n; i++) { 4289 a = new JXG.Complex(r * Math.cos(alpha1 * i + alpha0), r * Math.sin(alpha1 * i + alpha0)); 4290 init[i] = JXG.C.add(b, a); 4291 } 4292 4293 return init; 4294 }, 4295 4296 /** 4297 * Ehrlich-Aberth iteration. The stopping criterion is from 4298 * D.A. Bini, "Numerical computation of polynomial zeros 4299 * by means of Aberths's method", Numerical Algorithms (1996). 4300 * 4301 * @function 4302 * @param {Array} a Array of complex coefficients of the polynomial a[0] + a[1]*x+ a[2]*x**2... 4303 * @param {Number} mu Machine precision 4304 * @param {Number} max_it Maximum number of iterations 4305 * @param {Array} z Initial guess for the roots. Will be changed in place. 4306 * @returns {Number} Number of iterations 4307 * @ignore 4308 */ 4309 aberthIteration = function (cc, mu, max_it, z) { 4310 var k, i, j, 4311 done = [], 4312 cr = [], 4313 gamma, x, 4314 done_sum = 0, 4315 num, denom, s, pp, 4316 n = z.length; 4317 4318 for (i = 0; i < n; i++) { 4319 done.push(false); 4320 } 4321 for (i = 0; i < cc.length; i++) { 4322 cr.push(JXG.C.abs(cc[i]) * (4 * i + 1)); 4323 } 4324 for (k = 0; k < max_it && done_sum < n; k++) { 4325 for (i = 0; i < n; i++) { 4326 if (done[i]) { 4327 continue; 4328 } 4329 num = hornerComplex(cc, z[i]); 4330 x = JXG.C.abs(z[i]); 4331 4332 // Stopping criterion by D.A. Bini 4333 // "Numerical computation of polynomial zeros 4334 // by means of Aberths's method", Numerical Algorithms (1996). 4335 // 4336 if (JXG.C.abs(num) < mu * horner(cr, x)) { 4337 done[i] = true; 4338 done_sum++; 4339 if (done_sum === n) { 4340 break; 4341 } 4342 continue; 4343 } 4344 4345 // num = P(z_i) / P'(z_i) 4346 if (x > 1) { 4347 gamma = JXG.C.div(1, z[i]); 4348 pp = hornerRec(cc, gamma, true); 4349 pp.div(hornerRec(cc, gamma)); 4350 pp.mult(gamma); 4351 num = JXG.C.sub(n, pp); 4352 num = JXG.C.div(z[i], num); 4353 } else { 4354 num.div(hornerComplex(cc, z[i], true)); 4355 } 4356 4357 // denom = sum_{i\neq j} 1 / (z_i - z_j) 4358 denom = new JXG.Complex(0); 4359 for (j = 0; j < n; j++) { 4360 if (j === i) { 4361 continue; 4362 } 4363 s = JXG.C.sub(z[i], z[j]); 4364 s = JXG.C.div(1, s); 4365 denom.add(s); 4366 } 4367 4368 // num = num / 1 - num * sum_{i\neq j} 1 / (z_i - z_j) 4369 denom.mult(num); 4370 denom = JXG.C.sub(1, denom); 4371 num.div(denom); 4372 // z_i = z_i - num 4373 z[i].sub(num); 4374 } 4375 } 4376 4377 return k; 4378 }; 4379 4380 4381 tol = tol || Number.EPSILON; 4382 max_it = max_it || 30; 4383 4384 le = coeffs.length; 4385 if (JXG.isNumber(deg) && deg >= 0 && deg < le - 1) { 4386 le = deg + 1; 4387 } 4388 4389 // Convert coefficient array to complex numbers 4390 for (i = 0; i < le; i++) { 4391 cc.push(new JXG.Complex(coeffs[i])); 4392 } 4393 4394 // Search for (multiple) roots at x=0 4395 for (i = 0; i < le; i++) { 4396 if (cc[i].real !== 0 || cc[i].imaginary !== 0) { 4397 off = i; 4398 break; 4399 } 4400 } 4401 4402 // Deflate root x=0, store roots at x=0 in obvious 4403 for (i = 0; i < off; i++) { 4404 obvious.push(new JXG.Complex(0)); 4405 } 4406 cc = cc.slice(off); 4407 le = cc.length; 4408 4409 // Remove leading zeros from the coefficient array 4410 for (i = le - 1; i >= 0; i--) { 4411 if (cc[i].real !== 0 || cc[i].imaginary !== 0) { 4412 break; 4413 } 4414 cc.pop(); 4415 } 4416 le = cc.length; 4417 if (le === 0) { 4418 return []; 4419 } 4420 4421 // From now on we can assume that the 4422 // constant coefficient and the leading coefficient 4423 // are not zero. 4424 if (initial_values) { 4425 for (i = 0; i < le - 1; i++) { 4426 roots.push(new JXG.Complex(initial_values[i])); 4427 } 4428 } else { 4429 roots = initial_guess(cc); 4430 } 4431 it = aberthIteration(cc, tol, max_it, roots); 4432 4433 // Append the roots at x=0 4434 roots = obvious.concat(roots); 4435 4436 if (debug) { 4437 console.log("Iterations:", it); 4438 console.log('Roots:'); 4439 for (i = 0; i < roots.length; i++) { 4440 console.log(i, roots[i].toString(), JXG.C.abs(hornerComplex(cc, roots[i]))); 4441 } 4442 } 4443 4444 // Sort roots according to their real part 4445 roots.sort(function (a, b) { 4446 if (a.real < b.real) { 4447 return -1; 4448 } 4449 if (a.real > b.real) { 4450 return 1; 4451 } 4452 return 0; 4453 }); 4454 4455 return roots; 4456 }, 4457 4458 /** 4459 * Implements the Ramer-Douglas-Peucker algorithm. 4460 * It discards points which are not necessary from the polygonal line defined by the point array 4461 * pts. The computation is done in screen coordinates. 4462 * Average runtime is O(nlog(n)), worst case runtime is O(n^2), where n is the number of points. 4463 * @param {Array} pts Array of {@link JXG.Coords} 4464 * @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>. 4465 * @returns {Array} An array containing points which represent an apparently identical curve as the points of pts do, but contains fewer points. 4466 * @memberof JXG.Math.Numerics 4467 */ 4468 RamerDouglasPeucker: function (pts, eps) { 4469 var allPts = [], 4470 newPts = [], 4471 i, k, len, 4472 endless = true, 4473 4474 /** 4475 * findSplit() is a subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}. 4476 * It searches for the point between index i and j which 4477 * has the largest distance from the line between the points i and j. 4478 * @param {Array} pts Array of {@link JXG.Coords} 4479 * @param {Number} i Index of a point in pts 4480 * @param {Number} j Index of a point in pts 4481 * @ignore 4482 * @private 4483 */ 4484 findSplit = function (pts, i, j) { 4485 var d, k, ci, cj, ck, 4486 x0, y0, x1, y1, 4487 den, lbda, 4488 eps = Mat.eps * Mat.eps, 4489 huge = 10000, 4490 dist = 0, 4491 f = i; 4492 4493 if (j - i < 2) { 4494 return [-1.0, 0]; 4495 } 4496 4497 ci = pts[i].scrCoords; 4498 cj = pts[j].scrCoords; 4499 4500 if (isNaN(ci[1]) || isNaN(ci[2])) { 4501 return [NaN, i]; 4502 } 4503 if (isNaN(cj[1]) || isNaN(cj[2])) { 4504 return [NaN, j]; 4505 } 4506 4507 for (k = i + 1; k < j; k++) { 4508 ck = pts[k].scrCoords; 4509 if (isNaN(ck[1]) || isNaN(ck[2])) { 4510 return [NaN, k]; 4511 } 4512 4513 x0 = ck[1] - ci[1]; 4514 y0 = ck[2] - ci[2]; 4515 x1 = cj[1] - ci[1]; 4516 y1 = cj[2] - ci[2]; 4517 x0 = x0 === Infinity ? huge : x0; 4518 y0 = y0 === Infinity ? huge : y0; 4519 x1 = x1 === Infinity ? huge : x1; 4520 y1 = y1 === Infinity ? huge : y1; 4521 x0 = x0 === -Infinity ? -huge : x0; 4522 y0 = y0 === -Infinity ? -huge : y0; 4523 x1 = x1 === -Infinity ? -huge : x1; 4524 y1 = y1 === -Infinity ? -huge : y1; 4525 den = x1 * x1 + y1 * y1; 4526 4527 if (den > eps) { 4528 lbda = (x0 * x1 + y0 * y1) / den; 4529 4530 if (lbda < 0.0) { 4531 lbda = 0.0; 4532 } else if (lbda > 1.0) { 4533 lbda = 1.0; 4534 } 4535 4536 x0 = x0 - lbda * x1; 4537 y0 = y0 - lbda * y1; 4538 d = x0 * x0 + y0 * y0; 4539 } else { 4540 lbda = 0.0; 4541 d = x0 * x0 + y0 * y0; 4542 } 4543 4544 if (d > dist) { 4545 dist = d; 4546 f = k; 4547 } 4548 } 4549 return [Math.sqrt(dist), f]; 4550 }, 4551 /** 4552 * RDP() is a private subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}. 4553 * It runs recursively through the point set and searches the 4554 * point which has the largest distance from the line between the first point and 4555 * the last point. If the distance from the line is greater than eps, this point is 4556 * included in our new point set otherwise it is discarded. 4557 * If it is taken, we recursively apply the subroutine to the point set before 4558 * and after the chosen point. 4559 * @param {Array} pts Array of {@link JXG.Coords} 4560 * @param {Number} i Index of an element of pts 4561 * @param {Number} j Index of an element of pts 4562 * @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>. 4563 * @param {Array} newPts Array of {@link JXG.Coords} 4564 * @ignore 4565 * @private 4566 */ 4567 RDP = function (pts, i, j, eps, newPts) { 4568 var result = findSplit(pts, i, j), 4569 k = result[1]; 4570 4571 if (isNaN(result[0])) { 4572 RDP(pts, i, k - 1, eps, newPts); 4573 newPts.push(pts[k]); 4574 do { 4575 ++k; 4576 } while (k <= j && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])); 4577 if (k <= j) { 4578 newPts.push(pts[k]); 4579 } 4580 RDP(pts, k + 1, j, eps, newPts); 4581 } else if (result[0] > eps) { 4582 RDP(pts, i, k, eps, newPts); 4583 RDP(pts, k, j, eps, newPts); 4584 } else { 4585 newPts.push(pts[j]); 4586 } 4587 }; 4588 4589 len = pts.length; 4590 4591 i = 0; 4592 while (endless) { 4593 // Search for the next point without NaN coordinates 4594 while (i < len && isNaN(pts[i].scrCoords[1] + pts[i].scrCoords[2])) { 4595 i += 1; 4596 } 4597 // Search for the next position of a NaN point 4598 k = i + 1; 4599 while (k < len && !isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])) { 4600 k += 1; 4601 } 4602 k--; 4603 4604 // Only proceed if something is left 4605 if (i < len && k > i) { 4606 newPts = []; 4607 newPts[0] = pts[i]; 4608 RDP(pts, i, k, eps, newPts); 4609 allPts = allPts.concat(newPts); 4610 } 4611 if (i >= len) { 4612 break; 4613 } 4614 // Push the NaN point 4615 if (k < len - 1) { 4616 allPts.push(pts[k + 1]); 4617 } 4618 i = k + 1; 4619 } 4620 4621 return allPts; 4622 }, 4623 4624 /** 4625 * Old name for the implementation of the Ramer-Douglas-Peucker algorithm. 4626 * @deprecated Use {@link JXG.Math.Numerics.RamerDouglasPeucker} 4627 * @memberof JXG.Math.Numerics 4628 */ 4629 RamerDouglasPeuker: function (pts, eps) { 4630 JXG.deprecated("Numerics.RamerDouglasPeuker()", "Numerics.RamerDouglasPeucker()"); 4631 return this.RamerDouglasPeucker(pts, eps); 4632 }, 4633 4634 /** 4635 * Implements the Visvalingam-Whyatt algorithm. 4636 * See M. Visvalingam, J. D. Whyatt: 4637 * "Line generalisation by repeated elimination of the smallest area", C.I.S.R.G Discussion paper 10, July 1992 4638 * 4639 * The algorithm discards points which are not necessary from the polygonal line defined by the point array 4640 * pts (consisting of type JXG.Coords). 4641 * @param {Array} pts Array of {@link JXG.Coords} 4642 * @param {Number} numPoints Number of remaining intermediate points. The first and the last point of the original points will 4643 * be taken in any case. 4644 * @returns {Array} An array containing points which approximates the curve defined by pts. 4645 * @memberof JXG.Math.Numerics 4646 * 4647 * @example 4648 * var i, p = []; 4649 * for (i = 0; i < 5; ++i) { 4650 * p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6])); 4651 * } 4652 * 4653 * // Plot a cardinal spline curve 4654 * var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5); 4655 * var cu1 = board.create('curve', splineArr, {strokeColor: 'green'}); 4656 * 4657 * var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'}); 4658 * c.updateDataArray = function() { 4659 * var i, len, points; 4660 * 4661 * // Reduce number of intermediate points with Visvakingam-Whyatt to 6 4662 * points = JXG.Math.Numerics.Visvalingam(cu1.points, 6); 4663 * // Plot the remaining points 4664 * len = points.length; 4665 * this.dataX = []; 4666 * this.dataY = []; 4667 * for (i = 0; i < len; i++) { 4668 * this.dataX.push(points[i].usrCoords[1]); 4669 * this.dataY.push(points[i].usrCoords[2]); 4670 * } 4671 * }; 4672 * board.update(); 4673 * 4674 * </pre><div id="JXGce0cc55c-b592-11e6-8270-104a7d3be7eb" class="jxgbox" style="width: 300px; height: 300px;"></div> 4675 * <script type="text/javascript"> 4676 * (function() { 4677 * var board = JXG.JSXGraph.initBoard('JXGce0cc55c-b592-11e6-8270-104a7d3be7eb', 4678 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 4679 * 4680 * var i, p = []; 4681 * for (i = 0; i < 5; ++i) { 4682 * p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6])); 4683 * } 4684 * 4685 * // Plot a cardinal spline curve 4686 * var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5); 4687 * var cu1 = board.create('curve', splineArr, {strokeColor: 'green'}); 4688 * 4689 * var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'}); 4690 * c.updateDataArray = function() { 4691 * var i, len, points; 4692 * 4693 * // Reduce number of intermediate points with Visvakingam-Whyatt to 6 4694 * points = JXG.Math.Numerics.Visvalingam(cu1.points, 6); 4695 * // Plot the remaining points 4696 * len = points.length; 4697 * this.dataX = []; 4698 * this.dataY = []; 4699 * for (i = 0; i < len; i++) { 4700 * this.dataX.push(points[i].usrCoords[1]); 4701 * this.dataY.push(points[i].usrCoords[2]); 4702 * } 4703 * }; 4704 * board.update(); 4705 * 4706 * })(); 4707 * 4708 * </script><pre> 4709 * 4710 */ 4711 Visvalingam: function (pts, numPoints) { 4712 var i, 4713 len, 4714 vol, 4715 lastVol, 4716 linkedList = [], 4717 heap = [], 4718 points = [], 4719 lft, 4720 rt, 4721 lft2, 4722 rt2, 4723 obj; 4724 4725 len = pts.length; 4726 4727 // At least one intermediate point is needed 4728 if (len <= 2) { 4729 return pts; 4730 } 4731 4732 // Fill the linked list 4733 // Add first point to the linked list 4734 linkedList[0] = { 4735 used: true, 4736 lft: null, 4737 node: null 4738 }; 4739 4740 // Add all intermediate points to the linked list, 4741 // whose triangle area is nonzero. 4742 lft = 0; 4743 for (i = 1; i < len - 1; i++) { 4744 vol = Math.abs( 4745 JXG.Math.Numerics.det([ 4746 pts[i - 1].usrCoords, 4747 pts[i].usrCoords, 4748 pts[i + 1].usrCoords 4749 ]) 4750 ); 4751 if (!isNaN(vol)) { 4752 obj = { 4753 v: vol, 4754 idx: i 4755 }; 4756 heap.push(obj); 4757 linkedList[i] = { 4758 used: true, 4759 lft: lft, 4760 node: obj 4761 }; 4762 linkedList[lft].rt = i; 4763 lft = i; 4764 } 4765 } 4766 4767 // Add last point to the linked list 4768 linkedList[len - 1] = { 4769 used: true, 4770 rt: null, 4771 lft: lft, 4772 node: null 4773 }; 4774 linkedList[lft].rt = len - 1; 4775 4776 // Remove points until only numPoints intermediate points remain 4777 lastVol = -Infinity; 4778 while (heap.length > numPoints) { 4779 // Sort the heap with the updated volume values 4780 heap.sort(function (a, b) { 4781 // descending sort 4782 return b.v - a.v; 4783 }); 4784 4785 // Remove the point with the smallest triangle 4786 i = heap.pop().idx; 4787 linkedList[i].used = false; 4788 lastVol = linkedList[i].node.v; 4789 4790 // Update the pointers of the linked list 4791 lft = linkedList[i].lft; 4792 rt = linkedList[i].rt; 4793 linkedList[lft].rt = rt; 4794 linkedList[rt].lft = lft; 4795 4796 // Update the values for the volumes in the linked list 4797 lft2 = linkedList[lft].lft; 4798 if (lft2 !== null) { 4799 vol = Math.abs( 4800 JXG.Math.Numerics.det([ 4801 pts[lft2].usrCoords, 4802 pts[lft].usrCoords, 4803 pts[rt].usrCoords 4804 ]) 4805 ); 4806 4807 linkedList[lft].node.v = vol >= lastVol ? vol : lastVol; 4808 } 4809 rt2 = linkedList[rt].rt; 4810 if (rt2 !== null) { 4811 vol = Math.abs( 4812 JXG.Math.Numerics.det([ 4813 pts[lft].usrCoords, 4814 pts[rt].usrCoords, 4815 pts[rt2].usrCoords 4816 ]) 4817 ); 4818 4819 linkedList[rt].node.v = vol >= lastVol ? vol : lastVol; 4820 } 4821 } 4822 4823 // Return an array with the remaining points 4824 i = 0; 4825 points = [pts[i]]; 4826 do { 4827 i = linkedList[i].rt; 4828 points.push(pts[i]); 4829 } while (linkedList[i].rt !== null); 4830 4831 return points; 4832 } 4833 }; 4834 4835 export default Mat.Numerics; 4836