1 /* 2 Copyright 2008-2017 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', 'math/math'], function (JXG, Type, 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, 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, abserr, 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, resObj, 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. 1424 * @param {function} f We search for a solution of f(x)=0. 1425 * @param {Number} x initial guess for the root, i.e. starting value. 1426 * @param {Object} context optional object that is treated as "this" in the function body. This is useful if 1427 * the function is a method of an object and contains a reference to its parent object via "this". 1428 * @returns {Number} A root of the function f. 1429 * @memberof JXG.Math.Numerics 1430 */ 1431 root: function (f, x, context) { 1432 return this.fzero(f, x, context); 1433 }, 1434 1435 /** 1436 * Compute an intersection of the curves c1 and c2 1437 * with a generalized Newton method. 1438 * We want to find values t1, t2 such that 1439 * c1(t1) = c2(t2), i.e. 1440 * (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2)) = (0,0). 1441 * We set 1442 * (e,f) := (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2)) 1443 * 1444 * The Jacobian J is defined by 1445 * J = (a, b) 1446 * (c, d) 1447 * where 1448 * a = c1_x'(t1) 1449 * b = -c2_x'(t2) 1450 * c = c1_y'(t1) 1451 * d = -c2_y'(t2) 1452 * 1453 * The inverse J^(-1) of J is equal to 1454 * (d, -b)/ 1455 * (-c, a) / (ad-bc) 1456 * 1457 * Then, (t1new, t2new) := (t1,t2) - J^(-1)*(e,f). 1458 * If the function meetCurveCurve possesses the properties 1459 * t1memo and t2memo then these are taken as start values 1460 * for the Newton algorithm. 1461 * After stopping of the Newton algorithm the values of t1 and t2 are stored in 1462 * t1memo and t2memo. 1463 * 1464 * @param {JXG.Curve} c1 Curve, Line or Circle 1465 * @param {JXG.Curve} c2 Curve, Line or Circle 1466 * @param {Number} t1ini start value for t1 1467 * @param {Number} t2ini start value for t2 1468 * @returns {JXG.Coords} intersection point 1469 * @memberof JXG.Math.Numerics 1470 */ 1471 generalizedNewton: function (c1, c2, t1ini, t2ini) { 1472 var t1, t2, 1473 a, b, c, d, disc, 1474 e, f, F, 1475 D00, D01, 1476 D10, D11, 1477 count = 0; 1478 1479 if (this.generalizedNewton.t1memo) { 1480 t1 = this.generalizedNewton.t1memo; 1481 t2 = this.generalizedNewton.t2memo; 1482 } else { 1483 t1 = t1ini; 1484 t2 = t2ini; 1485 } 1486 1487 e = c1.X(t1) - c2.X(t2); 1488 f = c1.Y(t1) - c2.Y(t2); 1489 F = e * e + f * f; 1490 1491 D00 = this.D(c1.X, c1); 1492 D01 = this.D(c2.X, c2); 1493 D10 = this.D(c1.Y, c1); 1494 D11 = this.D(c2.Y, c2); 1495 1496 while (F > Mat.eps && count < 10) { 1497 a = D00(t1); 1498 b = -D01(t2); 1499 c = D10(t1); 1500 d = -D11(t2); 1501 disc = a * d - b * c; 1502 t1 -= (d * e - b * f) / disc; 1503 t2 -= (a * f - c * e) / disc; 1504 e = c1.X(t1) - c2.X(t2); 1505 f = c1.Y(t1) - c2.Y(t2); 1506 F = e * e + f * f; 1507 count += 1; 1508 } 1509 1510 this.generalizedNewton.t1memo = t1; 1511 this.generalizedNewton.t2memo = t2; 1512 1513 if (Math.abs(t1) < Math.abs(t2)) { 1514 return [c1.X(t1), c1.Y(t1)]; 1515 } 1516 1517 return [c2.X(t2), c2.Y(t2)]; 1518 }, 1519 1520 /** 1521 * Returns the Lagrange polynomials for curves with equidistant nodes, see 1522 * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation, 1523 * SIAM Review, Vol 46, No 3, (2004) 501-517. 1524 * The graph of the parametric curve [x(t),y(t)] runs through the given points. 1525 * @param {Array} p Array of JXG.Points 1526 * @returns {Array} An array consisting of two functions x(t), y(t) which define a parametric curve 1527 * f(t) = (x(t), y(t)) and two numbers x1 and x2 defining the curve's domain. x1 always equals zero. 1528 * @memberof JXG.Math.Numerics 1529 */ 1530 Neville: function (p) { 1531 var w = [], 1532 /** @ignore */ 1533 makeFct = function (fun) { 1534 return function (t, suspendedUpdate) { 1535 var i, d, s, 1536 bin = Mat.binomial, 1537 len = p.length, 1538 len1 = len - 1, 1539 num = 0.0, 1540 denom = 0.0; 1541 1542 if (!suspendedUpdate) { 1543 s = 1; 1544 for (i = 0; i < len; i++) { 1545 w[i] = bin(len1, i) * s; 1546 s *= (-1); 1547 } 1548 } 1549 1550 d = t; 1551 1552 for (i = 0; i < len; i++) { 1553 if (d === 0) { 1554 return p[i][fun](); 1555 } 1556 s = w[i] / d; 1557 d -= 1; 1558 num += p[i][fun]() * s; 1559 denom += s; 1560 } 1561 return num / denom; 1562 }; 1563 }, 1564 1565 xfct = makeFct('X'), 1566 yfct = makeFct('Y'); 1567 1568 return [xfct, yfct, 0, function () { 1569 return p.length - 1; 1570 }]; 1571 }, 1572 1573 /** 1574 * Calculates second derivatives at the given knots. 1575 * @param {Array} x x values of knots 1576 * @param {Array} y y values of knots 1577 * @returns {Array} Second derivatives of the interpolated function at the knots. 1578 * @see #splineEval 1579 * @memberof JXG.Math.Numerics 1580 */ 1581 splineDef: function (x, y) { 1582 var pair, i, l, 1583 n = Math.min(x.length, y.length), 1584 diag = [], 1585 z = [], 1586 data = [], 1587 dx = [], 1588 delta = [], 1589 F = []; 1590 1591 if (n === 2) { 1592 return [0, 0]; 1593 } 1594 1595 for (i = 0; i < n; i++) { 1596 pair = {X: x[i], Y: y[i]}; 1597 data.push(pair); 1598 } 1599 data.sort(function (a, b) { 1600 return a.X - b.X; 1601 }); 1602 for (i = 0; i < n; i++) { 1603 x[i] = data[i].X; 1604 y[i] = data[i].Y; 1605 } 1606 1607 for (i = 0; i < n - 1; i++) { 1608 dx.push(x[i + 1] - x[i]); 1609 } 1610 for (i = 0; i < n - 2; i++) { 1611 delta.push(6 * (y[i + 2] - y[i + 1]) / (dx[i + 1]) - 6 * (y[i + 1] - y[i]) / (dx[i])); 1612 } 1613 1614 // ForwardSolve 1615 diag.push(2 * (dx[0] + dx[1])); 1616 z.push(delta[0]); 1617 1618 for (i = 0; i < n - 3; i++) { 1619 l = dx[i + 1] / diag[i]; 1620 diag.push(2 * (dx[i + 1] + dx[i + 2]) - l * dx[i + 1]); 1621 z.push(delta[i + 1] - l * z[i]); 1622 } 1623 1624 // BackwardSolve 1625 F[n - 3] = z[n - 3] / diag[n - 3]; 1626 for (i = n - 4; i >= 0; i--) { 1627 F[i] = (z[i] - (dx[i + 1] * F[i + 1])) / diag[i]; 1628 } 1629 1630 // Generate f''-Vector 1631 for (i = n - 3; i >= 0; i--) { 1632 F[i + 1] = F[i]; 1633 } 1634 1635 // natural cubic spline 1636 F[0] = 0; 1637 F[n - 1] = 0; 1638 1639 return F; 1640 }, 1641 1642 /** 1643 * Evaluate points on spline. 1644 * @param {Number,Array} x0 A single float value or an array of values to evaluate 1645 * @param {Array} x x values of knots 1646 * @param {Array} y y values of knots 1647 * @param {Array} F Second derivatives at knots, calculated by {@link JXG.Math.Numerics.splineDef} 1648 * @see #splineDef 1649 * @returns {Number,Array} A single value or an array, depending on what is given as x0. 1650 * @memberof JXG.Math.Numerics 1651 */ 1652 splineEval: function (x0, x, y, F) { 1653 var i, j, a, b, c, d, x_, 1654 n = Math.min(x.length, y.length), 1655 l = 1, 1656 asArray = false, 1657 y0 = []; 1658 1659 // number of points to be evaluated 1660 if (Type.isArray(x0)) { 1661 l = x0.length; 1662 asArray = true; 1663 } else { 1664 x0 = [x0]; 1665 } 1666 1667 for (i = 0; i < l; i++) { 1668 // is x0 in defining interval? 1669 if ((x0[i] < x[0]) || (x[i] > x[n - 1])) { 1670 return NaN; 1671 } 1672 1673 // determine part of spline in which x0 lies 1674 for (j = 1; j < n; j++) { 1675 if (x0[i] <= x[j]) { 1676 break; 1677 } 1678 } 1679 1680 j -= 1; 1681 1682 // we're now in the j-th partial interval, i.e. x[j] < x0[i] <= x[j+1]; 1683 // determine the coefficients of the polynomial in this interval 1684 a = y[j]; 1685 b = (y[j + 1] - y[j]) / (x[j + 1] - x[j]) - (x[j + 1] - x[j]) / 6 * (F[j + 1] + 2 * F[j]); 1686 c = F[j] / 2; 1687 d = (F[j + 1] - F[j]) / (6 * (x[j + 1] - x[j])); 1688 // evaluate x0[i] 1689 x_ = x0[i] - x[j]; 1690 //y0.push(a + b*x_ + c*x_*x_ + d*x_*x_*x_); 1691 y0.push(a + (b + (c + d * x_) * x_) * x_); 1692 } 1693 1694 if (asArray) { 1695 return y0; 1696 } 1697 1698 return y0[0]; 1699 }, 1700 1701 /** 1702 * Generate a string containing the function term of a polynomial. 1703 * @param {Array} coeffs Coefficients of the polynomial. The position i belongs to x^i. 1704 * @param {Number} deg Degree of the polynomial 1705 * @param {String} varname Name of the variable (usually 'x') 1706 * @param {Number} prec Precision 1707 * @returns {String} A string containg the function term of the polynomial. 1708 * @memberof JXG.Math.Numerics 1709 */ 1710 generatePolynomialTerm: function (coeffs, deg, varname, prec) { 1711 var i, t = []; 1712 1713 for (i = deg; i >= 0; i--) { 1714 t = t.concat(['(', coeffs[i].toPrecision(prec), ')']); 1715 1716 if (i > 1) { 1717 t = t.concat(['*', varname, '<sup>', i, '<', '/sup> + ']); 1718 } else if (i === 1) { 1719 t = t.concat(['*', varname, ' + ']); 1720 } 1721 } 1722 1723 return t.join(''); 1724 }, 1725 1726 /** 1727 * Computes the polynomial through a given set of coordinates in Lagrange form. 1728 * Returns the Lagrange polynomials, see 1729 * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation, 1730 * SIAM Review, Vol 46, No 3, (2004) 501-517. 1731 * @param {Array} p Array of JXG.Points 1732 * @returns {function} A function of one parameter which returns the value of the polynomial, whose graph runs through the given points. 1733 * @memberof JXG.Math.Numerics 1734 */ 1735 lagrangePolynomial: function (p) { 1736 var w = [], 1737 /** @ignore */ 1738 fct = function (x, suspendedUpdate) { 1739 var i, j, k, xi, s, M, 1740 len = p.length, 1741 num = 0, 1742 denom = 0; 1743 1744 if (!suspendedUpdate) { 1745 for (i = 0; i < len; i++) { 1746 w[i] = 1.0; 1747 xi = p[i].X(); 1748 1749 for (k = 0; k < len; k++) { 1750 if (k !== i) { 1751 w[i] *= (xi - p[k].X()); 1752 } 1753 } 1754 1755 w[i] = 1 / w[i]; 1756 } 1757 1758 M = []; 1759 1760 for (j = 0; j < len; j++) { 1761 M.push([1]); 1762 } 1763 } 1764 1765 for (i = 0; i < len; i++) { 1766 xi = p[i].X(); 1767 1768 if (x === xi) { 1769 return p[i].Y(); 1770 } 1771 1772 s = w[i] / (x - xi); 1773 denom += s; 1774 num += s * p[i].Y(); 1775 } 1776 1777 return num / denom; 1778 }; 1779 1780 fct.getTerm = function () { 1781 return ''; 1782 }; 1783 1784 return fct; 1785 }, 1786 1787 /** 1788 * Determine the coefficients of a cardinal spline polynom, See 1789 * http://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections 1790 * @param {Number} x1 point 1 1791 * @param {Number} x2 point 2 1792 * @param {Number} t1 tangent slope 1 1793 * @param {Number} t2 tangent slope 2 1794 * @return {Array} coefficents array c for the polynomial t maps to 1795 * c[0] + c[1]*t + c[2]*t*t + c[3]*t*t*t 1796 */ 1797 _initCubicPoly: function(x1, x2, t1, t2) { 1798 return [ 1799 x1, 1800 t1, 1801 -3 * x1 + 3 * x2 - 2 * t1 - t2, 1802 2 * x1 - 2 * x2 + t1 + t2 1803 ]; 1804 }, 1805 1806 /** 1807 * Computes the cubic cardinal spline curve through a given set of points. The curve 1808 * is uniformly parametrized. 1809 * Two artificial control points at the beginning and the end are added. 1810 * 1811 * The implementation (especially the centripetal parametrization) is from 1812 * http://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections . 1813 * @param {Array} points Array consisting of JXG.Points. 1814 * @param {Number|Function} tau The tension parameter, either a constant number or a function returning a number. This number is between 0 and 1. 1815 * tau=1/2 give Catmull-Rom splines. 1816 * @param {String} type (Optional) parameter which allows to choose between "uniform" (default) and 1817 * "centripetal" parameterization. Thus the two possible values are "uniform" or "centripetal". 1818 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 1819 * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, 1820 * and a function simply returning the length of the points array 1821 * minus three. 1822 * @memberof JXG.Math.Numerics 1823 */ 1824 CardinalSpline: function (points, tau_param, type) { 1825 var p, 1826 coeffs = [], 1827 makeFct, 1828 tau, _tau, 1829 that = this; 1830 1831 if (Type.isFunction(tau_param)) { 1832 _tau = tau_param; 1833 } else { 1834 _tau = function () { return tau_param; }; 1835 } 1836 1837 if (type === undefined) { 1838 type = 'uniform'; 1839 } 1840 1841 /** @ignore */ 1842 makeFct = function (which) { 1843 return function (t, suspendedUpdate) { 1844 var s, c, 1845 // control point at the beginning and at the end 1846 first, last, 1847 t1, t2, dt0, dt1, dt2, 1848 dx, dy, 1849 len; 1850 1851 if (points.length < 2) { 1852 return NaN; 1853 } 1854 1855 if (!suspendedUpdate) { 1856 tau = _tau(); 1857 1858 // New point list p: [first, points ..., last] 1859 first = { 1860 X: function () { return 2 * points[0].X() - points[1].X(); }, 1861 Y: function () { return 2 * points[0].Y() - points[1].Y(); }, 1862 Dist: function(p) { 1863 var dx = this.X() - p.X(), 1864 dy = this.Y() - p.Y(); 1865 return Math.sqrt(dx * dx + dy * dy); 1866 } 1867 }; 1868 1869 last = { 1870 X: function () { return 2 * points[points.length - 1].X() - points[points.length - 2].X(); }, 1871 Y: function () { return 2 * points[points.length - 1].Y() - points[points.length - 2].Y(); }, 1872 Dist: function(p) { 1873 var dx = this.X() - p.X(), 1874 dy = this.Y() - p.Y(); 1875 return Math.sqrt(dx * dx + dy * dy); 1876 } 1877 }; 1878 1879 p = [first].concat(points, [last]); 1880 len = p.length; 1881 1882 coeffs[which] = []; 1883 1884 for (s = 0; s < len - 3; s++) { 1885 if (type === 'centripetal') { 1886 // The order is importeant, since p[0].coords === undefined 1887 dt0 = p[s].Dist(p[s + 1]); 1888 dt1 = p[s + 2].Dist(p[s + 1]); 1889 dt2 = p[s + 3].Dist(p[s + 2]); 1890 1891 dt0 = Math.sqrt(dt0); 1892 dt1 = Math.sqrt(dt1); 1893 dt2 = Math.sqrt(dt2); 1894 1895 if (dt1 < Mat.eps) dt1 = 1.0; 1896 if (dt0 < Mat.eps) dt0 = dt1; 1897 if (dt2 < Mat.eps) dt2 = dt1; 1898 1899 t1 = (p[s + 1][which]() - p[s][which]()) / dt0 - 1900 (p[s + 2][which]() - p[s][which]()) / (dt1 + dt0) + 1901 (p[s + 2][which]() - p[s + 1][which]()) / dt1; 1902 1903 t2 = (p[s + 2][which]() - p[s + 1][which]()) / dt1 - 1904 (p[s + 3][which]() - p[s + 1][which]()) / (dt2 + dt1) + 1905 (p[s + 3][which]() - p[s + 2][which]()) / dt2; 1906 1907 t1 *= dt1; 1908 t2 *= dt1; 1909 1910 coeffs[which][s] = that._initCubicPoly( 1911 p[s + 1][which](), 1912 p[s + 2][which](), 1913 tau * t1, 1914 tau * t2 1915 ); 1916 } else { 1917 coeffs[which][s] = that._initCubicPoly( 1918 p[s + 1][which](), 1919 p[s + 2][which](), 1920 tau * (p[s + 2][which]() - p[s][which]()), 1921 tau * (p[s + 3][which]() - p[s + 1][which]()) 1922 ); 1923 } 1924 } 1925 } 1926 1927 if (isNaN(t)) { 1928 return NaN; 1929 } 1930 1931 len = points.length; 1932 // This is necessary for our advanced plotting algorithm: 1933 if (t <= 0.0) { 1934 return points[0][which](); 1935 } else if (t >= len) { 1936 return points[len - 1][which](); 1937 } 1938 1939 s = Math.floor(t); 1940 if (s === t) { 1941 return points[s][which](); 1942 } 1943 1944 t -= s; 1945 c = coeffs[which][s]; 1946 1947 return (((c[3] * t + c[2]) * t + c[1]) * t + c[0]); 1948 }; 1949 }; 1950 1951 return [makeFct('X'), makeFct('Y'), 0, 1952 function () { 1953 return points.length - 1; 1954 }]; 1955 }, 1956 1957 /** 1958 * Computes the cubic Catmull-Rom spline curve through a given set of points. The curve 1959 * is uniformly parametrized. The curve is the cardinal spline curve for tau=0.5. 1960 * Two artificial control points at the beginning and the end are added. 1961 * @param {Array} points Array consisting of JXG.Points. 1962 * @param {String} type (Optional) parameter which allows to choose between "uniform" (default) and 1963 * "centripetal" parameterization. Thus the two possible values are "uniform" or "centripetal". 1964 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 1965 * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, and a function simply 1966 * returning the length of the points array minus three. 1967 * @memberof JXG.Math.Numerics 1968 */ 1969 CatmullRomSpline: function (points, type) { 1970 return this.CardinalSpline(points, 0.5, type); 1971 }, 1972 1973 /** 1974 * Computes the regression polynomial of a given degree through a given set of coordinates. 1975 * Returns the regression polynomial function. 1976 * @param {Number,function,Slider} degree number, function or slider. 1977 * Either 1978 * @param {Array} dataX Array containing either the x-coordinates of the data set or both coordinates in 1979 * an array of {@link JXG.Point}s or {@link JXG.Coords}. 1980 * In the latter case, the <tt>dataY</tt> parameter will be ignored. 1981 * @param {Array} dataY Array containing the y-coordinates of the data set, 1982 * @returns {function} A function of one parameter which returns the value of the regression polynomial of the given degree. 1983 * It possesses the method getTerm() which returns the string containing the function term of the polynomial. 1984 * @memberof JXG.Math.Numerics 1985 */ 1986 regressionPolynomial: function (degree, dataX, dataY) { 1987 var coeffs, deg, dX, dY, inputType, fct, 1988 term = ''; 1989 1990 // Slider 1991 if (Type.isPoint(degree) && Type.isFunction(degree.Value)) { 1992 /** @ignore */ 1993 deg = function () { 1994 return degree.Value(); 1995 }; 1996 // function 1997 } else if (Type.isFunction(degree)) { 1998 deg = degree; 1999 // number 2000 } else if (Type.isNumber(degree)) { 2001 /** @ignore */ 2002 deg = function () { 2003 return degree; 2004 }; 2005 } else { 2006 throw new Error("JSXGraph: Can't create regressionPolynomial from degree of type'" + (typeof degree) + "'."); 2007 } 2008 2009 // Parameters degree, dataX, dataY 2010 if (arguments.length === 3 && Type.isArray(dataX) && Type.isArray(dataY)) { 2011 inputType = 0; 2012 // Parameters degree, point array 2013 } else if (arguments.length === 2 && Type.isArray(dataX) && dataX.length > 0 && Type.isPoint(dataX[0])) { 2014 inputType = 1; 2015 } else if (arguments.length === 2 && Type.isArray(dataX) && dataX.length > 0 && dataX[0].usrCoords && dataX[0].scrCoords) { 2016 inputType = 2; 2017 } else { 2018 throw new Error("JSXGraph: Can't create regressionPolynomial. Wrong parameters."); 2019 } 2020 2021 /** @ignore */ 2022 fct = function (x, suspendedUpdate) { 2023 var i, j, M, MT, y, B, c, s, d, 2024 // input data 2025 len = dataX.length; 2026 2027 d = Math.floor(deg()); 2028 2029 if (!suspendedUpdate) { 2030 // point list as input 2031 if (inputType === 1) { 2032 dX = []; 2033 dY = []; 2034 2035 for (i = 0; i < len; i++) { 2036 dX[i] = dataX[i].X(); 2037 dY[i] = dataX[i].Y(); 2038 } 2039 } 2040 2041 if (inputType === 2) { 2042 dX = []; 2043 dY = []; 2044 2045 for (i = 0; i < len; i++) { 2046 dX[i] = dataX[i].usrCoords[1]; 2047 dY[i] = dataX[i].usrCoords[2]; 2048 } 2049 } 2050 2051 // check for functions 2052 if (inputType === 0) { 2053 dX = []; 2054 dY = []; 2055 2056 for (i = 0; i < len; i++) { 2057 if (Type.isFunction(dataX[i])) { 2058 dX.push(dataX[i]()); 2059 } else { 2060 dX.push(dataX[i]); 2061 } 2062 2063 if (Type.isFunction(dataY[i])) { 2064 dY.push(dataY[i]()); 2065 } else { 2066 dY.push(dataY[i]); 2067 } 2068 } 2069 } 2070 2071 M = []; 2072 2073 for (j = 0; j < len; j++) { 2074 M.push([1]); 2075 } 2076 2077 for (i = 1; i <= d; i++) { 2078 for (j = 0; j < len; j++) { 2079 M[j][i] = M[j][i - 1] * dX[j]; 2080 } 2081 } 2082 2083 y = dY; 2084 MT = Mat.transpose(M); 2085 B = Mat.matMatMult(MT, M); 2086 c = Mat.matVecMult(MT, y); 2087 coeffs = Mat.Numerics.Gauss(B, c); 2088 term = Mat.Numerics.generatePolynomialTerm(coeffs, d, 'x', 3); 2089 } 2090 2091 // Horner's scheme to evaluate polynomial 2092 s = coeffs[d]; 2093 2094 for (i = d - 1; i >= 0; i--) { 2095 s = (s * x + coeffs[i]); 2096 } 2097 2098 return s; 2099 }; 2100 2101 fct.getTerm = function () { 2102 return term; 2103 }; 2104 2105 return fct; 2106 }, 2107 2108 /** 2109 * Computes the cubic Bezier curve through a given set of points. 2110 * @param {Array} points Array consisting of 3*k+1 {@link JXG.Points}. 2111 * The points at position k with k mod 3 = 0 are the data points, 2112 * points at position k with k mod 3 = 1 or 2 are the control points. 2113 * @returns {Array} An array consisting of two functions of one parameter t which return the 2114 * x resp. y coordinates of the Bezier curve in t, one zero value, and a third function accepting 2115 * no parameters and returning one third of the length of the points. 2116 * @memberof JXG.Math.Numerics 2117 */ 2118 bezier: function (points) { 2119 var len, flen, 2120 /** @ignore */ 2121 makeFct = function (which) { 2122 return function (t, suspendedUpdate) { 2123 var z = Math.floor(t) * 3, 2124 t0 = t % 1, 2125 t1 = 1 - t0; 2126 2127 if (!suspendedUpdate) { 2128 flen = 3 * Math.floor((points.length - 1) / 3); 2129 len = Math.floor(flen / 3); 2130 } 2131 2132 if (t < 0) { 2133 return points[0][which](); 2134 } 2135 2136 if (t >= len) { 2137 return points[flen][which](); 2138 } 2139 2140 if (isNaN(t)) { 2141 return NaN; 2142 } 2143 2144 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; 2145 }; 2146 }; 2147 2148 return [makeFct('X'), makeFct('Y'), 0, 2149 function () { 2150 return Math.floor(points.length / 3); 2151 }]; 2152 }, 2153 2154 /** 2155 * Computes the B-spline curve of order k (order = degree+1) through a given set of points. 2156 * @param {Array} points Array consisting of JXG.Points. 2157 * @param {Number} order Order of the B-spline curve. 2158 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 2159 * which return the x resp. y coordinates of the B-spline curve in t, a zero value, and a function simply 2160 * returning the length of the points array minus one. 2161 * @memberof JXG.Math.Numerics 2162 */ 2163 bspline: function (points, order) { 2164 var knots, N = [], 2165 _knotVector = function (n, k) { 2166 var j, 2167 kn = []; 2168 2169 for (j = 0; j < n + k + 1; j++) { 2170 if (j < k) { 2171 kn[j] = 0.0; 2172 } else if (j <= n) { 2173 kn[j] = j - k + 1; 2174 } else { 2175 kn[j] = n - k + 2; 2176 } 2177 } 2178 2179 return kn; 2180 }, 2181 2182 _evalBasisFuncs = function (t, kn, n, k, s) { 2183 var i, j, a, b, den, 2184 N = []; 2185 2186 if (kn[s] <= t && t < kn[s + 1]) { 2187 N[s] = 1; 2188 } else { 2189 N[s] = 0; 2190 } 2191 2192 for (i = 2; i <= k; i++) { 2193 for (j = s - i + 1; j <= s; j++) { 2194 if (j <= s - i + 1 || j < 0) { 2195 a = 0.0; 2196 } else { 2197 a = N[j]; 2198 } 2199 2200 if (j >= s) { 2201 b = 0.0; 2202 } else { 2203 b = N[j + 1]; 2204 } 2205 2206 den = kn[j + i - 1] - kn[j]; 2207 2208 if (den === 0) { 2209 N[j] = 0; 2210 } else { 2211 N[j] = (t - kn[j]) / den * a; 2212 } 2213 2214 den = kn[j + i] - kn[j + 1]; 2215 2216 if (den !== 0) { 2217 N[j] += (kn[j + i] - t) / den * b; 2218 } 2219 } 2220 } 2221 return N; 2222 }, 2223 /** @ignore */ 2224 makeFct = function (which) { 2225 return function (t, suspendedUpdate) { 2226 var y, j, s, 2227 len = points.length, 2228 n = len - 1, 2229 k = order; 2230 2231 if (n <= 0) { 2232 return NaN; 2233 } 2234 2235 if (n + 2 <= k) { 2236 k = n + 1; 2237 } 2238 2239 if (t <= 0) { 2240 return points[0][which](); 2241 } 2242 2243 if (t >= n - k + 2) { 2244 return points[n][which](); 2245 } 2246 2247 s = Math.floor(t) + k - 1; 2248 knots = _knotVector(n, k); 2249 N = _evalBasisFuncs(t, knots, n, k, s); 2250 2251 y = 0.0; 2252 for (j = s - k + 1; j <= s; j++) { 2253 if (j < len && j >= 0) { 2254 y += points[j][which]() * N[j]; 2255 } 2256 } 2257 2258 return y; 2259 }; 2260 }; 2261 2262 return [makeFct('X'), makeFct('Y'), 0, 2263 function () { 2264 return points.length - 1; 2265 }]; 2266 }, 2267 2268 /** 2269 * Numerical (symmetric) approximation of derivative. suspendUpdate is piped through, 2270 * see {@link JXG.Curve#updateCurve} 2271 * and {@link JXG.Curve#hasPoint}. 2272 * @param {function} f Function in one variable to be differentiated. 2273 * @param {object} [obj] Optional object that is treated as "this" in the function body. This is useful, if the function is a 2274 * method of an object and contains a reference to its parent object via "this". 2275 * @returns {function} Derivative function of a given function f. 2276 * @memberof JXG.Math.Numerics 2277 */ 2278 D: function (f, obj) { 2279 if (!Type.exists(obj)) { 2280 return function (x, suspendUpdate) { 2281 var h = 0.00001, 2282 h2 = (h * 2.0); 2283 2284 // Experiments with Richardsons rule 2285 /* 2286 var phi = (f(x + h, suspendUpdate) - f(x - h, suspendUpdate)) / h2; 2287 var phi2; 2288 h *= 0.5; 2289 h2 *= 0.5; 2290 phi2 = (f(x + h, suspendUpdate) - f(x - h, suspendUpdate)) / h2; 2291 2292 return phi2 + (phi2 - phi) / 3.0; 2293 */ 2294 return (f(x + h, suspendUpdate) - f(x - h, suspendUpdate)) / h2; 2295 }; 2296 } 2297 2298 return function (x, suspendUpdate) { 2299 var h = 0.00001, 2300 h2 = (h * 2.0); 2301 2302 return (f.apply(obj, [x + h, suspendUpdate]) - f.apply(obj, [x - h, suspendUpdate])) / h2; 2303 }; 2304 }, 2305 2306 /** 2307 * Evaluate the function term for {@see #riemann}. 2308 * @private 2309 * @param {Number} x function argument 2310 * @param {function} f JavaScript function returning a number 2311 * @param {String} type Name of the Riemann sum type, e.g. 'lower', see {@see #riemann}. 2312 * @param {Number} delta Width of the bars in user coordinates 2313 * @returns {Number} Upper (delta > 0) or lower (delta < 0) value of the bar containing x of the Riemann sum. 2314 * 2315 * @memberof JXG.Math.Numerics 2316 */ 2317 _riemannValue: function (x, f, type, delta) { 2318 var y, y1, x1, delta1; 2319 2320 if (delta < 0) { // delta is negative if the lower function term is evaluated 2321 if (type !== 'trapezoidal') { 2322 x = x + delta; 2323 } 2324 delta *= -1; 2325 if (type === 'lower') { 2326 type = 'upper'; 2327 } else if (type === 'upper') { 2328 type = 'lower'; 2329 } 2330 } 2331 2332 delta1 = delta * 0.01; // for 'lower' and 'upper' 2333 2334 if (type === 'right') { 2335 y = f(x + delta); 2336 } else if (type === 'middle') { 2337 y = f(x + delta * 0.5); 2338 } else if (type === 'left' || type === 'trapezoidal') { 2339 y = f(x); 2340 } else if (type === 'lower') { 2341 y = f(x); 2342 2343 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 2344 y1 = f(x1); 2345 2346 if (y1 < y) { 2347 y = y1; 2348 } 2349 } 2350 2351 y1 = f(x + delta); 2352 if (y1 < y) { 2353 y = y1; 2354 } 2355 } else if (type === 'upper') { 2356 y = f(x); 2357 2358 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 2359 y1 = f(x1); 2360 if (y1 > y) { 2361 y = y1; 2362 } 2363 } 2364 2365 y1 = f(x + delta); 2366 if (y1 > y) { 2367 y = y1; 2368 } 2369 } else if (type === 'random') { 2370 y = f(x + delta * Math.random()); 2371 } else if (type === 'simpson') { 2372 y = (f(x) + 4 * f(x + delta * 0.5) + f(x + delta)) / 6.0; 2373 } else { 2374 y = f(x); // default is lower 2375 } 2376 2377 return y; 2378 }, 2379 2380 /** 2381 * Helper function to create curve which displays Riemann sums. 2382 * Compute coordinates for the rectangles showing the Riemann sum. 2383 * @param {Function,Array} f Function or array of two functions. 2384 * If f is a function the integral of this function is approximated by the Riemann sum. 2385 * If f is an array consisting of two functions the area between the two functions is filled 2386 * by the Riemann sum bars. 2387 * @param {Number} n number of rectangles. 2388 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson', or 'trapezoidal'. 2389 * @param {Number} start Left border of the approximation interval 2390 * @param {Number} end Right border of the approximation interval 2391 * @returns {Array} An array of two arrays containing the x and y coordinates for the rectangles showing the Riemann sum. This 2392 * 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 2393 * rectangles. 2394 * @memberof JXG.Math.Numerics 2395 */ 2396 riemann: function (gf, n, type, start, end) { 2397 var i, delta, 2398 xarr = [], 2399 yarr = [], 2400 j = 0, 2401 x = start, y, 2402 sum = 0, 2403 f, g, 2404 ylow, yup; 2405 2406 if (Type.isArray(gf)) { 2407 g = gf[0]; 2408 f = gf[1]; 2409 } else { 2410 f = gf; 2411 } 2412 2413 n = Math.floor(n); 2414 2415 if (n <= 0) { 2416 return [xarr, yarr, sum]; 2417 } 2418 2419 delta = (end - start) / n; 2420 2421 // Upper bar ends 2422 for (i = 0; i < n; i++) { 2423 y = this._riemannValue(x, f, type, delta); 2424 xarr[j] = x; 2425 yarr[j] = y; 2426 2427 j += 1; 2428 x += delta; 2429 if (type === 'trapezoidal') { 2430 y = f(x); 2431 } 2432 xarr[j] = x; 2433 yarr[j] = y; 2434 2435 j += 1; 2436 } 2437 2438 // Lower bar ends 2439 for (i = 0; i < n; i++) { 2440 if (g) { 2441 y = this._riemannValue(x, g, type, -delta); 2442 } else { 2443 y = 0.0; 2444 } 2445 xarr[j] = x; 2446 yarr[j] = y; 2447 2448 j += 1; 2449 x -= delta; 2450 if (type === 'trapezoidal' && g) { 2451 y = g(x); 2452 } 2453 xarr[j] = x; 2454 yarr[j] = y; 2455 2456 // Add the area of the bar to 'sum' 2457 if (type !== 'trapezoidal') { 2458 ylow = y; 2459 yup = yarr[2 * (n - 1) - 2 * i]; 2460 } else { 2461 yup = 0.5 * (f(x + delta) + f(x)); 2462 if (g) { 2463 ylow = 0.5 * (g(x + delta) + g(x)); 2464 } else { 2465 ylow = 0.0; 2466 } 2467 } 2468 sum += (yup - ylow) * delta; 2469 2470 // Draw the vertical lines 2471 j += 1; 2472 xarr[j] = x; 2473 yarr[j] = yarr[2 * (n - 1) - 2 * i]; 2474 2475 j += 1; 2476 } 2477 2478 return [xarr, yarr, sum]; 2479 }, 2480 2481 /** 2482 * Approximate the integral by Riemann sums. 2483 * Compute the area described by the riemann sum rectangles. 2484 * 2485 * If there is an element of type {@link Riemannsum}, then it is more efficient 2486 * to use the method JXG.Curve.Value() of this element instead. 2487 * 2488 * @param {Function_Array} f Function or array of two functions. 2489 * If f is a function the integral of this function is approximated by the Riemann sum. 2490 * If f is an array consisting of two functions the area between the two functions is approximated 2491 * by the Riemann sum. 2492 * @param {Number} n number of rectangles. 2493 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson' or 'trapezoidal'. 2494 * 2495 * @param {Number} start Left border of the approximation interval 2496 * @param {Number} end Right border of the approximation interval 2497 * @returns {Number} The sum of the areas of the rectangles. 2498 * @memberof JXG.Math.Numerics 2499 */ 2500 riemannsum: function (f, n, type, start, end) { 2501 JXG.deprecated('Numerics.riemannsum()', 'Numerics.riemann()'); 2502 return this.riemann(f, n, type, start, end)[2]; 2503 }, 2504 2505 /** 2506 * Solve initial value problems numerically using Runge-Kutta-methods. 2507 * See {@link http://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm. 2508 * @param {object,String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing 2509 * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure 2510 * <pre> 2511 * { 2512 * s: <Number>, 2513 * A: <matrix>, 2514 * b: <Array>, 2515 * c: <Array> 2516 * } 2517 * </pre> 2518 * 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 2519 * @param {Array} x0 Initial value vector. If the problem is of one-dimensional, the initial value also has to be given in an array. 2520 * @param {Array} I Interval on which to integrate. 2521 * @param {Number} N Number of evaluation points. 2522 * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode 2523 * 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 2524 * vector <tt>x</tt>, and has to return a vector of the same dimension as <tt>x</tt> has. 2525 * @returns {Array} An array of vectors describing the solution of the ode on the given interval I. 2526 * @example 2527 * // A very simple autonomous system dx(t)/dt = x(t); 2528 * function f(t, x) { 2529 * return x; 2530 * } 2531 * 2532 * // Solve it with initial value x(0) = 1 on the interval [0, 2] 2533 * // with 20 evaluation points. 2534 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 2535 * 2536 * // Prepare data for plotting the solution of the ode using a curve. 2537 * var dataX = []; 2538 * var dataY = []; 2539 * var h = 0.1; // (I[1] - I[0])/N = (2-0)/20 2540 * for(var i=0; i<data.length; i++) { 2541 * dataX[i] = i*h; 2542 * dataY[i] = data[i][0]; 2543 * } 2544 * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'}); 2545 * </pre><div class="jxgbox" id="d2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div> 2546 * <script type="text/javascript"> 2547 * var board = JXG.JSXGraph.initBoard('d2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false}); 2548 * function f(t, x) { 2549 * // we have to copy the value. 2550 * // return x; would just return the reference. 2551 * return [x[0]]; 2552 * } 2553 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 2554 * var dataX = []; 2555 * var dataY = []; 2556 * var h = 0.1; 2557 * for(var i=0; i<data.length; i++) { 2558 * dataX[i] = i*h; 2559 * dataY[i] = data[i][0]; 2560 * } 2561 * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'}); 2562 * </script><pre> 2563 * @memberof JXG.Math.Numerics 2564 */ 2565 rungeKutta: function (butcher, x0, I, N, f) { 2566 var e, i, j, k, l, s, 2567 x = [], 2568 y = [], 2569 h = (I[1] - I[0]) / N, 2570 t = I[0], 2571 dim = x0.length, 2572 result = [], 2573 r = 0; 2574 2575 if (Type.isString(butcher)) { 2576 butcher = predefinedButcher[butcher] || predefinedButcher.euler; 2577 } 2578 s = butcher.s; 2579 2580 // don't change x0, so copy it 2581 for (e = 0; e < dim; e++) { 2582 x[e] = x0[e]; 2583 } 2584 2585 for (i = 0; i < N; i++) { 2586 // Optimization doesn't work for ODEs plotted using time 2587 // if((i % quotient == 0) || (i == N-1)) { 2588 result[r] = []; 2589 for (e = 0; e < dim; e++) { 2590 result[r][e] = x[e]; 2591 } 2592 2593 r += 1; 2594 k = []; 2595 2596 for (j = 0; j < s; j++) { 2597 // init y = 0 2598 for (e = 0; e < dim; e++) { 2599 y[e] = 0.0; 2600 } 2601 2602 2603 // Calculate linear combination of former k's and save it in y 2604 for (l = 0; l < j; l++) { 2605 for (e = 0; e < dim; e++) { 2606 y[e] += (butcher.A[j][l]) * h * k[l][e]; 2607 } 2608 } 2609 2610 // add x(t) to y 2611 for (e = 0; e < dim; e++) { 2612 y[e] += x[e]; 2613 } 2614 2615 // calculate new k and add it to the k matrix 2616 k.push(f(t + butcher.c[j] * h, y)); 2617 } 2618 2619 // init y = 0 2620 for (e = 0; e < dim; e++) { 2621 y[e] = 0.0; 2622 } 2623 2624 for (l = 0; l < s; l++) { 2625 for (e = 0; e < dim; e++) { 2626 y[e] += butcher.b[l] * k[l][e]; 2627 } 2628 } 2629 2630 for (e = 0; e < dim; e++) { 2631 x[e] = x[e] + h * y[e]; 2632 } 2633 2634 t += h; 2635 } 2636 2637 return result; 2638 }, 2639 2640 /** 2641 * Maximum number of iterations in {@link JXG.Math.Numerics.fzero} 2642 * @type Number 2643 * @default 80 2644 * @memberof JXG.Math.Numerics 2645 */ 2646 maxIterationsRoot: 80, 2647 2648 /** 2649 * Maximum number of iterations in {@link JXG.Math.Numerics.fminbr} 2650 * @type Number 2651 * @default 500 2652 * @memberof JXG.Math.Numerics 2653 */ 2654 maxIterationsMinimize: 500, 2655 2656 /** 2657 * 2658 * Find zero of an univariate function f. 2659 * @param {function} f Function, whose root is to be found 2660 * @param {Array,Number} x0 Start value or start interval enclosing the root 2661 * @param {Object} object Parent object in case f is method of it 2662 * @returns {Number} the approximation of the root 2663 * Algorithm: 2664 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 2665 * computations. M., Mir, 1980, p.180 of the Russian edition 2666 * 2667 * If x0 is an array containing lower and upper bound for the zero 2668 * algorithm 748 is applied. Otherwise, if x0 is a number, 2669 * the algorithm tries to bracket a zero of f starting from x0. 2670 * If this fails, we fall back to Newton's method. 2671 * @memberof JXG.Math.Numerics 2672 */ 2673 fzero: function (f, x0, object) { 2674 var a, b, c, 2675 fa, fb, fc, 2676 aa, blist, i, len, u, fu, 2677 prev_step, t1, cb, t2, 2678 // Actual tolerance 2679 tol_act, 2680 // Interpolation step is calculated in the form p/q; division 2681 // operations is delayed until the last moment 2682 p, q, 2683 // Step at this iteration 2684 new_step, 2685 eps = Mat.eps, 2686 maxiter = this.maxIterationsRoot, 2687 niter = 0, 2688 nfev = 0; 2689 2690 if (Type.isArray(x0)) { 2691 if (x0.length < 2) { 2692 throw new Error("JXG.Math.Numerics.fzero: length of array x0 has to be at least two."); 2693 } 2694 2695 a = x0[0]; 2696 fa = f.call(object, a); 2697 nfev += 1; 2698 b = x0[1]; 2699 fb = f.call(object, b); 2700 nfev += 1; 2701 } else { 2702 a = x0; 2703 fa = f.call(object, a); 2704 nfev += 1; 2705 2706 // Try to get b. 2707 if (a === 0) { 2708 aa = 1; 2709 } else { 2710 aa = a; 2711 } 2712 2713 blist = [0.9 * aa, 1.1 * aa, aa - 1, aa + 1, 0.5 * aa, 1.5 * aa, -aa, 2 * aa, -10 * aa, 10 * aa]; 2714 len = blist.length; 2715 2716 for (i = 0; i < len; i++) { 2717 b = blist[i]; 2718 fb = f.call(object, b); 2719 nfev += 1; 2720 2721 if (fa * fb <= 0) { 2722 break; 2723 } 2724 } 2725 if (b < a) { 2726 u = a; 2727 a = b; 2728 b = u; 2729 2730 fu = fa; 2731 fa = fb; 2732 fb = fu; 2733 } 2734 } 2735 2736 if (fa * fb > 0) { 2737 // Bracketing not successful, fall back to Newton's method or to fminbr 2738 if (Type.isArray(x0)) { 2739 return this.fminbr(f, [a, b], object); 2740 } 2741 2742 return this.Newton(f, a, object); 2743 } 2744 2745 // OK, we have enclosed a zero of f. 2746 // Now we can start Brent's method 2747 2748 c = a; 2749 fc = fa; 2750 2751 // Main iteration loop 2752 while (niter < maxiter) { 2753 // Distance from the last but one to the last approximation 2754 prev_step = b - a; 2755 2756 // Swap data for b to be the best approximation 2757 if (Math.abs(fc) < Math.abs(fb)) { 2758 a = b; 2759 b = c; 2760 c = a; 2761 2762 fa = fb; 2763 fb = fc; 2764 fc = fa; 2765 } 2766 tol_act = 2 * eps * Math.abs(b) + eps * 0.5; 2767 new_step = (c - b) * 0.5; 2768 2769 if (Math.abs(new_step) <= tol_act && Math.abs(fb) <= eps) { 2770 // Acceptable approx. is found 2771 return b; 2772 } 2773 2774 // Decide if the interpolation can be tried 2775 // If prev_step was large enough and was in true direction Interpolatiom may be tried 2776 if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) { 2777 cb = c - b; 2778 2779 // If we have only two distinct points linear interpolation can only be applied 2780 if (a === c) { 2781 t1 = fb / fa; 2782 p = cb * t1; 2783 q = 1.0 - t1; 2784 // Quadric inverse interpolation 2785 } else { 2786 q = fa / fc; 2787 t1 = fb / fc; 2788 t2 = fb / fa; 2789 2790 p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0)); 2791 q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0); 2792 } 2793 2794 // p was calculated with the opposite sign; make p positive 2795 if (p > 0) { 2796 q = -q; 2797 // and assign possible minus to q 2798 } else { 2799 p = -p; 2800 } 2801 2802 // If b+p/q falls in [b,c] and isn't too large it is accepted 2803 // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent 2804 if (p < (0.75 * cb * q - Math.abs(tol_act * q) * 0.5) && 2805 p < Math.abs(prev_step * q * 0.5)) { 2806 new_step = p / q; 2807 } 2808 } 2809 2810 // Adjust the step to be not less than tolerance 2811 if (Math.abs(new_step) < tol_act) { 2812 if (new_step > 0) { 2813 new_step = tol_act; 2814 } else { 2815 new_step = -tol_act; 2816 } 2817 } 2818 2819 // Save the previous approx. 2820 a = b; 2821 fa = fb; 2822 b += new_step; 2823 fb = f.call(object, b); 2824 // Do step to a new approxim. 2825 nfev += 1; 2826 2827 // Adjust c for it to have a sign opposite to that of b 2828 if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) { 2829 c = a; 2830 fc = fa; 2831 } 2832 niter++; 2833 } // End while 2834 2835 return b; 2836 }, 2837 2838 /** 2839 * 2840 * Find minimum of an univariate function f. 2841 * <p> 2842 * Algorithm: 2843 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 2844 * computations. M., Mir, 1980, p.180 of the Russian edition 2845 * 2846 * @param {function} f Function, whose minimum is to be found 2847 * @param {Array} x0 Start interval enclosing the minimum 2848 * @param {Object} context Parent object in case f is method of it 2849 * @returns {Number} the approximation of the minimum value position 2850 * @memberof JXG.Math.Numerics 2851 **/ 2852 fminbr: function (f, x0, context) { 2853 var a, b, x, v, w, 2854 fx, fv, fw, 2855 range, middle_range, tol_act, new_step, 2856 p, q, t, ft, 2857 // Golden section ratio 2858 r = (3.0 - Math.sqrt(5.0)) * 0.5, 2859 tol = Mat.eps, 2860 sqrteps = Mat.eps, //Math.sqrt(Mat.eps), 2861 maxiter = this.maxIterationsMinimize, 2862 niter = 0, 2863 nfev = 0; 2864 2865 if (!Type.isArray(x0) || x0.length < 2) { 2866 throw new Error("JXG.Math.Numerics.fminbr: length of array x0 has to be at least two."); 2867 } 2868 2869 a = x0[0]; 2870 b = x0[1]; 2871 v = a + r * (b - a); 2872 fv = f.call(context, v); 2873 2874 // First step - always gold section 2875 nfev += 1; 2876 x = v; 2877 w = v; 2878 fx = fv; 2879 fw = fv; 2880 2881 while (niter < maxiter) { 2882 // Range over the interval in which we are looking for the minimum 2883 range = b - a; 2884 middle_range = (a + b) * 0.5; 2885 2886 // Actual tolerance 2887 tol_act = sqrteps * Math.abs(x) + tol / 3.0; 2888 2889 if (Math.abs(x - middle_range) + range * 0.5 <= 2.0 * tol_act) { 2890 // Acceptable approx. is found 2891 return x; 2892 } 2893 2894 // Obtain the golden section step 2895 new_step = r * (x < middle_range ? b - x : a - x); 2896 2897 // Decide if the interpolation can be tried. If x and w are distinct interpolatiom may be tried 2898 if (Math.abs(x - w) >= tol_act) { 2899 // Interpolation step is calculated as p/q; 2900 // division operation is delayed until last moment 2901 t = (x - w) * (fx - fv); 2902 q = (x - v) * (fx - fw); 2903 p = (x - v) * q - (x - w) * t; 2904 q = 2 * (q - t); 2905 2906 if (q > 0) { // q was calculated with the op- 2907 p = -p; // posite sign; make q positive 2908 } else { // and assign possible minus to 2909 q = -q; // p 2910 } 2911 if (Math.abs(p) < Math.abs(new_step * q) && // If x+p/q falls in [a,b] 2912 p > q * (a - x + 2 * tol_act) && // not too close to a and 2913 p < q * (b - x - 2 * tol_act)) { // b, and isn't too large 2914 new_step = p / q; // it is accepted 2915 } 2916 // If p/q is too large then the 2917 // golden section procedure can 2918 // reduce [a,b] range to more 2919 // extent 2920 } 2921 2922 // Adjust the step to be not less than tolerance 2923 if (Math.abs(new_step) < tol_act) { 2924 if (new_step > 0) { 2925 new_step = tol_act; 2926 } else { 2927 new_step = -tol_act; 2928 } 2929 } 2930 2931 // Obtain the next approximation to min 2932 // and reduce the enveloping range 2933 2934 // Tentative point for the min 2935 t = x + new_step; 2936 ft = f.call(context, t); 2937 nfev += 1; 2938 2939 // t is a better approximation 2940 if (ft <= fx) { 2941 // Reduce the range so that t would fall within it 2942 if (t < x) { 2943 b = x; 2944 } else { 2945 a = x; 2946 } 2947 2948 // Assign the best approx to x 2949 v = w; 2950 w = x; 2951 x = t; 2952 2953 fv = fw; 2954 fw = fx; 2955 fx = ft; 2956 // x remains the better approx 2957 } else { 2958 // Reduce the range enclosing x 2959 if (t < x) { 2960 a = t; 2961 } else { 2962 b = t; 2963 } 2964 2965 if (ft <= fw || w === x) { 2966 v = w; 2967 w = t; 2968 fv = fw; 2969 fw = ft; 2970 } else if (ft <= fv || v === x || v === w) { 2971 v = t; 2972 fv = ft; 2973 } 2974 } 2975 niter += 1; 2976 } 2977 2978 return x; 2979 }, 2980 2981 /** 2982 * Implements the Ramer-Douglas-Peucker algorithm. 2983 * It discards points which are not necessary from the polygonal line defined by the point array 2984 * pts. The computation is done in screen coordinates. 2985 * Average runtime is O(nlog(n)), worst case runtime is O(n^2), where n is the number of points. 2986 * @param {Array} pts Array of {@link JXG.Coords} 2987 * @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>. 2988 * @returns {Array} An array containing points which represent an apparently identical curve as the points of pts do, but contains fewer points. 2989 * @memberof JXG.Math.Numerics 2990 */ 2991 RamerDouglasPeucker: function (pts, eps) { 2992 var newPts = [], i, k, len, 2993 2994 /** 2995 * findSplit() is a subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}. 2996 * It searches for the point between index i and j which 2997 * has the largest distance from the line between the points i and j. 2998 * @param {Array} pts Array of {@link JXG.Coords} 2999 * @param {Number} i Index of a point in pts 3000 * @param {Number} j Index of a point in pts 3001 * @ignore 3002 * @private 3003 */ 3004 findSplit = function (pts, i, j) { 3005 var d, k, ci, cj, ck, 3006 x0, y0, x1, y1, 3007 den, lbda, 3008 dist = 0, 3009 f = i; 3010 3011 if (j - i < 2) { 3012 return [-1.0, 0]; 3013 } 3014 3015 ci = pts[i].scrCoords; 3016 cj = pts[j].scrCoords; 3017 3018 if (isNaN(ci[1] + ci[2])) { 3019 return [NaN, i]; 3020 } 3021 if (isNaN(cj[1] + cj[2])) { 3022 return [NaN, j]; 3023 } 3024 3025 for (k = i + 1; k < j; k++) { 3026 ck = pts[k].scrCoords; 3027 if (isNaN(ck[1] + ck[2])) { 3028 return [NaN, k]; 3029 } 3030 3031 x0 = ck[1] - ci[1]; 3032 y0 = ck[2] - ci[2]; 3033 x1 = cj[1] - ci[1]; 3034 y1 = cj[2] - ci[2]; 3035 den = x1 * x1 + y1 * y1; 3036 3037 if (den >= Mat.eps) { 3038 lbda = (x0 * x1 + y0 * y1) / den; 3039 3040 if (lbda < 0.0) { 3041 lbda = 0.0; 3042 } else if (lbda > 1.0) { 3043 lbda = 1.0; 3044 } 3045 3046 x0 = x0 - lbda * x1; 3047 y0 = y0 - lbda * y1; 3048 d = x0 * x0 + y0 * y0; 3049 } else { 3050 lbda = 0.0; 3051 d = x0 * x0 + y0 * y0; 3052 } 3053 3054 if (d > dist) { 3055 dist = d; 3056 f = k; 3057 } 3058 } 3059 return [Math.sqrt(dist), f]; 3060 }, 3061 3062 /** 3063 * RDP() is a private subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}. 3064 * It runs recursively through the point set and searches the 3065 * point which has the largest distance from the line between the first point and 3066 * the last point. If the distance from the line is greater than eps, this point is 3067 * included in our new point set otherwise it is discarded. 3068 * If it is taken, we recursively apply the subroutine to the point set before 3069 * and after the chosen point. 3070 * @param {Array} pts Array of {@link JXG.Coords} 3071 * @param {Number} i Index of an element of pts 3072 * @param {Number} j Index of an element of pts 3073 * @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>. 3074 * @param {Array} newPts Array of {@link JXG.Coords} 3075 * @ignore 3076 * @private 3077 */ 3078 RDP = function (pts, i, j, eps, newPts) { 3079 var result = findSplit(pts, i, j), 3080 k = result[1]; 3081 3082 if (isNaN(result[0])) { 3083 RDP(pts, i, k - 1, eps, newPts); 3084 newPts.push(pts[k]); 3085 do { 3086 ++k; 3087 } while (k <= j && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])); 3088 if (k <= j) { 3089 newPts.push(pts[k]); 3090 } 3091 RDP(pts, k + 1, j, eps, newPts); 3092 } else if (result[0] > eps) { 3093 RDP(pts, i, k, eps, newPts); 3094 RDP(pts, k, j, eps, newPts); 3095 } else { 3096 newPts.push(pts[j]); 3097 } 3098 }; 3099 3100 len = pts.length; 3101 3102 // Search for the left most point woithout NaN coordinates 3103 i = 0; 3104 while (i < len && isNaN(pts[i].scrCoords[1] + pts[i].scrCoords[2])) { 3105 i += 1; 3106 } 3107 // Search for the right most point woithout NaN coordinates 3108 k = len - 1; 3109 while (k > i && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])) { 3110 k -= 1; 3111 } 3112 3113 // Only proceed if something is left 3114 if (!(i > k || i === len)) { 3115 newPts[0] = pts[i]; 3116 RDP(pts, i, k, eps, newPts); 3117 } 3118 3119 return newPts; 3120 }, 3121 3122 /** 3123 * Old name for the implementation of the Ramer-Douglas-Peucker algorithm. 3124 * @deprecated Use {@link JXG.Math.Numerics.RamerDouglasPeucker} 3125 * @memberof JXG.Math.Numerics 3126 */ 3127 RamerDouglasPeuker: function (pts, eps) { 3128 JXG.deprecated('Numerics.RamerDouglasPeuker()', 'Numerics.RamerDouglasPeucker()'); 3129 return this.RamerDouglasPeucker(pts, eps); 3130 }, 3131 3132 /** 3133 * Implements the Visvalingam-Whyatt algorithm. 3134 * See M. Visvalingam, J. D. Whyatt: 3135 * "Line generalisation by repeated elimination of the smallest area", C.I.S.R.G Discussion paper 10, July 1992 3136 * 3137 * The algorithm discards points which are not necessary from the polygonal line defined by the point array 3138 * pts (consisting of type JXG.Coords). 3139 * @param {Array} pts Array of {@link JXG.Coords} 3140 * @param {Number} numPoints Number of remaining intermediate points. The first and the last point of the original points will 3141 * be taken in any case. 3142 * @returns {Array} An array containing points which approximates the curve defined by pts. 3143 * @memberof JXG.Math.Numerics 3144 * 3145 * @example 3146 * var i, p = []; 3147 * for (i = 0; i < 5; ++i) { 3148 * p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6])); 3149 * } 3150 * 3151 * // Plot a cardinal spline curve 3152 * var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5); 3153 * var cu1 = board.create('curve', splineArr, {strokeColor: 'green'}); 3154 * 3155 * var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'}); 3156 * c.updateDataArray = function() { 3157 * var i, len, points; 3158 * 3159 * // Reduce number of intermediate points with Visvakingam-Whyatt to 6 3160 * points = JXG.Math.Numerics.Visvalingam(cu1.points, 6); 3161 * // Plot the remaining points 3162 * len = points.length; 3163 * this.dataX = []; 3164 * this.dataY = []; 3165 * for (i = 0; i < len; i++) { 3166 * this.dataX.push(points[i].usrCoords[1]); 3167 * this.dataY.push(points[i].usrCoords[2]); 3168 * } 3169 * }; 3170 * board.update(); 3171 * 3172 * </pre><div id="ce0cc55c-b592-11e6-8270-104a7d3be7eb" class="jxgbox" style="width: 300px; height: 300px;"></div> 3173 * <script type="text/javascript"> 3174 * (function() { 3175 * var board = JXG.JSXGraph.initBoard('ce0cc55c-b592-11e6-8270-104a7d3be7eb', 3176 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 3177 * 3178 * var i, p = []; 3179 * for (i = 0; i < 5; ++i) { 3180 * p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6])); 3181 * } 3182 * 3183 * // Plot a cardinal spline curve 3184 * var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5); 3185 * var cu1 = board.create('curve', splineArr, {strokeColor: 'green'}); 3186 * 3187 * var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'}); 3188 * c.updateDataArray = function() { 3189 * var i, len, points; 3190 * 3191 * // Reduce number of intermediate points with Visvakingam-Whyatt to 6 3192 * points = JXG.Math.Numerics.Visvalingam(cu1.points, 6); 3193 * // Plot the remaining points 3194 * len = points.length; 3195 * this.dataX = []; 3196 * this.dataY = []; 3197 * for (i = 0; i < len; i++) { 3198 * this.dataX.push(points[i].usrCoords[1]); 3199 * this.dataY.push(points[i].usrCoords[2]); 3200 * } 3201 * }; 3202 * board.update(); 3203 * 3204 * })(); 3205 * 3206 * </script><pre> 3207 * 3208 */ 3209 Visvalingam: function(pts, numPoints) { 3210 var i, len, vol, lastVol, 3211 linkedList = [], 3212 heap = [], 3213 points = [], 3214 lft, rt, lft2, rt2, 3215 obj; 3216 3217 len = pts.length; 3218 3219 // At least one intermediate point is needed 3220 if (len <= 2) { 3221 return pts; 3222 } 3223 3224 // Fill the linked list 3225 // Add first point to the linked list 3226 linkedList[0] = { 3227 used: true, 3228 lft: null, 3229 node: null 3230 }; 3231 3232 // Add all intermediate points to the linked list, 3233 // whose triangle area is nonzero. 3234 lft = 0; 3235 for (i = 1; i < len - 1; i++) { 3236 vol = Math.abs(JXG.Math.Numerics.det([pts[i - 1].usrCoords, 3237 pts[i].usrCoords, 3238 pts[i + 1].usrCoords])); 3239 if (!isNaN(vol)) { 3240 obj = { 3241 v: vol, 3242 idx: i 3243 }; 3244 heap.push(obj); 3245 linkedList[i] = { 3246 used: true, 3247 lft: lft, 3248 node: obj 3249 }; 3250 linkedList[lft].rt = i; 3251 lft = i; 3252 } 3253 } 3254 3255 // Add last point to the linked list 3256 linkedList[len - 1] = { 3257 used: true, 3258 rt: null, 3259 lft: lft, 3260 node: null 3261 }; 3262 linkedList[lft].rt = len - 1; 3263 3264 // Remove points until only numPoints intermediate points remain 3265 lastVol = -Infinity; 3266 while (heap.length > numPoints) { 3267 // Sort the heap with the updated volume values 3268 heap.sort(function(a, b) { 3269 // descending sort 3270 return b.v - a.v; 3271 }); 3272 3273 // Remove the point with the smallest triangle 3274 i = heap.pop().idx; 3275 linkedList[i].used = false; 3276 lastVol = linkedList[i].node.v; 3277 3278 // Update the pointers of the linked list 3279 lft = linkedList[i].lft; 3280 rt = linkedList[i].rt; 3281 linkedList[lft].rt = rt; 3282 linkedList[rt].lft = lft; 3283 3284 // Update the values for the volumes in the linked list 3285 lft2 = linkedList[lft].lft; 3286 if (lft2 !== null) { 3287 vol = Math.abs(JXG.Math.Numerics.det( 3288 [pts[lft2].usrCoords, 3289 pts[lft].usrCoords, 3290 pts[rt].usrCoords])); 3291 3292 linkedList[lft].node.v = (vol >= lastVol) ? vol : lastVol; 3293 } 3294 rt2 = linkedList[rt].rt; 3295 if (rt2 !== null) { 3296 vol = Math.abs(JXG.Math.Numerics.det( 3297 [pts[lft].usrCoords, 3298 pts[rt].usrCoords, 3299 pts[rt2].usrCoords])); 3300 3301 linkedList[rt].node.v = (vol >= lastVol) ? vol : lastVol; 3302 } 3303 } 3304 3305 // Return an array with the remaining points 3306 i = 0; 3307 points = [pts[i]]; 3308 do { 3309 i = linkedList[i].rt; 3310 points.push(pts[i]); 3311 } while (linkedList[i].rt !== null); 3312 3313 return points; 3314 } 3315 3316 }; 3317 3318 return Mat.Numerics; 3319 }); 3320