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