I need to refine a level curve representing an implicit function that is not numerically solvable in sympy. The initial solution set generated by matplotlib's contour function has an accuracy of around ±0.02. The inaccuracy is due in part to the steep slope about the contour. However the inaccuracy of erf, catastrophic cancellation of exponential values, and floating point inaccuracies all play possible roles, which I hope to counter using arbitrarily-precise floats from mpmath. I'll know if the plot of this level curve and another intersect above zero.

To refine the curve I slice the function at regular intervals and solve for the roots of each slice using scipy's root_scalar function. Each slice has two horizontal asymptotes, see the example below. I've encountered several errors that I don't understand and I would like to automate the entire process.

slice points of interest

Right now I'm doing this manually, but it would be nice to slice the contour line generated by matplotlib for initial points. I have no idea how to do this.

Using derivative-based methods, initial points to the right of $B$ march up the right asymptote until the tolerance is exceeded. Providing the solver with the second derivative does not invalidate the solution; even if it did the solver would not search to the left of the initial point. Also iteration terminates on $A$ or $B$ where the derivative is zero. Instead I try bracketing the solution to within the green lines, $[g_1,g_2]$. Note that $f(x_1) \approx 10^{-17}$ but $f(g_2) \approx 10^{-100}$. And so every bracket solver returns $g_2$ every time. I don't understand why. The "bisect" solver should test the intervals $[g_1, (g_1+g_2)/2], [(g_1+g_2)/2, g_2]$ and the second interval should fail because the results do not differ in sign. The result of every iteration should differ in sign but brackets near $g_2$ cannot differ in sign.

How should I extend whichever process to find all roots within the blue region $[b_1, b_2]$? The minimum $A$ is not a fixed point across slices. Notice that $A$ is much closer to $x_1$ than $x_0$ so that a small error in the initial estimate of $x_1$ could easily find $x_0$ instead. Most slice also have no roots.

  • $\begingroup$ Re "However the inaccuracy of erf, catastrophic cancellation of exponential values": Could you clarify how these two items contribute to the overall problem? Have you looked at using erfc (complementary error function) and / or erfcx (exponentially scaled complementary error function) instead? $\endgroup$
    – njuffa
    Mar 8 at 21:37


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.