[FIXED] Create a curve in which the space between points is balanced

Issue

If I want to plot the curve of a function, say 1/x, between 0.1 and 10, I can do it like this:

xx = np.linspace(0.1,10,1000)
yy = 1.0/xx
plt.plot(xx,yy)

The problem is that the spacing between points is not balanced along the curve. Specifically, at the left of the curve, where x<y, the points are very sparse (only about 10% of the points are there), and at the right of the curve, where x>y, the points are much denser (about 90% of the points are there, even though the curve is symmetric in both parts).

How can I create the arrays xx and yy (in general, not only for this particular function) such that the spacing between adjacent points is similar throughout the entire curve?

Solution

If the solution in the comments is not what you're looking for, here are a few others that choose abscissae such that the distances between points along the curve are constant. The last one is probably the most direct of the lot, but it requires some private functions, and those provate functions aren't available unless you have SciPy 1.12.

They all generate plots that are visually similar to this.

Using splines:

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline, PchipInterpolator

def f(x):
    return 1/x

a, b = 0.1, 10
n = 100  # use more points for improved accuracy

# Generate interpolant for the curve and its derivative
# If you have the derivative, use it instead of dydx_spline
x0 = np.linspace(a, b, n)
y0 = f(x0)
y_spline = CubicSpline(x0, y0, extrapolate=False)
dydx_spline = y_spline.derivative()

# Generate interpolant for arc length `s` as function of `x`
dsdx = np.sqrt(1 + dydx_spline(x0)**2)
dsdx_spline = PchipInterpolator(x0, dsdx, extrapolate=False)  # Pchip is monotonic
s_spline = dsdx_spline.antiderivative()

# Solve for `x` at specified arc distances
s = np.linspace(0, s_spline(b), 10)
x = np.asarray([a] + [s_spline.solve(si)[0] for si in s[1:-1]] + [b])

# Plot the results
y = f(x)
plt.plot(x0, y0, '-')
plt.plot(x, y, '*')

---

By finding events of an initial value problem...

# After running the solution above...
from scipy.integrate import solve_ivp

# Assume we have the derivative of the function
# If not, we can use the spline approach above or use numerical differentiation (e.g. next solution)
dydx = dydx_spline

def dsdx(x, s):
  # derivative of arc length s with respect to x
  return np.sqrt(1 + dydx(x)**2)

def event(x, s):
  # event function crosses 0 when arc length is a multiple of ds
  return (s - ds/2) % ds - ds/2
event.direction = 1

ds = s[1] - s[0]
sol = solve_ivp(dsdx, y0=np.asarray([0]), t_span=(x[0], x[-1]),
                max_step=0.001, events=event)
x = np.ravel(sol.t_events)
x = np.append(x, b)

# Plot the results
y = f(x)
plt.plot(x0, y0, '-')
plt.plot(x, y, '*')

---

Use numerical differentiation to get derivative of function f, numerical integration to get arc length, and scalar rootfinding to get abscissa at desired arc lengths.

# After running the solutions above...
from scipy.integrate._tanhsinh import _tanhsinh
from scipy.optimize._zeros_py import _chandrupatla, _differentiate

def dydx(x):  # numerical differentiation of function
    return _differentiate(f, x).df

def dsdx(x):  # derivative of arc length
    return np.sqrt(1 + dydx(x) ** 2)

def s(x):  # numerical integration to get arc length
    return _tanhsinh(dsdx, a, x).integral

# solve for `x` s.t. `s(x) == s0`
s0 = np.linspace(0, s(b), 10)  # linearly spaced arc lengths

def g(x, s0):  
    return s(x) - s0

res = _chandrupatla(g, a, b, args=(s0,))
x = res.x

# Plot the results
y = f(x)
plt.plot(x0, y0, '-')
plt.plot(x, y, '*')

(The private functions are not available in SciPy 1.11, but they will be available in SciPy 1.12, which you can get now with the nightly wheels. Note that these shouldn't be relied on for production code; they can be changed or removed without notice. Substitutes for _differentiate, _tanhsinh, and _chandrupatla are scipy.optimize.approx_fprime, scipy.integrate.quad, and scipy.optimize.root_scalar, but you'd need to use loops or apply np.vectorize.)

Answered By - Matt Haberland