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