Dibujar una catenaria
Al hacer las cuerdas entre los globos de Balloon Platform Defense, tuve que averiguar cómo dibujar una cuerda que cuelga entre dos puntos, la forma conocida en general como catenaria. Todos los ejemplos que encontré en internet o imponían restricciones al problema, como que los dos puntos estuvieran a la misma altura, o daban por sabidos ciertos datos, como el ángulo con el que empieza la cuerda, que aquí no estaban disponibles: solo se conocen los puntos inicial y final de la cuerda y su longitud. Wikipedia ofrece mucha información sobre las ecuaciones de la catenaria, y de ahí partí. La ecuación de la catenaria es:y=a\cosh \left (\frac{x}{a} \right )=\frac{a\left ( e^{\frac{x}{a}}+e^{-\frac{x}{a}} \right )}{2}
Esto supone que el punto más bajo de la cuerda está donde cruza el eje y. En la práctica, a las variables x e y se les suman constantes, ambas aún por determinar, porque no sabemos dónde estará el punto más bajo; es uno de los factores que hacen esto más complicado que todos los ejemplos de los libros de texto. (¿Te has fijado en que los problemas matemáticos de la vida real siempre son mucho más complicados que los ejemplos de los libros? Todavía espero el problema real con una integral que se pueda resolver analíticamente.) Pero eso tendrá que esperar, porque primero hay que determinar el valor de a, la constante de la ecuación que decide básicamente lo estrecha o ancha que es la curva. Cuidado: si los puntos están en el mismo lugar del eje x, a será infinito, así que conviene haber tratado ese caso antes de llegar aquí, o incluso el caso en que están muy cerca y a es incalculablemente grande. (Ten en cuenta que «incalculablemente grande» significa aquí cualquier número que, usado como argumento de la función exponencial, la haría devolver NaN o infinito con la precisión de coma flotante que estés usando.) Y, por supuesto, evita que tu código llegue hasta aquí si los dos puntos están más separados que la longitud de la cuerda: el resultado será malo.
Calcular el factor de escala
Para calcular a, Wikipedia da la ecuación (basada en una propiedad de la catenaria: mientras que cosh da la posición de la curva, sinh da su longitud):\sqrt{s^{2}-v^{2}}=2a\sinh \left ( \frac{h}{2a} \right )
donde s es la longitud de la cuerda, y h y v son las distancias horizontal y vertical (en valor absoluto) entre los puntos inicial y final. Todos son valores conocidos, así que solo queda una incógnita, pero hay que hallarla numéricamente. Primero lo intenté sin preocuparme mucho del valor inicial, y que convergiera o no era cuestión de suerte. Por suerte, resulta fácil encontrar un buen valor inicial. La serie de Taylor del seno hiperbólico es:
\sinh x=x+\frac{x^{3}}{3!}+\frac{x^{5}}{5!}+\cdots
Si hacemos las sustituciones u = \frac{1}{4a^{2}} y c = \sqrt{s^{2}-v^{2}}, la ecuación que hay que resolver queda:
c=\frac{1}{\sqrt{u}}\sinh \left ( h\sqrt{u} \right )
que, tomando los tres primeros términos de la serie de Taylor de sinh y simplificando, da:
c=\frac{1}{\sqrt{u}}\left (h\sqrt{u}+\frac{\left ( h\sqrt{u} \right )^{3}}{3!} +\frac{\left ( h\sqrt{u} \right )^{5}}{5!} \right )
c=\frac{1}{\sqrt{u}}\left (h\sqrt{u}+\frac{h^{3}u\sqrt{u}}{3!} +\frac{h^{5}u^{2}\sqrt{u}}{5!} \right )
c=h+\frac{h^{3}u}{3!} +\frac{h^{5}u^{2}}{5!}
que, reordenada, queda:
\frac{h^{5}}{120}u^{2}+\frac{h^{3}}{6}u+\left ( h-c \right )=0
Es una sencilla ecuación de segundo grado en u, y al sustituir sus coeficientes en la fórmula cuadrática del colegio x = \frac{-b\pm \sqrt{b^{2}-4ac}}{2a}
obtenemos un buen valor inicial para u y, a partir de él, para a.
u = \frac{-\frac{1}{6}h^{3}+ \sqrt{\frac{1}{36}h^{6}-\frac{1}{30}h^{5}\left ( h-c \right )}}{\frac{1}{60}h^{5}}, donde a=\frac{1}{2\sqrt{u}}
Este valor inicial está lo bastante cerca como para hallar una solución con el método de Newton. Escribiéndolo en la forma f\left ( a \right )=0
tenemos f\left ( a \right )=2a \sinh \left ( \frac{h}{2a} \right )-c
y {f}'\left ( a \right )=2 \sinh \left ( \frac{h}{2a} \right )-\frac{h}{a}\cosh \left ( \frac{h}{2a} \right )
Como a_{n+1}=a_{n}-\frac{f\left ( a_{n} \right )}{{f}'\left ( a_{n} \right )}
,
a_{n+1}=a_{n}-\frac{2a \sinh \left ( \frac{h}{2a} \right )-c}{2 \sinh \left ( \frac{h}{2a} \right )-\frac{h}{a}\cosh \left ( \frac{h}{2a} \right )}=a_{n}-\frac{a \sinh \left ( \frac{h}{2a} \right )-0.5c}{ \sinh \left ( \frac{h}{2a} \right )-\frac{h}{2a}\cosh \left ( \frac{h}{2a} \right )}
Ahora mismo compruebo que la sucesión converja a menos de 0,001 de su valor anterior, lo que suele ocurrir en 2 a 4 iteraciones, aunque a veces hacen falta más de diez cuando las coordenadas x están más cerca. Probablemente se podrían usar métodos de Householder de orden superior sin mucho más cálculo, ya que la mayor parte del trabajo está en calcular las funciones hiperbólicas, y las derivadas sucesivas de la función solo deberían dar más múltiplos de ellas. Ahora que conocemos la escala, podemos calcular la traslación.
Calcular la traslación
La ecuación estándary=a\cosh \left (\frac{x}{a} \right )
supone que la catenaria es simétrica respecto al eje y y tiene su punto más bajo en (0,a), mientras que la que dibujamos puede estar en cualquier parte de la pantalla, así que hay que calcular una traslación (p,q) que la lleve a la posición correcta, con lo que la ecuación queda:
y-q=a\cosh \left (\frac{x-p}{a} \right )
(Recordemos que en este punto a ya es conocido y x e y son variables, así que solo p y q siguen siendo incógnitas.) Mi primera idea fue sustituir en la ecuación los valores conocidos de x e y, es decir, los dos extremos de la cuerda. Aunque en teoría funciona, las ecuaciones resultantes, a pesar de mis esfuerzos por simplificarlas, a menudo producían durante el cálculo valores que las variables de coma flotante de doble precisión no podían manejar. Al final tuve que buscar otro enfoque, que además resultó más eficiente.
Recordemos que, al cambiar cosh por sinh en la ecuación de la catenaria, se obtiene la longitud de la curva. Sustituyendo en esta forma las coordenadas x de los extremos izquierdo y derecho, junto con la longitud deseada de la cuerda (s), obtenemos la ecuación:
a\sinh \left (\frac{x_{right}-p}{a} \right )-a\sinh \left (\frac{x_{left}-p}{a} \right )=s
Una de las identidades de las funciones hiperbólicas dice que:
\sinh x - \sinh y=2\cosh\left ( \frac{x+y}{2} \right ) \sinh \left ( \frac{x-y}{2} \right )
Por tanto:
\sinh \left (\frac{x_{right}-p}{a} \right )-\sinh \left (\frac{x_{left}-p}{a} \right )=2\cosh\left ( \frac{x_{right}+x_{left}-2p}{2a} \right ) \sinh \left ( \frac{x_{right}-x_{left}}{2a} \right )
En el argumento de cosh, las coordenadas x derecha e izquierda se suman y se dividen entre 2, lo que da la coordenada del punto medio, así que podemos sustituirla directamente. La resta del argumento de sinh da el mismo valor que antes llamamos h, de modo que ese sinh también se calcula como una cantidad conocida, y la ecuación original se reordena así:
\cosh\left ( \frac{x_{middle}-p}{a} \right ) = \frac{s}{2a\sinh \left ( \frac{h}{2a}\right )}
Ahora podemos despejar p con el arcocoseno hiperbólico, pero con cuidado. Como cosh es una función par, siempre obtendremos una respuesta positiva, pero también hay una negativa, lo que da dos valores posibles para p.
p =x_{middle}\pm a\cosh^{-1}\left ( \frac{s}{2a\sinh \left ( \frac{h}{2a}\right )} \right )
Para saber cuál queremos, podemos mirar las coordenadas y de los extremos. Si experimentas con un trozo de cuerda verás que, para la misma longitud y las mismas coordenadas x, el punto más bajo de la curva (si es que está entre las coordenadas x) estará a la derecha del centro (y por tanto corresponderá al arcocoseno positivo) si la coordenada y derecha es más baja que la izquierda, y al revés.
Otro problema que puede surgir es que, cuando las coordenadas y son iguales, el punto más bajo de la curva debería estar justo en el centro, así que el argumento de cosh debería ser cero y el valor del que se toma el arcocoseno debería ser uno. Por las imprecisiones de la coma flotante, a menudo sale un poco menor que uno, lo que da un error al calcular el arcocoseno. Conviene comprobar si las coordenadas y son iguales (o casi iguales) antes de empezar este cálculo, para ahorrar cálculo y evitar este error, y en ese caso suponer simplemente que p está exactamente a medio camino entre los extremos.
Ahora que a y p son conocidos, q se calcula sustituyendo los valores de x e y de uno de los extremos en la ecuación
y-q=a\cosh \left (\frac{x-p}{a} \right )
y despejando:
q=y-a\cosh \left (\frac{x-p}{a} \right )
Calcula q con esto y ya tienes la ecuación final con la que dibujar la curva.
y=a\cosh \left (\frac{x-p}{a} \right )+q
Dibujar la curva ahora es sencillo. Durante la mayoría de los pasos del cálculo de los parámetros calculé las funciones hiperbólicas con la mayor precisión posible, porque era difícil saber cómo se propagaría cualquier imprecisión. Pero en este último paso, que requiere la mayoría de los cálculos, sabrás lo grande que es el argumento que usas y cuánta imprecisión puedes tolerar, y podrás prever si una serie de Taylor hasta cierta potencia será lo bastante precisa; en mi experiencia, casi siempre merece la pena usarla.