Why is the path of my parametric brachistochrone cycloid inverted?

94 views Asked by At

I am trying to parametrize the path of a brachistochrone cycloid. Most examples show the final position at the origin, but I want to generalize the path. I have tried adapting an approach similar to the one found on the scipython blog. While it is technically accurate to call my generated path a brachistochrone cycloid, it does not have the correct concavity.

My output path

I am not sure what I did wrong. Can someone help me identify my mistake?

The output path was obtained by running the following code for a minimal working example.

`import numpy as np
from scipy.optimize import newton
from scipy.integrate import quad
import matplotlib.pyplot as plt


class BrachistochroneCycloidPath():

    def __init__(self, initial_position, final_position, g):
        super().__init__()
        self.initial_position = initial_position
        self.final_position = final_position
        self.xi = initial_position[0]
        self.yi = initial_position[1]
        self.xf = final_position[0]
        self.yf = final_position[1]
        self.g = g
        self._initial_theta = None
        self._final_theta = None
        self._radius = None
        self._arc_length = None
        self._elapsed_duration = None

    @property
    def initial_theta(self):
        return self._initial_theta

    @property
    def final_theta(self):
        return self._final_theta

    @property
    def radius(self):
        return self._radius

    @property
    def arc_length(self):
        return self._arc_length

    @property
    def elapsed_duration(self):
        return self._elapsed_duration

    def initialize_theta_and_radius(self, x0):
        initial_theta = 0
        f = lambda theta : (self.yf - self.yi) / (self.xf - self.xi) - (1 - np.cos(theta)) / (theta - np.sin(theta))
        final_theta = newton(
            f,
            x0=x0)
        radius = (self.yi - self.yf) / (1 - np.cos(final_theta))
        self._initial_theta = initial_theta
        self._final_theta = final_theta
        self._radius = radius

    def initialize_arc_length(self):
        integral = quad(
            self.ds,
            self.initial_theta,
            self.final_theta)
        arc_length = abs(
            integral[0])
        self._arc_length = arc_length       

    def initialize_elapsed_duration(self):
        f = lambda theta : self.ds(theta) / self.v(theta)
        integral = quad(
            f,
            self.initial_theta,
            self.final_theta)
        elapsed_duration = abs(
            integral[0])
        self._elapsed_duration = elapsed_duration       

    def x(self, theta):
        x = self.radius * (theta - np.sin(theta)) + self.xf
        return x

    def y(self, theta):
        y = self.radius * (1 - np.cos(theta)) + self.yf
        return y

    def dx(self, theta):
        dx = self.radius * (1 - np.cos(theta))
        return dx

    def dy(self, theta):
        dy = self.radius * np.sin(theta)
        return dy

    def ds(self, theta):
        # dy_dx = self.dy(theta) / self.dx(theta)
        # ds = np.sqrt(
        #   1 + np.square(dy_dx))
        dx = self.dx(theta)
        dy = self.dy(theta)
        ds = np.sqrt(
            np.square(dx) + np.square(dy))
        return ds

    def v(self, theta):
        v = np.sqrt(2 * self.g * self.y(theta))
        return v


if __name__ == "__main__":

    ## initialize inputs
    initial_position = (1, 10)
    final_position = (10, 1)
    # x0 = np.pi / 2 # RuntimeError
    x0 = -1 * np.pi / 2
        g = 9.8

    ## initialize path
    brachistochrone_cycloid_path = BrachistochroneCycloidPath(
        initial_position=initial_position,
        final_position=final_position,
        g=g)
    brachistochrone_cycloid_path.initialize_theta_and_radius(
        x0=x0)
    brachistochrone_cycloid_path.initialize_arc_length()
    brachistochrone_cycloid_path.initialize_elapsed_duration()

    ## print-check arc-length and elapsed duration
    print("\n\tArc-Length [m]:\n{}\n".format(
        brachistochrone_cycloid_path.arc_length))
    print("\n\tElapsed Duration [s]:\n{}\n".format(
        brachistochrone_cycloid_path.elapsed_duration))

    ## view
    theta = np.linspace(
        brachistochrone_cycloid_path.initial_theta, ## 0
        brachistochrone_cycloid_path.final_theta, ## found by newton method
        100)

    path_x = brachistochrone_cycloid_path.x(
        theta=theta)
    path_y = brachistochrone_cycloid_path.y(
        theta=theta)
    path_label = "Path of Brachistochrone Cycloid (radius = {:,.4})".format(
        brachistochrone_cycloid_path.radius)

    positions_x = [
        initial_position[0],
        final_position[0]]
    positions_y = [
        initial_position[1],
        final_position[1]]
    positions_label = "Initial/Final Position"

    fig, ax = plt.subplots()
    ax.plot(
        path_x,
        path_y,
        color="darkorange",
        label=path_label)
    ax.scatter(
        positions_x,
        positions_y,
        color="black",
        label=positions_label)
    fig.subplots_adjust(
        bottom=0.2)
    fig.legend(
        mode="expand",
        loc="lower center",
        ncol=2)
    ax.grid(
        True)
    plt.show()
    plt.close()`
0

There are 0 answers