Consider the unit circle $x^2 + y^2 = 1$, it also represents the solution to the ODE

$$ x + y \frac{dy}{dx} = 0; y(0) = 1 $$

Suppose we don't know how to solve the above equation analytically (there are more complex cases where analytical solution is impossible), a finite difference discretization gives: $$ y^{n+1} = y^n - \frac{x^n}{y^n} \delta, x^0 = 0, y^0 = 1. $$ $\delta$ is the step size. This method suffers from the singularity at $(x=1,y=0)$, where the slope of the circle diverges to infinity.

I am wondering if there are better methods that can recover the entire circle. A key aspect appears to be the integration direction, i.e., the sign of $\delta$, must be reversed at some point.

  • 1
    $\begingroup$ If you look at the parametrization of the circle in the form $y=y(x)$, that singularity is there too, so it's not the fault of the finite difference method. It would be better to parameterize the ODE in a different way, for example by arc length: you'd get the coupled system $y \dot y + x\dot x = 0$, $\dot y^2+\dot x^2=0$. (It just says that $(\dot x, \dot y)$ is a unit vector perpendicular to $(x,y)$, and such a vector always exists so there's no singularity.) Does this help? $\endgroup$ – Kirill Jul 13 '17 at 17:48
  • $\begingroup$ @Kirill Arclength is just the time, if one starts from the original physical problem (like the SHO) where time is still there as the independent variable, then after obtaining x=x(t), y=y(t) using any ode solver, one just have to plot y against x to have the circle. In this sense, by removing t from the picture one shoots in his own foot. It is still interesting to know if the posted equation can be solved directly instead of tracing back to the original time dependent problem. $\endgroup$ – Taozi Jul 13 '17 at 20:31

I think that you can rewrite your equation as

$$x + \frac{1}{2}\frac{d u}{d x} = 0\, ,$$

with $u = y^2$. Then, the difference equation would be

$$u^{n+1} = u^{n} - 2x \Delta x\, .$$

One could recover the original variable by computing the square root of the result.

With these changes, I managed to solved it in the following snippet

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt

#%% Solution
npts = 100
x = np.linspace(0, 1, npts)
u = np.ones_like(x)
dx = x[1] - x[0]
for cont in range(1, npts):
    u[cont] = u[cont - 1] - 2*x[cont - 1]*dx

#%% Plotting
plt.figure(figsize=(2, 2))
plt.plot(x, np.sqrt(u))
plt.xticks([0, 0.5, 1])
plt.yticks([0, 0.5, 1])

and, this is the plot

enter image description here

  • $\begingroup$ Your solution neither recovers the whole unit circle nor treats the singularity at (1,0), which is the point of the question. $\endgroup$ – Endulum Aug 12 '17 at 17:33
  • $\begingroup$ @Endulum, if you solve the equation analytically you obtain $u = k - x^2$. And that's the solution that you recover numerically, when taking the positive root you are missing half a circle, but you can't take the negative root as well. Regarding the singularity, I don't see why the iteration would face a singularity. $\endgroup$ – nicoguaro Aug 12 '17 at 18:46

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.