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