How can I resolve issues in my chaotic billiards simulation using Python?

[carlosbgois]

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.

 

[DrClaude]

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.

 

[carlosbgois]

Here is the code: (I didn't know how to use code tags before, sorry)

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))

 

[DrClaude]

Maybe I can be of help after all. I'm not sure your approach is the best, but let's 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$?

 

[carlosbgois]

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.