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