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 * @author graphjs 45 */ 46 47 define(['utils/type', 'math/math'], function (Type, Mat) { 48 49 "use strict"; 50 51 // Predefined butcher tableaus for the common Runge-Kutta method (fourth order), Heun method (second order), and Euler method (first order). 52 var predefinedButcher = { 53 rk4: { 54 s: 4, 55 A: [ 56 [ 0, 0, 0, 0], 57 [0.5, 0, 0, 0], 58 [ 0, 0.5, 0, 0], 59 [ 0, 0, 1, 0] 60 ], 61 b: [1.0 / 6.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 6.0], 62 c: [0, 0.5, 0.5, 1] 63 }, 64 heun: { 65 s: 2, 66 A: [ 67 [0, 0], 68 [1, 0] 69 ], 70 b: [0.5, 0.5], 71 c: [0, 1] 72 }, 73 euler: { 74 s: 1, 75 A: [ 76 [0] 77 ], 78 b: [1], 79 c: [0] 80 } 81 }; 82 83 /** 84 * The JXG.Math.Numerics namespace holds numerical algorithms, constants, and variables. 85 * @name JXG.Math.Geometry 86 * @namespace 87 */ 88 Mat.Numerics = { 89 /** 90 * Solves a system of linear equations given by A and b using the Gauss-Jordan-elimination. 91 * The algorithm runs in-place. I.e. the entries of A and b are changed. 92 * @param {Array} A Square matrix represented by an array of rows, containing the coefficients of the lineare equation system. 93 * @param {Array} b A vector containing the linear equation system's right hand side. 94 * @throws {Error} If a non-square-matrix is given or if b has not the right length or A's rank is not full. 95 * @returns {Array} A vector that solves the linear equation system. 96 */ 97 Gauss: function (A, b) { 98 var i, j, k, 99 // copy the matrix to prevent changes in the original 100 Acopy, 101 // solution vector, to prevent changing b 102 x, 103 eps = Mat.eps, 104 // number of columns of A 105 n = A.length > 0 ? A[0].length : 0; 106 107 if ((n !== b.length) || (n !== A.length)) { 108 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."); 109 } 110 111 // initialize solution vector 112 Acopy = []; 113 x = b.slice(0, n); 114 115 for (i = 0; i < n; i++) { 116 Acopy[i] = A[i].slice(0, n); 117 } 118 119 // Gauss-Jordan-elimination 120 for (j = 0; j < n; j++) { 121 for (i = n - 1; i > j; i--) { 122 // Is the element which is to eliminate greater than zero? 123 if (Math.abs(Acopy[i][j]) > eps) { 124 // Equals pivot element zero? 125 if (Math.abs(Acopy[j][j]) < eps) { 126 // At least numerically, so we have to exchange the rows 127 Type.swap(Acopy, i, j); 128 Type.swap(x, i, j); 129 } else { 130 // Saves the L matrix of the LR-decomposition. unnecessary. 131 Acopy[i][j] /= Acopy[j][j]; 132 // Transform right-hand-side b 133 x[i] -= Acopy[i][j] * x[j]; 134 135 // subtract the multiple of A[i][j] / A[j][j] of the j-th row from the i-th. 136 for (k = j + 1; k < n; k++) { 137 Acopy[i][k] -= Acopy[i][j] * Acopy[j][k]; 138 } 139 } 140 } 141 } 142 143 // The absolute values of all coefficients below the j-th row in the j-th column are smaller than JXG.Math.eps. 144 if (Math.abs(Acopy[j][j]) < eps) { 145 throw new Error("JXG.Math.Numerics.Gauss(): The given matrix seems to be singular."); 146 } 147 } 148 149 this.backwardSolve(Acopy, x, true); 150 151 return x; 152 }, 153 154 /** 155 * Solves a system of linear equations given by the right triangular matrix R and vector b. 156 * @param {Array} R Right triangular matrix represented by an array of rows. All entries a_(i,j) with i < j are ignored. 157 * @param {Array} b Right hand side of the linear equation system. 158 * @param {Boolean} [canModify=false] If true, the right hand side vector is allowed to be changed by this method. 159 * @returns {Array} An array representing a vector that solves the system of linear equations. 160 */ 161 backwardSolve: function (R, b, canModify) { 162 var x, m, n, i, j; 163 164 if (canModify) { 165 x = b; 166 } else { 167 x = b.slice(0, b.length); 168 } 169 170 // m: number of rows of R 171 // n: number of columns of R 172 m = R.length; 173 n = R.length > 0 ? R[0].length : 0; 174 175 for (i = m - 1; i >= 0; i--) { 176 for (j = n - 1; j > i; j--) { 177 x[i] -= R[i][j] * x[j]; 178 } 179 x[i] /= R[i][i]; 180 } 181 182 return x; 183 }, 184 185 /** 186 * @private 187 * Gauss-Bareiss algorithm to compute the 188 * determinant of matrix without fractions. 189 * See H. Cohen, "A course in computational 190 * algebraic number thoery". 191 */ 192 gaussBareiss: function (mat) { 193 var k, c, s, i, j, p, n, M, t, 194 eps = Mat.eps; 195 196 n = mat.length; 197 198 if (n <= 0) { 199 return 0; 200 } 201 202 if (mat[0].length < n) { 203 n = mat[0].length; 204 } 205 206 // Copy the input matrix to M 207 M = []; 208 209 for (i = 0; i < n; i++) { 210 M[i] = mat[i].slice(0, n); 211 } 212 213 c = 1; 214 s = 1; 215 216 for (k = 0; k < n - 1; k++) { 217 p = M[k][k]; 218 219 // Pivot step 220 if (Math.abs(p) < eps) { 221 for (i = 0; i < n; i++) { 222 if (Math.abs(M[i][k]) >= eps) { 223 break; 224 } 225 } 226 227 // No nonzero entry found in column k -> det(M) = 0 228 if (i === n) { 229 return 0.0; 230 } 231 232 // swap row i and k partially 233 for (j = k; j < n; j++) { 234 t = M[i][j]; 235 M[i][j] = M[k][j]; 236 M[k][j] = t; 237 } 238 s = -s; 239 p = M[k][k]; 240 } 241 242 // Main step 243 for (i = k + 1; i < n; i++) { 244 for (j = k + 1; j < n; j++) { 245 t = p * M[i][j] - M[i][k] * M[k][j]; 246 M[i][j] = t / c; 247 } 248 } 249 250 c = p; 251 } 252 253 return s * M[n - 1][n - 1]; 254 }, 255 256 /** 257 * Computes the determinant of a square nxn matrix with the 258 * Gauss-Bareiss algorithm. 259 * @param {Array} mat Matrix. 260 * @returns {Number} The determinant pf the matrix mat. 261 * The empty matrix returns 0. 262 */ 263 det: function (mat) { 264 var n = mat.length; 265 266 if (n === 2 && mat[0].length === 2) { 267 return mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1]; 268 } 269 270 return this.gaussBareiss(mat); 271 }, 272 273 /** 274 * Compute the Eigenvalues and Eigenvectors of a symmetric 3x3 matrix with the Jacobi method 275 * Adaption of a FORTRAN program by Ed Wilson, Dec. 25, 1990 276 * @param {Array} Ain A symmetric 3x3 matrix. 277 * @returns {Array} [A,V] the matrices A and V. The diagonal of A contains the Eigenvalues, V contains the Eigenvectors. 278 */ 279 Jacobi: function (Ain) { 280 var i, j, k, aa, si, co, tt, ssum, amax, 281 eps = Mat.eps, 282 sum = 0.0, 283 n = Ain.length, 284 V = [ 285 [0, 0, 0], 286 [0, 0, 0], 287 [0, 0, 0] 288 ], 289 A = [ 290 [0, 0, 0], 291 [0, 0, 0], 292 [0, 0, 0] 293 ], 294 nloops = 0; 295 296 // Initialization. Set initial Eigenvectors. 297 for (i = 0; i < n; i++) { 298 for (j = 0; j < n; j++) { 299 V[i][j] = 0.0; 300 A[i][j] = Ain[i][j]; 301 sum += Math.abs(A[i][j]); 302 } 303 V[i][i] = 1.0; 304 } 305 306 // Trivial problems 307 if (n === 1) { 308 return [A, V]; 309 } 310 311 if (sum <= 0.0) { 312 return [A, V]; 313 } 314 315 sum /= (n * n); 316 317 // Reduce matrix to diagonal 318 do { 319 ssum = 0.0; 320 amax = 0.0; 321 for (j = 1; j < n; j++) { 322 for (i = 0; i < j; i++) { 323 // Check if A[i][j] is to be reduced 324 aa = Math.abs(A[i][j]); 325 326 if (aa > amax) { 327 amax = aa; 328 } 329 330 ssum += aa; 331 332 if (aa >= eps) { 333 // calculate rotation angle 334 aa = Math.atan2(2.0 * A[i][j], A[i][i] - A[j][j]) * 0.5; 335 si = Math.sin(aa); 336 co = Math.cos(aa); 337 338 // Modify 'i' and 'j' columns 339 for (k = 0; k < n; k++) { 340 tt = A[k][i]; 341 A[k][i] = co * tt + si * A[k][j]; 342 A[k][j] = -si * tt + co * A[k][j]; 343 tt = V[k][i]; 344 V[k][i] = co * tt + si * V[k][j]; 345 V[k][j] = -si * tt + co * V[k][j]; 346 } 347 348 // Modify diagonal terms 349 A[i][i] = co * A[i][i] + si * A[j][i]; 350 A[j][j] = -si * A[i][j] + co * A[j][j]; 351 A[i][j] = 0.0; 352 353 // Make 'A' matrix symmetrical 354 for (k = 0; k < n; k++) { 355 A[i][k] = A[k][i]; 356 A[j][k] = A[k][j]; 357 } 358 // A[i][j] made zero by rotation 359 } 360 } 361 } 362 nloops += 1; 363 } while (Math.abs(ssum) / sum > eps && nloops < 2000); 364 365 return [A, V]; 366 }, 367 368 /** 369 * Calculates the integral of function f over interval using Newton-Cotes-algorithm. 370 * @param {Array} interval The integration interval, e.g. [0, 3]. 371 * @param {function} f A function which takes one argument of type number and returns a number. 372 * @param {Object} [config] The algorithm setup. Accepted properties are number_of_nodes of type number and integration_type 373 * with value being either 'trapez', 'simpson', or 'milne'. 374 * @param {Number} [config.number_of_nodes=28] 375 * @param {String} [config.integration_type='milne'] Possible values are 'milne', 'simpson', 'trapez' 376 * @returns {Number} Integral value of f over interval 377 * @throws {Error} If config.number_of_nodes doesn't match config.integration_type an exception is thrown. If you want to use 378 * simpson rule respectively milne rule config.number_of_nodes must be dividable by 2 respectively 4. 379 * @example 380 * function f(x) { 381 * return x*x; 382 * } 383 * 384 * // calculates integral of <tt>f</tt> from 0 to 2. 385 * var area1 = JXG.Math.Numerics.NewtonCotes([0, 2], f); 386 * 387 * // the same with an anonymous function 388 * var area2 = JXG.Math.Numerics.NewtonCotes([0, 2], function (x) { return x*x; }); 389 * 390 * // use trapez rule with 16 nodes 391 * var area3 = JXG.Math.Numerics.NewtonCotes([0, 2], f, 392 * {number_of_nodes: 16, intergration_type: 'trapez'}); 393 */ 394 NewtonCotes: function (interval, f, config) { 395 var evaluation_point, i, number_of_intervals, 396 integral_value = 0.0, 397 number_of_nodes = config && typeof config.number_of_nodes === 'number' ? config.number_of_nodes : 28, 398 available_types = {trapez: true, simpson: true, milne: true}, 399 integration_type = config && config.integration_type && available_types.hasOwnProperty(config.integration_type) && available_types[config.integration_type] ? config.integration_type : 'milne', 400 step_size = (interval[1] - interval[0]) / number_of_nodes; 401 402 switch (integration_type) { 403 case 'trapez': 404 integral_value = (f(interval[0]) + f(interval[1])) * 0.5; 405 evaluation_point = interval[0]; 406 407 for (i = 0; i < number_of_nodes - 1; i++) { 408 evaluation_point += step_size; 409 integral_value += f(evaluation_point); 410 } 411 412 integral_value *= step_size; 413 break; 414 case 'simpson': 415 if (number_of_nodes % 2 > 0) { 416 throw new Error("JSXGraph: INT_SIMPSON requires config.number_of_nodes dividable by 2."); 417 } 418 419 number_of_intervals = number_of_nodes / 2.0; 420 integral_value = f(interval[0]) + f(interval[1]); 421 evaluation_point = interval[0]; 422 423 for (i = 0; i < number_of_intervals - 1; i++) { 424 evaluation_point += 2.0 * step_size; 425 integral_value += 2.0 * f(evaluation_point); 426 } 427 428 evaluation_point = interval[0] - step_size; 429 430 for (i = 0; i < number_of_intervals; i++) { 431 evaluation_point += 2.0 * step_size; 432 integral_value += 4.0 * f(evaluation_point); 433 } 434 435 integral_value *= step_size / 3.0; 436 break; 437 default: 438 if (number_of_nodes % 4 > 0) { 439 throw new Error("JSXGraph: Error in INT_MILNE: config.number_of_nodes must be a multiple of 4"); 440 } 441 442 number_of_intervals = number_of_nodes * 0.25; 443 integral_value = 7.0 * (f(interval[0]) + f(interval[1])); 444 evaluation_point = interval[0]; 445 446 for (i = 0; i < number_of_intervals - 1; i++) { 447 evaluation_point += 4.0 * step_size; 448 integral_value += 14.0 * f(evaluation_point); 449 } 450 451 evaluation_point = interval[0] - 3.0 * step_size; 452 453 for (i = 0; i < number_of_intervals; i++) { 454 evaluation_point += 4.0 * step_size; 455 integral_value += 32.0 * (f(evaluation_point) + f(evaluation_point + 2 * step_size)); 456 } 457 458 evaluation_point = interval[0] - 2.0 * step_size; 459 460 for (i = 0; i < number_of_intervals; i++) { 461 evaluation_point += 4.0 * step_size; 462 integral_value += 12.0 * f(evaluation_point); 463 } 464 465 integral_value *= 2.0 * step_size / 45.0; 466 } 467 return integral_value; 468 }, 469 470 /** 471 * Integral of function f over interval. 472 * @param {Array} interval The integration interval, e.g. [0, 3]. 473 * @param {function} f A function which takes one argument of type number and returns a number. 474 * @returns {Number} The value of the integral of f over interval 475 * @see JXG.Math.Numerics.NewtonCotes 476 */ 477 I: function (interval, f) { 478 return this.NewtonCotes(interval, f, {number_of_nodes: 16, integration_type: 'milne'}); 479 }, 480 481 /** 482 * Newton's method to find roots of a funtion in one variable. 483 * @param {function} f We search for a solution of f(x)=0. 484 * @param {Number} x initial guess for the root, i.e. start value. 485 * @param {Object} context optional object that is treated as "this" in the function body. This is useful if 486 * the function is a method of an object and contains a reference to its parent object via "this". 487 * @returns {Number} A root of the function f. 488 */ 489 Newton: function (f, x, context) { 490 var df, 491 i = 0, 492 h = Mat.eps, 493 newf = f.apply(context, [x]), 494 nfev = 1; 495 496 // For compatibility 497 if (Type.isArray(x)) { 498 x = x[0]; 499 } 500 501 while (i < 50 && Math.abs(newf) > h) { 502 df = this.D(f, context)(x); 503 nfev += 2; 504 505 if (Math.abs(df) > h) { 506 x -= newf / df; 507 } else { 508 x += (Math.random() * 0.2 - 1.0); 509 } 510 511 newf = f.apply(context, [x]); 512 nfev += 1; 513 i += 1; 514 } 515 516 return x; 517 }, 518 519 /** 520 * Abstract method to find roots of univariate functions. 521 * @param {function} f We search for a solution of f(x)=0. 522 * @param {Number} x initial guess for the root, i.e. starting value. 523 * @param {Object} context optional object that is treated as "this" in the function body. This is useful if 524 * the function is a method of an object and contains a reference to its parent object via "this". 525 * @returns {Number} A root of the function f. 526 */ 527 root: function (f, x, context) { 528 return this.fzero(f, x, context); 529 }, 530 531 532 /** 533 * Compute an intersection of the curves c1 and c2 534 * with a generalized Newton method. 535 * We want to find values t1, t2 such that 536 * c1(t1) = c2(t2), i.e. 537 * (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2)) = (0,0). 538 * We set 539 * (e,f) := (c1_x(t1)-c2_x(t2),c1_y(t1)-c2_y(t2)) 540 * 541 * The Jacobian J is defined by 542 * J = (a, b) 543 * (c, d) 544 * where 545 * a = c1_x'(t1) 546 * b = -c2_x'(t2) 547 * c = c1_y'(t1) 548 * d = -c2_y'(t2) 549 * 550 * The inverse J^(-1) of J is equal to 551 * (d, -b)/ 552 * (-c, a) / (ad-bc) 553 * 554 * Then, (t1new, t2new) := (t1,t2) - J^(-1)*(e,f). 555 * If the function meetCurveCurve possesses the properties 556 * t1memo and t2memo then these are taken as start values 557 * for the Newton algorithm. 558 * After stopping of the Newton algorithm the values of t1 and t2 are stored in 559 * t1memo and t2memo. 560 * 561 * @param {JXG.Curve} c1 Curve, Line or Circle 562 * @param {JXG.Curve} c2 Curve, Line or Circle 563 * @param {Number} t1ini start value for t1 564 * @param {Number} t2ini start value for t2 565 * @returns {JXG.Coords} intersection point 566 */ 567 generalizedNewton: function (c1, c2, t1ini, t2ini) { 568 var t1, t2, 569 a, b, c, d, disc, 570 e, f, F, 571 D00, D01, 572 D10, D11, 573 count = 0; 574 575 if (this.generalizedNewton.t1memo) { 576 t1 = this.generalizedNewton.t1memo; 577 t2 = this.generalizedNewton.t2memo; 578 } else { 579 t1 = t1ini; 580 t2 = t2ini; 581 } 582 583 e = c1.X(t1) - c2.X(t2); 584 f = c1.Y(t1) - c2.Y(t2); 585 F = e * e + f * f; 586 587 D00 = this.D(c1.X, c1); 588 D01 = this.D(c2.X, c2); 589 D10 = this.D(c1.Y, c1); 590 D11 = this.D(c2.Y, c2); 591 592 while (F > Mat.eps && count < 10) { 593 a = D00(t1); 594 b = -D01(t2); 595 c = D10(t1); 596 d = -D11(t2); 597 disc = a * d - b * c; 598 t1 -= (d * e - b * f) / disc; 599 t2 -= (a * f - c * e) / disc; 600 e = c1.X(t1) - c2.X(t2); 601 f = c1.Y(t1) - c2.Y(t2); 602 F = e * e + f * f; 603 count += 1; 604 } 605 606 this.generalizedNewton.t1memo = t1; 607 this.generalizedNewton.t2memo = t2; 608 609 if (Math.abs(t1) < Math.abs(t2)) { 610 return [c1.X(t1), c1.Y(t1)]; 611 } 612 613 return [c2.X(t2), c2.Y(t2)]; 614 }, 615 616 /** 617 * Returns the Lagrange polynomials for curves with equidistant nodes, see 618 * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation, 619 * SIAM Review, Vol 46, No 3, (2004) 501-517. 620 * The graph of the parametric curve [x(t),y(t)] runs through the given points. 621 * @param {Array} p Array of JXG.Points 622 * @returns {Array} An array consisting of two functions x(t), y(t) which define a parametric curve 623 * f(t) = (x(t), y(t)) and two numbers x1 and x2 defining the curve's domain. x1 always equals zero. 624 */ 625 Neville: function (p) { 626 var w = [], 627 /** @ignore */ 628 makeFct = function (fun) { 629 return function (t, suspendedUpdate) { 630 var i, d, s, 631 bin = Mat.binomial, 632 len = p.length, 633 len1 = len - 1, 634 num = 0.0, 635 denom = 0.0; 636 637 if (!suspendedUpdate) { 638 s = 1; 639 for (i = 0; i < len; i++) { 640 w[i] = bin(len1, i) * s; 641 s *= (-1); 642 } 643 } 644 645 d = t; 646 647 for (i = 0; i < len; i++) { 648 if (d === 0) { 649 return p[i][fun](); 650 } 651 s = w[i] / d; 652 d -= 1; 653 num += p[i][fun]() * s; 654 denom += s; 655 } 656 return num / denom; 657 }; 658 }, 659 660 xfct = makeFct('X'), 661 yfct = makeFct('Y'); 662 663 return [xfct, yfct, 0, function () { 664 return p.length - 1; 665 }]; 666 }, 667 668 /** 669 * Calculates second derivatives at the given knots. 670 * @param {Array} x x values of knots 671 * @param {Array} y y values of knots 672 * @returns {Array} Second derivatives of the interpolated function at the knots. 673 * @see #splineEval 674 */ 675 splineDef: function (x, y) { 676 var pair, i, l, 677 n = Math.min(x.length, y.length), 678 diag = [], 679 z = [], 680 data = [], 681 dx = [], 682 delta = [], 683 F = []; 684 685 if (n === 2) { 686 return [0, 0]; 687 } 688 689 for (i = 0; i < n; i++) { 690 pair = {X: x[i], Y: y[i]}; 691 data.push(pair); 692 } 693 data.sort(function (a, b) { 694 return a.X - b.X; 695 }); 696 for (i = 0; i < n; i++) { 697 x[i] = data[i].X; 698 y[i] = data[i].Y; 699 } 700 701 for (i = 0; i < n - 1; i++) { 702 dx.push(x[i + 1] - x[i]); 703 } 704 for (i = 0; i < n - 2; i++) { 705 delta.push(6 * (y[i + 2] - y[i + 1]) / (dx[i + 1]) - 6 * (y[i + 1] - y[i]) / (dx[i])); 706 } 707 708 // ForwardSolve 709 diag.push(2 * (dx[0] + dx[1])); 710 z.push(delta[0]); 711 712 for (i = 0; i < n - 3; i++) { 713 l = dx[i + 1] / diag[i]; 714 diag.push(2 * (dx[i + 1] + dx[i + 2]) - l * dx[i + 1]); 715 z.push(delta[i + 1] - l * z[i]); 716 } 717 718 // BackwardSolve 719 F[n - 3] = z[n - 3] / diag[n - 3]; 720 for (i = n - 4; i >= 0; i--) { 721 F[i] = (z[i] - (dx[i + 1] * F[i + 1])) / diag[i]; 722 } 723 724 // Generate f''-Vector 725 for (i = n - 3; i >= 0; i--) { 726 F[i + 1] = F[i]; 727 } 728 729 // natural cubic spline 730 F[0] = 0; 731 F[n - 1] = 0; 732 733 return F; 734 }, 735 736 /** 737 * Evaluate points on spline. 738 * @param {Number,Array} x0 A single float value or an array of values to evaluate 739 * @param {Array} x x values of knots 740 * @param {Array} y y values of knots 741 * @param {Array} F Second derivatives at knots, calculated by {@link #splineDef} 742 * @see #splineDef 743 * @returns {Number,Array} A single value or an array, depending on what is given as x0. 744 */ 745 splineEval: function (x0, x, y, F) { 746 var i, j, a, b, c, d, x_, 747 n = Math.min(x.length, y.length), 748 l = 1, 749 asArray = false, 750 y0 = []; 751 752 // number of points to be evaluated 753 if (Type.isArray(x0)) { 754 l = x0.length; 755 asArray = true; 756 } else { 757 x0 = [x0]; 758 } 759 760 for (i = 0; i < l; i++) { 761 // is x0 in defining interval? 762 if ((x0[i] < x[0]) || (x[i] > x[n - 1])) { 763 return NaN; 764 } 765 766 // determine part of spline in which x0 lies 767 for (j = 1; j < n; j++) { 768 if (x0[i] <= x[j]) { 769 break; 770 } 771 } 772 773 j -= 1; 774 775 // we're now in the j-th partial interval, i.e. x[j] < x0[i] <= x[j+1]; 776 // determine the coefficients of the polynomial in this interval 777 a = y[j]; 778 b = (y[j + 1] - y[j]) / (x[j + 1] - x[j]) - (x[j + 1] - x[j]) / 6 * (F[j + 1] + 2 * F[j]); 779 c = F[j] / 2; 780 d = (F[j + 1] - F[j]) / (6 * (x[j + 1] - x[j])); 781 // evaluate x0[i] 782 x_ = x0[i] - x[j]; 783 //y0.push(a + b*x_ + c*x_*x_ + d*x_*x_*x_); 784 y0.push(a + (b + (c + d * x_) * x_) * x_); 785 } 786 787 if (asArray) { 788 return y0; 789 } 790 791 return y0[0]; 792 }, 793 794 /** 795 * Generate a string containing the function term of a polynomial. 796 * @param {Array} coeffs Coefficients of the polynomial. The position i belongs to x^i. 797 * @param {Number} deg Degree of the polynomial 798 * @param {String} varname Name of the variable (usually 'x') 799 * @param {Number} prec Precision 800 * @returns {String} A string containg the function term of the polynomial. 801 */ 802 generatePolynomialTerm: function (coeffs, deg, varname, prec) { 803 var i, t = []; 804 805 for (i = deg; i >= 0; i--) { 806 t = t.concat(['(', coeffs[i].toPrecision(prec), ')']); 807 808 if (i > 1) { 809 t = t.concat(['*', varname, '<sup>', i, '<', '/sup> + ']); 810 } else if (i === 1) { 811 t = t.concat(['*', varname, ' + ']); 812 } 813 } 814 815 return t.join(''); 816 }, 817 818 /** 819 * Computes the polynomial through a given set of coordinates in Lagrange form. 820 * Returns the Lagrange polynomials, see 821 * Jean-Paul Berrut, Lloyd N. Trefethen: Barycentric Lagrange Interpolation, 822 * SIAM Review, Vol 46, No 3, (2004) 501-517. 823 * @param {Array} p Array of JXG.Points 824 * @returns {function} A function of one parameter which returns the value of the polynomial, whose graph runs through the given points. 825 */ 826 lagrangePolynomial: function (p) { 827 var w = [], 828 /** @ignore */ 829 fct = function (x, suspendedUpdate) { 830 var i, j, k, xi, s, M, 831 len = p.length, 832 num = 0, 833 denom = 0; 834 835 if (!suspendedUpdate) { 836 for (i = 0; i < len; i++) { 837 w[i] = 1.0; 838 xi = p[i].X(); 839 840 for (k = 0; k < len; k++) { 841 if (k !== i) { 842 w[i] *= (xi - p[k].X()); 843 } 844 } 845 846 w[i] = 1 / w[i]; 847 } 848 849 M = []; 850 851 for (j = 0; j < len; j++) { 852 M.push([1]); 853 } 854 } 855 856 for (i = 0; i < len; i++) { 857 xi = p[i].X(); 858 859 if (x === xi) { 860 return p[i].Y(); 861 } 862 863 s = w[i] / (x - xi); 864 denom += s; 865 num += s * p[i].Y(); 866 } 867 868 return num / denom; 869 }; 870 871 fct.getTerm = function () { 872 return ''; 873 }; 874 875 return fct; 876 }, 877 878 /** 879 * Computes the cubic Catmull-Rom spline curve through a given set of points. The curve 880 * is uniformly parametrized. 881 * Two artificial control points at the beginning and the end are added. 882 * @param {Array} points Array consisting of JXG.Points. 883 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 884 * which return the x resp. y coordinates of the Catmull-Rom-spline curve in t, a zero value, and a function simply 885 * returning the length of the points array minus three. 886 */ 887 CatmullRomSpline: function (points) { 888 var p, 889 coeffs = [], 890 // control point at the beginning and at the end 891 first = {}, 892 last = {}, 893 894 /** @ignore */ 895 makeFct = function (which) { 896 return function (t, suspendedUpdate) { 897 var s, c, 898 len = points.length; 899 900 if (len < 2) { 901 return NaN; 902 } 903 904 if (!suspendedUpdate) { 905 first[which] = function () { 906 return 2 * points[0][which]() - points[1][which](); 907 }; 908 909 last[which] = function () { 910 return 2 * points[len - 1][which]() - points[len - 2][which](); 911 }; 912 913 p = [first].concat(points, [last]); 914 coeffs[which] = []; 915 916 for (s = 0; s < len - 1; s++) { 917 coeffs[which][s] = [ 918 2 * p[s + 1][which](), 919 -p[s][which]() + p[s + 2][which](), 920 2 * p[s][which]() - 5 * p[s + 1][which]() + 4 * p[s + 2][which]() - p[s + 3][which](), 921 -p[s][which]() + 3 * p[s + 1][which]() - 3 * p[s + 2][which]() + p[s + 3][which]() 922 ]; 923 } 924 } 925 926 len += 2; // add the two control points 927 928 if (isNaN(t)) { 929 return NaN; 930 } 931 // This is necessary for our advanced plotting algorithm: 932 if (t < 0) { 933 return p[1][which](); 934 } 935 936 if (t >= len - 3) { 937 return p[len - 2][which](); 938 } 939 940 s = Math.floor(t); 941 942 if (s === t) { 943 return p[s][which](); 944 } 945 946 t -= s; 947 c = coeffs[which][s]; 948 949 return 0.5 * (((c[3] * t + c[2]) * t + c[1]) * t + c[0]); 950 }; 951 }; 952 953 return [makeFct('X'), makeFct('Y'), 0, 954 function () { 955 return points.length - 1; 956 }]; 957 }, 958 959 960 /** 961 * Computes the regression polynomial of a given degree through a given set of coordinates. 962 * Returns the regression polynomial function. 963 * @param {Number,function,Slider} degree number, function or slider. 964 * Either 965 * @param {Array} dataX Array containing either the x-coordinates of the data set or both coordinates in 966 * an array of {@link JXG.Point}s or {@link JXG.Coords}. In the latter case, the <tt>dataY</tt> parameter will be ignored. 967 * @param {Array} dataY Array containing the y-coordinates of the data set, 968 * @returns {function} A function of one parameter which returns the value of the regression polynomial of the given degree. 969 * It possesses the method getTerm() which returns the string containing the function term of the polynomial. 970 */ 971 regressionPolynomial: function (degree, dataX, dataY) { 972 var coeffs, deg, dX, dY, inputType, fct, 973 term = ''; 974 975 // Slider 976 if (Type.isPoint(degree) && typeof degree.Value === 'function') { 977 /** @ignore */ 978 deg = function () { 979 return degree.Value(); 980 }; 981 // function 982 } else if (Type.isFunction(degree)) { 983 deg = degree; 984 // number 985 } else if (Type.isNumber(degree)) { 986 /** @ignore */ 987 deg = function () { 988 return degree; 989 }; 990 } else { 991 throw new Error("JSXGraph: Can't create regressionPolynomial from degree of type'" + (typeof degree) + "'."); 992 } 993 994 // Parameters degree, dataX, dataY 995 if (arguments.length === 3 && Type.isArray(dataX) && Type.isArray(dataY)) { 996 inputType = 0; 997 // Parameters degree, point array 998 } else if (arguments.length === 2 && Type.isArray(dataX) && dataX.length > 0 && Type.isPoint(dataX[0])) { 999 inputType = 1; 1000 } else if (arguments.length === 2 && Type.isArray(dataX) && dataX.length > 0 && dataX[0].usrCoords && dataX[0].scrCoords) { 1001 inputType = 2; 1002 } else { 1003 throw new Error("JSXGraph: Can't create regressionPolynomial. Wrong parameters."); 1004 } 1005 1006 /** @ignore */ 1007 fct = function (x, suspendedUpdate) { 1008 var i, j, M, MT, y, B, c, s, d, 1009 // input data 1010 len = dataX.length; 1011 1012 d = Math.floor(deg()); 1013 1014 if (!suspendedUpdate) { 1015 // point list as input 1016 if (inputType === 1) { 1017 dX = []; 1018 dY = []; 1019 1020 for (i = 0; i < len; i++) { 1021 dX[i] = dataX[i].X(); 1022 dY[i] = dataX[i].Y(); 1023 } 1024 } 1025 1026 if (inputType === 2) { 1027 dX = []; 1028 dY = []; 1029 1030 for (i = 0; i < len; i++) { 1031 dX[i] = dataX[i].usrCoords[1]; 1032 dY[i] = dataX[i].usrCoords[2]; 1033 } 1034 } 1035 1036 // check for functions 1037 if (inputType === 0) { 1038 dX = []; 1039 dY = []; 1040 1041 for (i = 0; i < len; i++) { 1042 if (Type.isFunction(dataX[i])) { 1043 dX.push(dataX[i]()); 1044 } else { 1045 dX.push(dataX[i]); 1046 } 1047 1048 if (Type.isFunction(dataY[i])) { 1049 dY.push(dataY[i]()); 1050 } else { 1051 dY.push(dataY[i]); 1052 } 1053 } 1054 } 1055 1056 M = []; 1057 1058 for (j = 0; j < len; j++) { 1059 M.push([1]); 1060 } 1061 1062 for (i = 1; i <= d; i++) { 1063 for (j = 0; j < len; j++) { 1064 M[j][i] = M[j][i - 1] * dX[j]; 1065 } 1066 } 1067 1068 y = dY; 1069 MT = Mat.transpose(M); 1070 B = Mat.matMatMult(MT, M); 1071 c = Mat.matVecMult(MT, y); 1072 coeffs = Mat.Numerics.Gauss(B, c); 1073 term = Mat.Numerics.generatePolynomialTerm(coeffs, d, 'x', 3); 1074 } 1075 1076 // Horner's scheme to evaluate polynomial 1077 s = coeffs[d]; 1078 1079 for (i = d - 1; i >= 0; i--) { 1080 s = (s * x + coeffs[i]); 1081 } 1082 1083 return s; 1084 }; 1085 1086 fct.getTerm = function () { 1087 return term; 1088 }; 1089 1090 return fct; 1091 }, 1092 1093 /** 1094 * Computes the cubic Bezier curve through a given set of points. 1095 * @param {Array} points Array consisting of 3*k+1 {@link JXG.Points}. 1096 * The points at position k with k mod 3 = 0 are the data points, 1097 * points at position k with k mod 3 = 1 or 2 are the control points. 1098 * @returns {Array} An array consisting of two functions of one parameter t which return the 1099 * x resp. y coordinates of the Bezier curve in t, one zero value, and a third function accepting 1100 * no parameters and returning one third of the length of the points. 1101 */ 1102 bezier: function (points) { 1103 var len, flen, 1104 /** @ignore */ 1105 makeFct = function (which) { 1106 return function (t, suspendedUpdate) { 1107 var z = Math.floor(t) * 3, 1108 t0 = t % 1, 1109 t1 = 1 - t0; 1110 1111 if (!suspendedUpdate) { 1112 flen = 3 * Math.floor((points.length - 1) / 3); 1113 len = Math.floor(flen / 3); 1114 } 1115 1116 if (t < 0) { 1117 return points[0][which](); 1118 } 1119 1120 if (t >= len) { 1121 return points[flen][which](); 1122 } 1123 1124 if (isNaN(t)) { 1125 return NaN; 1126 } 1127 1128 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; 1129 }; 1130 }; 1131 1132 return [makeFct('X'), makeFct('Y'), 0, 1133 function () { 1134 return Math.floor(points.length / 3); 1135 }]; 1136 }, 1137 1138 /** 1139 * Computes the B-spline curve of order k (order = degree+1) through a given set of points. 1140 * @param {Array} points Array consisting of JXG.Points. 1141 * @param {Number} order Order of the B-spline curve. 1142 * @returns {Array} An Array consisting of four components: Two functions each of one parameter t 1143 * which return the x resp. y coordinates of the B-spline curve in t, a zero value, and a function simply 1144 * returning the length of the points array minus one. 1145 */ 1146 bspline: function (points, order) { 1147 var knots, N = [], 1148 _knotVector = function (n, k) { 1149 var j, 1150 kn = []; 1151 1152 for (j = 0; j < n + k + 1; j++) { 1153 if (j < k) { 1154 kn[j] = 0.0; 1155 } else if (j <= n) { 1156 kn[j] = j - k + 1; 1157 } else { 1158 kn[j] = n - k + 2; 1159 } 1160 } 1161 1162 return kn; 1163 }, 1164 1165 _evalBasisFuncs = function (t, kn, n, k, s) { 1166 var i, j, a, b, den, 1167 N = []; 1168 1169 if (kn[s] <= t && t < kn[s + 1]) { 1170 N[s] = 1; 1171 } else { 1172 N[s] = 0; 1173 } 1174 1175 for (i = 2; i <= k; i++) { 1176 for (j = s - i + 1; j <= s; j++) { 1177 if (j <= s - i + 1 || j < 0) { 1178 a = 0.0; 1179 } else { 1180 a = N[j]; 1181 } 1182 1183 if (j >= s) { 1184 b = 0.0; 1185 } else { 1186 b = N[j + 1]; 1187 } 1188 1189 den = kn[j + i - 1] - kn[j]; 1190 1191 if (den === 0) { 1192 N[j] = 0; 1193 } else { 1194 N[j] = (t - kn[j]) / den * a; 1195 } 1196 1197 den = kn[j + i] - kn[j + 1]; 1198 1199 if (den !== 0) { 1200 N[j] += (kn[j + i] - t) / den * b; 1201 } 1202 } 1203 } 1204 return N; 1205 }, 1206 /** @ignore */ 1207 makeFct = function (which) { 1208 return function (t, suspendedUpdate) { 1209 var y, j, s, 1210 len = points.length, 1211 n = len - 1, 1212 k = order; 1213 1214 if (n <= 0) { 1215 return NaN; 1216 } 1217 1218 if (n + 2 <= k) { 1219 k = n + 1; 1220 } 1221 1222 if (t <= 0) { 1223 return points[0][which](); 1224 } 1225 1226 if (t >= n - k + 2) { 1227 return points[n][which](); 1228 } 1229 1230 s = Math.floor(t) + k - 1; 1231 knots = _knotVector(n, k); 1232 N = _evalBasisFuncs(t, knots, n, k, s); 1233 1234 y = 0.0; 1235 for (j = s - k + 1; j <= s; j++) { 1236 if (j < len && j >= 0) { 1237 y += points[j][which]() * N[j]; 1238 } 1239 } 1240 1241 return y; 1242 }; 1243 }; 1244 1245 return [makeFct('X'), makeFct('Y'), 0, 1246 function () { 1247 return points.length - 1; 1248 }]; 1249 }, 1250 1251 /** 1252 * Numerical (symmetric) approximation of derivative. suspendUpdate is piped through, see {@link JXG.Curve#updateCurve} 1253 * and {@link JXG.Curve#hasPoint}. 1254 * @param {function} f Function in one variable to be differentiated. 1255 * @param {object} [obj] Optional object that is treated as "this" in the function body. This is useful, if the function is a 1256 * method of an object and contains a reference to its parent object via "this". 1257 * @returns {function} Derivative function of a given function f. 1258 */ 1259 D: function (f, obj) { 1260 var h = 0.00001, 1261 h2 = 1.0 / (h * 2.0); 1262 1263 if (!Type.exists(obj)) { 1264 return function (x, suspendUpdate) { 1265 return (f(x + h, suspendUpdate) - f(x - h, suspendUpdate)) * h2; 1266 }; 1267 } 1268 1269 return function (x, suspendUpdate) { 1270 return (f.apply(obj, [x + h, suspendUpdate]) - f.apply(obj, [x - h, suspendUpdate])) * h2; 1271 }; 1272 }, 1273 1274 /** 1275 * Helper function to create curve which displays Riemann sums. 1276 * Compute coordinates for the rectangles showing the Riemann sum. 1277 * @param {function} f Function, whose integral is approximated by the Riemann sum. 1278 * @param {Number} n number of rectangles. 1279 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson', or 'trapezoidal'. 1280 * @param {Number} start Left border of the approximation interval 1281 * @param {Number} end Right border of the approximation interval 1282 * @returns {Array} An array of two arrays containing the x and y coordinates for the rectangles showing the Riemann sum. This 1283 * 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 1284 * rectangles. 1285 */ 1286 riemann: function (f, n, type, start, end) { 1287 var i, x1, y1, delta1, delta, 1288 xarr = [], 1289 yarr = [], 1290 j = 0, 1291 x = start, y, 1292 sum = 0; 1293 1294 n = Math.round(n); 1295 1296 xarr[j] = x; 1297 yarr[j] = 0.0; 1298 1299 if (n > 0) { 1300 delta = (end - start) / n; 1301 // for 'lower' and 'upper' 1302 delta1 = delta * 0.01; 1303 1304 for (i = 0; i < n; i++) { 1305 if (type === 'right') { 1306 y = f(x + delta); 1307 } else if (type === 'middle') { 1308 y = f(x + delta * 0.5); 1309 } else if (type === 'left' || type === 'trapezoidal') { 1310 y = f(x); 1311 } else if (type === 'lower') { 1312 y = f(x); 1313 1314 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1315 y1 = f(x1); 1316 1317 if (y1 < y) { 1318 y = y1; 1319 } 1320 } 1321 } else if (type === 'upper') { 1322 y = f(x); 1323 1324 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1325 y1 = f(x1); 1326 1327 if (y1 > y) { 1328 y = y1; 1329 } 1330 } 1331 } else if (type === 'random') { 1332 y = f(x + delta * Math.random()); 1333 } else if (type === 'simpson') { 1334 y = (f(x) + 4 * f(x + delta * 0.5) + f(x + delta)) / 6.0; 1335 } else { 1336 y = f(x); // default is lower 1337 } 1338 1339 j += 1; 1340 xarr[j] = x; 1341 yarr[j] = y; 1342 j += 1; 1343 x += delta; 1344 1345 if (type === 'trapezoidal') { 1346 y = f(x); 1347 } 1348 1349 xarr[j] = x; 1350 yarr[j] = y; 1351 j += 1; 1352 xarr[j] = x; 1353 yarr[j] = 0.0; 1354 sum += y * delta; 1355 } 1356 } 1357 return [xarr, yarr, sum]; 1358 }, 1359 1360 /** 1361 * Approximate the integral by Riemann sums. 1362 * Compute the area described by the riemann sum rectangles. 1363 * @deprecated Replaced by JXG.Curve.Value(), see {@link JXG.Curve#riemannsum} 1364 * @param {function} f Function, whose integral is approximated by the Riemann sum. 1365 * @param {Number} n number of rectangles. 1366 * @param {String} type Type of approximation. Possible values are: 'left', 'right', 'middle', 'lower', 'upper', 'random', 'simpson' or 'trapezoidal'. 1367 * 1368 * @param {Number} start Left border of the approximation interval 1369 * @param {Number} end Right border of the approximation interval 1370 * @returns {Number} The sum of the areas of the rectangles. 1371 */ 1372 riemannsum: function (f, n, type, start, end) { 1373 var i, x1, y1, delta1, delta, y, 1374 sum = 0.0, 1375 x = start; 1376 1377 // this function looks very similar to this.riemann... maybe there is some merge potential? 1378 1379 n = Math.floor(n); 1380 1381 if (n > 0) { 1382 delta = (end - start) / n; 1383 // for 'lower' and 'upper' 1384 delta1 = delta * 0.01; 1385 1386 for (i = 0; i < n; i++) { 1387 if (type === 'right') { 1388 y = f(x + delta); 1389 } else if (type === 'middle') { 1390 y = f(x + delta * 0.5); 1391 } else if (type === 'trapezoidal') { 1392 y = 0.5 * (f(x + delta) + f(x)); 1393 } else if (type === 'left') { 1394 y = f(x); 1395 } else if (type === 'lower') { 1396 y = f(x); 1397 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1398 y1 = f(x1); 1399 1400 if (y1 < y) { 1401 y = y1; 1402 } 1403 } 1404 } else if (type === 'upper') { 1405 y = f(x); 1406 1407 for (x1 = x + delta1; x1 <= x + delta; x1 += delta1) { 1408 y1 = f(x1); 1409 1410 if (y1 > y) { 1411 y = y1; 1412 } 1413 } 1414 } else if (type === 'random') { 1415 y = f(x + delta * Math.random()); 1416 } else if (type === 'simpson') { 1417 y = (f(x) + 4 * f(x + delta * 0.5) + f(x + delta)) / 6.0; 1418 } else { 1419 y = f(x); // default is lower 1420 } 1421 1422 sum += delta * y; 1423 x += delta; 1424 } 1425 } 1426 1427 return sum; 1428 }, 1429 1430 /** 1431 * Solve initial value problems numerically using Runge-Kutta-methods. 1432 * See {@link http://en.wikipedia.org/wiki/Runge-Kutta_methods} for more information on the algorithm. 1433 * @param {object,String} butcher Butcher tableau describing the Runge-Kutta method to use. This can be either a string describing 1434 * a Runge-Kutta method with a Butcher tableau predefined in JSXGraph like 'euler', 'heun', 'rk4' or an object providing the structure 1435 * <pre> 1436 * { 1437 * s: <Number>, 1438 * A: <matrix>, 1439 * b: <Array>, 1440 * c: <Array> 1441 * } 1442 * </pre> 1443 * 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 1444 * @param {Array} x0 Initial value vector. If the problem is of one-dimensional, the initial value also has to be given in an array. 1445 * @param {Array} I Interval on which to integrate. 1446 * @param {Number} N Number of evaluation points. 1447 * @param {function} f Function describing the right hand side of the first order ordinary differential equation, i.e. if the ode 1448 * 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 1449 * vector <tt>x</tt>, and has to return a vector of the same dimension as <tt>x</tt> has. 1450 * @returns {Array} An array of vectors describing the solution of the ode on the given interval I. 1451 * @example 1452 * // A very simple autonomous system dx(t)/dt = x(t); 1453 * function f(t, x) { 1454 * return x; 1455 * } 1456 * 1457 * // Solve it with initial value x(0) = 1 on the interval [0, 2] 1458 * // with 20 evaluation points. 1459 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 1460 * 1461 * // Prepare data for plotting the solution of the ode using a curve. 1462 * var dataX = []; 1463 * var dataY = []; 1464 * var h = 0.1; // (I[1] - I[0])/N = (2-0)/20 1465 * for(var i=0; i<data.length; i++) { 1466 * dataX[i] = i*h; 1467 * dataY[i] = data[i][0]; 1468 * } 1469 * var g = board.create('curve', [dataX, dataY], {strokeWidth:'2px'}); 1470 * </pre><div id="d2432d04-4ef7-4159-a90b-a2eb8d38c4f6" style="width: 300px; height: 300px;"></div> 1471 * <script type="text/javascript"> 1472 * var board = JXG.JSXGraph.initBoard('d2432d04-4ef7-4159-a90b-a2eb8d38c4f6', {boundingbox: [-1, 5, 5, -1], axis: true, showcopyright: false, shownavigation: false}); 1473 * function f(t, x) { 1474 * // we have to copy the value. 1475 * // return x; would just return the reference. 1476 * return [x[0]]; 1477 * } 1478 * var data = JXG.Math.Numerics.rungeKutta('heun', [1], [0, 2], 20, f); 1479 * var dataX = []; 1480 * var dataY = []; 1481 * var h = 0.1; 1482 * for(var i=0; i<data.length; i++) { 1483 * dataX[i] = i*h; 1484 * dataY[i] = data[i][0]; 1485 * } 1486 * var g = board.create('curve', [dataX, dataY], {strokeColor:'red', strokeWidth:'2px'}); 1487 * </script><pre> 1488 */ 1489 rungeKutta: function (butcher, x0, I, N, f) { 1490 var e, i, j, k, l, s, 1491 x = [], 1492 y = [], 1493 h = (I[1] - I[0]) / N, 1494 t = I[0], 1495 dim = x0.length, 1496 result = [], 1497 r = 0; 1498 1499 if (Type.isString(butcher)) { 1500 butcher = predefinedButcher[butcher] || predefinedButcher.euler; 1501 } 1502 s = butcher.s; 1503 1504 // don't change x0, so copy it 1505 for (e = 0; e < dim; e++) { 1506 x[e] = x0[e]; 1507 } 1508 1509 for (i = 0; i < N; i++) { 1510 // Optimization doesn't work for ODEs plotted using time 1511 // if((i % quotient == 0) || (i == N-1)) { 1512 result[r] = []; 1513 for (e = 0; e < dim; e++) { 1514 result[r][e] = x[e]; 1515 } 1516 1517 r += 1; 1518 k = []; 1519 1520 for (j = 0; j < s; j++) { 1521 // init y = 0 1522 for (e = 0; e < dim; e++) { 1523 y[e] = 0.0; 1524 } 1525 1526 1527 // Calculate linear combination of former k's and save it in y 1528 for (l = 0; l < j; l++) { 1529 for (e = 0; e < dim; e++) { 1530 y[e] += (butcher.A[j][l]) * h * k[l][e]; 1531 } 1532 } 1533 1534 // add x(t) to y 1535 for (e = 0; e < dim; e++) { 1536 y[e] += x[e]; 1537 } 1538 1539 // calculate new k and add it to the k matrix 1540 k.push(f(t + butcher.c[j] * h, y)); 1541 } 1542 1543 // init y = 0 1544 for (e = 0; e < dim; e++) { 1545 y[e] = 0.0; 1546 } 1547 1548 for (l = 0; l < s; l++) { 1549 for (e = 0; e < dim; e++) { 1550 y[e] += butcher.b[l] * k[l][e]; 1551 } 1552 } 1553 1554 for (e = 0; e < dim; e++) { 1555 x[e] = x[e] + h * y[e]; 1556 } 1557 1558 t += h; 1559 } 1560 1561 return result; 1562 }, 1563 1564 1565 /** 1566 * Maximum number of iterations in {@link JXG.Math.Numerics#fzero} 1567 * @type Number 1568 * @default 80 1569 */ 1570 maxIterationsRoot: 80, 1571 1572 /** 1573 * Maximum number of iterations in {@link JXG.Math.Numerics#fminbr} 1574 * @type Number 1575 * @default 500 1576 */ 1577 maxIterationsMinimize: 500, 1578 1579 /** 1580 * 1581 * Find zero of an univariate function f. 1582 * @param {function} f Function, whose root is to be found 1583 * @param {Array,Number} x0 Start value or start interval enclosing the root 1584 * @param {Object} object Parent object in case f is method of it 1585 * @returns {Number} the approximation of the root 1586 * Algorithm: 1587 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 1588 * computations. M., Mir, 1980, p.180 of the Russian edition 1589 * 1590 * If x0 is an array containing lower and upper bound for the zero 1591 * algorithm 748 is applied. Otherwise, if x0 is a number, 1592 * the algorithm tries to bracket a zero of f starting from x0. 1593 * If this fails, we fall back to Newton's method. 1594 */ 1595 fzero: function (f, x0, object) { 1596 var a, b, c, 1597 fa, fb, fc, 1598 aa, blist, i, len, u, fu, 1599 prev_step, t1, cb, t2, 1600 // Actual tolerance 1601 tol_act, 1602 // Interpolation step is calculated in the form p/q; division 1603 // operations is delayed until the last moment 1604 p, q, 1605 // Step at this iteration 1606 new_step, 1607 eps = Mat.eps, 1608 maxiter = this.maxIterationsRoot, 1609 niter = 0, 1610 nfev = 0; 1611 1612 if (Type.isArray(x0)) { 1613 if (x0.length < 2) { 1614 throw new Error("JXG.Math.Numerics.fzero: length of array x0 has to be at least two."); 1615 } 1616 1617 a = x0[0]; 1618 fa = f.call(object, a); 1619 nfev += 1; 1620 b = x0[1]; 1621 fb = f.call(object, b); 1622 nfev += 1; 1623 } else { 1624 a = x0; 1625 fa = f.call(object, a); 1626 nfev += 1; 1627 1628 // Try to get b. 1629 if (a === 0) { 1630 aa = 1; 1631 } else { 1632 aa = a; 1633 } 1634 1635 blist = [0.9 * aa, 1.1 * aa, aa - 1, aa + 1, 0.5 * aa, 1.5 * aa, -aa, 2 * aa, -10 * aa, 10 * aa]; 1636 len = blist.length; 1637 1638 for (i = 0; i < len; i++) { 1639 b = blist[i]; 1640 fb = f.call(object, b); 1641 nfev += 1; 1642 1643 if (fa * fb <= 0) { 1644 break; 1645 } 1646 } 1647 if (b < a) { 1648 u = a; 1649 a = b; 1650 b = u; 1651 1652 fu = fa; 1653 fa = fb; 1654 fb = fu; 1655 } 1656 } 1657 1658 if (fa * fb > 0) { 1659 // Bracketing not successful, fall back to Newton's method or to fminbr 1660 if (Type.isArray(x0)) { 1661 return this.fminbr(f, [a, b], object); 1662 } 1663 1664 return this.Newton(f, a, object); 1665 } 1666 1667 // OK, we have enclosed a zero of f. 1668 // Now we can start Brent's method 1669 1670 c = a; 1671 fc = fa; 1672 1673 // Main iteration loop 1674 while (niter < maxiter) { 1675 // Distance from the last but one to the last approximation 1676 prev_step = b - a; 1677 1678 // Swap data for b to be the best approximation 1679 if (Math.abs(fc) < Math.abs(fb)) { 1680 a = b; 1681 b = c; 1682 c = a; 1683 1684 fa = fb; 1685 fb = fc; 1686 fc = fa; 1687 } 1688 tol_act = 2 * eps * Math.abs(b) + eps * 0.5; 1689 new_step = (c - b) * 0.5; 1690 1691 if (Math.abs(new_step) <= tol_act && Math.abs(fb) <= eps) { 1692 // Acceptable approx. is found 1693 return b; 1694 } 1695 1696 // Decide if the interpolation can be tried 1697 // If prev_step was large enough and was in true direction Interpolatiom may be tried 1698 if (Math.abs(prev_step) >= tol_act && Math.abs(fa) > Math.abs(fb)) { 1699 cb = c - b; 1700 1701 // If we have only two distinct points linear interpolation can only be applied 1702 if (a === c) { 1703 t1 = fb / fa; 1704 p = cb * t1; 1705 q = 1.0 - t1; 1706 // Quadric inverse interpolation 1707 } else { 1708 q = fa / fc; 1709 t1 = fb / fc; 1710 t2 = fb / fa; 1711 1712 p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0)); 1713 q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0); 1714 } 1715 1716 // p was calculated with the opposite sign; make p positive 1717 if (p > 0) { 1718 q = -q; 1719 // and assign possible minus to q 1720 } else { 1721 p = -p; 1722 } 1723 1724 // If b+p/q falls in [b,c] and isn't too large it is accepted 1725 // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent 1726 if (p < (0.75 * cb * q - Math.abs(tol_act * q) * 0.5) && 1727 p < Math.abs(prev_step * q * 0.5)) { 1728 new_step = p / q; 1729 } 1730 } 1731 1732 // Adjust the step to be not less than tolerance 1733 if (Math.abs(new_step) < tol_act) { 1734 if (new_step > 0) { 1735 new_step = tol_act; 1736 } else { 1737 new_step = -tol_act; 1738 } 1739 } 1740 1741 // Save the previous approx. 1742 a = b; 1743 fa = fb; 1744 b += new_step; 1745 fb = f.call(object, b); 1746 // Do step to a new approxim. 1747 nfev += 1; 1748 1749 // Adjust c for it to have a sign opposite to that of b 1750 if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) { 1751 c = a; 1752 fc = fa; 1753 } 1754 niter++; 1755 } // End while 1756 1757 return b; 1758 }, 1759 1760 1761 /** 1762 * 1763 * Find minimum of an univariate function f. 1764 * @param {function} f Function, whose minimum is to be found 1765 * @param {Array} x0 Start interval enclosing the minimum 1766 * @param {Object} context Parent object in case f is method of it 1767 * @returns {Number} the approximation of the minimum value position 1768 * Algorithm: 1769 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical 1770 * computations. M., Mir, 1980, p.180 of the Russian edition 1771 * x0 1772 **/ 1773 fminbr: function (f, x0, context) { 1774 var a, b, x, v, w, 1775 fx, fv, fw, 1776 range, middle_range, tol_act, new_step, 1777 p, q, t, ft, 1778 // Golden section ratio 1779 r = (3.0 - Math.sqrt(5.0)) * 0.5, 1780 tol = Mat.eps, 1781 sqrteps = Mat.eps, //Math.sqrt(Mat.eps), 1782 maxiter = this.maxIterationsMinimize, 1783 niter = 0, 1784 nfev = 0; 1785 1786 if (!Type.isArray(x0) || x0.length < 2) { 1787 throw new Error("JXG.Math.Numerics.fminbr: length of array x0 has to be at least two."); 1788 } 1789 1790 a = x0[0]; 1791 b = x0[1]; 1792 v = a + r * (b - a); 1793 fv = f.call(context, v); 1794 1795 // First step - always gold section 1796 nfev += 1; 1797 x = v; 1798 w = v; 1799 fx = fv; 1800 fw = fv; 1801 1802 while (niter < maxiter) { 1803 // Range over the interval in which we are looking for the minimum 1804 range = b - a; 1805 middle_range = (a + b) * 0.5; 1806 1807 // Actual tolerance 1808 tol_act = sqrteps * Math.abs(x) + tol / 3.0; 1809 1810 if (Math.abs(x - middle_range) + range * 0.5 <= 2.0 * tol_act) { 1811 // Acceptable approx. is found 1812 return x; 1813 } 1814 1815 // Obtain the golden section step 1816 new_step = r * (x < middle_range ? b - x : a - x); 1817 1818 // Decide if the interpolation can be tried. If x and w are distinct interpolatiom may be tried 1819 if (Math.abs(x - w) >= tol_act) { 1820 // Interpolation step is calculated as p/q; 1821 // division operation is delayed until last moment 1822 t = (x - w) * (fx - fv); 1823 q = (x - v) * (fx - fw); 1824 p = (x - v) * q - (x - w) * t; 1825 q = 2 * (q - t); 1826 1827 if (q > 0) { // q was calculated with the op- 1828 p = -p; // posite sign; make q positive 1829 } else { // and assign possible minus to 1830 q = -q; // p 1831 } 1832 if (Math.abs(p) < Math.abs(new_step * q) && // If x+p/q falls in [a,b] 1833 p > q * (a - x + 2 * tol_act) && // not too close to a and 1834 p < q * (b - x - 2 * tol_act)) { // b, and isn't too large 1835 new_step = p / q; // it is accepted 1836 } 1837 // If p/q is too large then the 1838 // golden section procedure can 1839 // reduce [a,b] range to more 1840 // extent 1841 } 1842 1843 // Adjust the step to be not less than tolerance 1844 if (Math.abs(new_step) < tol_act) { 1845 if (new_step > 0) { 1846 new_step = tol_act; 1847 } else { 1848 new_step = -tol_act; 1849 } 1850 } 1851 1852 // Obtain the next approximation to min 1853 // and reduce the enveloping range 1854 1855 // Tentative point for the min 1856 t = x + new_step; 1857 ft = f.call(context, t); 1858 nfev += 1; 1859 1860 // t is a better approximation 1861 if (ft <= fx) { 1862 // Reduce the range so that t would fall within it 1863 if (t < x) { 1864 b = x; 1865 } else { 1866 a = x; 1867 } 1868 1869 // Assign the best approx to x 1870 v = w; 1871 w = x; 1872 x = t; 1873 1874 fv = fw; 1875 fw = fx; 1876 fx = ft; 1877 // x remains the better approx 1878 } else { 1879 // Reduce the range enclosing x 1880 if (t < x) { 1881 a = t; 1882 } else { 1883 b = t; 1884 } 1885 1886 if (ft <= fw || w === x) { 1887 v = w; 1888 w = t; 1889 fv = fw; 1890 fw = ft; 1891 } else if (ft <= fv || v === x || v === w) { 1892 v = t; 1893 fv = ft; 1894 } 1895 } 1896 niter += 1; 1897 } 1898 1899 return x; 1900 }, 1901 1902 /** 1903 * Implements the Ramer-Douglas-Peuker algorithm. 1904 * It discards points which are not necessary from the polygonal line defined by the point array 1905 * pts. The computation is done in screen coordinates. 1906 * Average runtime is O(nlog(n)), worst case runtime is O(n^2), where n is the number of points. 1907 * @param {Array} pts Array of {@link JXG.Coords} 1908 * @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>. 1909 * @returns {Array} An array containing points which represent an apparently identical curve as the points of pts do, but contains fewer points. 1910 */ 1911 RamerDouglasPeuker: function (pts, eps) { 1912 var newPts = [], i, k, len, 1913 1914 /** 1915 * findSplit() is a subroutine of {@link JXG.Math.Numerics#RamerDouglasPeuker}. 1916 * It searches for the point between index i and j which 1917 * has the largest distance from the line between the points i and j. 1918 * @param {Array} pts Array of {@link JXG.Coords} 1919 * @param {Number} i Index of a point in pts 1920 * @param {Number} j Index of a point in pts 1921 * @ignore 1922 * @private 1923 */ 1924 findSplit = function (pts, i, j) { 1925 var d, k, ci, cj, ck, 1926 x0, y0, x1, y1, 1927 den, lbda, 1928 dist = 0, 1929 f = i; 1930 1931 if (j - i < 2) { 1932 return [-1.0, 0]; 1933 } 1934 1935 ci = pts[i].scrCoords; 1936 cj = pts[j].scrCoords; 1937 1938 if (isNaN(ci[1] + ci[2] + cj[1] + cj[2])) { 1939 return [NaN, j]; 1940 } 1941 1942 for (k = i + 1; k < j; k++) { 1943 ck = pts[k].scrCoords; 1944 x0 = ck[1] - ci[1]; 1945 y0 = ck[2] - ci[2]; 1946 x1 = cj[1] - ci[1]; 1947 y1 = cj[2] - ci[2]; 1948 den = x1 * x1 + y1 * y1; 1949 1950 if (den >= Mat.eps) { 1951 lbda = (x0 * x1 + y0 * y1) / den; 1952 1953 if (lbda < 0.0) { 1954 lbda = 0.0; 1955 } else if (lbda > 1.0) { 1956 lbda = 1.0; 1957 } 1958 1959 x0 = x0 - lbda * x1; 1960 y0 = y0 - lbda * y1; 1961 d = x0 * x0 + y0 * y0; 1962 } else { 1963 lbda = 0.0; 1964 d = x0 * x0 + y0 * y0; 1965 } 1966 1967 if (d > dist) { 1968 dist = d; 1969 f = k; 1970 } 1971 } 1972 return [Math.sqrt(dist), f]; 1973 }, 1974 1975 /** 1976 * RDP() is a private subroutine of {@link JXG.Math.Numerics#RamerDouglasPeuker}. 1977 * It runs recursively through the point set and searches the 1978 * point which has the largest distance from the line between the first point and 1979 * the last point. If the distance from the line is greater than eps, this point is 1980 * included in our new point set otherwise it is discarded. 1981 * If it is taken, we recursively apply the subroutine to the point set before 1982 * and after the chosen point. 1983 * @param {Array} pts Array of {@link JXG.Coords} 1984 * @param {Number} i Index of an element of pts 1985 * @param {Number} j Index of an element of pts 1986 * @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>. 1987 * @param {Array} newPts Array of {@link JXG.Coords} 1988 * @ignore 1989 * @private 1990 */ 1991 RDP = function (pts, i, j, eps, newPts) { 1992 var result = findSplit(pts, i, j); 1993 1994 if (result[0] > eps) { 1995 RDP(pts, i, result[1], eps, newPts); 1996 RDP(pts, result[1], j, eps, newPts); 1997 } else { 1998 newPts.push(pts[j]); 1999 } 2000 }; 2001 2002 len = pts.length; 2003 2004 // Search for the left most point woithout NaN coordinates 2005 i = 0; 2006 while (i < len && isNaN(pts[i].scrCoords[1] + pts[i].scrCoords[2])) { 2007 i += 1; 2008 } 2009 // Search for the right most point woithout NaN coordinates 2010 k = len - 1; 2011 while (k > i && isNaN(pts[k].scrCoords[1] + pts[k].scrCoords[2])) { 2012 k -= 1; 2013 } 2014 2015 // Only proceed if something is left 2016 if (!(i > k || i === len)) { 2017 newPts[0] = pts[i]; 2018 RDP(pts, i, k, eps, newPts); 2019 } 2020 2021 return newPts; 2022 } 2023 }; 2024 2025 return Mat.Numerics; 2026 }); 2027