# Chaotic billiards simulation

Tags:
1. Nov 30, 2014

### 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.

2. Dec 2, 2014

### 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.

3. Dec 2, 2014

### carlosbgois

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

4. Dec 2, 2014

### 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$?

5. Dec 3, 2014

### 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.