1 /* 2 Copyright 2008-2025 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 242 delta, 243 that = this, 244 245 fmi_x = function (t) { 246 return that.f(x, t); 247 }, 248 fma_x = function (t) { 249 return -that.f(x, t); 250 }, 251 fmi_y = function (t) { 252 return that.f(t, y); 253 }, 254 fma_y = function (t) { 255 return -that.f(t, y); 256 }; 257 258 // Vertical lines or circular search: 259 mi_x = Math.min(this.bbox[0], this.bbox[2]) - Mat.eps; 260 ma_x = Math.max(this.bbox[0], this.bbox[2]); 261 mi_y = Math.min(this.bbox[1], this.bbox[3]) + Mat.eps; 262 ma_y = Math.max(this.bbox[1], this.bbox[3]); 263 264 if (doVerticalSearch) { 265 delta = this.config.resolution_out / this.config.unitX; 266 delta *= (1 + Mat.eps); 267 // console.log("Outer delta x", delta) 268 269 for (x = mi_x; x < ma_x; x += delta) { 270 ret = this.searchLine( 271 fmi_x, fma_x, x, 272 [mi_y, ma_y], 'vertical', 273 num_components, dataX, dataY, 20); 274 275 if (ret !== false) { 276 dataX = ret[0]; 277 dataY = ret[1]; 278 num_components = ret[2]; 279 } 280 281 } 282 } 283 if (doHorizontalSearch) { 284 delta = this.config.resolution_out / this.config.unitY; 285 delta *= (1 + Mat.eps); 286 // console.log("Outer delta y", delta) 287 288 for (y = mi_y; y < ma_y; y += delta) { 289 ret = this.searchLine( 290 fmi_y, fma_y, y, 291 [mi_x, ma_x], 'horizontal', 292 num_components, dataX, dataY, 20); 293 294 if (ret !== false) { 295 dataX = ret[0]; 296 dataY = ret[1]; 297 num_components = ret[2]; 298 } 299 } 300 } 301 302 return [dataX, dataY, num_components]; 303 }, 304 305 /** 306 * Recursively search a horizontal or vertical line for points on the 307 * fulfilling the given equation. 308 * 309 * @param {Function} fmi Minimization function 310 * @param {Function} fma Maximization function 311 * @param {Number} fix Value of the fixed variable 312 * @param {Array} interval Search interval of the free variable 313 * @param {String} dir 'vertical' or 'horizontal' 314 * @param {Number} num_components Number of components before search 315 * @param {Array} dataX x-coordinates of points so far 316 * @param {Array} dataY y-coordinates of points so far 317 * @param {Number} level Recursion level 318 * @returns {Array} consisting of [dataX, dataY, number_of_components]- 319 * @private 320 */ 321 searchLine: function (fmi, fma, fix, interval, dir, 322 num_components, dataX, dataY, level) { 323 var t_mi, t_ma, t, 324 ft, 325 mi, ma, tmp, m, 326 is_in, 327 u0, i, le, 328 ret, 329 offset, 330 delta, 331 eps = this.config.tol_u0, 332 DEBUG = false, 333 b = interval[0], 334 e = interval[1]; 335 336 t_mi = Numerics.fminbr(fmi, [b, e]); 337 mi = fmi(t_mi); 338 t_ma = Numerics.fminbr(fma, [b, e]); 339 ma = fmi(t_ma); 340 341 if (mi < eps && ma > -eps) { 342 tmp = t_mi; 343 t_mi = Math.min(tmp, t_ma); 344 t_ma = Math.max(tmp, t_ma); 345 346 t = Numerics.fzero(fmi, [t_mi, t_ma]); 347 // t = Numerics.chandrupatla(fmi, [t_mi, t_ma]); 348 349 ft = fmi(t); 350 if (Math.abs(ft) > Math.max((ma - mi) * Mat.eps, 0.001)) { 351 //console.log("searchLine:", dir, fix, t, "no root " + ft); 352 return false; 353 // throw new Error("searchLine: no root " + ft); 354 } 355 if (dir === 'vertical') { 356 u0 = [1, fix, t]; 357 delta = this.config.resolution_in / this.config.unitY; 358 // console.log("Inner delta x", delta) 359 } else { 360 u0 = [1, t, fix]; 361 delta = this.config.resolution_in / this.config.unitX; 362 // console.log("Inner delta y", delta) 363 } 364 delta *= (1 + Mat.eps); 365 366 is_in = this.curveContainsPoint(u0, dataX, dataY, 367 delta * 2, // Allowed dist from segment 368 this.config.qdt_box // 0.5 of box size to search in qdt 369 ); 370 371 if (is_in) { 372 if (DEBUG) { 373 console.log("Found in quadtree", u0); 374 } 375 } else { 376 if (DEBUG) { 377 console.log("Not in quadtree", u0, dataX.length); 378 } 379 ret = this.traceComponent(u0, 1); 380 if (ret[0].length > 0) { 381 // Add jump in curve 382 if (num_components > 0) { 383 dataX.push(NaN); 384 dataY.push(NaN); 385 } 386 387 offset = dataX.length; 388 le = ret[0].length; 389 for (i = 1; i < le; i++) { 390 this.qdt.insertItem({ 391 xlb: Math.min(ret[0][i - 1], ret[0][i]), 392 xub: Math.max(ret[0][i - 1], ret[0][i]), 393 ylb: Math.min(ret[1][i - 1], ret[1][i]), 394 yub: Math.max(ret[1][i - 1], ret[1][i]), 395 idx1: offset + i - 1, 396 idx2: offset + i, 397 comp: num_components 398 }); 399 } 400 401 num_components++; 402 Type.concat(dataX, ret[0]); 403 Type.concat(dataY, ret[1]); 404 } 405 } 406 407 m = t - delta * 0.01; 408 if (m - b > delta && level > 0) { 409 ret = this.searchLine( 410 fmi, fma, fix, [b, m], dir, 411 num_components, dataX, dataY, level - 1); 412 if (ret !== false) { 413 dataX = ret[0]; 414 dataY = ret[1]; 415 num_components = ret[2]; 416 } 417 } 418 m = t + delta * 0.01; 419 if (e - m > delta && level > 0) { 420 ret = this.searchLine( 421 fmi, fma, fix, [m, e], dir, 422 num_components, dataX, dataY, level - 1); 423 if (ret !== false) { 424 dataX = ret[0]; 425 dataY = ret[1]; 426 num_components = ret[2]; 427 } 428 } 429 430 return [dataX, dataY, num_components]; 431 } 432 433 return false; 434 }, 435 436 /** 437 * Test if the data points contain a given coordinate, i.e. if the 438 * given coordinate is close enough to the polygonal chain 439 * through the data points. 440 * 441 * @param {Array} p Homogenous coordinates [1, x, y] of the coordinate point 442 * @param {Array} dataX x-coordinates of points so far 443 * @param {Array} dataY y-coordinates of points so far 444 * @param {Number} tol Maximal distance of p from the polygonal chain through the data points 445 * @param {Number} eps Helper tolerance used for the quadtree 446 * @returns Boolean 447 */ 448 curveContainsPoint: function (p, dataX, dataY, tol, eps) { 449 var i, le, hits, d, 450 x = p[1], 451 y = p[2]; 452 453 hits = this.qdt.find([x - eps, y + eps, x + eps, y - eps]); 454 455 le = hits.length; 456 for (i = 0; i < le; i++) { 457 d = Geometry.distPointSegment( 458 p, 459 [1, dataX[hits[i].idx1], dataY[hits[i].idx1]], 460 [1, dataX[hits[i].idx2], dataY[hits[i].idx2]] 461 ); 462 if (d < tol) { 463 return true; 464 } 465 } 466 return false; 467 }, 468 469 /** 470 * Starting at an initial point the curve is traced with a Euler-Newton method. 471 * After tracing in one direction the algorithm stops if the component is a closed loop. 472 * Otherwise, the curved is traced in the opposite direction, starting from 473 * the same initial point. Finally, the two components are glued together. 474 * 475 * @param {Array} u0 Initial point in homogenous coordinates [1, x, y]. 476 * @returns Array [dataX, dataY] containing a new component. 477 * @private 478 */ 479 traceComponent: function (u0) { 480 var dataX = [], 481 dataY = [], 482 arr = []; 483 484 // Trace in first direction 485 // console.log("---- Start tracing forward ---------") 486 arr = this.tracing(u0, 1); 487 488 if (arr.length === 0) { 489 // console.log("Could not start tracing due to singularity") 490 } else { 491 // console.log("Trace from", [arr[0][0], arr[1][0]], "to", [arr[0][arr[0].length - 1], arr[1][arr[1].length - 1]], 492 // "num points:", arr[0].length); 493 dataX = arr[0]; 494 dataY = arr[1]; 495 } 496 497 // Trace in the other direction 498 if (!arr[2]) { 499 // No loop in the first tracing step, 500 // now explore the other direction. 501 502 // console.log("---- Start tracing backward ---------") 503 arr = this.tracing(u0, -1); 504 505 if (arr.length === 0) { 506 // console.log("Could not start backward tracing due to singularity") 507 } else { 508 // console.log("Trace backwards from", [arr[0][0], arr[1][0]], "to", 509 // [arr[0][arr[0].length - 1], arr[1][arr[1].length - 1]], "num points:", arr[0].length); 510 dataX = arr[0].reverse().concat(dataX.slice(1)); 511 dataY = arr[1].reverse().concat(dataY.slice(1)); 512 } 513 } 514 515 if (dataX.length > 0 && dataX.length < 6) { 516 // Solitary point 517 dataX.push(dataX[dataX.length - 1]); 518 dataY.push(dataY[dataY.length - 1]); 519 } 520 521 return [dataX, dataY]; 522 }, 523 524 /** 525 * Starting at a point <i>u0</i>, this routine traces the curve <i>f(u)=0</i> until 526 * a loop is detected, a critical point is reached, the curve leaves the bounding box, 527 * or the maximum number of points is reached. 528 * <p> 529 * The method is a predictor / corrector method consisting of Euler and Newton steps 530 * together with step width adaption. 531 * <p> 532 * The algorithm is an adaption of the algorithm in 533 * Eugene L. Allgower, Kurt Georg: <i>Introduction to Numerical Continuation methods.</i> 534 * 535 * @param {Array} u0 Starting point in homogenous coordinates [1, x, y]. 536 * @param {Number} direction 1 or -1 537 * @returns Array [pathX, pathY, loop_closed] or [] 538 * @private 539 */ 540 tracing: function (u0, direction) { 541 var u = [], 542 ulast = [], 543 len, 544 v = [], 545 v_start = [], 546 w = [], 547 t_u, t_v, t_u_0, tloc, 548 A, 549 grad, 550 nrm, 551 dir, 552 steps = 0, 553 k = 0, 554 loop_closed = false, 555 k0, k1, denom, dist, progress, 556 kappa, delta, alpha, 557 factor, 558 point_added = false, 559 560 quasi = false, 561 cusp_or_bifurc = false, 562 kappa_0 = this.config.kappa_0, // Inverse of planned number of Newton steps 563 delta_0 = this.config.delta_0, // Distance of predictor point to curve 564 alpha_0 = this.config.alpha_0, // Angle between two successive tangents 565 h = this.config.h_initial, 566 max_steps = this.config.max_steps, 567 568 omega = direction, 569 pathX = [], 570 pathY = [], 571 572 T = [], // Gosper's loop detector table 573 n, m, i, e; 574 575 u = u0.slice(1); 576 pathX.push(u[0]); 577 pathY.push(u[1]); 578 579 t_u = this.tangent(u); 580 if (t_u === false) { 581 // We don't want to start at a singularity. 582 // Get out of here and search for another starting point. 583 return []; 584 } 585 A = [this.dfx(u[0], u[1]), this.dfy(u[0], u[1])]; 586 587 do { 588 589 if (quasi) { 590 t_u = this.tangent_A(A); 591 } else { 592 t_u = this.tangent(u); 593 } 594 if (t_u === false) { 595 u = v.slice(); 596 pathX.push(u[0]); 597 pathY.push(u[1]); 598 // console.log("-> Bail out: t_u undefined."); 599 break; 600 } 601 602 if (pathX.length === 1) { 603 // Store first point 604 t_u_0 = t_u.slice(); 605 } else if (pathX.length === 2) { 606 T.push(pathX.length - 1); // Put first point into Gosper table T 607 608 } else if (point_added && pathX.length > 2 && !cusp_or_bifurc) { 609 610 // Detect if loop has been closed 611 dist = Geometry.distPointSegment( 612 [1, u[0], u[1]], 613 [1, pathX[0], pathY[0]], 614 [1, pathX[1], pathY[1]] 615 ); 616 617 if (dist < this.config.loop_dist * h && 618 Mat.innerProduct(t_u, t_u_0, 2) > this.config.loop_dir 619 ) { 620 621 // console.log("Loop detected after", steps, "steps"); 622 // console.log("\t", "v", v, "u0:", u0) 623 // console.log("\t", "Dist(v, path0)", dist, config.loop_dist * h) 624 // console.log("\t", "t_u", t_u); 625 // console.log("\t", "inner:", Mat.innerProduct(t_u, t_u_0, 2)); 626 // console.log("\t", "h", h); 627 628 u = u0.slice(1); 629 pathX.push(u[0]); 630 pathY.push(u[1]); 631 632 loop_closed = true; 633 break; 634 } 635 636 // Gosper's loop detector 637 if (this.config.loop_detection) { 638 n = pathX.length - 1; 639 // console.log("Check Gosper", n); 640 m = Math.floor(Mat.log2(n)); 641 642 for (i = 0; i <= m; i++) { 643 dist = Geometry.distPointSegment( 644 [1, u[0], u[1]], 645 [1, pathX[T[i] - 1], pathY[T[i] - 1]], 646 [1, pathX[T[i]], pathY[T[i]]] 647 ); 648 649 if (dist < this.config.loop_dist * h) { 650 // console.log("!!!!!!!!!!!!!!! GOSPER LOOP CLOSED !!!!", i, n + 1, 651 // this.config.loop_dist * h 652 // ); 653 654 t_v = this.tangent([pathX[T[i]], pathY[T[i]]]); 655 if (Mat.innerProduct(t_u, t_v, 2) > this.config.loop_dir) { 656 // console.log("!!!!!!!!!!!!!!! angle is good enough"); 657 break; 658 } 659 } 660 } 661 if (i <= m) { 662 loop_closed = true; 663 break; 664 } 665 666 m = 1; 667 e = 0; 668 for (i = 0; i < 100; i++) { 669 if ((n + 1) % m !== 0) { 670 break; 671 } 672 m *= 2; 673 e++; 674 } 675 // console.log("Add at e", e); 676 T[e] = n; 677 } 678 679 } 680 681 // Predictor step 682 // if (true /*h < 2 * this.config.h_initial*/) { 683 // Euler 684 // console.log("euler") 685 v[0] = u[0] + h * omega * t_u[0]; 686 v[1] = u[1] + h * omega * t_u[1]; 687 // } else { 688 // // Heun 689 // // console.log("heun") 690 // v[0] = u[0] + h * omega * t_u[0]; 691 // v[1] = u[1] + h * omega * t_u[1]; 692 693 // t_v = this.tangent(v); 694 // v[0] = 0.5 * u[0] + 0.5 * (v[0] + h * omega * t_v[0]); 695 // v[1] = 0.5 * u[1] + 0.5 * (v[1] + h * omega * t_v[1]); 696 // } 697 if (quasi) { 698 A = this.updateA(A, u, v); 699 v_start = v.slice(); 700 } 701 702 // Corrector step: Newton 703 k = 0; 704 do { 705 if (quasi) { 706 grad = A; 707 } else { 708 grad = [this.dfx(v[0], v[1]), this.dfy(v[0], v[1])]; 709 } 710 711 // Compute w = v - A(v) * f(v), 712 // grad: row vector and A(v) is the Moore-Penrose inverse: 713 // grad^T * (grad * grad^T)^(-1) 714 denom = grad[0] * grad[0] + grad[1] * grad[1]; 715 nrm = this.f(v[0], v[1]) / denom; 716 717 w[0] = v[0] - grad[0] * nrm; 718 w[1] = v[1] - grad[1] * nrm; 719 if (k === 0) { 720 k0 = Math.abs(nrm) * Math.sqrt(denom); 721 } else if (k === 1) { 722 k1 = Math.abs(nrm) * Math.sqrt(denom); 723 } 724 725 v[0] = w[0]; 726 v[1] = w[1]; 727 k++; 728 } while (k < 20 && 729 Math.abs(this.f(v[0], v[1])) > this.config.tol_newton 730 ); 731 732 delta = k0; 733 if (k > 1) { 734 kappa = k1 / k0; 735 } else { 736 kappa = 0.0; 737 } 738 739 if (quasi) { 740 A = this.updateA(A, v_start, v); 741 t_v = this.tangent_A(A); 742 } else { 743 t_v = this.tangent(v); 744 } 745 746 dir = Mat.innerProduct(t_u, t_v, 2); 747 dir = Math.max(-1, Math.min(1, dir)); 748 alpha = Math.acos(dir); 749 750 // Look for simple bifurcation points and cusps 751 cusp_or_bifurc = false; 752 progress = Geometry.distance(u, v, 2); 753 if (progress < this.config.tol_progress) { 754 u = v.slice(); 755 pathX.push(u[0]); 756 pathY.push(u[1]); 757 // console.log("-> Bail out, no progress", progress, steps); 758 break; 759 760 } else if (dir < 0.0) { 761 if (h > this.config.h_critical) { 762 // 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); 763 764 } else { 765 766 cusp_or_bifurc = true; 767 if (this.isBifurcation(u, this.config.tol_cusp)) { 768 // console.log(steps, "bifurcation point between", u, "and", v, ":", dir, "h", h, "alpha", alpha); 769 // A = [dfx(v[0], v[1]), dfy(v[0], v[1])]; 770 omega *= (-1); 771 // If there is a bifurcation point, we 772 // ignore the angle alpha for subsequent step length 773 // adaption. Because then we might be able to 774 // "jump over the critical point" 775 alpha = 0; 776 } else { 777 // Cusp or something more weird 778 u = v.slice(); 779 pathX.push(u[0]); 780 pathY.push(u[1]); 781 // console.log("-> Bail out, cusp") 782 break; 783 } 784 } 785 } 786 787 // Adapt stepwidth 788 if (!cusp_or_bifurc) { 789 factor = Math.max( 790 Math.sqrt(kappa / kappa_0), 791 Math.sqrt(delta / delta_0), 792 alpha / alpha_0 793 ); 794 if (isNaN(factor)) { 795 factor = 1; 796 } 797 factor = Math.max(Math.min(factor, 2), 0.5); 798 h /= factor; 799 h = Math.min(this.config.h_max, h); 800 801 if (factor >= 2) { 802 steps++; 803 if (steps >= 3 * max_steps) { 804 break; 805 } 806 807 point_added = false; 808 continue; 809 } 810 } 811 812 u = v.slice(); 813 pathX.push(u[0]); 814 pathY.push(u[1]); 815 point_added = true; 816 817 steps++; 818 } while ( 819 steps < max_steps && 820 u[0] >= this.bbox[0] && 821 u[1] <= this.bbox[1] && 822 u[0] <= this.bbox[2] && 823 u[1] >= this.bbox[3] 824 ); 825 826 // Clipping to bounding box, last may be outside, interpolate between second last und last point 827 len = pathX.length; 828 ulast = [pathX[len - 2], pathY[len - 2]]; 829 830 // If u[0] is outside x-interval in bounding box, interpolate to the box. 831 if (u[0] < this.bbox[0]) { 832 if (u[0] !== ulast[0]) { 833 tloc = (this.bbox[0] - ulast[0]) / (u[0] - ulast[0]); 834 if (u[1] !== ulast[1]) { 835 u[1] = ulast[1] + tloc * (u[1] - ulast[1]); 836 } 837 } 838 u[0] = this.bbox[0]; 839 } 840 if (u[0] > this.bbox[2]) { 841 if (u[0] !== ulast[0]) { 842 tloc = (this.bbox[2] - ulast[0]) / (u[0] - ulast[0]); 843 if (u[1] !== ulast[1]) { 844 u[1] = ulast[1] + tloc * (u[1] - ulast[1]); 845 } 846 } 847 u[0] = this.bbox[2]; 848 } 849 850 // If u[1] is outside y-interval in bounding box, interpolate to the box. 851 if (u[1] < this.bbox[3]) { 852 if (u[1] !== ulast[1]) { 853 tloc = (this.bbox[3] - ulast[1]) / (u[1] - ulast[1]); 854 if (u[0] !== ulast[0]) { 855 u[0] = ulast[0] + tloc * (u[0] - ulast[0]); 856 } 857 } 858 u[1] = this.bbox[3]; 859 } 860 if (u[1] > this.bbox[1]) { 861 if (u[1] !== ulast[1]) { 862 tloc = (this.bbox[1] - ulast[1]) / (u[1] - ulast[1]); 863 if (u[0] !== ulast[0]) { 864 u[0] = ulast[0] + tloc * (u[0] - ulast[0]); 865 } 866 } 867 u[1] = this.bbox[1]; 868 } 869 870 // Update last point 871 pathX[len - 1] = u[0]; 872 pathY[len - 1] = u[1]; 873 874 // if (!loop_closed) { 875 // console.log("No loop", steps); 876 // } else { 877 // console.log("Loop", steps); 878 // } 879 880 return [pathX, pathY, loop_closed]; 881 }, 882 883 /** 884 * If both eigenvalues of the Hessian are different from zero, the critical point at u 885 * is a simple bifurcation point. 886 * 887 * @param {Array} u Critical point [x, y] 888 * @param {Number} tol Tolerance of the eigenvalues to be zero. 889 * @returns Boolean True if the point is a simple bifurcation point. 890 * @private 891 */ 892 isBifurcation: function (u, tol) { 893 // Former experiments: 894 // If the Hessian has exactly one zero eigenvalue, 895 // we claim that there is a cusp. 896 // Otherwise, we decide that there is a bifurcation point. 897 // In the latter case, if both eigenvalues are zero 898 // this is a somewhat crude decision. 899 // 900 var h = Mat.eps * Mat.eps * 100, 901 x, y, a, b, c, d, ad, 902 lbda1, lbda2, 903 dis; 904 905 x = u[0]; 906 y = u[1]; 907 a = 0.5 * (this.dfx(x + h, y) - this.dfx(x - h, y)) / h; 908 b = 0.5 * (this.dfx(x, y + h) - this.dfx(x, y - h)) / h; 909 c = 0.5 * (this.dfy(x + h, y) - this.dfy(x - h, y)) / h; 910 d = 0.5 * (this.dfy(x, y + h) - this.dfy(x, y - h)) / h; 911 912 // c = b 913 ad = a + d; 914 dis = ad * ad - 4 * (a * d - b * c); 915 lbda1 = 0.5 * (ad + Math.sqrt(dis)); 916 lbda2 = 0.5 * (ad - Math.sqrt(dis)); 917 918 // console.log(a, b, c, d) 919 // console.log("Eigenvals u:", lbda1, lbda2, tol); 920 921 if (Math.abs(lbda1) > tol && Math.abs(lbda2) > tol) { 922 // if (lbda1 * lbda2 > 0) { 923 // console.log("Seems to be isolated singularity at", u); 924 // } 925 return true; 926 } 927 928 return false; 929 }, 930 931 /** 932 * Search in an arc around a critical point for a further point on the curve. 933 * Unused for the moment. 934 * 935 * @param {Array} u Critical point [x, y] 936 * @param {Array} t_u Tangent at u 937 * @param {Number} r Radius 938 * @param {Number} omega angle 939 * @returns {Array} Coordinates [x, y] of a new point. 940 * @private 941 */ 942 handleCriticalPoint: function (u, t_u, r, omega) { 943 var a = Math.atan2(omega * t_u[1], omega * t_u[0]), 944 // s = a - 0.75 * Math.PI, 945 // e = a + 0.75 * Math.PI, 946 f_circ = function (t) { 947 var x = u[0] + r * Math.cos(t), 948 y = u[1] + r * Math.sin(t); 949 return this.f(x, y); 950 }, 951 x, y, t0; 952 953 // t0 = Numerics.fzero(f_circ, [s, e]); 954 t0 = Numerics.root(f_circ, a); 955 956 x = u[0] + r * Math.cos(t0); 957 y = u[1] + r * Math.sin(t0); 958 // console.log("\t", "result", x, y, "f", f(x, y)); 959 960 return [x, y]; 961 }, 962 963 /** 964 * Quasi-Newton update of the Moore-Penrose inverse. 965 * See (7.2.3) in Allgower, Georg. 966 * 967 * @param {Array} A 968 * @param {Array} u0 969 * @param {Array} u1 970 * @returns Array 971 * @private 972 */ 973 updateA: function (A, u0, u1) { 974 var s = [u1[0] - u0[0], u1[1] - u0[1]], 975 y = this.f(u1[0], u1[1]) - this.f(u0[0], u0[1]), 976 nom, denom; 977 978 denom = s[0] * s[0] + s[1] * s[1]; 979 nom = y - (A[0] * s[0] + A[1] * s[1]); 980 nom /= denom; 981 A[0] += nom * s[0]; 982 A[1] += nom * s[1]; 983 984 return A; 985 }, 986 987 /** 988 * Approximate tangent (of norm 1) with Quasi-Newton method 989 * @param {Array} A 990 * @returns Array 991 * @private 992 */ 993 tangent_A: function (A) { 994 var t = [-A[1], A[0]], 995 nrm = Mat.norm(t, 2); 996 997 if (nrm < Mat.eps) { 998 // console.log("Approx. Singularity", t, "is zero", nrm); 999 } 1000 return [t[0] / nrm, t[1] / nrm]; 1001 }, 1002 1003 /** 1004 * Tangent of norm 1 at point u. 1005 * @param {Array} u Point [x, y] 1006 * @returns Array 1007 * @private 1008 */ 1009 tangent: function (u) { 1010 var t = [-this.dfy(u[0], u[1]), this.dfx(u[0], u[1])], 1011 nrm = Mat.norm(t, 2); 1012 1013 if (nrm < Mat.eps * Mat.eps) { 1014 // console.log("Singularity", t, "is zero", "at", u, ":", nrm); 1015 return false; 1016 } 1017 return [t[0] / nrm, t[1] / nrm]; 1018 } 1019 } 1020 1021 ); 1022 1023 export default Mat.ImplicitPlot; 1024 1025