1 /* 2 Copyright 2008-2015 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 33 /*global JXG: true, define: true*/ 34 /*jslint nomen: true, plusplus: true*/ 35 36 /* depends: 37 utils/type 38 math/math 39 */ 40 41 /** 42 * @fileoverview In this file the namespace Math.Numerics is defined, which holds numerical 43 * algorithms for solving linear equations etc. 44 */ 45 46 define(['utils/type', 'math/math'], function (Type, Mat) { 47 48 "use strict"; 49 50 // Predefined butcher tableaus for the common Runge-Kutta method (fourth order), Heun method (second order), and Euler method (first order). 51 var predefinedButcher = { 52 rk4: { 53 s: 4, 54 A: [ 55 [ 0, 0, 0, 0], 56 [0.5, 0, 0, 0], 57 [ 0, 0.5, 0, 0], 58 [ 0, 0, 1, 0] 59 ], 60 b: [1.0 / 6.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 6.0], 61 c: [0, 0.5, 0.5, 1] 62 }, 63 heun: { 64 s: 2, 65 A: [ 66 [0, 0], 67 [1, 0] 68 ], 69 b: [0.5, 0.5], 70 c: [0, 1] 71 }, 72 euler: { 73 s: 1, 74 A: [ 75 [0] 76 ], 77 b: [1], 78 c: [0] 79 } 80 }; 81 82 /** 83 * The JXG.Math.Numerics namespace holds numerical algorithms, constants, and variables. 84 * @namespace 85 * @name JXG.Math.Numerics 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, intergration_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 && typeof config.number_of_nodes === 'number' ? 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 * Integral of function f over interval. 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 * @returns {Number} The value of the integral of f over interval 486 * @see JXG.Math.Numerics.NewtonCotes 487 * @memberof JXG.Math.Numerics 488 */ 489 I: function (interval, f) { 490 return this.NewtonCotes(interval, f, {number_of_nodes: 16, integration_type: 'milne'}); 491 }, 492 493 /** 494 * Newton's method to find roots of a funtion in one variable. 495 * @param {function} f We search for a solution of f(x)=0. 496 * @param {Number} x initial guess for the root, i.e. start value. 497 * @param {Object} context optional object that is treated as "this" in the function body. This is useful if 498 * the function is a method of an object and contains a reference to its parent object via "this". 499 * @returns {Number} A root of the function f. 500 * @memberof JXG.Math.Numerics 501 */ 502 Newton: function (f, x, context) { 503 var df, 504 i = 0, 505 h = Mat.eps, 506 newf = f.apply(context, [x]), 507 nfev = 1; 508 509 // For compatibility 510 if (Type.isArray(x)) { 511 x = x[0]; 512 } 513 514 while (i < 50 && Math.abs(newf) > h) { 515 df = this.D(f, context)(x); 516 nfev += 2; 517 518 if (Math.abs(df) > h) { 519 x -= newf / df; 520 } else { 521 x += (Math.random() * 0.2 - 1.0); 522 } 523 524 newf = f.apply(context, [x]); 525 nfev += 1; 526 i += 1; 527 } 528 529 return x; 530 }, 531 532 /** 533 * Abstract method to find roots of univariate functions. 534 * @param {function} f We search for a solution of f(x)=0. 535 * @param {Number} x initial guess for the root, i.e. starting value. 536 * @param {Object} context optional object that is treated as "this" in the function body. This is useful if 537 * the function is a method of an object and contains a reference to its parent object via "this". 538 * @returns {Number} A root of the function f. 539 * @memberof JXG.Math.Numerics 540 */ 541 root: function (f, x, context) { 542 return this.fzero(f, x, context); 543 }, 544 545 /** 546 * Compute an intersection of the curves c1 and c2 547 * with a generalized Newton method. 548 * We want to find values t1, t2 such that 549 * c1(t1) = c2(t2), i.e. 550 * (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2)) = (0,0). 551 * We set 552 * (e,f) := (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2)) 553 * 554 * The Jacobian J is defined by 555 * J = (a, b) 556 * (c, d) 557 * where 558 * a = c1_x'(t1) 559 * b = -c2_x'(t2) 560 * c = c1_y'(t1) 561 * d = -c2_y'(t2) 562 * 563 * The inverse J^(-1) of J is equal to 564 * (d, -b)/ 565 * (-c, a) / (ad-bc) 566 * 567 * Then, (t1new, t2new) := (t1,t2) - J^(-1)*(e,f). 568 * If the function meetCurveCurve possesses the properties 569 * t1memo and t2memo then these are taken as start values 570 * for the Newton algorithm. 571 * After stopping of the Newton algorithm the values of t1 and t2 are stored in 572 * t1memo and t2memo. 573 * 574 * @param {JXG.Curve} c1 Curve, Line or Circle 575 * @param {JXG.Curve} c2 Curve, Line or Circle 576 * @param {Number} t1ini start value for t1 577 * @param {Number} t2ini start value for t2 578 * @returns {JXG.Coords} intersection point 579 * @memberof JXG.Math.Numerics 580 */ 581 generalizedNewton: function (c1, c2, t1ini, t2ini) { 582 var t1, t2, 583 a, b, c, d, disc, 584 e, f, F, 585 D00, D01, 586 D10, D11, 587 count = 0; 588 589 if (this.generalizedNewton.t1memo) { 590 t1 = this.generalizedNewton.t1memo; 591 t2 = this.generalizedNewton.t2memo; 592 } else { 593 t1 = t1ini; 594 t2 = t2ini; 595 } 596 597 e = c1.X(t1) - c2.X(t2); 598 f = c1.Y(t1) - c2.Y(t2); 599 F = e * e + f * f; 600 601 D00 = this.D(c1.X, c1); 602 D01 = this.D(c2.X, c2); 603 D10 = this.D(c1.Y, c1); 604 D11 = this.D(c2.Y, c2); 605 606 while (F > Mat.eps && count < 10) { 607 a = D00(t1); 608 b = -D01(t2); 609 c = D10(t1); 610 d = -D11(t2); 611 disc = a * d - b * c; 612 t1 -= (d * e - b * f) / disc; 613 t2 -= (a * f - c * e) / disc; 614 e = c1.X(t1) - c2.X(t2); 615 f = c1.Y(t1) - c2.Y(t2); 616 F = e * e + f * f; 617 count += 1; 618 } 619 620 this.generalizedNewton.t1memo = t1; 621 this.generalizedNewton.t2memo = t2; 622 623 if (Math.abs(t1) < Math.abs(t2)) { 624 return [c1.X(t1), c1.Y(t1)]; 625 } 626 627 return [c2.X(t2), c2.Y(t2)]; 628 }, 629 630 /** 631 * Returns the Lagrange polynomials for curves with equidistant nodes, see 632 * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation, 633 * SIAM Review, Vol 46, No 3, (2004) 501-517. 634 * The graph of the parametric curve [x(t),y(t)] runs through the given points. 635 * @param {Array} p Array of JXG.Points 636 * @returns {Array} An array consisting of two functions x(t), y(t) which define a parametric curve 637 * f(t) = (x(t), y(t)) and two numbers x1 and x2 defining the curve's domain. x1 always equals zero. 638 * @memberof JXG.Math.Numerics 639 */ 640 Neville: function (p) { 641 var w = [], 642 /** @ignore */ 643 makeFct = function (fun) { 644 return function (t, suspendedUpdate) { 645 var i, d, s, 646 bin = Mat.binomial, 647 len = p.length, 648 len1 = len - 1, 649 num = 0.0, 650 denom = 0.0; 651 652 if (!suspendedUpdate) { 653 s = 1; 654 for (i = 0; i < len; i++) { 655 w[i] = bin(len1, i) * s; 656 s *= (-1); 657 } 658 } 659 660 d = t; 661 662 for (i = 0; i < len; i++) { 663 if (d === 0) { 664 return p[i][fun](); 665 } 666 s = w[i] / d; 667 d -= 1; 668 num += p[i][fun]() * s; 669 denom += s; 670 } 671 return num / denom; 672 }; 673 }, 674 675 xfct = makeFct('X'), 676 yfct = makeFct('Y'); 677 678 return [xfct, yfct, 0, function () { 679 return p.length - 1; 680 }]; 681 }, 682 683 /** 684 * Calculates second derivatives at the given knots. 685 * @param {Array} x x values of knots 686 * @param {Array} y y values of knots 687 * @returns {Array} Second derivatives of the interpolated function at the knots. 688 * @see #splineEval 689 * @memberof JXG.Math.Numerics 690 */ 691 splineDef: function (x, y) { 692 var pair, i, l, 693 n = Math.min(x.length, y.length), 694 diag = [], 695 z = [], 696 data = [], 697 dx = [], 698 delta = [], 699 F = []; 700 701 if (n === 2) { 702 return [0, 0]; 703 } 704 705 for (i = 0; i < n; i++) { 706 pair = {X: x[i], Y: y[i]}; 707 data.push(pair); 708 } 709 data.sort(function (a, b) { 710 return a.X - b.X; 711 }); 712 for (i = 0; i < n; i++) { 713 x[i] = data[i].X; 714 y[i] = data[i].Y; 715 } 716 717 for (i = 0; i < n - 1; i++) { 718 dx.push(x[i + 1] - x[i]); 719 } 720 for (i = 0; i < n - 2; i++) { 721 delta.push(6 * (y[i + 2] - y[i + 1]) / (dx[i + 1]) - 6 * (y[i + 1] - y[i]) / (dx[i])); 722 } 723 724 // ForwardSolve 725 diag.push(2 * (dx[0] + dx[1])); 726 z.push(delta[0]); 727 728 for (i = 0; i < n - 3; i++) { 729 l = dx[i + 1] / diag[i]; 730 diag.push(2 * (dx[i + 1] + dx[i + 2]) - l * dx[i + 1]); 731 z.push(delta[i + 1] - l * z[i]); 732 } 733 734 // BackwardSolve 735 F[n - 3] = z[n - 3] / diag[n - 3]; 736 for (i = n - 4; i >= 0; i--) { 737 F[i] = (z[i] - (dx[i + 1] * F[i + 1])) / diag[i]; 738 } 739 740 // Generate f''-Vector 741 for (i = n - 3; i >= 0; i--) { 742 F[i + 1] = F[i]; 743 } 744 745 // natural cubic spline 746 F[0] = 0; 747 F[n - 1] = 0; 748 749 return F; 750 }, 751 752 /** 753 * Evaluate points on spline. 754 * @param {Number,Array} x0 A single float value or an array of values to evaluate 755 * @param {Array} x x values of knots 756 * @param {Array} y y values of knots 757 * @param {Array} F Second derivatives at knots, calculated by {@link #splineDef} 758 * @see #splineDef 759 * @returns {Number,Array} A single value or an array, depending on what is given as x0. 760 * @memberof JXG.Math.Numerics 761 */ 762 splineEval: function (x0, x, y, F) { 763 var i, j, a, b, c, d, x_, 764 n = Math.min(x.length, y.length), 765 l = 1, 766 asArray = false, 767 y0 = []; 768 769 // number of points to be evaluated 770 if (Type.isArray(x0)) { 771 l = x0.length; 772 asArray = true; 773 } else { 774 x0 = [x0]; 775 } 776 777 for (i = 0; i < l; i++) { 778 // is x0 in defining interval? 779 if ((x0[i] < x[0]) || (x[i] > x[n - 1])) { 780 return NaN; 781 } 782 783 // determine part of spline in which x0 lies 784 for (j = 1; j < n; j++) { 785 if (x0[i] <= x[j]) { 786 break; 787 } 788 } 789 790 j -= 1; 791 792 // we're now in the j-th partial interval, i.e. x[j] < x0[i] <= x[j+1]; 793 // determine the coefficients of the polynomial in this interval 794 a = y[j]; 795 b = (y[j + 1] - y[j]) / (x[j + 1] - x[j]) - (x[j + 1] - x[j]) / 6 * (F[j + 1] + 2 * F[j]); 796 c = F[j] / 2; 797 d = (F[j + 1] - F[j]) / (6 * (x[j + 1] - x[j])); 798 // evaluate x0[i] 799 x_ = x0[i] - x[j]; 800 //y0.push(a + b*x_ + c*x_*x_ + d*x_*x_*x_); 801 y0.push(a + (b + (c + d * x_) * x_) * x_); 802 } 803 804 if (asArray) { 805 return y0; 806 } 807 808 return y0[0]; 809 }, 810 811 /** 812 * Generate a string containing the function term of a polynomial. 813 * @param {Array} coeffs Coefficients of the polynomial. The position i belongs to x^i. 814 * @param {Number} deg Degree of the polynomial 815 * @param {String} varname Name of the variable (usually 'x') 816 * @param {Number} prec Precision 817 * @returns {String} A string containg the function term of the polynomial. 818 * @memberof JXG.Math.Numerics 819 */ 820 generatePolynomialTerm: function (coeffs, deg, varname, prec) { 821 var i, t = []; 822 823 for (i = deg; i >= 0; i--) { 824 t = t.concat(['(', coeffs[i].toPrecision(prec), ')']); 825 826 if (i > 1) { 827 t = t.concat(['*', varname, '<sup>', i, '<', '/sup> + ']); 828 } else if (i === 1) { 829 t = t.concat(['*', varname, ' + ']); 830 } 831 } 832 833 return t.join(''); 834 }, 835 836 /** 837 * Computes the polynomial through a given set of coordinates in Lagrange form. 838 * Returns the Lagrange polynomials, see 839 * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation, 840 * SIAM Review, Vol 46, No 3, (2004) 501-517. 841 * @param {Array} p Array of JXG.Points 842 * @returns {function} A function of one parameter which returns the value of the polynomial, whose graph runs through the given points. 843 * @memberof JXG.Math.Numerics 844 */ 845 lagrangePolynomial: function (p) { 846 var w = [], 847 /** @ignore */ 848 fct = function (x, suspendedUpdate) { 849 var i, j, k, xi, s, M, 850 len = p.length, 851 num = 0, 852 denom = 0; 853 854 if (!suspendedUpdate) { 855 for (i = 0; i < len; i++) { 856 w[i] = 1.0; 857 xi = p[i].X(); 858 859 for (k = 0; k < len; k++) { 860 if (k !== i) { 861 w[i] *= (xi - p[k].X()); 862 } 863 } 864 865 w[i] = 1 / w[i]; 866 } 867 868 M = []; 869 870 for (j = 0; j < len; j++) { 871 M.push([1]); 872 } 873 } 874 875 for (i = 0; i < len; i++) { 876 xi = p[i].X(); 877 878 if (x === xi) { 879 return p[i].Y(); 880 } 881 882 s = w[i] / (x - xi); 883 denom += s; 884 num += s * p[i].Y(); 885 } 886 887 return num / denom; 888 }; 889 890 fct.getTerm = function () { 891 return ''; 892 }; 893 894 return fct; 895 }, 896 897 /** 898 * Computes the cubic cardinal spline curve through a given set of points. The curve 899 * is uniformly parametrized. 900 * Two artificial control points at the beginning and the end are added. 901 * @param {Array} points Array consisting of JXG.Points. 902 * @param {Number|Function} tau The tension parameter, either a constant number or a function returning a number. This number is between 0 and 1. 903 * tau=1/2 give Catmull-Rom splines. 904 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 905 * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, and a function simply 906 * returning the length of the points array minus three. 907 * @memberof JXG.Math.Numerics 908 */ 909 CardinalSpline: function (points, tau) { 910 var p, 911 coeffs = [], 912 // control point at the beginning and at the end 913 first = {}, 914 last = {}, 915 makeFct, 916 _tau; 917 918 if (Type.isFunction(tau)) { 919 _tau = tau; 920 } else { 921 _tau = function () { return tau; }; 922 } 923 924 /** @ignore */ 925 makeFct = function (which) { 926 return function (t, suspendedUpdate) { 927 var s, c, 928 len = points.length, 929 tau = _tau(); 930 931 if (len < 2) { 932 return NaN; 933 } 934 935 if (!suspendedUpdate) { 936 first[which] = function () { 937 return 2 * points[0][which]() - points[1][which](); 938 }; 939 940 last[which] = function () { 941 return 2 * points[len - 1][which]() - points[len - 2][which](); 942 }; 943 944 p = [first].concat(points, [last]); 945 coeffs[which] = []; 946 947 for (s = 0; s < len - 1; s++) { 948 coeffs[which][s] = [ 949 1 / tau * p[s + 1][which](), 950 -p[s][which]() + p[s + 2][which](), 951 2 * p[s][which]() + (-3 / tau + 1) * p[s + 1][which]() + (3 / tau - 2) * p[s + 2][which]() - p[s + 3][which](), 952 -p[s][which]() + (2 / tau - 1) * p[s + 1][which]() + (-2 / tau + 1) * p[s + 2][which]() + p[s + 3][which]() 953 ]; 954 } 955 } 956 957 len += 2; // add the two control points 958 959 if (isNaN(t)) { 960 return NaN; 961 } 962 963 // This is necessary for our advanced plotting algorithm: 964 if (t <= 0.0) { 965 return p[1][which](); 966 } 967 968 if (t >= len - 3) { 969 return p[len - 2][which](); 970 } 971 972 s = Math.floor(t); 973 974 if (s === t) { 975 return p[s][which](); 976 } 977 978 t -= s; 979 c = coeffs[which][s]; 980 981 return tau * (((c[3] * t + c[2]) * t + c[1]) * t + c[0]); 982 }; 983 }; 984 985 return [makeFct('X'), makeFct('Y'), 0, 986 function () { 987 return points.length - 1; 988 }]; 989 }, 990 991 /** 992 * Computes the cubic Catmull-Rom spline curve through a given set of points. The curve 993 * is uniformly parametrized. The curve is the cardinal spline curve for tau=0.5. 994 * Two artificial control points at the beginning and the end are added. 995 * @param {Array} points Array consisting of JXG.Points. 996 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 997 * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, and a function simply 998 * returning the length of the points array minus three. 999 * @memberof JXG.Math.Numerics 1000 */ 1001 CatmullRomSpline: function (points) { 1002 return this.CardinalSpline(points, 0.5); 1003 }, 1004 1005 /** 1006 * Computes the regression polynomial of a given degree through a given set of coordinates. 1007 * Returns the regression polynomial function. 1008 * @param {Number,function,Slider} degree number, function or slider. 1009 * Either 1010 * @param {Array} dataX Array containing either the x-coordinates of the data set or both coordinates in 1011 * an array of {@link JXG.Point}s or {@link JXG.Coords}. In the latter case, the <tt>dataY</tt> parameter will be ignored. 1012 * @param {Array} dataY Array containing the y-coordinates of the data set, 1013 * @returns {function} A function of one parameter which returns the value of the regression polynomial of the given degree. 1014 * It possesses the method getTerm() which returns the string containing the function term of the polynomial. 1015 * @memberof JXG.Math.Numerics 1016 */ 1017 regressionPolynomial: function (degree, dataX, dataY) { 1018 var coeffs, deg, dX, dY, inputType, fct, 1019 term = ''; 1020 1021 // Slider 1022 if (Type.isPoint(degree) && typeof degree.Value === 'function') { 1023 /** @ignore */ 1024 deg = function () { 1025 return degree.Value(); 1026 }; 1027 // function 1028 } else if (Type.isFunction(degree)) { 1029 deg = degree; 1030 // number 1031 } else if (Type.isNumber(degree)) { 1032 /** @ignore */ 1033 deg = function () { 1034 return degree; 1035 }; 1036 } else { 1037 throw new Error("JSXGraph: Can't create regressionPolynomial from degree of type'" + (typeof degree) + "'."); 1038 } 1039 1040 // Parameters degree, dataX, dataY 1041 if (arguments.length === 3 && Type.isArray(dataX) && Type.isArray(dataY)) { 1042 inputType = 0; 1043 // Parameters degree, point array 1044 } else if (arguments.length === 2 && Type.isArray(dataX) && dataX.length > 0 && Type.isPoint(dataX[0])) { 1045 inputType = 1; 1046 } else if (arguments.length === 2 && Type.isArray(dataX) && dataX.length > 0 && dataX[0].usrCoords && dataX[0].scrCoords) { 1047 inputType = 2; 1048 } else { 1049 throw new Error("JSXGraph: Can't create regressionPolynomial. Wrong parameters."); 1050 } 1051 1052 /** @ignore */ 1053 fct = function (x, suspendedUpdate) { 1054 var i, j, M, MT, y, B, c, s, d, 1055 // input data 1056 len = dataX.length; 1057 1058 d = Math.floor(deg()); 1059 1060 if (!suspendedUpdate) { 1061 // point list as input 1062 if (inputType === 1) { 1063 dX = []; 1064 dY = []; 1065 1066 for (i = 0; i < len; i++) { 1067 dX[i] = dataX[i].X(); 1068 dY[i] = dataX[i].Y(); 1069 } 1070 } 1071 1072 if (inputType === 2) { 1073 dX = []; 1074 dY = []; 1075 1076 for (i = 0; i < len; i++) { 1077 dX[i] = dataX[i].usrCoords[1]; 1078 dY[i] = dataX[i].usrCoords[2]; 1079 } 1080 } 1081 1082 // check for functions 1083 if (inputType === 0) { 1084 dX = []; 1085 dY = []; 1086 1087 for (i = 0; i < len; i++) { 1088 if (Type.isFunction(dataX[i])) { 1089 dX.push(dataX[i]()); 1090 } else { 1091 dX.push(dataX[i]); 1092 } 1093 1094 if (Type.isFunction(dataY[i])) { 1095 dY.push(dataY[i]()); 1096 } else { 1097 dY.push(dataY[i]); 1098 } 1099 } 1100 } 1101 1102 M = []; 1103 1104 for (j = 0; j < len; j++) { 1105 M.push([1]); 1106 } 1107 1108 for (i = 1; i <= d; i++) { 1109 for (j = 0; j < len; j++) { 1110 M[j][i] = M[j][i - 1] * dX[j]; 1111 } 1112 } 1113 1114 y = dY; 1115 MT = Mat.transpose(M); 1116 B = Mat.matMatMult(MT, M); 1117 c = Mat.matVecMult(MT, y); 1118 coeffs = Mat.Numerics.Gauss(B, c); 1119 term = Mat.Numerics.generatePolynomialTerm(coeffs, d, 'x', 3); 1120 } 1121 1122 // Horner's scheme to evaluate polynomial 1123 s = coeffs[d]; 1124 1125 for (i = d - 1; i >= 0; i--) { 1126 s = (s * x + coeffs[i]); 1127 } 1128 1129 return s; 1130 }; 1131 1132 fct.getTerm = function () { 1133 return term; 1134 }; 1135 1136 return fct; 1137 }, 1138 1139 /** 1140 * Computes the cubic Bezier curve through a given set of points. 1141 * @param {Array} points Array consisting of 3*k+1 {@link JXG.Points}. 1142 * The points at position k with k mod 3 = 0 are the data points, 1143 * points at position k with k mod 3 = 1 or 2 are the control points. 1144 * @returns {Array} An array consisting of two functions of one parameter t which return the 1145 * x resp. y coordinates of the Bezier curve in t, one zero value, and a third function accepting 1146 * no parameters and returning one third of the length of the points. 1147 * @memberof JXG.Math.Numerics 1148 */ 1149 bezier: function (points) { 1150 var len, flen, 1151 /** @ignore */ 1152 makeFct = function (which) { 1153 return function (t, suspendedUpdate) { 1154 var z = Math.floor(t) * 3, 1155 t0 = t % 1, 1156 t1 = 1 - t0; 1157 1158 if (!suspendedUpdate) { 1159 flen = 3 * Math.floor((points.length - 1) / 3); 1160 len = Math.floor(flen / 3); 1161 } 1162 1163 if (t < 0) { 1164 return points[0][which](); 1165 } 1166 1167 if (t >= len) { 1168 return points[flen][which](); 1169 } 1170 1171 if (isNaN(t)) { 1172 return NaN; 1173 } 1174 1175 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; 1176 }; 1177 }; 1178 1179 return [makeFct('X'), makeFct('Y'), 0, 1180 function () { 1181 return Math.floor(points.length / 3); 1182 }]; 1183 }, 1184 1185 /** 1186 * Computes the B-spline curve of order k (order = degree+1) through a given set of points. 1187 * @param {Array} points Array consisting of JXG.Points. 1188 * @param {Number} order Order of the B-spline curve. 1189 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 1190 * which return the x resp. y coordinates of the B-spline curve in t, a zero value, and a function simply 1191 * returning the length of the points array minus one. 1192 * @memberof JXG.Math.Numerics 1193 */ 1194 bspline: function (points, order) { 1195 var knots, N = [], 1196 _knotVector = function (n, k) { 1197 var j, 1198 kn = []; 1199 1200 for (j = 0; j < n + k + 1; j++) { 1201 if (j < k) { 1202 kn[j] = 0.0; 1203 } else if (j <= n) { 1204 kn[j] = j - k + 1; 1205 } else { 1206 kn[j] = n - k + 2; 1207 } 1208 } 1209 1210 return kn; 1211 }, 1212 1213 _evalBasisFuncs = function (t, kn, n, k, s) { 1214 var i, j, a, b, den, 1215 N = []; 1216 1217 if (kn[s] <= t && t < kn[s + 1]) { 1218 N[s] = 1; 1219 } else { 1220 N[s] = 0; 1221 } 1222 1223 for (i = 2; i <= k; i++) { 1224 for (j = s - i + 1; j <= s; j++) { 1225 if (j <= s - i + 1 || j < 0) { 1226 a = 0.0; 1227 } else { 1228 a = N[j]; 1229 } 1230 1231 if (j >= s) { 1232 b = 0.0; 1233 } else { 1234 b = N[j + 1]; 1235 } 1236 1237 den = kn[j + i - 1] - kn[j]; 1238 1239 if (den === 0) { 1240 N[j] = 0; 1241 } else { 1242 N[j] = (t - kn[j]) / den * a; 1243 } 1244 1245 den = kn[j + i] - kn[j + 1]; 1246 1247 if (den !== 0) { 1248 N[j] += (kn[j + i] - t) / den * b; 1249 } 1250 } 1251 } 1252 return N; 1253 }, 1254 /** @ignore */ 1255 makeFct = function (which) { 1256 return function (t, suspendedUpdate) { 1257 var y, j, s, 1258 len = points.length, 1259 n = len - 1, 1260 k = order; 1261 1262 if (n <= 0) { 1263 return NaN; 1264 } 1265 1266 if (n + 2 <= k) { 1267 k = n + 1; 1268 } 1269 1270 if (t <= 0) { 1271 return points[0][which](); 1272 } 1273 1274 if (t >= n - k + 2) { 1275 return points[n][which](); 1276 } 1277 1278 s = Math.floor(t) + k - 1; 1279 knots = _knotVector(n, k); 1280 N = _evalBasisFuncs(t, knots, n, k, s); 1281 1282 y = 0.0; 1283 for (j = s - k + 1; j <= s; j++) { 1284 if (j < len && j >= 0) { 1285 y += points[j][which]() * N[j]; 1286 } 1287 } 1288 1289 return y; 1290 }; 1291 }; 1292 1293 return [makeFct('X'), makeFct('Y'), 0, 1294 function () { 1295 return points.length - 1; 1296 }]; 1297 }, 1298 1299 /** 1300 * Numerical (symmetric) approximation of derivative. suspendUpdate is piped through, see {@link JXG.Curve#updateCurve} 1301 * and {@link JXG.Curve#hasPoint}. 1302 * @param {function} f Function in one variable to be differentiated. 1303 * @param {object} [obj] Optional object that is treated as "this" in the function body. This is useful, if the function is a 1304 * method of an object and contains a reference to its parent object via "this". 1305 * @returns {function} Derivative function of a given function f. 1306 * @memberof JXG.Math.Numerics 1307 */ 1308 D: function (f, obj) { 1309 var h = 0.00001, 1310 h2 = 1.0 / (h * 2.0); 1311 1312 if (!Type.exists(obj)) { 1313 return function (x, suspendUpdate) { 1314 return (f(x + h, suspendUpdate) - f(x - h, suspendUpdate)) * h2; 1315 }; 1316 } 1317 1318 return function (x, suspendUpdate) { 1319 return (f.apply(obj, [x + h, suspendUpdate]) - f.apply(obj, [x - h, suspendUpdate])) * h2; 1320 }; 1321 }, 1322 1323 /** 1324 * Helper function to create curve which displays Riemann sums. 1325 * Compute coordinates for the rectangles showing the Riemann sum. 1326 * @param {function} f Function, whose integral is approximated by the Riemann sum. 1327 * @param {Number} n number of rectangles. 1328 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson', or 'trapezoidal'. 1329 * @param {Number} start Left border of the approximation interval 1330 * @param {Number} end Right border of the approximation interval 1331 * @returns {Array} An array of two arrays containing the x and y coordinates for the rectangles showing the Riemann sum. This 1332 * 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 1333 * rectangles. 1334 * @memberof JXG.Math.Numerics 1335 */ 1336 riemann: function (f, n, type, start, end) { 1337 var i, x1, y1, delta1, delta, 1338 xarr = [], 1339 yarr = [], 1340 j = 0, 1341 x = start, y, 1342 sum = 0; 1343 1344 n = Math.round(n); 1345 1346 xarr[j] = x; 1347 yarr[j] = 0.0; 1348 1349 if (n > 0) { 1350 delta = (end - start) / n; 1351 // for 'lower' and 'upper' 1352 delta1 = delta * 0.01; 1353 1354 for (i = 0; i < n; i++) { 1355 if (type === 'right') { 1356 y = f(x + delta); 1357 } else if (type === 'middle') { 1358 y = f(x + delta * 0.5); 1359 } else if (type === 'left' || type === 'trapezoidal') { 1360 y = f(x); 1361 } else if (type === 'lower') { 1362 y = f(x); 1363 1364 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1365 y1 = f(x1); 1366 1367 if (y1 < y) { 1368 y = y1; 1369 } 1370 } 1371 1372 y1 = f(x + delta); 1373 if (y1 < y) { 1374 y = y1; 1375 } 1376 } else if (type === 'upper') { 1377 y = f(x); 1378 1379 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1380 y1 = f(x1); 1381 1382 if (y1 > y) { 1383 y = y1; 1384 } 1385 } 1386 1387 y1 = f(x + delta); 1388 if (y1 > y) { 1389 y = y1; 1390 } 1391 } else if (type === 'random') { 1392 y = f(x + delta * Math.random()); 1393 } else if (type === 'simpson') { 1394 y = (f(x) + 4 * f(x + delta * 0.5) + f(x + delta)) / 6.0; 1395 } else { 1396 y = f(x); // default is lower 1397 } 1398 1399 j += 1; 1400 xarr[j] = x; 1401 yarr[j] = y; 1402 j += 1; 1403 x += delta; 1404 1405 if (type === 'trapezoidal') { 1406 y = f(x); 1407 } 1408 1409 xarr[j] = x; 1410 yarr[j] = y; 1411 j += 1; 1412 xarr[j] = x; 1413 yarr[j] = 0.0; 1414 sum += y * delta; 1415 } 1416 } 1417 return [xarr, yarr, sum]; 1418 }, 1419 1420 /** 1421 * Approximate the integral by Riemann sums. 1422 * Compute the area described by the riemann sum rectangles. 1423 * @deprecated Replaced by JXG.Curve.Value(), see {@link JXG.Curve#riemannsum} 1424 * @param {function} f Function, whose integral is approximated by the Riemann sum. 1425 * @param {Number} n number of rectangles. 1426 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson' or 'trapezoidal'. 1427 * 1428 * @param {Number} start Left border of the approximation interval 1429 * @param {Number} end Right border of the approximation interval 1430 * @returns {Number} The sum of the areas of the rectangles. 1431 * @memberof JXG.Math.Numerics 1432 */ 1433 riemannsum: function (f, n, type, start, end) { 1434 var i, x1, y1, delta1, delta, y, 1435 sum = 0.0, 1436 x = start; 1437 1438 // this function looks very similar to this.riemann... maybe there is some merge potential? 1439 1440 n = Math.floor(n); 1441 1442 if (n > 0) { 1443 delta = (end - start) / n; 1444 // for 'lower' and 'upper' 1445 delta1 = delta * 0.01; 1446 1447 for (i = 0; i < n; i++) { 1448 if (type === 'right') { 1449 y = f(x + delta); 1450 } else if (type === 'middle') { 1451 y = f(x + delta * 0.5); 1452 } else if (type === 'trapezoidal') { 1453 y = 0.5 * (f(x + delta) + f(x)); 1454 } else if (type === 'left') { 1455 y = f(x); 1456 } else if (type === 'lower') { 1457 y = f(x); 1458 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1459 y1 = f(x1); 1460 1461 if (y1 < y) { 1462 y = y1; 1463 } 1464 } 1465 1466 y1 = f(x + delta); 1467 if (y1 < y) { 1468 y = y1; 1469 } 1470 } else if (type === 'upper') { 1471 y = f(x); 1472 1473 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1474 y1 = f(x1); 1475 1476 if (y1 > y) { 1477 y = y1; 1478 } 1479 } 1480 1481 y1 = f(x + delta); 1482 if (y1 > y) { 1483 y = y1; 1484 } 1485 } else if (type === 'random') { 1486 y = f(x + delta * Math.random()); 1487 } else if (type === 'simpson') { 1488 y = (f(x) + 4 * f(x + delta * 0.5) + f(x + delta)) / 6.0; 1489 } else { 1490 y = f(x); // default is lower 1491 } 1492 1493 sum += delta * y; 1494 x += delta; 1495 } 1496 } 1497 1498 return sum; 1499 }, 1500 1501 /** 1502 * Solve initial value problems numerically using Runge-Kutta-methods. 1503 * See {@link http://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm. 1504 * @param {object,String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing 1505 * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure 1506 * <pre> 1507 * { 1508 * s: <Number>, 1509 * A: <matrix>, 1510 * b: <Array>, 1511 * c: <Array> 1512 * } 1513 * </pre> 1514 * 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 1515 * @param {Array} x0 Initial value vector. If the problem is of one-dimensional, the initial value also has to be given in an array. 1516 * @param {Array} I Interval on which to integrate. 1517 * @param {Number} N Number of evaluation points. 1518 * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode 1519 * 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 1520 * vector <tt>x</tt>, and has to return a vector of the same dimension as <tt>x</tt> has. 1521 * @returns {Array} An array of vectors describing the solution of the ode on the given interval I. 1522 * @example 1523 * // A very simple autonomous system dx(t)/dt = x(t); 1524 * function f(t, x) { 1525 * return x; 1526 * } 1527 * 1528 * // Solve it with initial value x(0) = 1 on the interval [0, 2] 1529 * // with 20 evaluation points. 1530 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 1531 * 1532 * // Prepare data for plotting the solution of the ode using a curve. 1533 * var dataX = []; 1534 * var dataY = []; 1535 * var h = 0.1; // (I[1] - I[0])/N = (2-0)/20 1536 * for(var i=0; i<data.length; i++) { 1537 * dataX[i] = i*h; 1538 * dataY[i] = data[i][0]; 1539 * } 1540 * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'}); 1541 * </pre><div id="d2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div> 1542 * <script type="text/javascript"> 1543 * var board = JXG.JSXGraph.initBoard('d2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false}); 1544 * function f(t, x) { 1545 * // we have to copy the value. 1546 * // return x; would just return the reference. 1547 * return [x[0]]; 1548 * } 1549 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 1550 * var dataX = []; 1551 * var dataY = []; 1552 * var h = 0.1; 1553 * for(var i=0; i<data.length; i++) { 1554 * dataX[i] = i*h; 1555 * dataY[i] = data[i][0]; 1556 * } 1557 * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'}); 1558 * </script><pre> 1559 * @memberof JXG.Math.Numerics 1560 */ 1561 rungeKutta: function (butcher, x0, I, N, f) { 1562 var e, i, j, k, l, s, 1563 x = [], 1564 y = [], 1565 h = (I[1] - I[0]) / N, 1566 t = I[0], 1567 dim = x0.length, 1568 result = [], 1569 r = 0; 1570 1571 if (Type.isString(butcher)) { 1572 butcher = predefinedButcher[butcher] || predefinedButcher.euler; 1573 } 1574 s = butcher.s; 1575 1576 // don't change x0, so copy it 1577 for (e = 0; e < dim; e++) { 1578 x[e] = x0[e]; 1579 } 1580 1581 for (i = 0; i < N; i++) { 1582 // Optimization doesn't work for ODEs plotted using time 1583 // if((i % quotient == 0) || (i == N-1)) { 1584 result[r] = []; 1585 for (e = 0; e < dim; e++) { 1586 result[r][e] = x[e]; 1587 } 1588 1589 r += 1; 1590 k = []; 1591 1592 for (j = 0; j < s; j++) { 1593 // init y = 0 1594 for (e = 0; e < dim; e++) { 1595 y[e] = 0.0; 1596 } 1597 1598 1599 // Calculate linear combination of former k's and save it in y 1600 for (l = 0; l < j; l++) { 1601 for (e = 0; e < dim; e++) { 1602 y[e] += (butcher.A[j][l]) * h * k[l][e]; 1603 } 1604 } 1605 1606 // add x(t) to y 1607 for (e = 0; e < dim; e++) { 1608 y[e] += x[e]; 1609 } 1610 1611 // calculate new k and add it to the k matrix 1612 k.push(f(t + butcher.c[j] * h, y)); 1613 } 1614 1615 // init y = 0 1616 for (e = 0; e < dim; e++) { 1617 y[e] = 0.0; 1618 } 1619 1620 for (l = 0; l < s; l++) { 1621 for (e = 0; e < dim; e++) { 1622 y[e] += butcher.b[l] * k[l][e]; 1623 } 1624 } 1625 1626 for (e = 0; e < dim; e++) { 1627 x[e] = x[e] + h * y[e]; 1628 } 1629 1630 t += h; 1631 } 1632 1633 return result; 1634 }, 1635 1636 1637 /** 1638 * Maximum number of iterations in {@link JXG.Math.Numerics#fzero} 1639 * @type Number 1640 * @default 80 1641 * @memberof JXG.Math.Numerics 1642 */ 1643 maxIterationsRoot: 80, 1644 1645 /** 1646 * Maximum number of iterations in {@link JXG.Math.Numerics#fminbr} 1647 * @type Number 1648 * @default 500 1649 * @memberof JXG.Math.Numerics 1650 */ 1651 maxIterationsMinimize: 500, 1652 1653 /** 1654 * 1655 * Find zero of an univariate function f. 1656 * @param {function} f Function, whose root is to be found 1657 * @param {Array,Number} x0 Start value or start interval enclosing the root 1658 * @param {Object} object Parent object in case f is method of it 1659 * @returns {Number} the approximation of the root 1660 * Algorithm: 1661 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 1662 * computations. M., Mir, 1980, p.180 of the Russian edition 1663 * 1664 * If x0 is an array containing lower and upper bound for the zero 1665 * algorithm 748 is applied. Otherwise, if x0 is a number, 1666 * the algorithm tries to bracket a zero of f starting from x0. 1667 * If this fails, we fall back to Newton's method. 1668 * @memberof JXG.Math.Numerics 1669 */ 1670 fzero: function (f, x0, object) { 1671 var a, b, c, 1672 fa, fb, fc, 1673 aa, blist, i, len, u, fu, 1674 prev_step, t1, cb, t2, 1675 // Actual tolerance 1676 tol_act, 1677 // Interpolation step is calculated in the form p/q; division 1678 // operations is delayed until the last moment 1679 p, q, 1680 // Step at this iteration 1681 new_step, 1682 eps = Mat.eps, 1683 maxiter = this.maxIterationsRoot, 1684 niter = 0, 1685 nfev = 0; 1686 1687 if (Type.isArray(x0)) { 1688 if (x0.length < 2) { 1689 throw new Error("JXG.Math.Numerics.fzero: length of array x0 has to be at least two."); 1690 } 1691 1692 a = x0[0]; 1693 fa = f.call(object, a); 1694 nfev += 1; 1695 b = x0[1]; 1696 fb = f.call(object, b); 1697 nfev += 1; 1698 } else { 1699 a = x0; 1700 fa = f.call(object, a); 1701 nfev += 1; 1702 1703 // Try to get b. 1704 if (a === 0) { 1705 aa = 1; 1706 } else { 1707 aa = a; 1708 } 1709 1710 blist = [0.9 * aa, 1.1 * aa, aa - 1, aa + 1, 0.5 * aa, 1.5 * aa, -aa, 2 * aa, -10 * aa, 10 * aa]; 1711 len = blist.length; 1712 1713 for (i = 0; i < len; i++) { 1714 b = blist[i]; 1715 fb = f.call(object, b); 1716 nfev += 1; 1717 1718 if (fa * fb <= 0) { 1719 break; 1720 } 1721 } 1722 if (b < a) { 1723 u = a; 1724 a = b; 1725 b = u; 1726 1727 fu = fa; 1728 fa = fb; 1729 fb = fu; 1730 } 1731 } 1732 1733 if (fa * fb > 0) { 1734 // Bracketing not successful, fall back to Newton's method or to fminbr 1735 if (Type.isArray(x0)) { 1736 return this.fminbr(f, [a, b], object); 1737 } 1738 1739 return this.Newton(f, a, object); 1740 } 1741 1742 // OK, we have enclosed a zero of f. 1743 // Now we can start Brent's method 1744 1745 c = a; 1746 fc = fa; 1747 1748 // Main iteration loop 1749 while (niter < maxiter) { 1750 // Distance from the last but one to the last approximation 1751 prev_step = b - a; 1752 1753 // Swap data for b to be the best approximation 1754 if (Math.abs(fc) < Math.abs(fb)) { 1755 a = b; 1756 b = c; 1757 c = a; 1758 1759 fa = fb; 1760 fb = fc; 1761 fc = fa; 1762 } 1763 tol_act = 2 * eps * Math.abs(b) + eps * 0.5; 1764 new_step = (c - b) * 0.5; 1765 1766 if (Math.abs(new_step) <= tol_act && Math.abs(fb) <= eps) { 1767 // Acceptable approx. is found 1768 return b; 1769 } 1770 1771 // Decide if the interpolation can be tried 1772 // If prev_step was large enough and was in true direction Interpolatiom may be tried 1773 if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) { 1774 cb = c - b; 1775 1776 // If we have only two distinct points linear interpolation can only be applied 1777 if (a === c) { 1778 t1 = fb / fa; 1779 p = cb * t1; 1780 q = 1.0 - t1; 1781 // Quadric inverse interpolation 1782 } else { 1783 q = fa / fc; 1784 t1 = fb / fc; 1785 t2 = fb / fa; 1786 1787 p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0)); 1788 q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0); 1789 } 1790 1791 // p was calculated with the opposite sign; make p positive 1792 if (p > 0) { 1793 q = -q; 1794 // and assign possible minus to q 1795 } else { 1796 p = -p; 1797 } 1798 1799 // If b+p/q falls in [b,c] and isn't too large it is accepted 1800 // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent 1801 if (p < (0.75 * cb * q - Math.abs(tol_act * q) * 0.5) && 1802 p < Math.abs(prev_step * q * 0.5)) { 1803 new_step = p / q; 1804 } 1805 } 1806 1807 // Adjust the step to be not less than tolerance 1808 if (Math.abs(new_step) < tol_act) { 1809 if (new_step > 0) { 1810 new_step = tol_act; 1811 } else { 1812 new_step = -tol_act; 1813 } 1814 } 1815 1816 // Save the previous approx. 1817 a = b; 1818 fa = fb; 1819 b += new_step; 1820 fb = f.call(object, b); 1821 // Do step to a new approxim. 1822 nfev += 1; 1823 1824 // Adjust c for it to have a sign opposite to that of b 1825 if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) { 1826 c = a; 1827 fc = fa; 1828 } 1829 niter++; 1830 } // End while 1831 1832 return b; 1833 }, 1834 1835 1836 /** 1837 * 1838 * Find minimum of an univariate function f. 1839 * @param {function} f Function, whose minimum is to be found 1840 * @param {Array} x0 Start interval enclosing the minimum 1841 * @param {Object} context Parent object in case f is method of it 1842 * @returns {Number} the approximation of the minimum value position 1843 * Algorithm: 1844 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 1845 * computations. M., Mir, 1980, p.180 of the Russian edition 1846 * x0 1847 * @memberof JXG.Math.Numerics 1848 **/ 1849 fminbr: function (f, x0, context) { 1850 var a, b, x, v, w, 1851 fx, fv, fw, 1852 range, middle_range, tol_act, new_step, 1853 p, q, t, ft, 1854 // Golden section ratio 1855 r = (3.0 - Math.sqrt(5.0)) * 0.5, 1856 tol = Mat.eps, 1857 sqrteps = Mat.eps, //Math.sqrt(Mat.eps), 1858 maxiter = this.maxIterationsMinimize, 1859 niter = 0, 1860 nfev = 0; 1861 1862 if (!Type.isArray(x0) || x0.length < 2) { 1863 throw new Error("JXG.Math.Numerics.fminbr: length of array x0 has to be at least two."); 1864 } 1865 1866 a = x0[0]; 1867 b = x0[1]; 1868 v = a + r * (b - a); 1869 fv = f.call(context, v); 1870 1871 // First step - always gold section 1872 nfev += 1; 1873 x = v; 1874 w = v; 1875 fx = fv; 1876 fw = fv; 1877 1878 while (niter < maxiter) { 1879 // Range over the interval in which we are looking for the minimum 1880 range = b - a; 1881 middle_range = (a + b) * 0.5; 1882 1883 // Actual tolerance 1884 tol_act = sqrteps * Math.abs(x) + tol / 3.0; 1885 1886 if (Math.abs(x - middle_range) + range * 0.5 <= 2.0 * tol_act) { 1887 // Acceptable approx. is found 1888 return x; 1889 } 1890 1891 // Obtain the golden section step 1892 new_step = r * (x < middle_range ? b - x : a - x); 1893 1894 // Decide if the interpolation can be tried. If x and w are distinct interpolatiom may be tried 1895 if (Math.abs(x - w) >= tol_act) { 1896 // Interpolation step is calculated as p/q; 1897 // division operation is delayed until last moment 1898 t = (x - w) * (fx - fv); 1899 q = (x - v) * (fx - fw); 1900 p = (x - v) * q - (x - w) * t; 1901 q = 2 * (q - t); 1902 1903 if (q > 0) { // q was calculated with the op- 1904 p = -p; // posite sign; make q positive 1905 } else { // and assign possible minus to 1906 q = -q; // p 1907 } 1908 if (Math.abs(p) < Math.abs(new_step * q) && // If x+p/q falls in [a,b] 1909 p > q * (a - x + 2 * tol_act) && // not too close to a and 1910 p < q * (b - x - 2 * tol_act)) { // b, and isn't too large 1911 new_step = p / q; // it is accepted 1912 } 1913 // If p/q is too large then the 1914 // golden section procedure can 1915 // reduce [a,b] range to more 1916 // extent 1917 } 1918 1919 // Adjust the step to be not less than tolerance 1920 if (Math.abs(new_step) < tol_act) { 1921 if (new_step > 0) { 1922 new_step = tol_act; 1923 } else { 1924 new_step = -tol_act; 1925 } 1926 } 1927 1928 // Obtain the next approximation to min 1929 // and reduce the enveloping range 1930 1931 // Tentative point for the min 1932 t = x + new_step; 1933 ft = f.call(context, t); 1934 nfev += 1; 1935 1936 // t is a better approximation 1937 if (ft <= fx) { 1938 // Reduce the range so that t would fall within it 1939 if (t < x) { 1940 b = x; 1941 } else { 1942 a = x; 1943 } 1944 1945 // Assign the best approx to x 1946 v = w; 1947 w = x; 1948 x = t; 1949 1950 fv = fw; 1951 fw = fx; 1952 fx = ft; 1953 // x remains the better approx 1954 } else { 1955 // Reduce the range enclosing x 1956 if (t < x) { 1957 a = t; 1958 } else { 1959 b = t; 1960 } 1961 1962 if (ft <= fw || w === x) { 1963 v = w; 1964 w = t; 1965 fv = fw; 1966 fw = ft; 1967 } else if (ft <= fv || v === x || v === w) { 1968 v = t; 1969 fv = ft; 1970 } 1971 } 1972 niter += 1; 1973 } 1974 1975 return x; 1976 }, 1977 1978 /** 1979 * Implements the Ramer-Douglas-Peucker algorithm. 1980 * It discards points which are not necessary from the polygonal line defined by the point array 1981 * pts. The computation is done in screen coordinates. 1982 * Average runtime is O(nlog(n)), worst case runtime is O(n^2), where n is the number of points. 1983 * @param {Array} pts Array of {@link JXG.Coords} 1984 * @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>. 1985 * @returns {Array} An array containing points which represent an apparently identical curve as the points of pts do, but contains fewer points. 1986 * @memberof JXG.Math.Numerics 1987 */ 1988 RamerDouglasPeucker: function (pts, eps) { 1989 var newPts = [], i, k, len, 1990 1991 /** 1992 * findSplit() is a subroutine of {@link JXG.Math.Numerics#RamerDouglasPeucker}. 1993 * It searches for the point between index i and j which 1994 * has the largest distance from the line between the points i and j. 1995 * @param {Array} pts Array of {@link JXG.Coords} 1996 * @param {Number} i Index of a point in pts 1997 * @param {Number} j Index of a point in pts 1998 * @ignore 1999 * @private 2000 */ 2001 findSplit = function (pts, i, j) { 2002 var d, k, ci, cj, ck, 2003 x0, y0, x1, y1, 2004 den, lbda, 2005 dist = 0, 2006 f = i; 2007 2008 if (j - i < 2) { 2009 return [-1.0, 0]; 2010 } 2011 2012 ci = pts[i].scrCoords; 2013 cj = pts[j].scrCoords; 2014 2015 if (isNaN(ci[1] + ci[2])) { 2016 return [NaN, i]; 2017 } 2018 if (isNaN(cj[1] + cj[2])) { 2019 return [NaN, j]; 2020 } 2021 2022 for (k = i + 1; k < j; k++) { 2023 ck = pts[k].scrCoords; 2024 if (isNaN(ck[1] + ck[2])) { 2025 return [NaN, k]; 2026 } 2027 2028 x0 = ck[1] - ci[1]; 2029 y0 = ck[2] - ci[2]; 2030 x1 = cj[1] - ci[1]; 2031 y1 = cj[2] - ci[2]; 2032 den = x1 * x1 + y1 * y1; 2033 2034 if (den >= Mat.eps) { 2035 lbda = (x0 * x1 + y0 * y1) / den; 2036 2037 if (lbda < 0.0) { 2038 lbda = 0.0; 2039 } else if (lbda > 1.0) { 2040 lbda = 1.0; 2041 } 2042 2043 x0 = x0 - lbda * x1; 2044 y0 = y0 - lbda * y1; 2045 d = x0 * x0 + y0 * y0; 2046 } else { 2047 lbda = 0.0; 2048 d = x0 * x0 + y0 * y0; 2049 } 2050 2051 if (d > dist) { 2052 dist = d; 2053 f = k; 2054 } 2055 } 2056 return [Math.sqrt(dist), f]; 2057 }, 2058 2059 /** 2060 * RDP() is a private subroutine of {@link JXG.Math.Numerics#RamerDouglasPeucker}. 2061 * It runs recursively through the point set and searches the 2062 * point which has the largest distance from the line between the first point and 2063 * the last point. If the distance from the line is greater than eps, this point is 2064 * included in our new point set otherwise it is discarded. 2065 * If it is taken, we recursively apply the subroutine to the point set before 2066 * and after the chosen point. 2067 * @param {Array} pts Array of {@link JXG.Coords} 2068 * @param {Number} i Index of an element of pts 2069 * @param {Number} j Index of an element of pts 2070 * @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>. 2071 * @param {Array} newPts Array of {@link JXG.Coords} 2072 * @ignore 2073 * @private 2074 */ 2075 RDP = function (pts, i, j, eps, newPts) { 2076 var result = findSplit(pts, i, j), 2077 k = result[1]; 2078 2079 if (isNaN(result[0])) { 2080 RDP(pts, i, k - 1, eps, newPts); 2081 newPts.push(pts[k]); 2082 do { 2083 ++k; 2084 } while (k <= j && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])); 2085 if (k <= j) { 2086 newPts.push(pts[k]); 2087 } 2088 RDP(pts, k + 1, j, eps, newPts); 2089 } else if (result[0] > eps) { 2090 RDP(pts, i, k, eps, newPts); 2091 RDP(pts, k, j, eps, newPts); 2092 } else { 2093 newPts.push(pts[j]); 2094 } 2095 }; 2096 2097 len = pts.length; 2098 2099 // Search for the left most point woithout NaN coordinates 2100 i = 0; 2101 while (i < len && isNaN(pts[i].scrCoords[1] + pts[i].scrCoords[2])) { 2102 i += 1; 2103 } 2104 // Search for the right most point woithout NaN coordinates 2105 k = len - 1; 2106 while (k > i && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])) { 2107 k -= 1; 2108 } 2109 2110 // Only proceed if something is left 2111 if (!(i > k || i === len)) { 2112 newPts[0] = pts[i]; 2113 RDP(pts, i, k, eps, newPts); 2114 } 2115 2116 return newPts; 2117 }, 2118 2119 /** 2120 * Old name for the implementation of the Ramer-Douglas-Peucker algorithm. 2121 * @deprecated Use {@link JXG.Math.Numerics#RamerDouglasPeucker} 2122 * @memberof JXG.Math.Numerics 2123 */ 2124 RamerDouglasPeuker: function (pts, eps) { 2125 return this.RamerDouglasPeucker(pts, eps); 2126 } 2127 }; 2128 2129 return Mat.Numerics; 2130 }); 2131