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