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