1 /* 2 Copyright 2008-2023 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"; 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