1 /* 2 Copyright 2008-2018 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, 1740 k, xi, s, //M, 1741 len = p.length, 1742 num = 0, 1743 denom = 0; 1744 1745 if (!suspendedUpdate) { 1746 for (i = 0; i < len; i++) { 1747 w[i] = 1.0; 1748 xi = p[i].X(); 1749 1750 for (k = 0; k < len; k++) { 1751 if (k !== i) { 1752 w[i] *= (xi - p[k].X()); 1753 } 1754 } 1755 1756 w[i] = 1 / w[i]; 1757 } 1758 1759 // M = []; 1760 // for (k = 0; k < len; k++) { 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 } 1936 if (t >= len) { 1937 return points[len - 1][which](); 1938 } 1939 1940 s = Math.floor(t); 1941 if (s === t) { 1942 return points[s][which](); 1943 } 1944 1945 t -= s; 1946 c = coeffs[which][s]; 1947 if (c === undefined) { 1948 return NaN; 1949 } 1950 1951 return (((c[3] * t + c[2]) * t + c[1]) * t + c[0]); 1952 }; 1953 }; 1954 1955 return [makeFct('X'), makeFct('Y'), 0, 1956 function () { 1957 return points.length - 1; 1958 }]; 1959 }, 1960 1961 /** 1962 * Computes the cubic Catmull-Rom spline curve through a given set of points. The curve 1963 * is uniformly parametrized. The curve is the cardinal spline curve for tau=0.5. 1964 * Two artificial control points at the beginning and the end are added. 1965 * @param {Array} points Array consisting of JXG.Points. 1966 * @param {String} type (Optional) parameter which allows to choose between "uniform" (default) and 1967 * "centripetal" parameterization. Thus the two possible values are "uniform" or "centripetal". 1968 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 1969 * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, and a function simply 1970 * returning the length of the points array minus three. 1971 * @memberof JXG.Math.Numerics 1972 */ 1973 CatmullRomSpline: function (points, type) { 1974 return this.CardinalSpline(points, 0.5, type); 1975 }, 1976 1977 /** 1978 * Computes the regression polynomial of a given degree through a given set of coordinates. 1979 * Returns the regression polynomial function. 1980 * @param {Number,function,Slider} degree number, function or slider. 1981 * Either 1982 * @param {Array} dataX Array containing either the x-coordinates of the data set or both coordinates in 1983 * an array of {@link JXG.Point}s or {@link JXG.Coords}. 1984 * In the latter case, the <tt>dataY</tt> parameter will be ignored. 1985 * @param {Array} dataY Array containing the y-coordinates of the data set, 1986 * @returns {function} A function of one parameter which returns the value of the regression polynomial of the given degree. 1987 * It possesses the method getTerm() which returns the string containing the function term of the polynomial. 1988 * @memberof JXG.Math.Numerics 1989 */ 1990 regressionPolynomial: function (degree, dataX, dataY) { 1991 var coeffs, deg, dX, dY, inputType, fct, 1992 term = ''; 1993 1994 // Slider 1995 if (Type.isPoint(degree) && Type.isFunction(degree.Value)) { 1996 /** @ignore */ 1997 deg = function () { 1998 return degree.Value(); 1999 }; 2000 // function 2001 } else if (Type.isFunction(degree)) { 2002 deg = degree; 2003 // number 2004 } else if (Type.isNumber(degree)) { 2005 /** @ignore */ 2006 deg = function () { 2007 return degree; 2008 }; 2009 } else { 2010 throw new Error("JSXGraph: Can't create regressionPolynomial from degree of type'" + (typeof degree) + "'."); 2011 } 2012 2013 // Parameters degree, dataX, dataY 2014 if (arguments.length === 3 && Type.isArray(dataX) && Type.isArray(dataY)) { 2015 inputType = 0; 2016 // Parameters degree, point array 2017 } else if (arguments.length === 2 && Type.isArray(dataX) && dataX.length > 0 && Type.isPoint(dataX[0])) { 2018 inputType = 1; 2019 } else if (arguments.length === 2 && Type.isArray(dataX) && dataX.length > 0 && dataX[0].usrCoords && dataX[0].scrCoords) { 2020 inputType = 2; 2021 } else { 2022 throw new Error("JSXGraph: Can't create regressionPolynomial. Wrong parameters."); 2023 } 2024 2025 /** @ignore */ 2026 fct = function (x, suspendedUpdate) { 2027 var i, j, M, MT, y, B, c, s, d, 2028 // input data 2029 len = dataX.length; 2030 2031 d = Math.floor(deg()); 2032 2033 if (!suspendedUpdate) { 2034 // point list as input 2035 if (inputType === 1) { 2036 dX = []; 2037 dY = []; 2038 2039 for (i = 0; i < len; i++) { 2040 dX[i] = dataX[i].X(); 2041 dY[i] = dataX[i].Y(); 2042 } 2043 } 2044 2045 if (inputType === 2) { 2046 dX = []; 2047 dY = []; 2048 2049 for (i = 0; i < len; i++) { 2050 dX[i] = dataX[i].usrCoords[1]; 2051 dY[i] = dataX[i].usrCoords[2]; 2052 } 2053 } 2054 2055 // check for functions 2056 if (inputType === 0) { 2057 dX = []; 2058 dY = []; 2059 2060 for (i = 0; i < len; i++) { 2061 if (Type.isFunction(dataX[i])) { 2062 dX.push(dataX[i]()); 2063 } else { 2064 dX.push(dataX[i]); 2065 } 2066 2067 if (Type.isFunction(dataY[i])) { 2068 dY.push(dataY[i]()); 2069 } else { 2070 dY.push(dataY[i]); 2071 } 2072 } 2073 } 2074 2075 M = []; 2076 2077 for (j = 0; j < len; j++) { 2078 M.push([1]); 2079 } 2080 2081 for (i = 1; i <= d; i++) { 2082 for (j = 0; j < len; j++) { 2083 M[j][i] = M[j][i - 1] * dX[j]; 2084 } 2085 } 2086 2087 y = dY; 2088 MT = Mat.transpose(M); 2089 B = Mat.matMatMult(MT, M); 2090 c = Mat.matVecMult(MT, y); 2091 coeffs = Mat.Numerics.Gauss(B, c); 2092 term = Mat.Numerics.generatePolynomialTerm(coeffs, d, 'x', 3); 2093 } 2094 2095 // Horner's scheme to evaluate polynomial 2096 s = coeffs[d]; 2097 2098 for (i = d - 1; i >= 0; i--) { 2099 s = (s * x + coeffs[i]); 2100 } 2101 2102 return s; 2103 }; 2104 2105 fct.getTerm = function () { 2106 return term; 2107 }; 2108 2109 return fct; 2110 }, 2111 2112 /** 2113 * Computes the cubic Bezier curve through a given set of points. 2114 * @param {Array} points Array consisting of 3*k+1 {@link JXG.Points}. 2115 * The points at position k with k mod 3 = 0 are the data points, 2116 * points at position k with k mod 3 = 1 or 2 are the control points. 2117 * @returns {Array} An array consisting of two functions of one parameter t which return the 2118 * x resp. y coordinates of the Bezier curve in t, one zero value, and a third function accepting 2119 * no parameters and returning one third of the length of the points. 2120 * @memberof JXG.Math.Numerics 2121 */ 2122 bezier: function (points) { 2123 var len, flen, 2124 /** @ignore */ 2125 makeFct = function (which) { 2126 return function (t, suspendedUpdate) { 2127 var z = Math.floor(t) * 3, 2128 t0 = t % 1, 2129 t1 = 1 - t0; 2130 2131 if (!suspendedUpdate) { 2132 flen = 3 * Math.floor((points.length - 1) / 3); 2133 len = Math.floor(flen / 3); 2134 } 2135 2136 if (t < 0) { 2137 return points[0][which](); 2138 } 2139 2140 if (t >= len) { 2141 return points[flen][which](); 2142 } 2143 2144 if (isNaN(t)) { 2145 return NaN; 2146 } 2147 2148 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; 2149 }; 2150 }; 2151 2152 return [makeFct('X'), makeFct('Y'), 0, 2153 function () { 2154 return Math.floor(points.length / 3); 2155 }]; 2156 }, 2157 2158 /** 2159 * Computes the B-spline curve of order k (order = degree+1) through a given set of points. 2160 * @param {Array} points Array consisting of JXG.Points. 2161 * @param {Number} order Order of the B-spline curve. 2162 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 2163 * which return the x resp. y coordinates of the B-spline curve in t, a zero value, and a function simply 2164 * returning the length of the points array minus one. 2165 * @memberof JXG.Math.Numerics 2166 */ 2167 bspline: function (points, order) { 2168 var knots, N = [], 2169 _knotVector = function (n, k) { 2170 var j, 2171 kn = []; 2172 2173 for (j = 0; j < n + k + 1; j++) { 2174 if (j < k) { 2175 kn[j] = 0.0; 2176 } else if (j <= n) { 2177 kn[j] = j - k + 1; 2178 } else { 2179 kn[j] = n - k + 2; 2180 } 2181 } 2182 2183 return kn; 2184 }, 2185 2186 _evalBasisFuncs = function (t, kn, n, k, s) { 2187 var i, j, a, b, den, 2188 N = []; 2189 2190 if (kn[s] <= t && t < kn[s + 1]) { 2191 N[s] = 1; 2192 } else { 2193 N[s] = 0; 2194 } 2195 2196 for (i = 2; i <= k; i++) { 2197 for (j = s - i + 1; j <= s; j++) { 2198 if (j <= s - i + 1 || j < 0) { 2199 a = 0.0; 2200 } else { 2201 a = N[j]; 2202 } 2203 2204 if (j >= s) { 2205 b = 0.0; 2206 } else { 2207 b = N[j + 1]; 2208 } 2209 2210 den = kn[j + i - 1] - kn[j]; 2211 2212 if (den === 0) { 2213 N[j] = 0; 2214 } else { 2215 N[j] = (t - kn[j]) / den * a; 2216 } 2217 2218 den = kn[j + i] - kn[j + 1]; 2219 2220 if (den !== 0) { 2221 N[j] += (kn[j + i] - t) / den * b; 2222 } 2223 } 2224 } 2225 return N; 2226 }, 2227 /** @ignore */ 2228 makeFct = function (which) { 2229 return function (t, suspendedUpdate) { 2230 var y, j, s, 2231 len = points.length, 2232 n = len - 1, 2233 k = order; 2234 2235 if (n <= 0) { 2236 return NaN; 2237 } 2238 2239 if (n + 2 <= k) { 2240 k = n + 1; 2241 } 2242 2243 if (t <= 0) { 2244 return points[0][which](); 2245 } 2246 2247 if (t >= n - k + 2) { 2248 return points[n][which](); 2249 } 2250 2251 s = Math.floor(t) + k - 1; 2252 knots = _knotVector(n, k); 2253 N = _evalBasisFuncs(t, knots, n, k, s); 2254 2255 y = 0.0; 2256 for (j = s - k + 1; j <= s; j++) { 2257 if (j < len && j >= 0) { 2258 y += points[j][which]() * N[j]; 2259 } 2260 } 2261 2262 return y; 2263 }; 2264 }; 2265 2266 return [makeFct('X'), makeFct('Y'), 0, 2267 function () { 2268 return points.length - 1; 2269 }]; 2270 }, 2271 2272 /** 2273 * Numerical (symmetric) approximation of derivative. suspendUpdate is piped through, 2274 * see {@link JXG.Curve#updateCurve} 2275 * and {@link JXG.Curve#hasPoint}. 2276 * @param {function} f Function in one variable to be differentiated. 2277 * @param {object} [obj] Optional object that is treated as "this" in the function body. This is useful, if the function is a 2278 * method of an object and contains a reference to its parent object via "this". 2279 * @returns {function} Derivative function of a given function f. 2280 * @memberof JXG.Math.Numerics 2281 */ 2282 D: function (f, obj) { 2283 if (!Type.exists(obj)) { 2284 return function (x, suspendUpdate) { 2285 var h = 0.00001, 2286 h2 = (h * 2.0); 2287 2288 // Experiments with Richardsons rule 2289 /* 2290 var phi = (f(x + h, suspendUpdate) - f(x - h, suspendUpdate)) / h2; 2291 var phi2; 2292 h *= 0.5; 2293 h2 *= 0.5; 2294 phi2 = (f(x + h, suspendUpdate) - f(x - h, suspendUpdate)) / h2; 2295 2296 return phi2 + (phi2 - phi) / 3.0; 2297 */ 2298 return (f(x + h, suspendUpdate) - f(x - h, suspendUpdate)) / h2; 2299 }; 2300 } 2301 2302 return function (x, suspendUpdate) { 2303 var h = 0.00001, 2304 h2 = (h * 2.0); 2305 2306 return (f.apply(obj, [x + h, suspendUpdate]) - f.apply(obj, [x - h, suspendUpdate])) / h2; 2307 }; 2308 }, 2309 2310 /** 2311 * Evaluate the function term for {@see #riemann}. 2312 * @private 2313 * @param {Number} x function argument 2314 * @param {function} f JavaScript function returning a number 2315 * @param {String} type Name of the Riemann sum type, e.g. 'lower', see {@see #riemann}. 2316 * @param {Number} delta Width of the bars in user coordinates 2317 * @returns {Number} Upper (delta > 0) or lower (delta < 0) value of the bar containing x of the Riemann sum. 2318 * 2319 * @memberof JXG.Math.Numerics 2320 */ 2321 _riemannValue: function (x, f, type, delta) { 2322 var y, y1, x1, delta1; 2323 2324 if (delta < 0) { // delta is negative if the lower function term is evaluated 2325 if (type !== 'trapezoidal') { 2326 x = x + delta; 2327 } 2328 delta *= -1; 2329 if (type === 'lower') { 2330 type = 'upper'; 2331 } else if (type === 'upper') { 2332 type = 'lower'; 2333 } 2334 } 2335 2336 delta1 = delta * 0.01; // for 'lower' and 'upper' 2337 2338 if (type === 'right') { 2339 y = f(x + delta); 2340 } else if (type === 'middle') { 2341 y = f(x + delta * 0.5); 2342 } else if (type === 'left' || type === 'trapezoidal') { 2343 y = f(x); 2344 } else if (type === 'lower') { 2345 y = f(x); 2346 2347 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 2348 y1 = f(x1); 2349 2350 if (y1 < y) { 2351 y = y1; 2352 } 2353 } 2354 2355 y1 = f(x + delta); 2356 if (y1 < y) { 2357 y = y1; 2358 } 2359 } else if (type === 'upper') { 2360 y = f(x); 2361 2362 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 2363 y1 = f(x1); 2364 if (y1 > y) { 2365 y = y1; 2366 } 2367 } 2368 2369 y1 = f(x + delta); 2370 if (y1 > y) { 2371 y = y1; 2372 } 2373 } else if (type === 'random') { 2374 y = f(x + delta * Math.random()); 2375 } else if (type === 'simpson') { 2376 y = (f(x) + 4 * f(x + delta * 0.5) + f(x + delta)) / 6.0; 2377 } else { 2378 y = f(x); // default is lower 2379 } 2380 2381 return y; 2382 }, 2383 2384 /** 2385 * Helper function to create curve which displays Riemann sums. 2386 * Compute coordinates for the rectangles showing the Riemann sum. 2387 * @param {Function,Array} f Function or array of two functions. 2388 * If f is a function the integral of this function is approximated by the Riemann sum. 2389 * If f is an array consisting of two functions the area between the two functions is filled 2390 * by the Riemann sum bars. 2391 * @param {Number} n number of rectangles. 2392 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson', or 'trapezoidal'. 2393 * @param {Number} start Left border of the approximation interval 2394 * @param {Number} end Right border of the approximation interval 2395 * @returns {Array} An array of two arrays containing the x and y coordinates for the rectangles showing the Riemann sum. This 2396 * 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 2397 * rectangles. 2398 * @memberof JXG.Math.Numerics 2399 */ 2400 riemann: function (gf, n, type, start, end) { 2401 var i, delta, 2402 xarr = [], 2403 yarr = [], 2404 j = 0, 2405 x = start, y, 2406 sum = 0, 2407 f, g, 2408 ylow, yup; 2409 2410 if (Type.isArray(gf)) { 2411 g = gf[0]; 2412 f = gf[1]; 2413 } else { 2414 f = gf; 2415 } 2416 2417 n = Math.floor(n); 2418 2419 if (n <= 0) { 2420 return [xarr, yarr, sum]; 2421 } 2422 2423 delta = (end - start) / n; 2424 2425 // Upper bar ends 2426 for (i = 0; i < n; i++) { 2427 y = this._riemannValue(x, f, type, delta); 2428 xarr[j] = x; 2429 yarr[j] = y; 2430 2431 j += 1; 2432 x += delta; 2433 if (type === 'trapezoidal') { 2434 y = f(x); 2435 } 2436 xarr[j] = x; 2437 yarr[j] = y; 2438 2439 j += 1; 2440 } 2441 2442 // Lower bar ends 2443 for (i = 0; i < n; i++) { 2444 if (g) { 2445 y = this._riemannValue(x, g, type, -delta); 2446 } else { 2447 y = 0.0; 2448 } 2449 xarr[j] = x; 2450 yarr[j] = y; 2451 2452 j += 1; 2453 x -= delta; 2454 if (type === 'trapezoidal' && g) { 2455 y = g(x); 2456 } 2457 xarr[j] = x; 2458 yarr[j] = y; 2459 2460 // Add the area of the bar to 'sum' 2461 if (type !== 'trapezoidal') { 2462 ylow = y; 2463 yup = yarr[2 * (n - 1) - 2 * i]; 2464 } else { 2465 yup = 0.5 * (f(x + delta) + f(x)); 2466 if (g) { 2467 ylow = 0.5 * (g(x + delta) + g(x)); 2468 } else { 2469 ylow = 0.0; 2470 } 2471 } 2472 sum += (yup - ylow) * delta; 2473 2474 // Draw the vertical lines 2475 j += 1; 2476 xarr[j] = x; 2477 yarr[j] = yarr[2 * (n - 1) - 2 * i]; 2478 2479 j += 1; 2480 } 2481 2482 return [xarr, yarr, sum]; 2483 }, 2484 2485 /** 2486 * Approximate the integral by Riemann sums. 2487 * Compute the area described by the riemann sum rectangles. 2488 * 2489 * If there is an element of type {@link Riemannsum}, then it is more efficient 2490 * to use the method JXG.Curve.Value() of this element instead. 2491 * 2492 * @param {Function_Array} f Function or array of two functions. 2493 * If f is a function the integral of this function is approximated by the Riemann sum. 2494 * If f is an array consisting of two functions the area between the two functions is approximated 2495 * by the Riemann sum. 2496 * @param {Number} n number of rectangles. 2497 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson' or 'trapezoidal'. 2498 * 2499 * @param {Number} start Left border of the approximation interval 2500 * @param {Number} end Right border of the approximation interval 2501 * @returns {Number} The sum of the areas of the rectangles. 2502 * @memberof JXG.Math.Numerics 2503 */ 2504 riemannsum: function (f, n, type, start, end) { 2505 JXG.deprecated('Numerics.riemannsum()', 'Numerics.riemann()'); 2506 return this.riemann(f, n, type, start, end)[2]; 2507 }, 2508 2509 /** 2510 * Solve initial value problems numerically using Runge-Kutta-methods. 2511 * See {@link http://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm. 2512 * @param {object,String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing 2513 * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure 2514 * <pre> 2515 * { 2516 * s: <Number>, 2517 * A: <matrix>, 2518 * b: <Array>, 2519 * c: <Array> 2520 * } 2521 * </pre> 2522 * 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 2523 * @param {Array} x0 Initial value vector. If the problem is of one-dimensional, the initial value also has to be given in an array. 2524 * @param {Array} I Interval on which to integrate. 2525 * @param {Number} N Number of evaluation points. 2526 * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode 2527 * 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 2528 * vector <tt>x</tt>, and has to return a vector of the same dimension as <tt>x</tt> has. 2529 * @returns {Array} An array of vectors describing the solution of the ode on the given interval I. 2530 * @example 2531 * // A very simple autonomous system dx(t)/dt = x(t); 2532 * function f(t, x) { 2533 * return x; 2534 * } 2535 * 2536 * // Solve it with initial value x(0) = 1 on the interval [0, 2] 2537 * // with 20 evaluation points. 2538 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 2539 * 2540 * // Prepare data for plotting the solution of the ode using a curve. 2541 * var dataX = []; 2542 * var dataY = []; 2543 * var h = 0.1; // (I[1] - I[0])/N = (2-0)/20 2544 * for(var i=0; i<data.length; i++) { 2545 * dataX[i] = i*h; 2546 * dataY[i] = data[i][0]; 2547 * } 2548 * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'}); 2549 * </pre><div class="jxgbox" id="d2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div> 2550 * <script type="text/javascript"> 2551 * var board = JXG.JSXGraph.initBoard('d2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false}); 2552 * function f(t, x) { 2553 * // we have to copy the value. 2554 * // return x; would just return the reference. 2555 * return [x[0]]; 2556 * } 2557 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 2558 * var dataX = []; 2559 * var dataY = []; 2560 * var h = 0.1; 2561 * for(var i=0; i<data.length; i++) { 2562 * dataX[i] = i*h; 2563 * dataY[i] = data[i][0]; 2564 * } 2565 * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'}); 2566 * </script><pre> 2567 * @memberof JXG.Math.Numerics 2568 */ 2569 rungeKutta: function (butcher, x0, I, N, f) { 2570 var e, i, j, k, l, s, 2571 x = [], 2572 y = [], 2573 h = (I[1] - I[0]) / N, 2574 t = I[0], 2575 dim = x0.length, 2576 result = [], 2577 r = 0; 2578 2579 if (Type.isString(butcher)) { 2580 butcher = predefinedButcher[butcher] || predefinedButcher.euler; 2581 } 2582 s = butcher.s; 2583 2584 // don't change x0, so copy it 2585 for (e = 0; e < dim; e++) { 2586 x[e] = x0[e]; 2587 } 2588 2589 for (i = 0; i < N; i++) { 2590 // Optimization doesn't work for ODEs plotted using time 2591 // if((i % quotient == 0) || (i == N-1)) { 2592 result[r] = []; 2593 for (e = 0; e < dim; e++) { 2594 result[r][e] = x[e]; 2595 } 2596 2597 r += 1; 2598 k = []; 2599 2600 for (j = 0; j < s; j++) { 2601 // init y = 0 2602 for (e = 0; e < dim; e++) { 2603 y[e] = 0.0; 2604 } 2605 2606 2607 // Calculate linear combination of former k's and save it in y 2608 for (l = 0; l < j; l++) { 2609 for (e = 0; e < dim; e++) { 2610 y[e] += (butcher.A[j][l]) * h * k[l][e]; 2611 } 2612 } 2613 2614 // add x(t) to y 2615 for (e = 0; e < dim; e++) { 2616 y[e] += x[e]; 2617 } 2618 2619 // calculate new k and add it to the k matrix 2620 k.push(f(t + butcher.c[j] * h, y)); 2621 } 2622 2623 // init y = 0 2624 for (e = 0; e < dim; e++) { 2625 y[e] = 0.0; 2626 } 2627 2628 for (l = 0; l < s; l++) { 2629 for (e = 0; e < dim; e++) { 2630 y[e] += butcher.b[l] * k[l][e]; 2631 } 2632 } 2633 2634 for (e = 0; e < dim; e++) { 2635 x[e] = x[e] + h * y[e]; 2636 } 2637 2638 t += h; 2639 } 2640 2641 return result; 2642 }, 2643 2644 /** 2645 * Maximum number of iterations in {@link JXG.Math.Numerics.fzero} 2646 * @type Number 2647 * @default 80 2648 * @memberof JXG.Math.Numerics 2649 */ 2650 maxIterationsRoot: 80, 2651 2652 /** 2653 * Maximum number of iterations in {@link JXG.Math.Numerics.fminbr} 2654 * @type Number 2655 * @default 500 2656 * @memberof JXG.Math.Numerics 2657 */ 2658 maxIterationsMinimize: 500, 2659 2660 /** 2661 * 2662 * Find zero of an univariate function f. 2663 * @param {function} f Function, whose root is to be found 2664 * @param {Array,Number} x0 Start value or start interval enclosing the root 2665 * @param {Object} object Parent object in case f is method of it 2666 * @returns {Number} the approximation of the root 2667 * Algorithm: 2668 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 2669 * computations. M., Mir, 1980, p.180 of the Russian edition 2670 * 2671 * If x0 is an array containing lower and upper bound for the zero 2672 * algorithm 748 is applied. Otherwise, if x0 is a number, 2673 * the algorithm tries to bracket a zero of f starting from x0. 2674 * If this fails, we fall back to Newton's method. 2675 * @memberof JXG.Math.Numerics 2676 */ 2677 fzero: function (f, x0, object) { 2678 var a, b, c, 2679 fa, fb, fc, 2680 aa, blist, i, len, u, fu, 2681 prev_step, t1, cb, t2, 2682 // Actual tolerance 2683 tol_act, 2684 // Interpolation step is calculated in the form p/q; division 2685 // operations is delayed until the last moment 2686 p, q, 2687 // Step at this iteration 2688 new_step, 2689 eps = Mat.eps, 2690 maxiter = this.maxIterationsRoot, 2691 niter = 0, 2692 nfev = 0; 2693 2694 if (Type.isArray(x0)) { 2695 if (x0.length < 2) { 2696 throw new Error("JXG.Math.Numerics.fzero: length of array x0 has to be at least two."); 2697 } 2698 2699 a = x0[0]; 2700 fa = f.call(object, a); 2701 nfev += 1; 2702 b = x0[1]; 2703 fb = f.call(object, b); 2704 nfev += 1; 2705 } else { 2706 a = x0; 2707 fa = f.call(object, a); 2708 nfev += 1; 2709 2710 // Try to get b. 2711 if (a === 0) { 2712 aa = 1; 2713 } else { 2714 aa = a; 2715 } 2716 2717 blist = [0.9 * aa, 1.1 * aa, aa - 1, aa + 1, 0.5 * aa, 1.5 * aa, -aa, 2 * aa, -10 * aa, 10 * aa]; 2718 len = blist.length; 2719 2720 for (i = 0; i < len; i++) { 2721 b = blist[i]; 2722 fb = f.call(object, b); 2723 nfev += 1; 2724 2725 if (fa * fb <= 0) { 2726 break; 2727 } 2728 } 2729 if (b < a) { 2730 u = a; 2731 a = b; 2732 b = u; 2733 2734 fu = fa; 2735 fa = fb; 2736 fb = fu; 2737 } 2738 } 2739 2740 if (fa * fb > 0) { 2741 // Bracketing not successful, fall back to Newton's method or to fminbr 2742 if (Type.isArray(x0)) { 2743 return this.fminbr(f, [a, b], object); 2744 } 2745 2746 return this.Newton(f, a, object); 2747 } 2748 2749 // OK, we have enclosed a zero of f. 2750 // Now we can start Brent's method 2751 2752 c = a; 2753 fc = fa; 2754 2755 // Main iteration loop 2756 while (niter < maxiter) { 2757 // Distance from the last but one to the last approximation 2758 prev_step = b - a; 2759 2760 // Swap data for b to be the best approximation 2761 if (Math.abs(fc) < Math.abs(fb)) { 2762 a = b; 2763 b = c; 2764 c = a; 2765 2766 fa = fb; 2767 fb = fc; 2768 fc = fa; 2769 } 2770 tol_act = 2 * eps * Math.abs(b) + eps * 0.5; 2771 new_step = (c - b) * 0.5; 2772 2773 if (Math.abs(new_step) <= tol_act && Math.abs(fb) <= eps) { 2774 // Acceptable approx. is found 2775 return b; 2776 } 2777 2778 // Decide if the interpolation can be tried 2779 // If prev_step was large enough and was in true direction Interpolatiom may be tried 2780 if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) { 2781 cb = c - b; 2782 2783 // If we have only two distinct points linear interpolation can only be applied 2784 if (a === c) { 2785 t1 = fb / fa; 2786 p = cb * t1; 2787 q = 1.0 - t1; 2788 // Quadric inverse interpolation 2789 } else { 2790 q = fa / fc; 2791 t1 = fb / fc; 2792 t2 = fb / fa; 2793 2794 p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0)); 2795 q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0); 2796 } 2797 2798 // p was calculated with the opposite sign; make p positive 2799 if (p > 0) { 2800 q = -q; 2801 // and assign possible minus to q 2802 } else { 2803 p = -p; 2804 } 2805 2806 // If b+p/q falls in [b,c] and isn't too large it is accepted 2807 // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent 2808 if (p < (0.75 * cb * q - Math.abs(tol_act * q) * 0.5) && 2809 p < Math.abs(prev_step * q * 0.5)) { 2810 new_step = p / q; 2811 } 2812 } 2813 2814 // Adjust the step to be not less than tolerance 2815 if (Math.abs(new_step) < tol_act) { 2816 if (new_step > 0) { 2817 new_step = tol_act; 2818 } else { 2819 new_step = -tol_act; 2820 } 2821 } 2822 2823 // Save the previous approx. 2824 a = b; 2825 fa = fb; 2826 b += new_step; 2827 fb = f.call(object, b); 2828 // Do step to a new approxim. 2829 nfev += 1; 2830 2831 // Adjust c for it to have a sign opposite to that of b 2832 if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) { 2833 c = a; 2834 fc = fa; 2835 } 2836 niter++; 2837 } // End while 2838 2839 return b; 2840 }, 2841 2842 /** 2843 * 2844 * Find minimum of an univariate function f. 2845 * <p> 2846 * Algorithm: 2847 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 2848 * computations. M., Mir, 1980, p.180 of the Russian edition 2849 * 2850 * @param {function} f Function, whose minimum is to be found 2851 * @param {Array} x0 Start interval enclosing the minimum 2852 * @param {Object} context Parent object in case f is method of it 2853 * @returns {Number} the approximation of the minimum value position 2854 * @memberof JXG.Math.Numerics 2855 **/ 2856 fminbr: function (f, x0, context) { 2857 var a, b, x, v, w, 2858 fx, fv, fw, 2859 range, middle_range, tol_act, new_step, 2860 p, q, t, ft, 2861 // Golden section ratio 2862 r = (3.0 - Math.sqrt(5.0)) * 0.5, 2863 tol = Mat.eps, 2864 sqrteps = Mat.eps, //Math.sqrt(Mat.eps), 2865 maxiter = this.maxIterationsMinimize, 2866 niter = 0, 2867 nfev = 0; 2868 2869 if (!Type.isArray(x0) || x0.length < 2) { 2870 throw new Error("JXG.Math.Numerics.fminbr: length of array x0 has to be at least two."); 2871 } 2872 2873 a = x0[0]; 2874 b = x0[1]; 2875 v = a + r * (b - a); 2876 fv = f.call(context, v); 2877 2878 // First step - always gold section 2879 nfev += 1; 2880 x = v; 2881 w = v; 2882 fx = fv; 2883 fw = fv; 2884 2885 while (niter < maxiter) { 2886 // Range over the interval in which we are looking for the minimum 2887 range = b - a; 2888 middle_range = (a + b) * 0.5; 2889 2890 // Actual tolerance 2891 tol_act = sqrteps * Math.abs(x) + tol / 3.0; 2892 2893 if (Math.abs(x - middle_range) + range * 0.5 <= 2.0 * tol_act) { 2894 // Acceptable approx. is found 2895 return x; 2896 } 2897 2898 // Obtain the golden section step 2899 new_step = r * (x < middle_range ? b - x : a - x); 2900 2901 // Decide if the interpolation can be tried. If x and w are distinct interpolatiom may be tried 2902 if (Math.abs(x - w) >= tol_act) { 2903 // Interpolation step is calculated as p/q; 2904 // division operation is delayed until last moment 2905 t = (x - w) * (fx - fv); 2906 q = (x - v) * (fx - fw); 2907 p = (x - v) * q - (x - w) * t; 2908 q = 2 * (q - t); 2909 2910 if (q > 0) { // q was calculated with the op- 2911 p = -p; // posite sign; make q positive 2912 } else { // and assign possible minus to 2913 q = -q; // p 2914 } 2915 if (Math.abs(p) < Math.abs(new_step * q) && // If x+p/q falls in [a,b] 2916 p > q * (a - x + 2 * tol_act) && // not too close to a and 2917 p < q * (b - x - 2 * tol_act)) { // b, and isn't too large 2918 new_step = p / q; // it is accepted 2919 } 2920 // If p/q is too large then the 2921 // golden section procedure can 2922 // reduce [a,b] range to more 2923 // extent 2924 } 2925 2926 // Adjust the step to be not less than tolerance 2927 if (Math.abs(new_step) < tol_act) { 2928 if (new_step > 0) { 2929 new_step = tol_act; 2930 } else { 2931 new_step = -tol_act; 2932 } 2933 } 2934 2935 // Obtain the next approximation to min 2936 // and reduce the enveloping range 2937 2938 // Tentative point for the min 2939 t = x + new_step; 2940 ft = f.call(context, t); 2941 nfev += 1; 2942 2943 // t is a better approximation 2944 if (ft <= fx) { 2945 // Reduce the range so that t would fall within it 2946 if (t < x) { 2947 b = x; 2948 } else { 2949 a = x; 2950 } 2951 2952 // Assign the best approx to x 2953 v = w; 2954 w = x; 2955 x = t; 2956 2957 fv = fw; 2958 fw = fx; 2959 fx = ft; 2960 // x remains the better approx 2961 } else { 2962 // Reduce the range enclosing x 2963 if (t < x) { 2964 a = t; 2965 } else { 2966 b = t; 2967 } 2968 2969 if (ft <= fw || w === x) { 2970 v = w; 2971 w = t; 2972 fv = fw; 2973 fw = ft; 2974 } else if (ft <= fv || v === x || v === w) { 2975 v = t; 2976 fv = ft; 2977 } 2978 } 2979 niter += 1; 2980 } 2981 2982 return x; 2983 }, 2984 2985 /** 2986 * Implements the Ramer-Douglas-Peucker algorithm. 2987 * It discards points which are not necessary from the polygonal line defined by the point array 2988 * pts. The computation is done in screen coordinates. 2989 * Average runtime is O(nlog(n)), worst case runtime is O(n^2), where n is the number of points. 2990 * @param {Array} pts Array of {@link JXG.Coords} 2991 * @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>. 2992 * @returns {Array} An array containing points which represent an apparently identical curve as the points of pts do, but contains fewer points. 2993 * @memberof JXG.Math.Numerics 2994 */ 2995 RamerDouglasPeucker: function (pts, eps) { 2996 var newPts = [], i, k, len, 2997 2998 /** 2999 * findSplit() is a subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}. 3000 * It searches for the point between index i and j which 3001 * has the largest distance from the line between the points i and j. 3002 * @param {Array} pts Array of {@link JXG.Coords} 3003 * @param {Number} i Index of a point in pts 3004 * @param {Number} j Index of a point in pts 3005 * @ignore 3006 * @private 3007 */ 3008 findSplit = function (pts, i, j) { 3009 var d, k, ci, cj, ck, 3010 x0, y0, x1, y1, 3011 den, lbda, 3012 dist = 0, 3013 f = i; 3014 3015 if (j - i < 2) { 3016 return [-1.0, 0]; 3017 } 3018 3019 ci = pts[i].scrCoords; 3020 cj = pts[j].scrCoords; 3021 3022 if (isNaN(ci[1] + ci[2])) { 3023 return [NaN, i]; 3024 } 3025 if (isNaN(cj[1] + cj[2])) { 3026 return [NaN, j]; 3027 } 3028 3029 for (k = i + 1; k < j; k++) { 3030 ck = pts[k].scrCoords; 3031 if (isNaN(ck[1] + ck[2])) { 3032 return [NaN, k]; 3033 } 3034 3035 x0 = ck[1] - ci[1]; 3036 y0 = ck[2] - ci[2]; 3037 x1 = cj[1] - ci[1]; 3038 y1 = cj[2] - ci[2]; 3039 den = x1 * x1 + y1 * y1; 3040 3041 if (den >= Mat.eps) { 3042 lbda = (x0 * x1 + y0 * y1) / den; 3043 3044 if (lbda < 0.0) { 3045 lbda = 0.0; 3046 } else if (lbda > 1.0) { 3047 lbda = 1.0; 3048 } 3049 3050 x0 = x0 - lbda * x1; 3051 y0 = y0 - lbda * y1; 3052 d = x0 * x0 + y0 * y0; 3053 } else { 3054 lbda = 0.0; 3055 d = x0 * x0 + y0 * y0; 3056 } 3057 3058 if (d > dist) { 3059 dist = d; 3060 f = k; 3061 } 3062 } 3063 return [Math.sqrt(dist), f]; 3064 }, 3065 3066 /** 3067 * RDP() is a private subroutine of {@link JXG.Math.Numerics.RamerDouglasPeucker}. 3068 * It runs recursively through the point set and searches the 3069 * point which has the largest distance from the line between the first point and 3070 * the last point. If the distance from the line is greater than eps, this point is 3071 * included in our new point set otherwise it is discarded. 3072 * If it is taken, we recursively apply the subroutine to the point set before 3073 * and after the chosen point. 3074 * @param {Array} pts Array of {@link JXG.Coords} 3075 * @param {Number} i Index of an element of pts 3076 * @param {Number} j Index of an element of pts 3077 * @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>. 3078 * @param {Array} newPts Array of {@link JXG.Coords} 3079 * @ignore 3080 * @private 3081 */ 3082 RDP = function (pts, i, j, eps, newPts) { 3083 var result = findSplit(pts, i, j), 3084 k = result[1]; 3085 3086 if (isNaN(result[0])) { 3087 RDP(pts, i, k - 1, eps, newPts); 3088 newPts.push(pts[k]); 3089 do { 3090 ++k; 3091 } while (k <= j && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])); 3092 if (k <= j) { 3093 newPts.push(pts[k]); 3094 } 3095 RDP(pts, k + 1, j, eps, newPts); 3096 } else if (result[0] > eps) { 3097 RDP(pts, i, k, eps, newPts); 3098 RDP(pts, k, j, eps, newPts); 3099 } else { 3100 newPts.push(pts[j]); 3101 } 3102 }; 3103 3104 len = pts.length; 3105 3106 // Search for the left most point woithout NaN coordinates 3107 i = 0; 3108 while (i < len && isNaN(pts[i].scrCoords[1] + pts[i].scrCoords[2])) { 3109 i += 1; 3110 } 3111 // Search for the right most point woithout NaN coordinates 3112 k = len - 1; 3113 while (k > i && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])) { 3114 k -= 1; 3115 } 3116 3117 // Only proceed if something is left 3118 if (!(i > k || i === len)) { 3119 newPts[0] = pts[i]; 3120 RDP(pts, i, k, eps, newPts); 3121 } 3122 3123 return newPts; 3124 }, 3125 3126 /** 3127 * Old name for the implementation of the Ramer-Douglas-Peucker algorithm. 3128 * @deprecated Use {@link JXG.Math.Numerics.RamerDouglasPeucker} 3129 * @memberof JXG.Math.Numerics 3130 */ 3131 RamerDouglasPeuker: function (pts, eps) { 3132 JXG.deprecated('Numerics.RamerDouglasPeuker()', 'Numerics.RamerDouglasPeucker()'); 3133 return this.RamerDouglasPeucker(pts, eps); 3134 }, 3135 3136 /** 3137 * Implements the Visvalingam-Whyatt algorithm. 3138 * See M. Visvalingam, J. D. Whyatt: 3139 * "Line generalisation by repeated elimination of the smallest area", C.I.S.R.G Discussion paper 10, July 1992 3140 * 3141 * The algorithm discards points which are not necessary from the polygonal line defined by the point array 3142 * pts (consisting of type JXG.Coords). 3143 * @param {Array} pts Array of {@link JXG.Coords} 3144 * @param {Number} numPoints Number of remaining intermediate points. The first and the last point of the original points will 3145 * be taken in any case. 3146 * @returns {Array} An array containing points which approximates the curve defined by pts. 3147 * @memberof JXG.Math.Numerics 3148 * 3149 * @example 3150 * var i, p = []; 3151 * for (i = 0; i < 5; ++i) { 3152 * p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6])); 3153 * } 3154 * 3155 * // Plot a cardinal spline curve 3156 * var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5); 3157 * var cu1 = board.create('curve', splineArr, {strokeColor: 'green'}); 3158 * 3159 * var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'}); 3160 * c.updateDataArray = function() { 3161 * var i, len, points; 3162 * 3163 * // Reduce number of intermediate points with Visvakingam-Whyatt to 6 3164 * points = JXG.Math.Numerics.Visvalingam(cu1.points, 6); 3165 * // Plot the remaining points 3166 * len = points.length; 3167 * this.dataX = []; 3168 * this.dataY = []; 3169 * for (i = 0; i < len; i++) { 3170 * this.dataX.push(points[i].usrCoords[1]); 3171 * this.dataY.push(points[i].usrCoords[2]); 3172 * } 3173 * }; 3174 * board.update(); 3175 * 3176 * </pre><div id="ce0cc55c-b592-11e6-8270-104a7d3be7eb" class="jxgbox" style="width: 300px; height: 300px;"></div> 3177 * <script type="text/javascript"> 3178 * (function() { 3179 * var board = JXG.JSXGraph.initBoard('ce0cc55c-b592-11e6-8270-104a7d3be7eb', 3180 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 3181 * 3182 * var i, p = []; 3183 * for (i = 0; i < 5; ++i) { 3184 * p.push(board.create('point', [Math.random() * 12 - 6, Math.random() * 12 - 6])); 3185 * } 3186 * 3187 * // Plot a cardinal spline curve 3188 * var splineArr = JXG.Math.Numerics.CardinalSpline(p, 0.5); 3189 * var cu1 = board.create('curve', splineArr, {strokeColor: 'green'}); 3190 * 3191 * var c = board.create('curve', [[0],[0]], {strokeWidth: 2, strokeColor: 'black'}); 3192 * c.updateDataArray = function() { 3193 * var i, len, points; 3194 * 3195 * // Reduce number of intermediate points with Visvakingam-Whyatt to 6 3196 * points = JXG.Math.Numerics.Visvalingam(cu1.points, 6); 3197 * // Plot the remaining points 3198 * len = points.length; 3199 * this.dataX = []; 3200 * this.dataY = []; 3201 * for (i = 0; i < len; i++) { 3202 * this.dataX.push(points[i].usrCoords[1]); 3203 * this.dataY.push(points[i].usrCoords[2]); 3204 * } 3205 * }; 3206 * board.update(); 3207 * 3208 * })(); 3209 * 3210 * </script><pre> 3211 * 3212 */ 3213 Visvalingam: function(pts, numPoints) { 3214 var i, len, vol, lastVol, 3215 linkedList = [], 3216 heap = [], 3217 points = [], 3218 lft, rt, lft2, rt2, 3219 obj; 3220 3221 len = pts.length; 3222 3223 // At least one intermediate point is needed 3224 if (len <= 2) { 3225 return pts; 3226 } 3227 3228 // Fill the linked list 3229 // Add first point to the linked list 3230 linkedList[0] = { 3231 used: true, 3232 lft: null, 3233 node: null 3234 }; 3235 3236 // Add all intermediate points to the linked list, 3237 // whose triangle area is nonzero. 3238 lft = 0; 3239 for (i = 1; i < len - 1; i++) { 3240 vol = Math.abs(JXG.Math.Numerics.det([pts[i - 1].usrCoords, 3241 pts[i].usrCoords, 3242 pts[i + 1].usrCoords])); 3243 if (!isNaN(vol)) { 3244 obj = { 3245 v: vol, 3246 idx: i 3247 }; 3248 heap.push(obj); 3249 linkedList[i] = { 3250 used: true, 3251 lft: lft, 3252 node: obj 3253 }; 3254 linkedList[lft].rt = i; 3255 lft = i; 3256 } 3257 } 3258 3259 // Add last point to the linked list 3260 linkedList[len - 1] = { 3261 used: true, 3262 rt: null, 3263 lft: lft, 3264 node: null 3265 }; 3266 linkedList[lft].rt = len - 1; 3267 3268 // Remove points until only numPoints intermediate points remain 3269 lastVol = -Infinity; 3270 while (heap.length > numPoints) { 3271 // Sort the heap with the updated volume values 3272 heap.sort(function(a, b) { 3273 // descending sort 3274 return b.v - a.v; 3275 }); 3276 3277 // Remove the point with the smallest triangle 3278 i = heap.pop().idx; 3279 linkedList[i].used = false; 3280 lastVol = linkedList[i].node.v; 3281 3282 // Update the pointers of the linked list 3283 lft = linkedList[i].lft; 3284 rt = linkedList[i].rt; 3285 linkedList[lft].rt = rt; 3286 linkedList[rt].lft = lft; 3287 3288 // Update the values for the volumes in the linked list 3289 lft2 = linkedList[lft].lft; 3290 if (lft2 !== null) { 3291 vol = Math.abs(JXG.Math.Numerics.det( 3292 [pts[lft2].usrCoords, 3293 pts[lft].usrCoords, 3294 pts[rt].usrCoords])); 3295 3296 linkedList[lft].node.v = (vol >= lastVol) ? vol : lastVol; 3297 } 3298 rt2 = linkedList[rt].rt; 3299 if (rt2 !== null) { 3300 vol = Math.abs(JXG.Math.Numerics.det( 3301 [pts[lft].usrCoords, 3302 pts[rt].usrCoords, 3303 pts[rt2].usrCoords])); 3304 3305 linkedList[rt].node.v = (vol >= lastVol) ? vol : lastVol; 3306 } 3307 } 3308 3309 // Return an array with the remaining points 3310 i = 0; 3311 points = [pts[i]]; 3312 do { 3313 i = linkedList[i].rt; 3314 points.push(pts[i]); 3315 } while (linkedList[i].rt !== null); 3316 3317 return points; 3318 } 3319 3320 }; 3321 3322 return Mat.Numerics; 3323 }); 3324