1 /*
  2     Copyright 2008-2021
  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