Ejemplo de búsqueda de raíces: Billar circular

Para jugar al billar en una mesa circular, hemos de goplear la bola $Q$ de la figura con la bola $P$, tras un impacto $I$ en la banda. Conocido el radio $R$ de la mesa y las posiciones en coordenadas cartesianas (cuyo origen es el centro de la mesa) de los puntos $P(x_P, y_P)$ y $Q(x_Q, y_Q)$, el punto de impacto viene definido por el ángulo central $\theta$.

Tras un análisis geométrico del problema se prueba que los valores de $\theta$ que proporcionan los puntos del impacto posibles son los ceros de la función:

$$f(\theta) = \frac{x_P \sin\theta - y_P\cos\theta}{\sqrt{(R\cos\theta - x_P)^2+ (R \sin\theta - y_P)^2}} + \frac{x_Q \sin\theta - y_Q\cos\theta}{\sqrt{(R\cos\theta - x_Q)^2 + (R \sin\theta - y_Q)^2}}\, .$$

Consideremos un billar de radio unitario, $P=(0.6, -0.3)$ y $Q=(-0.6, -0.3)$. Determinar los puntos de impacto utilizando los métodos de bisección y Newton. Comparar los resultados obtenidos.



In [1]:

    
%matplotlib notebook



In [2]:

    
import numpy as np
from numpy import sqrt, sin, cos, pi
from scipy.optimize import bisect, newton
import matplotlib.pyplot as plt
import ipywidgets as widgets



In [3]:

    
# Configurar las gráficas
repo = "https://raw.githubusercontent.com/nicoguaro/matplotlib_styles/master"
style = repo + "/styles/minimalist.mplstyle"
plt.style.use(style)
plt.rcParams["figure.figsize"] = 6, 4



In [4]:

    
def fun(theta, R, P, Q):
    xP, yP = P
    xQ, yQ = Q
    RP = sqrt((R*cos(theta) - xP)**2 + (R*sin(theta) - yP)**2)
    RQ = sqrt((R*cos(theta) - xQ)**2 + (R*sin(theta) - yQ)**2)
    termP = (xP*sin(theta) - yP*cos(theta))/RP
    termQ = (xQ*sin(theta) - yQ*cos(theta))/RQ
    return termP + termQ



In [5]:

    
def grad(theta, R, P, Q):
    xP, yP = P
    xQ, yQ = Q
    RP = sqrt((R*cos(theta) - xP)**2 + (R*sin(theta) - yP)**2)
    RQ = sqrt((R*cos(theta) - xQ)**2 + (R*sin(theta) - yQ)**2)
    termP = (xP*cos(theta) + yP*sin(theta))/RP
    termQ = (xQ*cos(theta) + yQ*sin(theta))/RQ
    termP2 = -R*(xP*sin(theta) - yP*cos(theta))**2/RP**3
    termQ2 = -R*(xQ*sin(theta) - yQ*cos(theta))**2/RQ**3
    return termP + termQ + termP2 + termQ2



In [6]:

    
P = (0.5, -0.3)
Q = (-0.7, -0.3)
R = 1



In [7]:

    
theta = np.linspace(0, 2*pi, 201)
plt.figure()
plt.plot(theta, fun(theta, R, P, Q))
plt.hlines(0, 0, 2*pi, alpha=0.3)
plt.xlabel(r"$\theta$")
plt.ylabel(r"$f(\theta)$");



In [8]:

    
bisect(fun, 0, pi, args=(R, P, Q))









    Out[8]:





1.3971569215409811



In [9]:

    
newton(fun, 1, fprime=grad, args=(R, P, Q))









    Out[9]:





1.3971569215401547

Referencias

Julio Mulero (2019). Ecuaciones, métodos y un billar circular. Blog: El último verso de Fermat.