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