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

Question

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?

Answer 1

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.

Answer 2

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.