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