Continuation procedure to solve for a 2D curve that satisfies f(x,y) = 0

I have some function of $R^2$, that must be numerically computed. For instance, I might be interested in a real-valued contour integral that begins from (x,y) = 0.

$$f(x,y) = \Re\left[\int_0^{x + iy} t^3 + t^5 \, dt \right]$$

where $\Re$ is for the real part. (In practice I don't have the explicit form for f and it must be numerically computed). The contours $f = 0$ look like: I know that there exists some continuous and smooth curve that leaves the origin, and for all points along this curve $x = x(s)$, $y = y(s)$, that $f = 0$ where $s$ is is some parameterization.

I also know a good initial guess to use for the initial direction of this line. For instance, I might know that initially, the curve leaves along some direction, $\theta$.

How do I numerically solve for the curve?

1. One way is to just compute $f$ over a grid of values in $x$ and in $y$. Then create a contour plot and interpolate for all constant-value contours satisfying f = 0. This is not a computationally friendly way.

2. Another way might be to start at (x,y) = (0,0) and take a small step in the direction of your guess of where the curve lies. Then from this new point, take another small step, but in a direction that minimizes $f$.

The problem with 2. is that I'm not sure how to code something adaptive (I'd like to take variable step sizes, particularly around regions where the curve may be highly-curved).

I'm sure there must be well-known numerical methods for these problems. This is really just root-finding but with the advantage that you know the solution set (that you are interested in) must be continuous and smooth. Can someone lend a hand?

• This question has some notational trouble. $z$ can't both be a limit of integration and the dummy variable of integration. Also, $x$ and $y$ don't appear in the right hand side of your equation. I have some guesses as to what you might mean, but can you clarify? Sep 8 '13 at 4:09
• @BillBarth: Thank you. I've edited to give a simple example.
– TSGM
Sep 8 '13 at 8:35

Let us parameterize the curve you are looking for by $(x(t),y(t))$ and let us for a moment assume that at all points on this curve $\nabla f(x(t),y(t)) \neq 0$, i.e., the curve never intersects another isocontour line. Then you know that at each point the tangent to the curve $(x(t),y(t))$ is parallel to the isocontour levels of $f$, i.e., that it is perpendicular to the gradient of $f$. That means that one possible parameterization would be to define the curve as $$\frac d{dt} \left({x(t) \atop y(t)}\right) = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \nabla f(x(t),y(t))$$ which is an ODE that you can integrate up from an arbitrarily chosen point $(x(0),y(0))$ that you know lies on the curve. Another possible choice for the time parameterization would be to use the form $$\frac d{dt} \left({x(t) \atop y(t)}\right) = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \frac{\nabla f(x(t),y(t))}{\|\nabla f(x(t),y(t))\|}.$$

In either case, you will get into trouble once you get to a point where the gradient of $f$ is zero, i.e., where two zero lines cross. You'll have to treat this case separately.

• Two quick questions: (i) what is the strategy at points where the gradient is zero? (ii) what is the advantage of forcing 't' to be measured via arclength? Is there any advantage?
– TSGM
Sep 9 '13 at 14:57
• (i) It's difficult. I don't have much experience with this, but I'd try making a large step over the critical point and then projecting back to the line where $f=0$. (ii) Not really, other than the fact that the choice of step length $\Delta t$ in your integrator is trivial. Sep 10 '13 at 3:59

I answered a similar question here. That method works if you have a good way of solving the segment intersection query (easiest if your $f$ returns a signed value so that you can use sign-based bisection).

Otherwise, I would try analysis. You can compute the gradient of $f$ with respect to $x$ and $y$ presumably. It's best to carry out the differentiation by hand, but you could approximate it with a finite difference. You could possibly compute higher order derivatives to estimate curvature, but this rapidly gets annoying. The curve goes in the directions of zero gradient. You will have to decide on a large-ish step size (which you can then refine) and then compute the correction to get back onto the curve (either using the line segment query from above (in the direction perpendicular to the step) or using something like Newton's method).

Note that in order to perform adaptive sampling, you have to have some idea of the maximum curvature that is possible or that you'll be happy with. If you know the spatial extent of your curve, then you can take a steps of more reasonable size.

• This was a helpful answer that led me to the right path. Your first suggestion related to the older post is unfortunately too complicated for me. I ended up learning a little bit about pseudo-arclength continuation methods (basically, a Newton scheme that walks along the arclength of the desired curve). This works decently.
– TSGM
Sep 9 '13 at 9:07
• Currently, my scheme does not work particularly well near points where there are multiple level curves intersecting. You mentioned that one could do the corrector step by either solving using Newton's method (which I am doing now) or by solving along a perpendicular line segment. I'm trying to find some more information about the differences between these two approaches and whether one is better than the other. But it's hard to find any details. Can you advise?
– TSGM
Sep 10 '13 at 9:28
• It's an extremely difficult case, so there probably isn't much to say about it. Points of near degeneracy are notoriously difficult in theory and in practice. The "solve along perpendicular line" approach is just some ad hoc thing I came up with, since it's much easier to doing repeated Delaunay triangulations. I think the triangulation method will properly sample your curves as long as you specify the minimum spacing between curves you expect (it's a provably good sampling algorithm). Sep 10 '13 at 17:54
• Unfortunately, the triangulation algorithm you linked to is a bit too complicated for me to decipher. You wouldn't happen to know a more readable version, would you? (particular for 2D applications and not 3D applications).
– TSGM
Sep 10 '13 at 18:54
• @TGSM: I found the original paper by Boissonat to be pretty clear, but I've read a lot of computational geometry papers. I don't think there's a better explanation available, but let me try to summarize the main ideas. Basically, if you can find points on the curve, you can triangulate the points, and ask whether the curve intersects the dual edges (perpendicular bisectors) of the triangulation. You can use any triangulation which is "nice", of which Delaunay is the main one. You could also sample on an adaptive rectangular grid, but that's actually harder, just thinking about it quickly. Sep 11 '13 at 7:09