1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Chaotic billiards simulation

  1. Nov 30, 2014 #1
    What am I trying to do? I'm trying to implement a simulation of a chaotic billiard system, following the algorithm in this excerpt.

    How am I trying it? Using numpy and matplotlib, I implemented the following code


    What is the problem? When calculating phi_new, the equation has two solutions (assuming the boundary is convex, which is.) I must enforce that phi_new is the solution which is different from phi but I don't know how to do that. Are there more issues with the code?

    What should the output be? A phase space diagram of S x Alpha, looking like this.

    Any help is very appreciated! Thanks in advance.
  2. jcsd
  3. Dec 2, 2014 #2


    User Avatar

    Staff: Mentor

    I don't think I'll be the one to help you, but posting code as an image can almost guarantee that no one will. Repost it as text inside code tags.
  4. Dec 2, 2014 #3
    Here is the code: (I didn't know how to use code tags before, sorry)

    Code (Text):

    def boundaryFunction(parameter):
        return 1 + 0.1 * np.cos(parameter)

    def boundaryDerivative(parameter):
        return -0.1 * np.sin(parameter)

    def trajectoryFunction(parameter):
        aux = np.sin(beta - phi) / np.sin(beta - parameter)
        return boundaryFunction(phi) * aux

    def difference(parameter):
        return trajectoryFunction(parameter) - boundaryFunction(parameter)

    def integrand(parameter):
        rr = boundaryFunction(parameter)
        dd = boundaryDerivative (parameter)
        return np.sqrt(rr ** 2 + dd ** 2)

    ##### Main #####

    length_vals = np.array([], dtype=np.float64)
    alpha_vals = np.array([], dtype=np.float64)

    # nof initial phi angles, alpha angles, and nof collisions for each.
    n_phi, n_alpha, n_cols, count = 10, 10, 10, 0
    # Length of the boundary
    total_length, err = integrate.quad(integrand, 0, 2 * np.pi)
    setPlot('Polar Billiard')

    #fig = plt.figure()

    for phi in np.linspace(0, 2 * np.pi, n_phi):
        for alpha in np.linspace(0, 2 * np.pi, n_alpha):
            for n in np.arange(1, n_cols):

                nu = np.arctan(boundaryFunction(phi) / boundaryDerivative(phi))
                beta = np.pi + phi + alpha - nu

                # Determines next impact coordinate.
                bnds = (0, 2 * np.pi)
                phi_new = optimize.minimize_scalar(difference, bounds=bnds, method='bounded').x
                if phi_new == phi:
                    print 'phi_new = phi'

                nu_new =  np.arctan(boundaryFunction(phi_new) / boundaryDerivative(phi_new))
                # Reflection angle with relation to tangent.
                alpha_new = phi_new - phi + nu - nu_new - alpha
                # Arc length for current phi value.
                arc_length, err = integrate.quad(integrand, 0, phi_new)

                # Append values to list
                length_vals = np.append(length_vals, arc_length / total_length)
                alpha_vals = np.append(alpha_vals, alpha)

            count += 1
        print  "{}%" .format(100 * count / (n_phi * n_alpha))
  5. Dec 2, 2014 #4


    User Avatar

    Staff: Mentor

    Maybe I can be of help after all. I'm not sure your approach is the best, but lets assume that your sticking with optimize.minimize_scalar.

    Why don't you define the bounds as ##[\phi + \epsilon, \phi + 2 \pi - \epsilon]##, where ##\epsilon## is a small number, to avoid the solution at ##\phi##?
  6. Dec 3, 2014 #5
    That's a nice solution. Thanks you. What would be a better approach?

    Unfortunately, this alone didn't fix the problem. There must be other mistakes, though I did exactly what was in the book excerpt. I'm wondering if varying only Φ and α parameters spans the whole possible starting parameters for the problem, but am not sure about that.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook