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