<p>As mentioned in the comments on <a href="https://math.stackexchange.com/q/4666525">your Math.SE question</a> there <em>is</em> an analytical function of <span class="math-container">$x$</span> in terms of <span class="math-container">$y$</span>, but it's enormous, and totally impractical.</p>
<p>As James K mentions, it's <em>much</em> easier to calculate <em>C</em> values for various <span class="math-container">$(x, y)$</span> pairs in the region of interest. Then we just need to interpolate between those points to plot the desired contours.</p>
<p>In regions where the function behaves nicely, it may be sufficient to use linear interpolation between the known points. For higher accuracy, we can use bisection. Or we could use the equation's partial derivatives with Newton's method to find a point with an exact <span class="math-container">$C$</span> value from a nearby starting point. Also, once we have a point on a particular contour we can &quot;step&quot; along the contour.</p>
<p>Modern contour plotting packages tend to use bisection, since it's simple and reasonably fast for the accuracy needed in plotting. Bisection gives 1 more bit of precision per iteration loop. Newton's method doubles the number of correct bits per loop, if the initial approximation is close enough to the true value, but of course it needs the derivatives, or a numerical approximation of the derivatives.</p>
<hr />
<p>I assume <span class="math-container">$m$</span> is a dimensionless mass parameter, in the range <span class="math-container">$0 \le m \le 1$</span>.</p>
<p>Here's a typical plot, created using <a href="http://www.sagemath.org/" rel="nofollow noreferrer">SageMath</a>, a mathematics software system based on Python. It uses Matplotlib for 2D plotting. This plot is for <span class="math-container">$m=0.3, x_1= 1.2, x_2=0.8$</span>.</p>
<p><a href="https://i.sstatic.net/3vWL0.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/3vWL0.png" alt="Contour plot" /></a></p>
<p>The SageMath <code>contour_plot</code> function is quite versatile, you can read about it <a href="https://doc.sagemath.org/html/en/reference/plotting/sage/plot/contour_plot.html#sage.plot.contour_plot.contour_plot" rel="nofollow noreferrer">here</a>.</p>
<p>Here's the <a href="https://sagecell.sagemath.org/?z=eJxNUt-L1DAQfl_Y_2HIy6Wa7W27inIQEAXxRVhUEHy4JW3TvUja1HS613r4vztJ22X7kM588_ubYYzBL-3d7qKtKw1OULoW3eB7MC3gk4ZxN0FnVatXwOsevSlRV1AaXw5WeTjsCleRn3eF1c12A_T99AZRt1BMcPwK-TfTniHf54d0f0jzd9sNY2y7MU3nPEKjsLMOrSnSbgoSqJ6q4nZzUZ7fjWK6S7abD6ZF7VVJcKVrOPFG7tODgDGTWZrTPyf9vYC-U618I0IbIdepcxTYy2y_J5v5q0mKxrJRnfyurS7ReU7l0tJZ5wnteSKAaqjBomRembZwz0yAVYW2kkU3arpjSUyEI8ov2NiPbuTs00LgYutReZQZVUbXybfhrzuZ3eezVV20_KxsrwWoAd1p6CqFC5Q8bDczl54PCUjo_3jkfIQdDMljDq9hesyT2aMmMx-vINxDEPNXPCPvJuiej1kSsWbWQugcbJ2AJ0MZdoG6mcDZcr0GKu5Ve9Y8zjNPM88i6DBKO1T6pNsqUi1_-EEvfaFBqym6Zi-NfEgP9T8BL7Swq5jPIlt7OZLzUvUUtsdrQYOJpUdaC59ulLU_uQrLjvrYg5hThjFoYYXyK9Vx8-ERUBtrV_j2Wm7ka5o4jIxvCDzHYwqPgLM3lTWt7iVr1G_n6Vhqrxp9y8Ux7Z_cM180U8ftP6zJo50AzpZZ0q49s-Q_jNEYQw==&amp;lang=sage" rel="nofollow noreferrer">plotting script</a> running on the SageMathCell server. You can see the available colormaps <a href="https://matplotlib.org/stable/gallery/color/colormap_reference.html#colormap-reference" rel="nofollow noreferrer">here</a>.</p>
<hr />
<p>Here's a short script that demonstrates <a href="https://sagecell.sagemath.org/?z=eJx9klFLwzAQx98L_Q5HfFi6dbXNBjpxT3vyQRERhCGTzaQzaNOZpTPpp_fSbq448aBw97vrP_9eSgiBm2LzIV-lgTvxZUrV20IhzFvJgYuiTMIAMJ60NEYoWDm4vwX2INUaWMpGSTpK2EUYEELCYLfUtGdj14vCIAy4yEHupV9UI03zGGwag8NnHl210jy3MIU84TLPqY0O0B2ha_Q8zksNEqQCvVRrQSeTg4iPjZbKUPKsHspKcRKDjI7N2svR9vAO5p7PYQj1iY4lrVniyK-X9gM18hr73Ce87gzIHJarLUUG15CJYTbu2PSx0mL5fvgoHxXawEWcGtw1DfeHc-47NZwDrRYMBrBbsE7bpjCYomwfB4_UNXTXoVqYSqv9rXhHZ_AotiYMCpRPk1EY2Az3wLDKEhYjuvRTmlYRou2nNpRa3F8VNSZcYwIXjXJ-gNofjE59yvo0w_nC15raLGpY0VasuenWjD8xBhYGM8zGHfrPTzWLvgHxOKbo&amp;lang=sage" rel="nofollow noreferrer">Newton's method for implicit functions</a>. Sage is handy for this because it can do symbolic differentiation. The algorithm uses the partial derivatives to move the point along the gradient towards the desired contour.</p>
<p>There are other ways to use Newton's method for this problem. For example, you can fix one of <span class="math-container">$x$</span> or <span class="math-container">$y$</span> and solve for the other. Or you can step towards the desired contour by alternately solving for <span class="math-container">$x$</span> and <span class="math-container">$y$</span>. That is, hold <span class="math-container">$x$</span> fixed and  do one iteration of Newton's to find <span class="math-container">$y$</span>, then hold that <span class="math-container">$y$</span> fixed and do an iteration to find a new <span class="math-container">$x$</span>. Eventually, the point will lie on the desired contour.</p>
<p>However, in some situations, Newton's method can fail to converge, as explained on <a href="https://en.wikipedia.org/wiki/Newton%27s_method#Failure_analysis" rel="nofollow noreferrer">Wikipedia</a>.</p>

CC BY-SA 4.0
2023-03-26T14:30:19.807Z
2023-03-26T16:04:47.570Z
2023-03-26T16:04:47.570Z
PM 2Ring
<blockquote>
<p>My second intention behind this post is that how Moulton did this?</p>
</blockquote>
<p>To do this the way someone would in 1914, I would start with graph paper and a lattice of points. (x,y)  Presumably I would have some sensible values of m, x_1 and x_2 to use.</p>
<p>For each point, I would measure (with dividers, for example) the distance of the point from (x1,0) and (x2,0) and the distance from the origin (which I would then square with my log tables/slide rule).  that gives me r_1,r_2 and d², which I would record in a table.</p>
<p>Then I would calculate 0.5×d2 +2(1-m)/r_1 +2m/r_2  (again with log tables or slide rule to speed up the process.  And on my graph paper I would write the resulting value at the point (x,y) using my sharpest pencil.</p>
<p>Then repeat - If this was getting boring I'd hire a computer (perhaps a graduate student) to do the calculations  (computers did exist, they just ran differently, and wanted breaks for lunch)</p>
<p>Once I had a piece of graph paper with values, I'd start plotting contours. I'd then trace these onto plain paper ready for the engravers to print.</p>
<p>The whole process might be a day's work.</p>

CC BY-SA 4.0
2023-03-26T13:22:20.540Z
Doesn't ask about finding zeros, but about contour plots but if you have a better way that would work in 1914, when animal fat candles and pens made out of feathers were still much used, I happily upvote.
James K
CC BY-SA 4.0
2023-03-26T13:25:14.053Z
But the general theme would remain the same - the contour plots would have been calculated, by hand, and plotted by hand.  Perhaps by the principle researcher, perhaps by a PhD student, perhaps by a paid computer.
James K
CC BY-SA 4.0
2023-03-26T13:29:21.973Z
I've made an edit and reverse my vote, how does it look?
uhoh
CC BY-SA 4.0
2023-03-26T13:38:14.970Z
@PM2Ring at the very end of the OP's question, that is described explicitly as a "second intention behind this post" which pretty much says it's not the primary intention. I've edited this answer slightly at the beginning to make it clearer that that's what this answer is addressing, and reversed my vote, and deleted my original comment.
uhoh
CC BY-SA 4.0
2023-03-26T13:41:57.633Z
@JamesK I've made an edit and reverse my vote, how does it look? Per Wikipedia's [Zero-velocity surface (of the Three-body problem)](https://en.wikipedia.org/wiki/Zero-velocity_surface#Three-body_problem) to find the zero-velocity surface, one would look for the funny-looking reentrant path in (x, y) where $$\omega^2 (x^2 + y^2) + 2 \left( \frac{\mu_1}{r_1} + \frac{\mu_2}{r_2} \right) - C = 0$$. In general it's common to look for solutions to equalities by setting LHS - RHS = 0 and then going "looking for zeros".
uhoh
CC BY-SA 4.0
2023-03-26T08:34:37.317Z
2023-03-26T13:29:01.347Z
2023-03-26T13:29:01.347Z
James K
<p>(I am having physics background; I am not an astronomer/astrophysicist.)</p>
<p>How to solve a <strong>single</strong> <strong>nonlinear</strong> algebraic equation in <strong>two variables</strong>, <span class="math-container">$x$</span> and <span class="math-container">$y$</span>?</p>
<p>(I know that - if there are two variables, one needs two different equations to solve them.)</p>
<p>Equation is this (It represents zero-velocity contours in x-y plane in restricted circular 3-body problem):</p>
<p><span class="math-container">$\frac{1}{2} (x^2+y^2) + \frac{2(1-m)}{r_1} + \frac{2m}{r_2} = C$</span></p>
<p>where
<span class="math-container">$m$</span> and <span class="math-container">$C$</span> are constants and <span class="math-container">$x,y$</span> are two variables and</p>
<p>where</p>
<p><span class="math-container">$r_1 = \sqrt {(x-x_1)^2+y^2}$</span>  and</p>
<p><span class="math-container">$r_2 = \sqrt {(x-x_2)^2+y^2}$</span></p>
<p>and where <span class="math-container">$x_1$</span> and <span class="math-container">$x_2$</span> are constants - i.e. two fixed values of variable <span class="math-container">$x$</span>.</p>
<p>Is it possible to solve this equation analytically; if not, how to solve it numerically?</p>
<p>Basically my intention is to find all those <span class="math-container">$(x,y)$</span> points in <span class="math-container">$x-y$</span> plane such that LHS of this equation becomes RHS i.e. constant <span class="math-container">$C$</span>;  and that is what is the solution is.</p>
<p>IMPORTANT : This equation is from Forest Ray Moulton's a classic and a great book on celestial mechanics, which was published in 1914. That time there were no computers. Still author has drawn the contours (i.e. solution of this equation) in his book ! Its really amazing. My second intention behind this post is that how Moulton did this ?</p>

CC BY-SA 4.0
2023-03-26T08:38:43.977Z
Please do not cross-post. https://math.stackexchange.com/q/4666525
PM 2Ring
CC BY-SA 4.0
2023-03-26T11:37:59.087Z
*Welcome to Astronomy SE!* It's a great question, but I recommend you edit this one and the one in Math SE to be somewhat different from each other, in each one add a link to the other. The Math question can ask about the analytical solution, this one could focus on more practical numerical solutions. Can you use a bit of Python? Only slightly related: [Why shouldn't we illustrate spacecraft trajectories on top of static zero velocity pseudo-potential surfaces?](https://space.stackexchange.com/q/49002/12102) The ugly way would be to use Python to plot the LHS and ask it to make a contour plot
uhoh
CC BY-SA 4.0
2023-03-26T11:43:41.893Z
with $C$ as one of the contours (or something like that). Then you can query the contour line to get the points on it. But much better to use simplex or steepest decent or some minimizer to find one zero velocity point, then find a nearby point a small distance away then keep stepping around until you end up near where you started. If Python is helpful I could post a scrip to do it, or more likely @PM2Ring will beat me to it :-)
uhoh
CC BY-SA 4.0
2023-03-26T13:43:53.043Z
@PM2Ring excellent! I suspected you might ;-)
uhoh
<celestial-mechanics><n-body-simulations>
How to find zero velocity contours for circular restricted 3-body problem in x-y plane?
atom
# How to find zero velocity contours for circular restricted 3-body problem in x-y plane?

(Asked by: username_0 on 2023-03-26)

(I am having physics background; I am not an astronomer/astrophysicist.)

How to solve a **single** **nonlinear** algebraic equation in **two variables**, $x$ and $y$?

(I know that - if there are two variables, one needs two different equations to solve them.)

Equation is this (It represents zero-velocity contours in x-y plane in restricted circular 3-body problem):

$\frac{1}{2} (x^2+y^2) + \frac{2(1-m)}{r\_1} + \frac{2m}{r\_2} = C$

where
$m$ and $C$ are constants and $x,y$ are two variables and

where

$r\_1 = \sqrt {(x-x\_1)^2+y^2}$ and

$r\_2 = \sqrt {(x-x\_2)^2+y^2}$

and where $x\_1$ and $x\_2$ are constants - i.e. two fixed values of variable $x$.

Is it possible to solve this equation analytically; if not, how to solve it numerically?

Basically my intention is to find all those $(x,y)$ points in $x-y$ plane such that LHS of this equation becomes RHS i.e. constant $C$; and that is what is the solution is.

IMPORTANT : This equation is from Forest Ray Moulton's a classic and a great book on celestial mechanics, which was published in 1914. That time there were no computers. Still author has drawn the contours (i.e. solution of this equation) in his book ! Its really amazing. My second intention behind this post is that how Moulton did this ?

**Tags:** celestial-mechanics, n-body-simulations
> **Comments:**
> **username_1**: Please do not cross-post. https://math.stackexchange.com/q/4666525  
> **username_2**: *Welcome to Astronomy SE!* It's a great question, but I recommend you edit this one and the one in Math SE to be somewhat different from each other, in each one add a link to the other. The Math question can ask about the analytical solution, this one could focus on more practical numerical solutions. Can you use a bit of Python? Only slightly related: [Why shouldn't we illustrate spacecraft trajectories on top of static zero velocity pseudo-potential surfaces?](https://space.stackexchange.com/q/49002/12102) The ugly way would be to use Python to plot the LHS and ask it to make a contour plot  
> **username_2**: with $C$ as one of the contours (or something like that). Then you can query the contour line to get the points on it. But much better to use simplex or steepest decent or some minimizer to find one zero velocity point, then find a nearby point a small distance away then keep stepping around until you end up near where you started. If Python is helpful I could post a scrip to do it, or more likely @PM2Ring will beat me to it :-)  
> **username_2**: @PM2Ring excellent! I suspected you might ;-)  

---

## 2 answers

### Answer by username_1 on 2023-03-26. Score: 3 | Accepted answer

As mentioned in the comments on [your Math.SE question](https://math.stackexchange.com/q/4666525) there *is* an analytical function of $x$ in terms of $y$, but it's enormous, and totally impractical.

As James K mentions, it's *much* easier to calculate *C* values for various $(x, y)$ pairs in the region of interest. Then we just need to interpolate between those points to plot the desired contours.

In regions where the function behaves nicely, it may be sufficient to use linear interpolation between the known points. For higher accuracy, we can use bisection. Or we could use the equation's partial derivatives with Newton's method to find a point with an exact $C$ value from a nearby starting point. Also, once we have a point on a particular contour we can "step" along the contour.

Modern contour plotting packages tend to use bisection, since it's simple and reasonably fast for the accuracy needed in plotting. Bisection gives 1 more bit of precision per iteration loop. Newton's method doubles the number of correct bits per loop, if the initial approximation is close enough to the true value, but of course it needs the derivatives, or a numerical approximation of the derivatives.

---

I assume $m$ is a dimensionless mass parameter, in the range $0 \le m \le 1$.

Here's a typical plot, created using [SageMath](http://www.sagemath.org/), a mathematics software system based on Python. It uses Matplotlib for 2D plotting. This plot is for $m=0.3, x\_1= 1.2, x\_2=0.8$.

[![Contour plot](https://i.sstatic.net/3vWL0.png)](https://i.sstatic.net/3vWL0.png)

The SageMath `contour_plot` function is quite versatile, you can read about it [here](https://doc.sagemath.org/html/en/reference/plotting/sage/plot/contour_plot.html#sage.plot.contour_plot.contour_plot).

Here's the [plotting script](https://sagecell.sagemath.org/?z=eJxNUt-L1DAQfl_Y_2HIy6Wa7W27inIQEAXxRVhUEHy4JW3TvUja1HS613r4vztJ22X7kM588_ubYYzBL-3d7qKtKw1OULoW3eB7MC3gk4ZxN0FnVatXwOsevSlRV1AaXw5WeTjsCleRn3eF1c12A_T99AZRt1BMcPwK-TfTniHf54d0f0jzd9sNY2y7MU3nPEKjsLMOrSnSbgoSqJ6q4nZzUZ7fjWK6S7abD6ZF7VVJcKVrOPFG7tODgDGTWZrTPyf9vYC-U618I0IbIdepcxTYy2y_J5v5q0mKxrJRnfyurS7ReU7l0tJZ5wnteSKAaqjBomRembZwz0yAVYW2kkU3arpjSUyEI8ov2NiPbuTs00LgYutReZQZVUbXybfhrzuZ3eezVV20_KxsrwWoAd1p6CqFC5Q8bDczl54PCUjo_3jkfIQdDMljDq9hesyT2aMmMx-vINxDEPNXPCPvJuiej1kSsWbWQugcbJ2AJ0MZdoG6mcDZcr0GKu5Ve9Y8zjNPM88i6DBKO1T6pNsqUi1_-EEvfaFBqym6Zi-NfEgP9T8BL7Swq5jPIlt7OZLzUvUUtsdrQYOJpUdaC59ulLU_uQrLjvrYg5hThjFoYYXyK9Vx8-ERUBtrV_j2Wm7ka5o4jIxvCDzHYwqPgLM3lTWt7iVr1G_n6Vhqrxp9y8Ux7Z_cM180U8ftP6zJo50AzpZZ0q49s-Q_jNEYQw==&lang=sage) running on the SageMathCell server. You can see the available colormaps [here](https://matplotlib.org/stable/gallery/color/colormap_reference.html#colormap-reference).

---

Here's a short script that demonstrates [Newton's method for implicit functions](https://sagecell.sagemath.org/?z=eJx9klFLwzAQx98L_Q5HfFi6dbXNBjpxT3vyQRERhCGTzaQzaNOZpTPpp_fSbq448aBw97vrP_9eSgiBm2LzIV-lgTvxZUrV20IhzFvJgYuiTMIAMJ60NEYoWDm4vwX2INUaWMpGSTpK2EUYEELCYLfUtGdj14vCIAy4yEHupV9UI03zGGwag8NnHl210jy3MIU84TLPqY0O0B2ha_Q8zksNEqQCvVRrQSeTg4iPjZbKUPKsHspKcRKDjI7N2svR9vAO5p7PYQj1iY4lrVniyK-X9gM18hr73Ce87gzIHJarLUUG15CJYTbu2PSx0mL5fvgoHxXawEWcGtw1DfeHc-47NZwDrRYMBrBbsE7bpjCYomwfB4_UNXTXoVqYSqv9rXhHZ_AotiYMCpRPk1EY2Az3wLDKEhYjuvRTmlYRou2nNpRa3F8VNSZcYwIXjXJ-gNofjE59yvo0w_nC15raLGpY0VasuenWjD8xBhYGM8zGHfrPTzWLvgHxOKbo&lang=sage). Sage is handy for this because it can do symbolic differentiation. The algorithm uses the partial derivatives to move the point along the gradient towards the desired contour.

There are other ways to use Newton's method for this problem. For example, you can fix one of $x$ or $y$ and solve for the other. Or you can step towards the desired contour by alternately solving for $x$ and $y$. That is, hold $x$ fixed and do one iteration of Newton's to find $y$, then hold that $y$ fixed and do an iteration to find a new $x$. Eventually, the point will lie on the desired contour.

However, in some situations, Newton's method can fail to converge, as explained on [Wikipedia](https://en.wikipedia.org/wiki/Newton%27s_method#Failure_analysis).


---

### Answer by username_3 on 2023-03-26. Score: 5

> My second intention behind this post is that how Moulton did this?

To do this the way someone would in 1914, I would start with graph paper and a lattice of points. (x,y) Presumably I would have some sensible values of m, x\_1 and x\_2 to use.

For each point, I would measure (with dividers, for example) the distance of the point from (x1,0) and (x2,0) and the distance from the origin (which I would then square with my log tables/slide rule). that gives me r\_1,r\_2 and d², which I would record in a table.

Then I would calculate 0.5×d2 +2(1-m)/r\_1 +2m/r\_2 (again with log tables or slide rule to speed up the process. And on my graph paper I would write the resulting value at the point (x,y) using my sharpest pencil.

Then repeat - If this was getting boring I'd hire a computer (perhaps a graduate student) to do the calculations (computers did exist, they just ran differently, and wanted breaks for lunch)

Once I had a piece of graph paper with values, I'd start plotting contours. I'd then trace these onto plain paper ready for the engravers to print.

The whole process might be a day's work.

> **Comments:**
> **username_3**: Doesn't ask about finding zeros, but about contour plots but if you have a better way that would work in 1914, when animal fat candles and pens made out of feathers were still much used, I happily upvote.  
> **username_3**: But the general theme would remain the same - the contour plots would have been calculated, by hand, and plotted by hand.  Perhaps by the principle researcher, perhaps by a PhD student, perhaps by a paid computer.  
> **username_2**: I've made an edit and reverse my vote, how does it look?  
> **username_2**: @PM2Ring at the very end of the OP's question, that is described explicitly as a "second intention behind this post" which pretty much says it's not the primary intention. I've edited this answer slightly at the beginning to make it clearer that that's what this answer is addressing, and reversed my vote, and deleted my original comment.  
> **username_2**: @JamesK I've made an edit and reverse my vote, how does it look? Per Wikipedia's [Zero-velocity surface (of the Three-body problem)](https://en.wikipedia.org/wiki/Zero-velocity_surface#Three-body_problem) to find the zero-velocity surface, one would look for the funny-looking reentrant path in (x, y) where $$\omega^2 (x^2 + y^2) + 2 \left( \frac{\mu_1}{r_1} + \frac{\mu_2}{r_2} \right) - C = 0$$. In general it's common to look for solutions to equalities by setting LHS - RHS = 0 and then going "looking for zeros".  

---

