      1 /*
      2     Copyright 2008-2022
      3         Matthias Ehmann,
      4         Michael Gerhaeuser,
      5         Carsten Miller,
      6         Alfred Wassermann
      7
      8     This file is part of JSXGraph.
      9
     10     JSXGraph is free software dual licensed under the GNU LGPL or MIT License.
     11
     12     You can redistribute it and/or modify it under the terms of the
     13
     14       * GNU Lesser General Public License as published by
     15         the Free Software Foundation, either version 3 of the License, or
     16         (at your option) any later version
     17       OR
     18       * MIT License: https://github.com/jsxgraph/jsxgraph/blob/master/LICENSE.MIT
     19
     20     JSXGraph is distributed in the hope that it will be useful,
     21     but WITHOUT ANY WARRANTY; without even the implied warranty of
     22     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     23     GNU Lesser General Public License for more details.
     24
     25     You should have received a copy of the GNU Lesser General Public License and
     26     the MIT License along with JSXGraph. If not, see <http://www.gnu.org/licenses/>
     27     and <http://opensource.org/licenses/MIT/>.
     28  */
     29
     30
     31 /*global JXG: true, define: true*/
     32 /*jslint nomen: true, plusplus: true*/
     33
     34 /* depends:
     35  jxg
     36  math/math
     37  utils/type
     38  */
     39
     40 define(['math/math'], function (Mat) {
     41
     42     "use strict";
     43
     44     /**
     45      * Functions for extrapolation of sequences. Used for finding limits of sequences which is used for curve plotting.
     46      * @name JXG.Math.Extrapolate
     47      * @exports Mat.Extrapolate as JXG.Math.Extrapolate
     48      * @namespace
     49      */
     50     Mat.Extrapolate = {
     51         upper: 15,
     52         infty: 1.e+4,
     53
     54         /**
     55          * Wynn's epsilon algorithm. Ported from the FORTRAN version in
     56          * Ernst Joachim Weniger, "Nonlinear sequence transformations for the acceleration of convergence
     57          * and the summation of divergent series", Computer Physics Reports Vol. 10, 189-371 (1989).
     58          *
     59          * @param {Number} s_n next value of sequence, i.e. n-th element of sequence
     60          * @param {Number} n index of s_n in the sequence
     61          * @param {Array} e One-dimensional array containing the extrapolation data. Has to be supplied by the calling routine.
     62          * @returns {Number} New estimate of the limit of the sequence.
     63          *
     64          * @memberof JXG.Math.Extrapolate
     65          */
     66         wynnEps: function(s_n, n, e) {
     67             var HUGE = 1.e+20,
     68                 TINY = 1.e-15,
     69                 f0 = 1, // f0 may be changed to other values, see vanden Broeck, Schwartz (1979)
     70                 f, j, aux1, aux2, diff, estlim;
     71
     72             e[n] = s_n;
     73             if (n === 0) {
     74                 estlim = s_n;
     75             } else {
     76                 aux2 = 0.0;
     77                 for (j = n; j > 0; j--) {
     78                     aux1 = aux2;
     79                     aux2 = e[j - 1];
     80                     diff = e[j] - aux2;
     81                     if (Math.abs(diff) <= TINY) {
     82                         e[j - 1] = HUGE;
     83                     } else {
     84                         f = ((n - j + 1) % 2 === 1) ? f0 : 1;
     85                         e[j - 1] = aux1 * f + 1 / diff;
     86                     }
     87                 }
     88                 estlim = e[n % 2];
     89             }
     90
     91             return estlim;
     92         },
     93
     94         // wynnRho: function(s_n, n, e) {
     95         //     var HUGE = 1.e+20,
     96         //         TINY = 1.e-15,
     97         //         j, f,
     98         //         aux1, aux2, diff, estlim;
     99
    100         //     e[n] = s_n;
    101         //     if (n === 0) {
    102         //         estlim = s_n;
    103         //     } else {
    104         //         aux2 = 0.0;
    105         //         for (j = n; j >= 1; j--) {
    106         //             aux1 = aux2;
    107         //             aux2 = e[j - 1];
    108         //             diff = e[j] - aux2;
    109         //             if (Math.abs(diff) <= TINY) {
    110         //                 e[j - 1] = HUGE;
    111         //             } else {
    112         //                 f = ((n - j + 1) % 2 === 1) ? n - j + 1  : 1;
    113         //                 e[j - 1] = aux1 + f / diff;
    114         //             }
    115         //         }
    116         //         estlim = e[n % 2];
    117         //     }
    118
    119         //     return estlim;
    120         // },
    121
    122         /**
    123          * Aitken transformation. Ported from the FORTRAN version in
    124          * Ernst Joachim Weniger, "Nonlinear sequence transformations for the acceleration of convergence
    125          * and the summation of divergent series", Computer Physics Reports Vol. 10, 189-371 (1989).
    126          *
    127          * @param {Number} s_n next value of sequence, i.e. n-th element of sequence
    128          * @param {Number} n index of s_n in the sequence
    129          * @param {Array} a One-dimensional array containing the extrapolation data. Has to be supplied by the calling routine.
    130          * @returns {Number} New estimate of the limit of the sequence.
    131          *
    132          * @memberof JXG.Math.Extrapolate
    133          */
    134         aitken: function(s_n, n, a) {
    135             var estlim,
    136                 HUGE = 1.e+20,
    137                 TINY = 1.e-15,
    138                 denom, v,
    139                 lowmax, j, m;
    140
    141             a[n] = s_n;
    142             if (n < 2) {
    143                 estlim = s_n;
    144             } else {
    145                 lowmax = n / 2;
    146                 for (j = 1; j <= lowmax; j++) {
    147                     m = n - 2 * j;
    148                     denom = a[m + 2] - 2 * a[m + 1] + a[m];
    149                     if (Math.abs(denom) < TINY) {
    150                          a[m] = HUGE;
    151                     } else {
    152                         v = a[m] - a[m + 1];
    153                         a[m] -= v * v / denom;
    154                     }
    155                 }
    156                 estlim = a[n % 2];
    157             }
    158             return estlim;
    159         },
    160
    161         /**
    162          * Iterated Brezinski transformation. Ported from the FORTRAN version in
    163          * Ernst Joachim Weniger, "Nonlinear sequence transformations for the acceleration of convergence
    164          * and the summation of divergent series", Computer Physics Reports Vol. 10, 189-371 (1989).
    165          *
    166          * @param {Number} s_n next value of sequence, i.e. n-th element of sequence
    167          * @param {Number} n index of s_n in the sequence
    168          * @param {Array} a One-dimensional array containing the extrapolation data. Has to be supplied by the calling routine.
    169          * @returns {Number} New estimate of the limit of the sequence.
    170          *
    171          * @memberof JXG.Math.Extrapolate
    172          */
    173         brezinski: function(s_n, n, a) {
    174             var estlim,
    175                 HUGE = 1.e+20,
    176                 TINY = 1.e-15,
    177                 denom,
    178                 d0, d1, d2,
    179                 lowmax, j, m;
    180
    181             a[n] = s_n;
    182             if (n < 3) {
    183                 estlim = s_n;
    184             } else {
    185                 lowmax = n / 3;
    186                 m = n;
    187                 for (j = 1; j <= lowmax; j++) {
    188                     m -= 3;
    189                     d0 = a[m + 1] - a[m];
    190                     d1 = a[m + 2] - a[m + 1];
    191                     d2 = a[m + 3] - a[m + 2];
    192                     denom = d2 * (d1 - d0) - d0 * (d2 - d1);
    193                     if (Math.abs(denom) < TINY) {
    194                         a[m] = HUGE;
    195                     } else {
    196                         a[m] = a[m + 1] - d0 * d1 * (d2 - d1) / denom;
    197                     }
    198                 }
    199                 estlim = a[n % 3];
    200             }
    201             return estlim;
    202         },
    203
    204         /**
    205          * Extrapolated iteration to approximate the value f(x_0).
    206          *
    207          * @param {Number} x0 Value for which the limit of f is to be determined. f(x0) may or may not exist.
    208          * @param {Number} h0 Initial (signed) distance from x0.
    209          * @param {Function} f Function for which the limit at x0 is to be determined
    210          * @param {String} method String to choose the method. Available values: "wynnEps", "aitken", "brezinski"
    211          * @param {Number} step_type Approximation method. step_type = 0 uses the sequence x0 + h0/n; step_type = 1 uses the sequence x0 + h0 * 2^(-n)
    212          *
    213          * @returns {Array} Array of length 3. Position 0: estimated value for f(x0), position 1: 'finite', 'infinite', or 'NaN'.
    214          * Position 2: value between 0 and 1 judging the reliability of the result (1: high, 0: not successful).
    215          *
    216          * @memberof JXG.Math.Extrapolate
    217          * @see JXG.Math.Extrapolate.limit
    218          * @see JXG.Math.Extrapolate.wynnEps
    219          * @see JXG.Math.Extrapolate.aitken
    220          * @see JXG.Math.Extrapolate.brezinski
    221          */
    222         iteration: function(x0, h0, f, method, step_type) {
    223             var n, v, w,
    224                 estlim = NaN,
    225                 diff,
    226                 r = 0.5,
    227                 E = [],
    228                 result = 'finite',
    229                 h = h0;
    230
    231             step_type = step_type || 0;
    232
    233             for (n = 1; n <= this.upper; n++) {
    234                 h = (step_type === 0) ?  h0 / (n + 1) : h * r;
    235                 v = f(x0 + h, true);
    236
    237                 w = this[method](v, n - 1, E);
    238 //console.log(n, x0 + h, v, w);
    239                 if (isNaN(w)) {
    240                     result = 'NaN';
    241                     break;
    242                 }
    243                 if (v !== 0 && w / v > this.infty) {
    244                     estlim = w;
    245                     result = 'infinite';
    246                     break;
    247                 }
    248                 diff = w - estlim;
    249                 if (Math.abs(diff) < 1.e-7) {
    250                     break;
    251                 }
    252                 estlim = w;
    253             }
    254             return [estlim, result, 1 - (n - 1) / this.upper];
    255         },
    256
    257         /**
    258          * Levin transformation. See Numerical Recipes, ed. 3.
    259          * Not yet ready for use.
    260          *
    261          * @param {Number} s_n next value of sequence, i.e. n-th element of sequence
    262          * @param {Number} n index of s_n in the sequence
    263          * @param {Array} numer One-dimensional array containing the extrapolation data for the numerator. Has to be supplied by the calling routine.
    264          * @param {Array} denom One-dimensional array containing the extrapolation data for the denominator. Has to be supplied by the calling routine.
    265          *
    266          * @memberof JXG.Math.Extrapolate
    267         */
    268         levin: function(s_n, n, omega, beta, numer, denom) {
    269             var HUGE = 1.e+20,
    270                 TINY = 1.e-15,
    271                 j,
    272                 fact, ratio, term, estlim;
    273
    274             term = 1.0 / (beta + n);
    275             numer[n] = s_n / omega;
    276             denom[n] = 1 / omega;
    277             if (n > 0) {
    278                 numer[n - 1] = numer[n] - numer[n - 1];
    279                 denom[n - 1] = denom[n] - denom[n - 1];
    280                 if (n > 1) {
    281                     ratio = (beta + n - 1) * term;
    282                     for (j = 2; j <= n; j++) {
    283                         fact = (beta + n - j) * Math.pow(ratio, j - 2) * term;
    284                         numer[n - j] = numer[n - j + 1] - fact * numer[n - j];
    285                         denom[n - j] = denom[n - j + 1] - fact * denom[n - j];
    286                         term *= ratio;
    287                     }
    288                 }
    289             }
    290             if (Math.abs(denom[0]) < TINY) {
    291                 estlim = HUGE;
    292             } else {
    293                 estlim = numer[0] / denom[0];
    294             }
    295             return estlim;
    296         },
    297
    298         iteration_levin: function(x0, h0, f, step_type) {
    299             var n, v, w,
    300                 estlim = NaN,
    301                 v_prev,
    302                 delta, diff, omega,
    303                 beta = 1,
    304                 r = 0.5,
    305                 numer = [],
    306                 denom = [],
    307                 result = 'finite',
    308                 h = h0, transform = 'u';
    309
    310             step_type = step_type || 0;
    311
    312             v_prev = f(x0 + h0, true);
    313             for (n = 1; n <= this.upper; n++) {
    314                 h = (step_type === 0) ?  h0 / (n + 1) : h * r;
    315                 v = f(x0 + h, true);
    316                 delta = v - v_prev;
    317                 if (Math.abs(delta) < 1) {
    318                     transform = 'u';
    319                 } else {
    320                     transform = 't';
    321                 }
    322                 if (transform === 'u') {
    323                     omega = (beta + n) * delta; // u transformation
    324                 } else {
    325                     omega = delta;              // t transformation
    326                 }
    327
    328                 v_prev = v;
    329                 w = this.levin(v, n - 1, omega, beta, numer, denom);
    330                 diff = w - estlim;
    331 // console.log(n, delta, transform, x0 + h, v, w, diff);
    332
    333                 if (isNaN(w)) {
    334                     result = 'NaN';
    335                     break;
    336                 }
    337                 if (v !== 0 && w / v > this.infty) {
    338                     estlim = w;
    339                     result = 'infinite';
    340                     break;
    341                 }
    342                 if (Math.abs(diff) < 1.e-7) {
    343                     break;
    344                 }
    345                 estlim = w;
    346             }
    347             return [estlim, result, 1 - (n - 1) / this.upper];
    348         },
    349
    350         /**
    351          *
    352          * @param {Number} x0 Value for which the limit of f is to be determined. f(x0) may or may not exist.
    353          * @param {Number} h0 Initial (signed) distance from x0.
    354          * @param {Function} f Function for which the limit at x0 is to be determined
    355          *
    356          * @returns {Array} Array of length 3. Position 0: estimated value for f(x0), position 1: 'finite', 'infinite', or 'NaN'.
    357          * Position 2: value between 0 and 1 judging the reliability of the result (1: high, 0: not successful).
    358          * In case that the extrapolation fails, position 1 and 2 contain 'direct' and 0.
    359          *
    360          * @example
    361          * var f1 = (x) => Math.log(x),
    362          *     f2 = (x) => Math.tan(x - Math.PI * 0.5),
    363          *     f3 = (x) => 4 / x;
    364          *
    365          * var x0 = 0.0000001;
    366          * var h = 0.1;
    367          * for (let f of [f1, f2, f3]) {
    368          *     console.log("x0=", x0, f.toString());
    369          *     console.log(JXG.Math.Extrapolate.limit(x0, h, f));
    370          *  }
    371          *
    372          * </pre><div id="JXG5e8c6a7e-eeae-43fb-a669-26b5c9e40cab" class="jxgbox" style="width: 300px; height: 300px;"></div>
    373          * <script type="text/javascript">
    374          *     (function() {
    375          *         var board = JXG.JSXGraph.initBoard('JXG5e8c6a7e-eeae-43fb-a669-26b5c9e40cab',
    376          *             {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false});
    377          *     var f1 = (x) => Math.log(x),
    378          *         f2 = (x) => Math.tan(x - Math.PI * 0.5),
    379          *         f3 = (x) => 4 / x;
    380          *
    381          *     var x0 = 0.0000001;
    382          *     var h = 0.1;
    383          *     for (let f of [f1, f2, f3]) {
    384          *         console.log("x0=", x0, f.toString());
    385          *         console.log(JXG.Math.Extrapolate.limit(x0, h, f));
    386          *      }
    387          *
    388          *     })();
    389          *
    390          * </script><pre>
    391          *
    392          *
    393          * @see JXG.Math.Extrapolate.iteration
    394          * @memberof JXG.Math.Extrapolate
    395          */
    396         limit: function(x0, h0, f) {
    397             return this.iteration_levin(x0, h0, f, 0);
    398             //return this.iteration(x0, h0, f, 'wynnEps', 1);
    399
    400             // var algs = ['wynnEps', 'levin'], //, 'wynnEps', 'levin', 'aitken', 'brezinski'],
    401             //     le = algs.length,
    402             //     i, t, res;
    403             // for (i = 0; i < le; i++) {
    404             //     for (t = 0; t < 1; t++) {
    405             //         if (algs[i] === 'levin') {
    406             //             res = this.iteration_levin(x0, h0, f, t);
    407             //         } else {
    408             //             res = this.iteration(x0, h0, f, algs[i], t);
    409             //         }
    410             //         if (res[2] > 0.6) {
    411             //             return res;
    412             //         }
    413             //         console.log(algs[i], t, res)
    414             //     }
    415             // }
    416             // return [f(x0 + Math.sign(h0) * Math.sqrt(Mat.eps)), 'direct', 0];
    417         }
    418     };
    419
    420     return Mat.Extrapolate;
    421 });
    422