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 v = [], 543 v_start = [], 544 w = [], 545 t_u, t_v, t_u_0, 546 A, 547 grad, 548 nrm, 549 dir, 550 steps = 0, 551 k = 0, 552 loop_closed = false, 553 k0, k1, denom, dist, progress, 554 kappa, delta, alpha, 555 factor, 556 point_added = false, 557 558 quasi = false, 559 cusp_or_bifurc = false, 560 kappa_0 = this.config.kappa_0, // Inverse of planned number of Newton steps 561 delta_0 = this.config.delta_0, // Distance of predictor point to curve 562 alpha_0 = this.config.alpha_0, // Angle between two successive tangents 563 h = this.config.h_initial, 564 max_steps = this.config.max_steps, 565 566 omega = direction, 567 pathX = [], 568 pathY = [], 569 570 T = [], // Gosper's loop detector table 571 n, m, i, e; 572 573 u = u0.slice(1); 574 pathX.push(u[0]); 575 pathY.push(u[1]); 576 577 t_u = this.tangent(u); 578 if (t_u === false) { 579 // We don't want to start at a singularity. 580 // Get out of here and search for another starting point. 581 return []; 582 } 583 A = [this.dfx(u[0], u[1]), this.dfy(u[0], u[1])]; 584 585 do { 586 587 if (quasi) { 588 t_u = this.tangent_A(A); 589 } else { 590 t_u = this.tangent(u); 591 } 592 if (t_u === false) { 593 u = v.slice(); 594 pathX.push(u[0]); 595 pathY.push(u[1]); 596 // console.log("-> Bail out: t_u undefined."); 597 break; 598 } 599 600 if (pathX.length === 1) { 601 // Store first point 602 t_u_0 = t_u.slice(); 603 } else if (pathX.length === 2) { 604 T.push(pathX.length - 1); // Put first point into Gosper table T 605 606 } else if (point_added && pathX.length > 2 && !cusp_or_bifurc) { 607 608 // Detect if loop has been closed 609 dist = Geometry.distPointSegment( 610 [1, u[0], u[1]], 611 [1, pathX[0], pathY[0]], 612 [1, pathX[1], pathY[1]] 613 ); 614 615 if (dist < this.config.loop_dist * h && 616 Mat.innerProduct(t_u, t_u_0, 2) > this.config.loop_dir 617 ) { 618 619 // console.log("Loop detected after", steps, "steps"); 620 // console.log("\t", "v", v, "u0:", u0) 621 // console.log("\t", "Dist(v, path0)", dist, config.loop_dist * h) 622 // console.log("\t", "t_u", t_u); 623 // console.log("\t", "inner:", Mat.innerProduct(t_u, t_u_0, 2)); 624 // console.log("\t", "h", h); 625 626 u = u0.slice(1); 627 pathX.push(u[0]); 628 pathY.push(u[1]); 629 630 loop_closed = true; 631 break; 632 } 633 634 // Gosper's loop detector 635 if (this.config.loop_detection) { 636 n = pathX.length - 1; 637 // console.log("Check Gosper", n); 638 m = Math.floor(Mat.log2(n)); 639 640 for (i = 0; i <= m; i++) { 641 dist = Geometry.distPointSegment( 642 [1, u[0], u[1]], 643 [1, pathX[T[i] - 1], pathY[T[i] - 1]], 644 [1, pathX[T[i]], pathY[T[i]]] 645 ); 646 647 if (dist < this.config.loop_dist * h) { 648 // console.log("!!!!!!!!!!!!!!! GOSPER LOOP CLOSED !!!!", i, n + 1, 649 // this.config.loop_dist * h 650 // ); 651 652 t_v = this.tangent([pathX[T[i]], pathY[T[i]]]); 653 if (Mat.innerProduct(t_u, t_v, 2) > this.config.loop_dir) { 654 // console.log("!!!!!!!!!!!!!!! angle is good enough"); 655 break; 656 } 657 } 658 } 659 if (i <= m) { 660 loop_closed = true; 661 break; 662 } 663 664 m = 1; 665 e = 0; 666 for (i = 0; i < 100; i++) { 667 if ((n + 1) % m !== 0) { 668 break; 669 } 670 m *= 2; 671 e++; 672 } 673 // console.log("Add at e", e); 674 T[e] = n; 675 } 676 677 } 678 679 // Predictor step 680 // if (true /*h < 2 * this.config.h_initial*/) { 681 // Euler 682 // console.log("euler") 683 v[0] = u[0] + h * omega * t_u[0]; 684 v[1] = u[1] + h * omega * t_u[1]; 685 // } else { 686 // // Heun 687 // // console.log("heun") 688 // v[0] = u[0] + h * omega * t_u[0]; 689 // v[1] = u[1] + h * omega * t_u[1]; 690 691 // t_v = this.tangent(v); 692 // v[0] = 0.5 * u[0] + 0.5 * (v[0] + h * omega * t_v[0]); 693 // v[1] = 0.5 * u[1] + 0.5 * (v[1] + h * omega * t_v[1]); 694 // } 695 if (quasi) { 696 A = this.updateA(A, u, v); 697 v_start = v.slice(); 698 } 699 700 // Corrector step: Newton 701 k = 0; 702 do { 703 if (quasi) { 704 grad = A; 705 } else { 706 grad = [this.dfx(v[0], v[1]), this.dfy(v[0], v[1])]; 707 } 708 709 // Compute w = v - A(v) * f(v), 710 // grad: row vector and A(v) is the Moore-Penrose inverse: 711 // grad^T * (grad * grad^T)^(-1) 712 denom = grad[0] * grad[0] + grad[1] * grad[1]; 713 nrm = this.f(v[0], v[1]) / denom; 714 715 w[0] = v[0] - grad[0] * nrm; 716 w[1] = v[1] - grad[1] * nrm; 717 if (k === 0) { 718 k0 = Math.abs(nrm) * Math.sqrt(denom); 719 } else if (k === 1) { 720 k1 = Math.abs(nrm) * Math.sqrt(denom); 721 } 722 723 v[0] = w[0]; 724 v[1] = w[1]; 725 k++; 726 } while (k < 20 && 727 Math.abs(this.f(v[0], v[1])) > this.config.tol_newton 728 ); 729 730 delta = k0; 731 if (k > 1) { 732 kappa = k1 / k0; 733 } else { 734 kappa = 0.0; 735 } 736 737 if (quasi) { 738 A = this.updateA(A, v_start, v); 739 t_v = this.tangent_A(A); 740 } else { 741 t_v = this.tangent(v); 742 } 743 744 dir = Mat.innerProduct(t_u, t_v, 2); 745 dir = Math.max(-1, Math.min(1, dir)); 746 alpha = Math.acos(dir); 747 748 // Look for simple bifurcation points and cusps 749 cusp_or_bifurc = false; 750 progress = Geometry.distance(u, v, 2); 751 if (progress < this.config.tol_progress) { 752 u = v.slice(); 753 pathX.push(u[0]); 754 pathY.push(u[1]); 755 // console.log("-> Bail out, no progress", progress, steps); 756 break; 757 758 } else if (dir < 0.0) { 759 if (h > this.config.h_critical) { 760 // 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); 761 762 } else { 763 764 cusp_or_bifurc = true; 765 if (this.isBifurcation(u, this.config.tol_cusp)) { 766 // console.log(steps, "bifurcation point between", u, "and", v, ":", dir, "h", h, "alpha", alpha); 767 // A = [dfx(v[0], v[1]), dfy(v[0], v[1])]; 768 omega *= (-1); 769 // If there is a bifurcation point, we 770 // ignore the angle alpha for subsequent step length 771 // adaption. Because then we might be able to 772 // "jump over the critical point" 773 alpha = 0; 774 } else { 775 // Cusp or something more weird 776 u = v.slice(); 777 pathX.push(u[0]); 778 pathY.push(u[1]); 779 // console.log("-> Bail out, cusp") 780 break; 781 } 782 } 783 } 784 785 // Adapt stepwidth 786 if (!cusp_or_bifurc) { 787 factor = Math.max( 788 Math.sqrt(kappa / kappa_0), 789 Math.sqrt(delta / delta_0), 790 alpha / alpha_0 791 ); 792 if (isNaN(factor)) { 793 factor = 1; 794 } 795 factor = Math.max(Math.min(factor, 2), 0.5); 796 h /= factor; 797 h = Math.min(this.config.h_max, h); 798 799 if (factor >= 2) { 800 steps++; 801 if (steps >= 3 * max_steps) { 802 break; 803 } 804 805 point_added = false; 806 continue; 807 } 808 } 809 810 u = v.slice(); 811 pathX.push(u[0]); 812 pathY.push(u[1]); 813 point_added = true; 814 815 steps++; 816 } while ( 817 steps < max_steps && 818 u[0] >= this.bbox[0] && 819 u[1] <= this.bbox[1] && 820 u[0] <= this.bbox[2] && 821 u[1] >= this.bbox[3] 822 ); 823 824 // if (!loop_closed) { 825 // console.log("No loop", steps); 826 // } else { 827 // console.log("Loop", steps); 828 // } 829 830 return [pathX, pathY, loop_closed]; 831 }, 832 833 /** 834 * If both eigenvalues of the Hessian are different from zero, the critical point at u 835 * is a simple bifurcation point. 836 * 837 * @param {Array} u Critical point [x, y] 838 * @param {Number} tol Tolerance of the eigenvalues to be zero. 839 * @returns Boolean True if the point is a simple bifurcation point. 840 * @private 841 */ 842 isBifurcation: function (u, tol) { 843 // Former experiments: 844 // If the Hessian has exactly one zero eigenvalue, 845 // we claim that there is a cusp. 846 // Otherwise, we decide that there is a bifurcation point. 847 // In the latter case, if both eigenvalues are zero 848 // this is a somewhat crude decision. 849 // 850 var h = Mat.eps * Mat.eps * 100, 851 x, y, a, b, c, d, ad, 852 lbda1, lbda2, 853 dis; 854 855 x = u[0]; 856 y = u[1]; 857 a = 0.5 * (this.dfx(x + h, y) - this.dfx(x - h, y)) / h; 858 b = 0.5 * (this.dfx(x, y + h) - this.dfx(x, y - h)) / h; 859 c = 0.5 * (this.dfy(x + h, y) - this.dfy(x - h, y)) / h; 860 d = 0.5 * (this.dfy(x, y + h) - this.dfy(x, y - h)) / h; 861 862 // c = b 863 ad = a + d; 864 dis = ad * ad - 4 * (a * d - b * c); 865 lbda1 = 0.5 * (ad + Math.sqrt(dis)); 866 lbda2 = 0.5 * (ad - Math.sqrt(dis)); 867 868 // console.log(a, b, c, d) 869 // console.log("Eigenvals u:", lbda1, lbda2, tol); 870 871 if (Math.abs(lbda1) > tol && Math.abs(lbda2) > tol) { 872 // if (lbda1 * lbda2 > 0) { 873 // console.log("Seems to be isolated singularity at", u); 874 // } 875 return true; 876 } 877 878 return false; 879 }, 880 881 /** 882 * Search in an arc around a critical point for a further point on the curve. 883 * Unused for the moment. 884 * 885 * @param {Array} u Critical point [x, y] 886 * @param {Array} t_u Tangent at u 887 * @param {Number} r Radius 888 * @param {Number} omega angle 889 * @returns {Array} Coordinates [x, y] of a new point. 890 * @private 891 */ 892 handleCriticalPoint: function (u, t_u, r, omega) { 893 var a = Math.atan2(omega * t_u[1], omega * t_u[0]), 894 // s = a - 0.75 * Math.PI, 895 // e = a + 0.75 * Math.PI, 896 f_circ = function (t) { 897 var x = u[0] + r * Math.cos(t), 898 y = u[1] + r * Math.sin(t); 899 return this.f(x, y); 900 }, 901 x, y, t0; 902 903 // t0 = Numerics.fzero(f_circ, [s, e]); 904 t0 = Numerics.root(f_circ, a); 905 906 x = u[0] + r * Math.cos(t0); 907 y = u[1] + r * Math.sin(t0); 908 // console.log("\t", "result", x, y, "f", f(x, y)); 909 910 return [x, y]; 911 }, 912 913 /** 914 * Quasi-Newton update of the Moore-Penrose inverse. 915 * See (7.2.3) in Allgower, Georg. 916 * 917 * @param {Array} A 918 * @param {Array} u0 919 * @param {Array} u1 920 * @returns Array 921 * @private 922 */ 923 updateA: function (A, u0, u1) { 924 var s = [u1[0] - u0[0], u1[1] - u0[1]], 925 y = this.f(u1[0], u1[1]) - this.f(u0[0], u0[1]), 926 nom, denom; 927 928 denom = s[0] * s[0] + s[1] * s[1]; 929 nom = y - (A[0] * s[0] + A[1] * s[1]); 930 nom /= denom; 931 A[0] += nom * s[0]; 932 A[1] += nom * s[1]; 933 934 return A; 935 }, 936 937 /** 938 * Approximate tangent (of norm 1) with Quasi-Newton method 939 * @param {Array} A 940 * @returns Array 941 * @private 942 */ 943 tangent_A: function (A) { 944 var t = [-A[1], A[0]], 945 nrm = Mat.norm(t, 2); 946 947 if (nrm < Mat.eps) { 948 // console.log("Approx. Singularity", t, "is zero", nrm); 949 } 950 return [t[0] / nrm, t[1] / nrm]; 951 }, 952 953 /** 954 * Tangent of norm 1 at point u. 955 * @param {Array} u Point [x, y] 956 * @returns Array 957 * @private 958 */ 959 tangent: function (u) { 960 var t = [-this.dfy(u[0], u[1]), this.dfx(u[0], u[1])], 961 nrm = Mat.norm(t, 2); 962 963 if (nrm < Mat.eps * Mat.eps) { 964 // console.log("Singularity", t, "is zero", "at", u, ":", nrm); 965 return false; 966 } 967 return [t[0] / nrm, t[1] / nrm]; 968 } 969 } 970 971 ); 972 973 export default Mat.ImplicitPlot; 974 975