1 /* 2 Copyright 2008-2013 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 } else if (type === 'upper') { 1372 y = f(x); 1373 1374 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1375 y1 = f(x1); 1376 1377 if (y1 > y) { 1378 y = y1; 1379 } 1380 } 1381 } else if (type === 'random') { 1382 y = f(x + delta * Math.random()); 1383 } else if (type === 'simpson') { 1384 y = (f(x) + 4 * f(x + delta * 0.5) + f(x + delta)) / 6.0; 1385 } else { 1386 y = f(x); // default is lower 1387 } 1388 1389 j += 1; 1390 xarr[j] = x; 1391 yarr[j] = y; 1392 j += 1; 1393 x += delta; 1394 1395 if (type === 'trapezoidal') { 1396 y = f(x); 1397 } 1398 1399 xarr[j] = x; 1400 yarr[j] = y; 1401 j += 1; 1402 xarr[j] = x; 1403 yarr[j] = 0.0; 1404 sum += y * delta; 1405 } 1406 } 1407 return [xarr, yarr, sum]; 1408 }, 1409 1410 /** 1411 * Approximate the integral by Riemann sums. 1412 * Compute the area described by the riemann sum rectangles. 1413 * @deprecated Replaced by JXG.Curve.Value(), see {@link JXG.Curve#riemannsum} 1414 * @param {function} f Function, whose integral is approximated by the Riemann sum. 1415 * @param {Number} n number of rectangles. 1416 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson' or 'trapezoidal'. 1417 * 1418 * @param {Number} start Left border of the approximation interval 1419 * @param {Number} end Right border of the approximation interval 1420 * @returns {Number} The sum of the areas of the rectangles. 1421 * @memberof JXG.Math.Numerics 1422 */ 1423 riemannsum: function (f, n, type, start, end) { 1424 var i, x1, y1, delta1, delta, y, 1425 sum = 0.0, 1426 x = start; 1427 1428 // this function looks very similar to this.riemann... maybe there is some merge potential? 1429 1430 n = Math.floor(n); 1431 1432 if (n > 0) { 1433 delta = (end - start) / n; 1434 // for 'lower' and 'upper' 1435 delta1 = delta * 0.01; 1436 1437 for (i = 0; i < n; i++) { 1438 if (type === 'right') { 1439 y = f(x + delta); 1440 } else if (type === 'middle') { 1441 y = f(x + delta * 0.5); 1442 } else if (type === 'trapezoidal') { 1443 y = 0.5 * (f(x + delta) + f(x)); 1444 } else if (type === 'left') { 1445 y = f(x); 1446 } else if (type === 'lower') { 1447 y = f(x); 1448 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1449 y1 = f(x1); 1450 1451 if (y1 < y) { 1452 y = y1; 1453 } 1454 } 1455 } else if (type === 'upper') { 1456 y = f(x); 1457 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 } else if (type === 'random') { 1466 y = f(x + delta * Math.random()); 1467 } else if (type === 'simpson') { 1468 y = (f(x) + 4 * f(x + delta * 0.5) + f(x + delta)) / 6.0; 1469 } else { 1470 y = f(x); // default is lower 1471 } 1472 1473 sum += delta * y; 1474 x += delta; 1475 } 1476 } 1477 1478 return sum; 1479 }, 1480 1481 /** 1482 * Solve initial value problems numerically using Runge-Kutta-methods. 1483 * See {@link http://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm. 1484 * @param {object,String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing 1485 * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure 1486 * <pre> 1487 * { 1488 * s: <Number>, 1489 * A: <matrix>, 1490 * b: <Array>, 1491 * c: <Array> 1492 * } 1493 * </pre> 1494 * 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 1495 * @param {Array} x0 Initial value vector. If the problem is of one-dimensional, the initial value also has to be given in an array. 1496 * @param {Array} I Interval on which to integrate. 1497 * @param {Number} N Number of evaluation points. 1498 * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode 1499 * 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 1500 * vector <tt>x</tt>, and has to return a vector of the same dimension as <tt>x</tt> has. 1501 * @returns {Array} An array of vectors describing the solution of the ode on the given interval I. 1502 * @example 1503 * // A very simple autonomous system dx(t)/dt = x(t); 1504 * function f(t, x) { 1505 * return x; 1506 * } 1507 * 1508 * // Solve it with initial value x(0) = 1 on the interval [0, 2] 1509 * // with 20 evaluation points. 1510 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 1511 * 1512 * // Prepare data for plotting the solution of the ode using a curve. 1513 * var dataX = []; 1514 * var dataY = []; 1515 * var h = 0.1; // (I[1] - I[0])/N = (2-0)/20 1516 * for(var i=0; i<data.length; i++) { 1517 * dataX[i] = i*h; 1518 * dataY[i] = data[i][0]; 1519 * } 1520 * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'}); 1521 * </pre><div id="d2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div> 1522 * <script type="text/javascript"> 1523 * var board = JXG.JSXGraph.initBoard('d2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false}); 1524 * function f(t, x) { 1525 * // we have to copy the value. 1526 * // return x; would just return the reference. 1527 * return [x[0]]; 1528 * } 1529 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 1530 * var dataX = []; 1531 * var dataY = []; 1532 * var h = 0.1; 1533 * for(var i=0; i<data.length; i++) { 1534 * dataX[i] = i*h; 1535 * dataY[i] = data[i][0]; 1536 * } 1537 * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'}); 1538 * </script><pre> 1539 * @memberof JXG.Math.Numerics 1540 */ 1541 rungeKutta: function (butcher, x0, I, N, f) { 1542 var e, i, j, k, l, s, 1543 x = [], 1544 y = [], 1545 h = (I[1] - I[0]) / N, 1546 t = I[0], 1547 dim = x0.length, 1548 result = [], 1549 r = 0; 1550 1551 if (Type.isString(butcher)) { 1552 butcher = predefinedButcher[butcher] || predefinedButcher.euler; 1553 } 1554 s = butcher.s; 1555 1556 // don't change x0, so copy it 1557 for (e = 0; e < dim; e++) { 1558 x[e] = x0[e]; 1559 } 1560 1561 for (i = 0; i < N; i++) { 1562 // Optimization doesn't work for ODEs plotted using time 1563 // if((i % quotient == 0) || (i == N-1)) { 1564 result[r] = []; 1565 for (e = 0; e < dim; e++) { 1566 result[r][e] = x[e]; 1567 } 1568 1569 r += 1; 1570 k = []; 1571 1572 for (j = 0; j < s; j++) { 1573 // init y = 0 1574 for (e = 0; e < dim; e++) { 1575 y[e] = 0.0; 1576 } 1577 1578 1579 // Calculate linear combination of former k's and save it in y 1580 for (l = 0; l < j; l++) { 1581 for (e = 0; e < dim; e++) { 1582 y[e] += (butcher.A[j][l]) * h * k[l][e]; 1583 } 1584 } 1585 1586 // add x(t) to y 1587 for (e = 0; e < dim; e++) { 1588 y[e] += x[e]; 1589 } 1590 1591 // calculate new k and add it to the k matrix 1592 k.push(f(t + butcher.c[j] * h, y)); 1593 } 1594 1595 // init y = 0 1596 for (e = 0; e < dim; e++) { 1597 y[e] = 0.0; 1598 } 1599 1600 for (l = 0; l < s; l++) { 1601 for (e = 0; e < dim; e++) { 1602 y[e] += butcher.b[l] * k[l][e]; 1603 } 1604 } 1605 1606 for (e = 0; e < dim; e++) { 1607 x[e] = x[e] + h * y[e]; 1608 } 1609 1610 t += h; 1611 } 1612 1613 return result; 1614 }, 1615 1616 1617 /** 1618 * Maximum number of iterations in {@link JXG.Math.Numerics#fzero} 1619 * @type Number 1620 * @default 80 1621 * @memberof JXG.Math.Numerics 1622 */ 1623 maxIterationsRoot: 80, 1624 1625 /** 1626 * Maximum number of iterations in {@link JXG.Math.Numerics#fminbr} 1627 * @type Number 1628 * @default 500 1629 * @memberof JXG.Math.Numerics 1630 */ 1631 maxIterationsMinimize: 500, 1632 1633 /** 1634 * 1635 * Find zero of an univariate function f. 1636 * @param {function} f Function, whose root is to be found 1637 * @param {Array,Number} x0 Start value or start interval enclosing the root 1638 * @param {Object} object Parent object in case f is method of it 1639 * @returns {Number} the approximation of the root 1640 * Algorithm: 1641 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 1642 * computations. M., Mir, 1980, p.180 of the Russian edition 1643 * 1644 * If x0 is an array containing lower and upper bound for the zero 1645 * algorithm 748 is applied. Otherwise, if x0 is a number, 1646 * the algorithm tries to bracket a zero of f starting from x0. 1647 * If this fails, we fall back to Newton's method. 1648 * @memberof JXG.Math.Numerics 1649 */ 1650 fzero: function (f, x0, object) { 1651 var a, b, c, 1652 fa, fb, fc, 1653 aa, blist, i, len, u, fu, 1654 prev_step, t1, cb, t2, 1655 // Actual tolerance 1656 tol_act, 1657 // Interpolation step is calculated in the form p/q; division 1658 // operations is delayed until the last moment 1659 p, q, 1660 // Step at this iteration 1661 new_step, 1662 eps = Mat.eps, 1663 maxiter = this.maxIterationsRoot, 1664 niter = 0, 1665 nfev = 0; 1666 1667 if (Type.isArray(x0)) { 1668 if (x0.length < 2) { 1669 throw new Error("JXG.Math.Numerics.fzero: length of array x0 has to be at least two."); 1670 } 1671 1672 a = x0[0]; 1673 fa = f.call(object, a); 1674 nfev += 1; 1675 b = x0[1]; 1676 fb = f.call(object, b); 1677 nfev += 1; 1678 } else { 1679 a = x0; 1680 fa = f.call(object, a); 1681 nfev += 1; 1682 1683 // Try to get b. 1684 if (a === 0) { 1685 aa = 1; 1686 } else { 1687 aa = a; 1688 } 1689 1690 blist = [0.9 * aa, 1.1 * aa, aa - 1, aa + 1, 0.5 * aa, 1.5 * aa, -aa, 2 * aa, -10 * aa, 10 * aa]; 1691 len = blist.length; 1692 1693 for (i = 0; i < len; i++) { 1694 b = blist[i]; 1695 fb = f.call(object, b); 1696 nfev += 1; 1697 1698 if (fa * fb <= 0) { 1699 break; 1700 } 1701 } 1702 if (b < a) { 1703 u = a; 1704 a = b; 1705 b = u; 1706 1707 fu = fa; 1708 fa = fb; 1709 fb = fu; 1710 } 1711 } 1712 1713 if (fa * fb > 0) { 1714 // Bracketing not successful, fall back to Newton's method or to fminbr 1715 if (Type.isArray(x0)) { 1716 return this.fminbr(f, [a, b], object); 1717 } 1718 1719 return this.Newton(f, a, object); 1720 } 1721 1722 // OK, we have enclosed a zero of f. 1723 // Now we can start Brent's method 1724 1725 c = a; 1726 fc = fa; 1727 1728 // Main iteration loop 1729 while (niter < maxiter) { 1730 // Distance from the last but one to the last approximation 1731 prev_step = b - a; 1732 1733 // Swap data for b to be the best approximation 1734 if (Math.abs(fc) < Math.abs(fb)) { 1735 a = b; 1736 b = c; 1737 c = a; 1738 1739 fa = fb; 1740 fb = fc; 1741 fc = fa; 1742 } 1743 tol_act = 2 * eps * Math.abs(b) + eps * 0.5; 1744 new_step = (c - b) * 0.5; 1745 1746 if (Math.abs(new_step) <= tol_act && Math.abs(fb) <= eps) { 1747 // Acceptable approx. is found 1748 return b; 1749 } 1750 1751 // Decide if the interpolation can be tried 1752 // If prev_step was large enough and was in true direction Interpolatiom may be tried 1753 if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) { 1754 cb = c - b; 1755 1756 // If we have only two distinct points linear interpolation can only be applied 1757 if (a === c) { 1758 t1 = fb / fa; 1759 p = cb * t1; 1760 q = 1.0 - t1; 1761 // Quadric inverse interpolation 1762 } else { 1763 q = fa / fc; 1764 t1 = fb / fc; 1765 t2 = fb / fa; 1766 1767 p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0)); 1768 q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0); 1769 } 1770 1771 // p was calculated with the opposite sign; make p positive 1772 if (p > 0) { 1773 q = -q; 1774 // and assign possible minus to q 1775 } else { 1776 p = -p; 1777 } 1778 1779 // If b+p/q falls in [b,c] and isn't too large it is accepted 1780 // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent 1781 if (p < (0.75 * cb * q - Math.abs(tol_act * q) * 0.5) && 1782 p < Math.abs(prev_step * q * 0.5)) { 1783 new_step = p / q; 1784 } 1785 } 1786 1787 // Adjust the step to be not less than tolerance 1788 if (Math.abs(new_step) < tol_act) { 1789 if (new_step > 0) { 1790 new_step = tol_act; 1791 } else { 1792 new_step = -tol_act; 1793 } 1794 } 1795 1796 // Save the previous approx. 1797 a = b; 1798 fa = fb; 1799 b += new_step; 1800 fb = f.call(object, b); 1801 // Do step to a new approxim. 1802 nfev += 1; 1803 1804 // Adjust c for it to have a sign opposite to that of b 1805 if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) { 1806 c = a; 1807 fc = fa; 1808 } 1809 niter++; 1810 } // End while 1811 1812 return b; 1813 }, 1814 1815 1816 /** 1817 * 1818 * Find minimum of an univariate function f. 1819 * @param {function} f Function, whose minimum is to be found 1820 * @param {Array} x0 Start interval enclosing the minimum 1821 * @param {Object} context Parent object in case f is method of it 1822 * @returns {Number} the approximation of the minimum value position 1823 * Algorithm: 1824 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 1825 * computations. M., Mir, 1980, p.180 of the Russian edition 1826 * x0 1827 * @memberof JXG.Math.Numerics 1828 **/ 1829 fminbr: function (f, x0, context) { 1830 var a, b, x, v, w, 1831 fx, fv, fw, 1832 range, middle_range, tol_act, new_step, 1833 p, q, t, ft, 1834 // Golden section ratio 1835 r = (3.0 - Math.sqrt(5.0)) * 0.5, 1836 tol = Mat.eps, 1837 sqrteps = Mat.eps, //Math.sqrt(Mat.eps), 1838 maxiter = this.maxIterationsMinimize, 1839 niter = 0, 1840 nfev = 0; 1841 1842 if (!Type.isArray(x0) || x0.length < 2) { 1843 throw new Error("JXG.Math.Numerics.fminbr: length of array x0 has to be at least two."); 1844 } 1845 1846 a = x0[0]; 1847 b = x0[1]; 1848 v = a + r * (b - a); 1849 fv = f.call(context, v); 1850 1851 // First step - always gold section 1852 nfev += 1; 1853 x = v; 1854 w = v; 1855 fx = fv; 1856 fw = fv; 1857 1858 while (niter < maxiter) { 1859 // Range over the interval in which we are looking for the minimum 1860 range = b - a; 1861 middle_range = (a + b) * 0.5; 1862 1863 // Actual tolerance 1864 tol_act = sqrteps * Math.abs(x) + tol / 3.0; 1865 1866 if (Math.abs(x - middle_range) + range * 0.5 <= 2.0 * tol_act) { 1867 // Acceptable approx. is found 1868 return x; 1869 } 1870 1871 // Obtain the golden section step 1872 new_step = r * (x < middle_range ? b - x : a - x); 1873 1874 // Decide if the interpolation can be tried. If x and w are distinct interpolatiom may be tried 1875 if (Math.abs(x - w) >= tol_act) { 1876 // Interpolation step is calculated as p/q; 1877 // division operation is delayed until last moment 1878 t = (x - w) * (fx - fv); 1879 q = (x - v) * (fx - fw); 1880 p = (x - v) * q - (x - w) * t; 1881 q = 2 * (q - t); 1882 1883 if (q > 0) { // q was calculated with the op- 1884 p = -p; // posite sign; make q positive 1885 } else { // and assign possible minus to 1886 q = -q; // p 1887 } 1888 if (Math.abs(p) < Math.abs(new_step * q) && // If x+p/q falls in [a,b] 1889 p > q * (a - x + 2 * tol_act) && // not too close to a and 1890 p < q * (b - x - 2 * tol_act)) { // b, and isn't too large 1891 new_step = p / q; // it is accepted 1892 } 1893 // If p/q is too large then the 1894 // golden section procedure can 1895 // reduce [a,b] range to more 1896 // extent 1897 } 1898 1899 // Adjust the step to be not less than tolerance 1900 if (Math.abs(new_step) < tol_act) { 1901 if (new_step > 0) { 1902 new_step = tol_act; 1903 } else { 1904 new_step = -tol_act; 1905 } 1906 } 1907 1908 // Obtain the next approximation to min 1909 // and reduce the enveloping range 1910 1911 // Tentative point for the min 1912 t = x + new_step; 1913 ft = f.call(context, t); 1914 nfev += 1; 1915 1916 // t is a better approximation 1917 if (ft <= fx) { 1918 // Reduce the range so that t would fall within it 1919 if (t < x) { 1920 b = x; 1921 } else { 1922 a = x; 1923 } 1924 1925 // Assign the best approx to x 1926 v = w; 1927 w = x; 1928 x = t; 1929 1930 fv = fw; 1931 fw = fx; 1932 fx = ft; 1933 // x remains the better approx 1934 } else { 1935 // Reduce the range enclosing x 1936 if (t < x) { 1937 a = t; 1938 } else { 1939 b = t; 1940 } 1941 1942 if (ft <= fw || w === x) { 1943 v = w; 1944 w = t; 1945 fv = fw; 1946 fw = ft; 1947 } else if (ft <= fv || v === x || v === w) { 1948 v = t; 1949 fv = ft; 1950 } 1951 } 1952 niter += 1; 1953 } 1954 1955 return x; 1956 }, 1957 1958 /** 1959 * Implements the Ramer-Douglas-Peuker algorithm. 1960 * It discards points which are not necessary from the polygonal line defined by the point array 1961 * pts. The computation is done in screen coordinates. 1962 * Average runtime is O(nlog(n)), worst case runtime is O(n^2), where n is the number of points. 1963 * @param {Array} pts Array of {@link JXG.Coords} 1964 * @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>. 1965 * @returns {Array} An array containing points which represent an apparently identical curve as the points of pts do, but contains fewer points. 1966 * @memberof JXG.Math.Numerics 1967 */ 1968 RamerDouglasPeuker: function (pts, eps) { 1969 var newPts = [], i, k, len, 1970 1971 /** 1972 * findSplit() is a subroutine of {@link JXG.Math.Numerics#RamerDouglasPeuker}. 1973 * It searches for the point between index i and j which 1974 * has the largest distance from the line between the points i and j. 1975 * @param {Array} pts Array of {@link JXG.Coords} 1976 * @param {Number} i Index of a point in pts 1977 * @param {Number} j Index of a point in pts 1978 * @ignore 1979 * @private 1980 */ 1981 findSplit = function (pts, i, j) { 1982 var d, k, ci, cj, ck, 1983 x0, y0, x1, y1, 1984 den, lbda, 1985 dist = 0, 1986 f = i; 1987 1988 if (j - i < 2) { 1989 return [-1.0, 0]; 1990 } 1991 1992 ci = pts[i].scrCoords; 1993 cj = pts[j].scrCoords; 1994 1995 if (isNaN(ci[1] + ci[2] + cj[1] + cj[2])) { 1996 return [NaN, j]; 1997 } 1998 1999 for (k = i + 1; k < j; k++) { 2000 ck = pts[k].scrCoords; 2001 x0 = ck[1] - ci[1]; 2002 y0 = ck[2] - ci[2]; 2003 x1 = cj[1] - ci[1]; 2004 y1 = cj[2] - ci[2]; 2005 den = x1 * x1 + y1 * y1; 2006 2007 if (den >= Mat.eps) { 2008 lbda = (x0 * x1 + y0 * y1) / den; 2009 2010 if (lbda < 0.0) { 2011 lbda = 0.0; 2012 } else if (lbda > 1.0) { 2013 lbda = 1.0; 2014 } 2015 2016 x0 = x0 - lbda * x1; 2017 y0 = y0 - lbda * y1; 2018 d = x0 * x0 + y0 * y0; 2019 } else { 2020 lbda = 0.0; 2021 d = x0 * x0 + y0 * y0; 2022 } 2023 2024 if (d > dist) { 2025 dist = d; 2026 f = k; 2027 } 2028 } 2029 return [Math.sqrt(dist), f]; 2030 }, 2031 2032 /** 2033 * RDP() is a private subroutine of {@link JXG.Math.Numerics#RamerDouglasPeuker}. 2034 * It runs recursively through the point set and searches the 2035 * point which has the largest distance from the line between the first point and 2036 * the last point. If the distance from the line is greater than eps, this point is 2037 * included in our new point set otherwise it is discarded. 2038 * If it is taken, we recursively apply the subroutine to the point set before 2039 * and after the chosen point. 2040 * @param {Array} pts Array of {@link JXG.Coords} 2041 * @param {Number} i Index of an element of pts 2042 * @param {Number} j Index of an element of pts 2043 * @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>. 2044 * @param {Array} newPts Array of {@link JXG.Coords} 2045 * @ignore 2046 * @private 2047 */ 2048 RDP = function (pts, i, j, eps, newPts) { 2049 var result = findSplit(pts, i, j); 2050 2051 if (result[0] > eps) { 2052 RDP(pts, i, result[1], eps, newPts); 2053 RDP(pts, result[1], j, eps, newPts); 2054 } else { 2055 newPts.push(pts[j]); 2056 } 2057 }; 2058 2059 len = pts.length; 2060 2061 // Search for the left most point woithout NaN coordinates 2062 i = 0; 2063 while (i < len && isNaN(pts[i].scrCoords[1] + pts[i].scrCoords[2])) { 2064 i += 1; 2065 } 2066 // Search for the right most point woithout NaN coordinates 2067 k = len - 1; 2068 while (k > i && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])) { 2069 k -= 1; 2070 } 2071 2072 // Only proceed if something is left 2073 if (!(i > k || i === len)) { 2074 newPts[0] = pts[i]; 2075 RDP(pts, i, k, eps, newPts); 2076 } 2077 2078 return newPts; 2079 } 2080 }; 2081 2082 return Mat.Numerics; 2083 }); 2084