1 /* 2 Copyright 2008-2026 3 Matthias Ehmann, 4 Carsten Miller, 5 Alfred Wassermann 6 7 This file is part of JSXGraph. 8 9 JSXGraph is free software dual licensed under the GNU LGPL or MIT License. 10 11 You can redistribute it and/or modify it under the terms of the 12 13 * GNU Lesser General Public License as published by 14 the Free Software Foundation, either version 3 of the License, or 15 (at your option) any later version 16 OR 17 * MIT License: https://github.com/jsxgraph/jsxgraph/blob/master/LICENSE.MIT 18 19 JSXGraph is distributed in the hope that it will be useful, 20 but WITHOUT ANY WARRANTY; without even the implied warranty of 21 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 22 GNU Lesser General Public License for more details. 23 24 You should have received a copy of the GNU Lesser General Public License and 25 the MIT License along with JSXGraph. If not, see <https://www.gnu.org/licenses/> 26 and <https://opensource.org/licenses/MIT/>. 27 */ 28 29 "use strict"; 30 31 import Type from "../utils/type.js"; 32 import Mat from "./math.js"; 33 import Geometry from "./geometry.js"; 34 import Numerics from "./numerics.js"; 35 import Quadtree from "./bqdt.js"; 36 37 /** 38 * Plotting of curves which are given implicitly as the set of points solving an equation 39 * <i>f(x,y) = 0</i>. 40 * <p> 41 * The main class initializes a new implicit plot instance. 42 * <p> 43 * The algorithm should be able to plot most implicit curves as long as the equations 44 * are not too complex. We are aware of the paper by Oliver Labs, 45 * <a href="https://link.springer.com/chapter/10.1007/978-1-4419-0999-2_6">A List of Challenges for Real Algebraic Plane Curve Visualization Software</a> 46 * which contains many equations where this algorithm may fail. 47 * For example, at the time being there is no attempt to detect <i>solitary points</i>. 48 * Also, it is always a trade off to find all components of the curve and 49 * keep the construction responsive. 50 * 51 * @name JXG.Math.ImplicitPlot 52 * @exports Mat.ImplicitPlot as JXG.Math.ImplicitPlot 53 * @param {Array} bbox Bounding box of the area in which solutions of the equation 54 * are determined. 55 * @param {Object} config Configuration object. Default: 56 * <pre> 57 * { 58 * resolution_out: 5, // Horizontal resolution: distance between vertical lines to search for components 59 * resolution_in: 5, // Vertical resolution to search for components 60 * max_steps: 1024, // Max number of points in one call of tracing 61 * alpha_0: 0.05, // Angle between two successive tangents: smoothness of curve 62 * 63 * tol_u0: Mat.eps, // Tolerance to find starting points for tracing. 64 * tol_newton: 1.0e-7, // Tolerance for Newton steps. 65 * tol_cusp: 0.05, // Tolerance for cusp / bifurcation detection 66 * tol_progress: 0.0001, // If two points are closer than this value, we bail out 67 * qdt_box: 0.2, // half of box size to search in qdt 68 * kappa_0: 0.2, // Inverse of planned number of Newton steps 69 * delta_0: 0.05, // Distance of predictor point to curve 70 * 71 * h_initial: 0.1, // Initial stepwidth 72 * h_critical: 0.001, // If h is below this threshold we bail out 73 * h_max: 1, // Maximal value of h (user units) 74 * loop_dist: 0.09, // Allowed distance (multiplied by actual stepwidth) to detect loop 75 * loop_dir: 0.99, // Should be > 0.95 76 * loop_detection: true, // Use Gosper's loop detector 77 * unitX: 10, // unitX of board 78 * unitY: 10 // unitX of board 79 * }; 80 * </pre> 81 * @param {function} f function from <b>R</b><sup>2</sup> to <b>R</b> 82 * @param {function} [dfx] Optional partial derivative of <i>f</i> with regard to <i>x</i> 83 * @param {function} [dfy] Optional partial derivative of <i>f</i> with regard to <i>y</i> 84 * 85 * @constructor 86 * @example 87 * var f = (x, y) => x**3 - 2 * x * y + y**3; 88 * var c = board.create('curve', [[], []], { 89 * strokeWidth: 3, 90 * strokeColor: JXG.palette.red 91 * }); 92 * 93 * c.updateDataArray = function () { 94 * var bbox = this.board.getBoundingBox(), 95 * ip, cfg, 96 * ret = [], 97 * mgn = 1; 98 * 99 * bbox[0] -= mgn; 100 * bbox[1] += mgn; 101 * bbox[2] += mgn; 102 * bbox[3] -= mgn; 103 * 104 * cfg = { 105 * resolution_out: 5, 106 * resolution_in: 5, 107 * unitX: this.board.unitX, 108 * unitY: this.board.unitX 109 * }; 110 * 111 * this.dataX = []; 112 * this.dataY = []; 113 * ip = new JXG.Math.ImplicitPlot(bbox, cfg, f, null, null); 114 * ret = ip.plot(); 115 * this.dataX = ret[0]; 116 * this.dataY = ret[1]; 117 * }; 118 * board.update(); 119 * </pre><div id="JXGf3e8cd82-2b67-4efb-900a-471eb92b3b96" class="jxgbox" style="width: 300px; height: 300px;"></div> 120 * <script type="text/javascript"> 121 * (function() { 122 * var board = JXG.JSXGraph.initBoard('JXGf3e8cd82-2b67-4efb-900a-471eb92b3b96', 123 * {boundingbox: [-8, 8, 8,-8], axis: true, showcopyright: false, shownavigation: false}); 124 * var f = (x, y) => x**3 - 2 * x * y + y**3; 125 * var c = board.create('curve', [[], []], { 126 * strokeWidth: 3, 127 * strokeColor: JXG.palette.red 128 * }); 129 * 130 * c.updateDataArray = function () { 131 * var bbox = this.board.getBoundingBox(), 132 * ip, cfg, 133 * ret = [], 134 * mgn = 1; 135 * 136 * bbox[0] -= mgn; 137 * bbox[1] += mgn; 138 * bbox[2] += mgn; 139 * bbox[3] -= mgn; 140 * 141 * cfg = { 142 * resolution_out: 5, 143 * resolution_in: 5, 144 * unitX: this.board.unitX, 145 * unitY: this.board.unitX 146 * }; 147 * 148 * this.dataX = []; 149 * this.dataY = []; 150 * 151 * ip = new JXG.Math.ImplicitPlot(bbox, cfg, f, null, null); 152 * ret = ip.plot(); 153 * 154 * this.dataX = ret[0]; 155 * this.dataY = ret[1]; 156 * }; 157 * board.update(); 158 * 159 * })(); 160 * 161 * </script><pre> 162 * 163 */ 164 Mat.ImplicitPlot = function (bbox, config, f, dfx, dfy) { 165 166 // Default values 167 var cfg_default = { 168 resolution_out: 5, // Distance between vertical lines to search for components 169 resolution_in: 5, // Distance between vertical lines to search for components 170 max_steps: 1024, // Max number of points in one call of tracing 171 alpha_0: 0.05, // Angle between two successive tangents: smoothness of curve 172 173 tol_u0: Mat.eps, // Tolerance to find starting points for tracing. 174 tol_newton: 1.0e-7, // Tolerance for Newton steps. 175 tol_cusp: 0.05, // Tolerance for cusp / bifurcation detection 176 tol_progress: 0.0001, // If two points are closer than this value, we bail out 177 qdt_box: 0.2, // half of box size to search in qdt 178 kappa_0: 0.2, // Inverse of planned number of Newton steps 179 delta_0: 0.05, // Distance of predictor point to curve 180 181 h_initial: 0.1, // Initial step width 182 h_critical: 0.001, // If h is below this threshold we bail out 183 h_max: 1, // Maximum value of h (user units) 184 loop_dist: 0.09, // Allowed distance (multiplied by actual step width) to detect loop 185 loop_dir: 0.99, // Should be > 0.95 186 loop_detection: true, // Use Gosper's loop detector 187 unitX: 10, // unitX of board 188 unitY: 10 // unitX of board 189 }; 190 191 this.config = Type.merge(cfg_default, config); 192 193 this.f = f; 194 195 this.dfx = null; 196 this.dfy = null; 197 198 if (Type.isFunction(dfx)) { 199 this.dfx = dfx; 200 } else { 201 this.dfx = function (x, y) { 202 var h = Mat.eps * Mat.eps; 203 return (this.f(x + h, y) - this.f(x - h, y)) * 0.5 / h; 204 }; 205 } 206 207 if (Type.isFunction(dfy)) { 208 this.dfy = dfy; 209 } else { 210 this.dfy = function (x, y) { 211 var h = Mat.eps * Mat.eps; 212 return (this.f(x, y + h) - this.f(x, y - h)) * 0.5 / h; 213 }; 214 } 215 216 this.bbox = bbox; 217 this.qdt = new Quadtree(20, 5, bbox); 218 219 this.components = []; 220 }; 221 222 Type.extend( 223 Mat.ImplicitPlot.prototype, 224 /** @lends JXG.Math.ImplicitPlot.prototype */ { 225 226 /** 227 * Implicit plotting method. 228 * 229 * @returns {Array} consisting of [dataX, dataY, number_of_components] 230 */ 231 plot: function () { 232 var // components = [], 233 doVerticalSearch = true, 234 doHorizontalSearch = true, 235 x, y, 236 mi_x, ma_x, mi_y, ma_y, 237 dataX = [], 238 dataY = [], 239 ret = [], 240 num_components = 0, 241 max_level = 8, 242 243 delta, 244 that = this, 245 246 fmi_x = function (t) { 247 return that.f(x, t); 248 }, 249 fma_x = function (t) { 250 return -that.f(x, t); 251 }, 252 fmi_y = function (t) { 253 return that.f(t, y); 254 }, 255 fma_y = function (t) { 256 return -that.f(t, y); 257 }; 258 259 // Vertical lines or circular search: 260 mi_x = Math.min(this.bbox[0], this.bbox[2]) - Mat.eps; 261 ma_x = Math.max(this.bbox[0], this.bbox[2]); 262 mi_y = Math.min(this.bbox[1], this.bbox[3]) + Mat.eps; 263 ma_y = Math.max(this.bbox[1], this.bbox[3]); 264 265 if (doVerticalSearch) { 266 delta = this.config.resolution_out / this.config.unitX; 267 delta *= (1 + Mat.eps); 268 // console.log("Outer delta x", delta) 269 270 for (x = mi_x; x < ma_x; x += delta) { 271 ret = this.searchLine( 272 fmi_x, fma_x, x, 273 [mi_y, ma_y], 'vertical', 274 num_components, dataX, dataY, max_level); 275 276 if (ret !== false) { 277 dataX = ret[0]; 278 dataY = ret[1]; 279 num_components = ret[2]; 280 } 281 282 } 283 } 284 if (doHorizontalSearch) { 285 delta = this.config.resolution_out / this.config.unitY; 286 delta *= (1 + Mat.eps); 287 // console.log("Outer delta y", delta) 288 289 for (y = mi_y; y < ma_y; y += delta) { 290 ret = this.searchLine( 291 fmi_y, fma_y, y, 292 [mi_x, ma_x], 'horizontal', 293 num_components, dataX, dataY, max_level); 294 295 if (ret !== false) { 296 dataX = ret[0]; 297 dataY = ret[1]; 298 num_components = ret[2]; 299 } 300 } 301 } 302 303 return [dataX, dataY, num_components]; 304 }, 305 306 /** 307 * Recursively search a horizontal or vertical line for points on the 308 * fulfilling the given equation. 309 * 310 * @param {Function} fmi Minimization function 311 * @param {Function} fma Maximization function 312 * @param {Number} fix Value of the fixed variable 313 * @param {Array} interval Search interval of the free variable 314 * @param {String} dir 'vertical' or 'horizontal' 315 * @param {Number} num_components Number of components before search 316 * @param {Array} dataX x-coordinates of points so far 317 * @param {Array} dataY y-coordinates of points so far 318 * @param {Number} level Recursion level 319 * @returns {Array} consisting of [dataX, dataY, number_of_components]- 320 * @private 321 */ 322 searchLine: function (fmi, fma, fix, interval, dir, 323 num_components, dataX, dataY, level) { 324 var t_mi, t_ma, t, 325 ft, 326 mi, ma, tmp, m, 327 is_in, 328 u0, i, le, 329 ret, 330 offset, 331 delta, 332 eps = this.config.tol_u0, 333 DEBUG = false, 334 b = interval[0], 335 e = interval[1]; 336 337 t_mi = Numerics.fminbr(fmi, [b, e]); 338 mi = fmi(t_mi); 339 t_ma = Numerics.fminbr(fma, [b, e]); 340 ma = fmi(t_ma); 341 342 if (mi < eps && ma > -eps) { 343 tmp = t_mi; 344 t_mi = Math.min(tmp, t_ma); 345 t_ma = Math.max(tmp, t_ma); 346 347 t = Numerics.fzero(fmi, [t_mi, t_ma]); 348 // t = Numerics.chandrupatla(fmi, [t_mi, t_ma]); 349 350 ft = fmi(t); 351 if (Math.abs(ft) > Math.max((ma - mi) * Mat.eps, 0.001)) { 352 //console.log("searchLine:", dir, fix, t, "no root " + ft); 353 return false; 354 // throw new Error("searchLine: no root " + ft); 355 } 356 if (dir === 'vertical') { 357 u0 = [1, fix, t]; 358 delta = this.config.resolution_in / this.config.unitY; 359 // console.log("Inner delta x", delta) 360 } else { 361 u0 = [1, t, fix]; 362 delta = this.config.resolution_in / this.config.unitX; 363 // console.log("Inner delta y", delta) 364 } 365 delta *= (1 + Mat.eps); 366 367 is_in = this.curveContainsPoint(u0, dataX, dataY, 368 delta * 2, // Allowed dist from segment 369 this.config.qdt_box // 0.5 of box size to search in qdt 370 ); 371 372 if (is_in) { 373 if (DEBUG) { 374 console.log("Found in quadtree", u0, "level:", level); 375 } 376 } else { 377 if (DEBUG) { 378 console.log("Not in quadtree", u0, dataX.length); 379 } 380 ret = this.traceComponent(u0, 1); 381 if (ret[0].length > 0) { 382 // Add jump in curve 383 if (num_components > 0) { 384 dataX.push(NaN); 385 dataY.push(NaN); 386 } 387 388 offset = dataX.length; 389 le = ret[0].length; 390 for (i = 1; i < le; i++) { 391 this.qdt.insertItem({ 392 xlb: Math.min(ret[0][i - 1], ret[0][i]), 393 xub: Math.max(ret[0][i - 1], ret[0][i]), 394 ylb: Math.min(ret[1][i - 1], ret[1][i]), 395 yub: Math.max(ret[1][i - 1], ret[1][i]), 396 idx1: offset + i - 1, 397 idx2: offset + i, 398 comp: num_components 399 }); 400 } 401 402 num_components++; 403 Type.concat(dataX, ret[0]); 404 Type.concat(dataY, ret[1]); 405 } 406 } 407 408 m = t - delta * 0.01; 409 if (m - b > delta && level > 0) { 410 ret = this.searchLine( 411 fmi, fma, fix, [b, m], dir, 412 num_components, dataX, dataY, level - 1); 413 if (ret !== false) { 414 dataX = ret[0]; 415 dataY = ret[1]; 416 num_components = ret[2]; 417 } 418 } 419 m = t + delta * 0.01; 420 if (e - m > delta && level > 0) { 421 ret = this.searchLine( 422 fmi, fma, fix, [m, e], dir, 423 num_components, dataX, dataY, level - 1); 424 if (ret !== false) { 425 dataX = ret[0]; 426 dataY = ret[1]; 427 num_components = ret[2]; 428 } 429 } 430 431 return [dataX, dataY, num_components]; 432 } 433 434 return false; 435 }, 436 437 /** 438 * Test if the data points contain a given coordinate, i.e. if the 439 * given coordinate is close enough to the polygonal chain 440 * through the data points. 441 * 442 * @param {Array} p Homogenous coordinates [1, x, y] of the coordinate point 443 * @param {Array} dataX x-coordinates of points so far 444 * @param {Array} dataY y-coordinates of points so far 445 * @param {Number} tol Maximal distance of p from the polygonal chain through the data points 446 * @param {Number} eps Helper tolerance used for the quadtree 447 * @returns Boolean 448 */ 449 curveContainsPoint: function (p, dataX, dataY, tol, eps) { 450 var i, le, hits, d, 451 x = p[1], 452 y = p[2]; 453 454 hits = this.qdt.find([x - eps, y + eps, x + eps, y - eps]); 455 456 le = hits.length; 457 for (i = 0; i < le; i++) { 458 d = Geometry.distPointSegment( 459 p, 460 [1, dataX[hits[i].idx1], dataY[hits[i].idx1]], 461 [1, dataX[hits[i].idx2], dataY[hits[i].idx2]] 462 ); 463 if (d < tol) { 464 return true; 465 } 466 } 467 return false; 468 }, 469 470 /** 471 * Starting at an initial point the curve is traced with a Euler-Newton method. 472 * After tracing in one direction the algorithm stops if the component is a closed loop. 473 * Otherwise, the curved is traced in the opposite direction, starting from 474 * the same initial point. Finally, the two components are glued together. 475 * 476 * @param {Array} u0 Initial point in homogenous coordinates [1, x, y]. 477 * @returns Array [dataX, dataY] containing a new component. 478 * @private 479 */ 480 traceComponent: function (u0) { 481 var dataX = [], 482 dataY = [], 483 arr = []; 484 485 // Trace in first direction 486 // console.log("---- Start tracing forward ---------") 487 arr = this.tracing(u0, 1); 488 489 if (arr.length === 0) { 490 // console.log("Could not start tracing due to singularity") 491 } else { 492 // console.log("Trace from", [arr[0][0], arr[1][0]], "to", [arr[0][arr[0].length - 1], arr[1][arr[1].length - 1]], 493 // "num points:", arr[0].length); 494 dataX = arr[0]; 495 dataY = arr[1]; 496 } 497 498 // Trace in the other direction 499 if (!arr[2]) { 500 // No loop in the first tracing step, 501 // now explore the other direction. 502 503 // console.log("---- Start tracing backward ---------") 504 arr = this.tracing(u0, -1); 505 506 if (arr.length === 0) { 507 // console.log("Could not start backward tracing due to singularity") 508 } else { 509 // console.log("Trace backwards from", [arr[0][0], arr[1][0]], "to", 510 // [arr[0][arr[0].length - 1], arr[1][arr[1].length - 1]], "num points:", arr[0].length); 511 dataX = arr[0].reverse().concat(dataX.slice(1)); 512 dataY = arr[1].reverse().concat(dataY.slice(1)); 513 } 514 } 515 516 if (dataX.length > 0 && dataX.length < 6) { 517 // Solitary point 518 dataX.push(dataX[dataX.length - 1]); 519 dataY.push(dataY[dataY.length - 1]); 520 } 521 522 return [dataX, dataY]; 523 }, 524 525 /** 526 * Starting at a point <i>u0</i>, this routine traces the curve <i>f(u)=0</i> until 527 * a loop is detected, a critical point is reached, the curve leaves the bounding box, 528 * or the maximum number of points is reached. 529 * <p> 530 * The method is a predictor / corrector method consisting of Euler and Newton steps 531 * together with step width adaption. 532 * <p> 533 * The algorithm is an adaption of the algorithm in 534 * Eugene L. Allgower, Kurt Georg: <i>Introduction to Numerical Continuation methods.</i> 535 * 536 * @param {Array} u0 Starting point in homogenous coordinates [1, x, y]. 537 * @param {Number} direction 1 or -1 538 * @returns Array [pathX, pathY, loop_closed] or [] 539 * @private 540 */ 541 tracing: function (u0, direction) { 542 var u = [], 543 ulast = [], 544 len, 545 v = [], 546 v_start = [], 547 w = [], 548 t_u, t_v, t_u_0, tloc, 549 A, 550 grad, 551 nrm, 552 dir, 553 steps = 0, 554 k = 0, 555 loop_closed = false, 556 k0, k1, denom, dist, progress, 557 kappa, delta, alpha, 558 factor, 559 point_added = false, 560 561 quasi = false, 562 cusp_or_bifurc = false, 563 kappa_0 = this.config.kappa_0, // Inverse of planned number of Newton steps 564 delta_0 = this.config.delta_0, // Distance of predictor point to curve 565 alpha_0 = this.config.alpha_0, // Angle between two successive tangents 566 h = this.config.h_initial, 567 max_steps = this.config.max_steps, 568 569 omega = direction, 570 pathX = [], 571 pathY = [], 572 573 T = [], // Gosper's loop detector table 574 n, m, i, e; 575 576 u = u0.slice(1); 577 pathX.push(u[0]); 578 pathY.push(u[1]); 579 580 t_u = this.tangent(u); 581 if (t_u === false) { 582 // We don't want to start at a singularity. 583 // Get out of here and search for another starting point. 584 return []; 585 } 586 A = [this.dfx(u[0], u[1]), this.dfy(u[0], u[1])]; 587 588 do { 589 590 if (quasi) { 591 t_u = this.tangent_A(A); 592 } else { 593 t_u = this.tangent(u); 594 } 595 if (t_u === false) { 596 u = v.slice(); 597 pathX.push(u[0]); 598 pathY.push(u[1]); 599 // console.log("-> Bail out: t_u undefined."); 600 break; 601 } 602 603 if (pathX.length === 1) { 604 // Store first point 605 t_u_0 = t_u.slice(); 606 } else if (pathX.length === 2) { 607 T.push(pathX.length - 1); // Put first point into Gosper table T 608 609 } else if (point_added && pathX.length > 2 && !cusp_or_bifurc) { 610 611 // Detect if loop has been closed 612 dist = Geometry.distPointSegment( 613 [1, u[0], u[1]], 614 [1, pathX[0], pathY[0]], 615 [1, pathX[1], pathY[1]] 616 ); 617 618 if (dist < this.config.loop_dist * h && 619 Mat.innerProduct(t_u, t_u_0, 2) > this.config.loop_dir 620 ) { 621 622 // console.log("Loop detected after", steps, 'steps'); 623 // console.log("\t", "v", v, "u0:", u0) 624 // console.log("\t", "Dist(v, path0)", dist, config.loop_dist * h) 625 // console.log("\t", "t_u", t_u); 626 // console.log("\t", "inner:", Mat.innerProduct(t_u, t_u_0, 2)); 627 // console.log("\t", "h", h); 628 629 u = u0.slice(1); 630 pathX.push(u[0]); 631 pathY.push(u[1]); 632 633 loop_closed = true; 634 break; 635 } 636 637 // Gosper's loop detector 638 if (this.config.loop_detection) { 639 n = pathX.length - 1; 640 // console.log("Check Gosper", n); 641 m = Math.floor(Mat.log2(n)); 642 643 for (i = 0; i <= m; i++) { 644 dist = Geometry.distPointSegment( 645 [1, u[0], u[1]], 646 [1, pathX[T[i] - 1], pathY[T[i] - 1]], 647 [1, pathX[T[i]], pathY[T[i]]] 648 ); 649 650 if (dist < this.config.loop_dist * h) { 651 // console.log("!!!!!!!!!!!!!!! GOSPER LOOP CLOSED !!!!", i, n + 1, 652 // this.config.loop_dist * h 653 // ); 654 655 t_v = this.tangent([pathX[T[i]], pathY[T[i]]]); 656 if (Mat.innerProduct(t_u, t_v, 2) > this.config.loop_dir) { 657 // console.log("!!!!!!!!!!!!!!! angle is good enough"); 658 break; 659 } 660 } 661 } 662 if (i <= m) { 663 loop_closed = true; 664 break; 665 } 666 667 m = 1; 668 e = 0; 669 for (i = 0; i < 100; i++) { 670 if ((n + 1) % m !== 0) { 671 break; 672 } 673 m *= 2; 674 e++; 675 } 676 // console.log("Add at e", e); 677 T[e] = n; 678 } 679 680 } 681 682 // Predictor step 683 // if (true /*h < 2 * this.config.h_initial*/) { 684 // Euler 685 // console.log('euler') 686 v[0] = u[0] + h * omega * t_u[0]; 687 v[1] = u[1] + h * omega * t_u[1]; 688 // } else { 689 // // Heun 690 // // console.log('heun') 691 // v[0] = u[0] + h * omega * t_u[0]; 692 // v[1] = u[1] + h * omega * t_u[1]; 693 694 // t_v = this.tangent(v); 695 // v[0] = 0.5 * u[0] + 0.5 * (v[0] + h * omega * t_v[0]); 696 // v[1] = 0.5 * u[1] + 0.5 * (v[1] + h * omega * t_v[1]); 697 // } 698 if (quasi) { 699 A = this.updateA(A, u, v); 700 v_start = v.slice(); 701 } 702 703 // Corrector step: Newton 704 k = 0; 705 do { 706 if (quasi) { 707 grad = A; 708 } else { 709 grad = [this.dfx(v[0], v[1]), this.dfy(v[0], v[1])]; 710 } 711 712 // Compute w = v - A(v) * f(v), 713 // grad: row vector and A(v) is the Moore-Penrose inverse: 714 // grad^T * (grad * grad^T)^(-1) 715 denom = grad[0] * grad[0] + grad[1] * grad[1]; 716 nrm = this.f(v[0], v[1]) / denom; 717 718 w[0] = v[0] - grad[0] * nrm; 719 w[1] = v[1] - grad[1] * nrm; 720 if (k === 0) { 721 k0 = Math.abs(nrm) * Math.sqrt(denom); 722 } else if (k === 1) { 723 k1 = Math.abs(nrm) * Math.sqrt(denom); 724 } 725 726 v[0] = w[0]; 727 v[1] = w[1]; 728 k++; 729 } while (k < 20 && 730 Math.abs(this.f(v[0], v[1])) > this.config.tol_newton 731 ); 732 733 delta = k0; 734 if (k > 1) { 735 kappa = k1 / k0; 736 } else { 737 kappa = 0.0; 738 } 739 740 if (quasi) { 741 A = this.updateA(A, v_start, v); 742 t_v = this.tangent_A(A); 743 } else { 744 t_v = this.tangent(v); 745 } 746 747 dir = Mat.innerProduct(t_u, t_v, 2); 748 dir = Math.max(-1, Math.min(1, dir)); 749 alpha = Math.acos(dir); 750 751 // Look for simple bifurcation points and cusps 752 cusp_or_bifurc = false; 753 progress = Geometry.distance(u, v, 2); 754 if (progress < this.config.tol_progress) { 755 u = v.slice(); 756 pathX.push(u[0]); 757 pathY.push(u[1]); 758 // console.log("-> Bail out, no progress", progress, steps); 759 break; 760 761 } else if (dir < 0.0) { 762 if (h > this.config.h_critical) { 763 // console.log("Critical point at [", u[0].toFixed(4), u[1].toFixed(4), "], v: [", v[0].toFixed(4), v[1].toFixed(4), "], but large h:", h); 764 765 } else { 766 767 cusp_or_bifurc = true; 768 if (this.isBifurcation(u, this.config.tol_cusp)) { 769 // console.log(steps, "bifurcation point between", u, "and", v, ":", dir, "h", h, "alpha", alpha); 770 // A = [dfx(v[0], v[1]), dfy(v[0], v[1])]; 771 omega *= (-1); 772 // If there is a bifurcation point, we 773 // ignore the angle alpha for subsequent step length 774 // adaption. Because then we might be able to 775 // "jump over the critical point" 776 alpha = 0; 777 } else { 778 // Cusp or something more weird 779 u = v.slice(); 780 pathX.push(u[0]); 781 pathY.push(u[1]); 782 // console.log("-> Bail out, cusp") 783 break; 784 } 785 } 786 } 787 788 // Adapt stepwidth 789 if (!cusp_or_bifurc) { 790 factor = Math.max( 791 Math.sqrt(kappa / kappa_0), 792 Math.sqrt(delta / delta_0), 793 alpha / alpha_0 794 ); 795 if (isNaN(factor)) { 796 factor = 1; 797 } 798 factor = Math.max(Math.min(factor, 2), 0.5); 799 h /= factor; 800 h = Math.min(this.config.h_max, h); 801 802 if (factor >= 2) { 803 steps++; 804 if (steps >= 3 * max_steps) { 805 break; 806 } 807 808 point_added = false; 809 continue; 810 } 811 } 812 813 u = v.slice(); 814 pathX.push(u[0]); 815 pathY.push(u[1]); 816 point_added = true; 817 818 steps++; 819 } while ( 820 steps < max_steps && 821 u[0] >= this.bbox[0] && 822 u[1] <= this.bbox[1] && 823 u[0] <= this.bbox[2] && 824 u[1] >= this.bbox[3] 825 ); 826 827 // Clipping to bounding box, last may be outside, interpolate between second last und last point 828 len = pathX.length; 829 ulast = [pathX[len - 2], pathY[len - 2]]; 830 831 // If u[0] is outside x-interval in bounding box, interpolate to the box. 832 if (u[0] < this.bbox[0]) { 833 if (u[0] !== ulast[0]) { 834 tloc = (this.bbox[0] - ulast[0]) / (u[0] - ulast[0]); 835 if (u[1] !== ulast[1]) { 836 u[1] = ulast[1] + tloc * (u[1] - ulast[1]); 837 } 838 } 839 u[0] = this.bbox[0]; 840 } 841 if (u[0] > this.bbox[2]) { 842 if (u[0] !== ulast[0]) { 843 tloc = (this.bbox[2] - ulast[0]) / (u[0] - ulast[0]); 844 if (u[1] !== ulast[1]) { 845 u[1] = ulast[1] + tloc * (u[1] - ulast[1]); 846 } 847 } 848 u[0] = this.bbox[2]; 849 } 850 851 // If u[1] is outside y-interval in bounding box, interpolate to the box. 852 if (u[1] < this.bbox[3]) { 853 if (u[1] !== ulast[1]) { 854 tloc = (this.bbox[3] - ulast[1]) / (u[1] - ulast[1]); 855 if (u[0] !== ulast[0]) { 856 u[0] = ulast[0] + tloc * (u[0] - ulast[0]); 857 } 858 } 859 u[1] = this.bbox[3]; 860 } 861 if (u[1] > this.bbox[1]) { 862 if (u[1] !== ulast[1]) { 863 tloc = (this.bbox[1] - ulast[1]) / (u[1] - ulast[1]); 864 if (u[0] !== ulast[0]) { 865 u[0] = ulast[0] + tloc * (u[0] - ulast[0]); 866 } 867 } 868 u[1] = this.bbox[1]; 869 } 870 871 // Update last point 872 pathX[len - 1] = u[0]; 873 pathY[len - 1] = u[1]; 874 875 // if (!loop_closed) { 876 // console.log("No loop", steps); 877 // } else { 878 // console.log("Loop", steps); 879 // } 880 881 return [pathX, pathY, loop_closed]; 882 }, 883 884 /** 885 * If both eigenvalues of the Hessian are different from zero, the critical point at u 886 * is a simple bifurcation point. 887 * 888 * @param {Array} u Critical point [x, y] 889 * @param {Number} tol Tolerance of the eigenvalues to be zero. 890 * @returns Boolean True if the point is a simple bifurcation point. 891 * @private 892 */ 893 isBifurcation: function (u, tol) { 894 // Former experiments: 895 // If the Hessian has exactly one zero eigenvalue, 896 // we claim that there is a cusp. 897 // Otherwise, we decide that there is a bifurcation point. 898 // In the latter case, if both eigenvalues are zero 899 // this is a somewhat crude decision. 900 // 901 var h = Mat.eps * Mat.eps * 100, 902 x, y, a, b, c, d, ad, 903 lbda1, lbda2, 904 dis; 905 906 x = u[0]; 907 y = u[1]; 908 a = 0.5 * (this.dfx(x + h, y) - this.dfx(x - h, y)) / h; 909 b = 0.5 * (this.dfx(x, y + h) - this.dfx(x, y - h)) / h; 910 c = 0.5 * (this.dfy(x + h, y) - this.dfy(x - h, y)) / h; 911 d = 0.5 * (this.dfy(x, y + h) - this.dfy(x, y - h)) / h; 912 913 // c = b 914 ad = a + d; 915 dis = ad * ad - 4 * (a * d - b * c); 916 lbda1 = 0.5 * (ad + Math.sqrt(dis)); 917 lbda2 = 0.5 * (ad - Math.sqrt(dis)); 918 919 // console.log(a, b, c, d) 920 // console.log("Eigenvals u:", lbda1, lbda2, tol); 921 922 if (Math.abs(lbda1) > tol && Math.abs(lbda2) > tol) { 923 // if (lbda1 * lbda2 > 0) { 924 // console.log("Seems to be isolated singularity at", u); 925 // } 926 return true; 927 } 928 929 return false; 930 }, 931 932 /** 933 * Search in an arc around a critical point for a further point on the curve. 934 * Unused for the moment. 935 * 936 * @param {Array} u Critical point [x, y] 937 * @param {Array} t_u Tangent at u 938 * @param {Number} r Radius 939 * @param {Number} omega angle 940 * @returns {Array} Coordinates [x, y] of a new point. 941 * @private 942 */ 943 handleCriticalPoint: function (u, t_u, r, omega) { 944 var a = Math.atan2(omega * t_u[1], omega * t_u[0]), 945 // s = a - 0.75 * Math.PI, 946 // e = a + 0.75 * Math.PI, 947 f_circ = function (t) { 948 var x = u[0] + r * Math.cos(t), 949 y = u[1] + r * Math.sin(t); 950 return this.f(x, y); 951 }, 952 x, y, t0; 953 954 // t0 = Numerics.fzero(f_circ, [s, e]); 955 t0 = Numerics.root(f_circ, a); 956 957 x = u[0] + r * Math.cos(t0); 958 y = u[1] + r * Math.sin(t0); 959 // console.log("\t", "result", x, y, "f", f(x, y)); 960 961 return [x, y]; 962 }, 963 964 /** 965 * Quasi-Newton update of the Moore-Penrose inverse. 966 * See (7.2.3) in Allgower, Georg. 967 * 968 * @param {Array} A 969 * @param {Array} u0 970 * @param {Array} u1 971 * @returns Array 972 * @private 973 */ 974 updateA: function (A, u0, u1) { 975 var s = [u1[0] - u0[0], u1[1] - u0[1]], 976 y = this.f(u1[0], u1[1]) - this.f(u0[0], u0[1]), 977 nom, denom; 978 979 denom = s[0] * s[0] + s[1] * s[1]; 980 nom = y - (A[0] * s[0] + A[1] * s[1]); 981 nom /= denom; 982 A[0] += nom * s[0]; 983 A[1] += nom * s[1]; 984 985 return A; 986 }, 987 988 /** 989 * Approximate tangent (of norm 1) with Quasi-Newton method 990 * @param {Array} A 991 * @returns Array 992 * @private 993 */ 994 tangent_A: function (A) { 995 var t = [-A[1], A[0]], 996 nrm = Mat.norm(t, 2); 997 998 if (nrm < Mat.eps) { 999 // console.log("Approx. Singularity", t, "is zero", nrm); 1000 } 1001 return [t[0] / nrm, t[1] / nrm]; 1002 }, 1003 1004 /** 1005 * Tangent of norm 1 at point u. 1006 * @param {Array} u Point [x, y] 1007 * @returns Array 1008 * @private 1009 */ 1010 tangent: function (u) { 1011 var t = [-this.dfy(u[0], u[1]), this.dfx(u[0], u[1])], 1012 nrm = Mat.norm(t, 2); 1013 1014 if (nrm < Mat.eps * Mat.eps) { 1015 // console.log("Singularity", t, "is zero", "at", u, ":", nrm); 1016 return false; 1017 } 1018 return [t[0] / nrm, t[1] / nrm]; 1019 } 1020 } 1021 1022 ); 1023 1024 export default Mat.ImplicitPlot; 1025 1026