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