Eine Kettenlinie zeichnen

Als ich in Balloon Platform Defense Schnüre zwischen Ballons zeichnen wollte, musste ich herausfinden, wie man eine Schnur zeichnet, die zwischen zwei Punkten hängt – die Form, die allgemein als Kettenlinie bekannt ist. Alle Beispiele, die ich im Internet fand, schränkten das Problem entweder ein, etwa indem beide Punkte auf gleicher Höhe lagen, oder setzten Wissen voraus, das hier nicht verfügbar war, etwa den Winkel, unter dem die Schnur beginnt. Bekannt sind nur der Anfangs- und der Endpunkt der Schnur und ihre Länge. Wikipedia bietet viele Informationen über die Gleichungen der Kettenlinie, und davon bin ich ausgegangen. Die Gleichung der Kettenlinie lautet:
y=a\cosh \left (\frac{x}{a} \right )=\frac{a\left ( e^{\frac{x}{a}}+e^{-\frac{x}{a}} \right )}{2}

Dabei wird angenommen, dass der tiefste Punkt der Schnur dort liegt, wo sie die y-Achse schneidet. In der Praxis kommen zu x und y noch Konstanten hinzu, die beide erst bestimmt werden müssen, denn wir wissen nicht, wo der tiefste Punkt liegen wird – einer der Gründe, warum das hier komplizierter ist als alle Lehrbuchbeispiele. (Ist dir schon aufgefallen, dass mathematische Probleme aus dem echten Leben immer viel komplizierter sind als die Beispiele im Lehrbuch? Ich warte immer noch auf das Problem aus dem echten Leben mit einem Integral, das sich analytisch lösen lässt.) Das muss aber warten, denn zuerst müssen wir den Wert von a bestimmen, der Konstante in der Gleichung, die im Wesentlichen festlegt, wie schmal oder breit die Kurve ist. Achtung: Liegen die Punkte an derselben Stelle der x-Achse, ist a unendlich. Diesen Fall solltest du also behandelt haben, bevor du so weit kommst – ebenso den Fall, dass sie sehr nahe beieinander liegen und a unberechenbar groß wird. (Bedenke, dass „unberechenbar groß“ hier jede Zahl ist, die als Argument der Exponentialfunktion bei der verwendeten Gleitkommagenauigkeit NaN oder Unendlich liefern würde.) Und natürlich sollte dein Code gar nicht erst so weit kommen, wenn die beiden Punkte weiter auseinander liegen, als die Schnur lang ist – das Ergebnis wäre unbrauchbar.

Den Skalierungsfaktor berechnen

Um a zu berechnen, gibt Wikipedia folgende Gleichung an (sie beruht auf einer Eigenschaft der Kettenlinie: Während cosh die Lage der Kurve angibt, gibt sinh ihre Länge an):
\sqrt{s^{2}-v^{2}}=2a\sinh \left ( \frac{h}{2a} \right )

Dabei ist s die Länge der Schnur, und h und v sind der (betragsmäßige) waagerechte und senkrechte Abstand zwischen Anfangs- und Endpunkt. Das sind alles bekannte Werte, es bleibt also nur eine Unbekannte, doch die muss numerisch bestimmt werden. Zuerst habe ich das versucht, ohne auf den Startwert zu achten, und ob das Verfahren konvergierte, war Glückssache. Zum Glück lässt sich leicht ein guter Startwert finden. Die Taylorreihe des Sinus hyperbolicus lautet:
\sinh x=x+\frac{x^{3}}{3!}+\frac{x^{5}}{5!}+\cdots

Mit den Substitutionen u = \frac{1}{4a^{2}} und c = \sqrt{s^{2}-v^{2}} wird die zu lösende Gleichung zu:
c=\frac{1}{\sqrt{u}}\sinh \left ( h\sqrt{u} \right )

Nimmt man die ersten drei Glieder der Taylorreihe für sinh und vereinfacht, ergibt sich:
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!}

und umgestellt:

\frac{h^{5}}{120}u^{2}+\frac{h^{3}}{6}u+\left ( h-c \right )=0

Das ist eine einfache quadratische Gleichung in u, und setzt man ihre Koeffizienten in die aus der Schule bekannte a-b-c-Formel
x = \frac{-b\pm \sqrt{b^{2}-4ac}}{2a}
ein, erhält man einen guten Startwert für u und damit für 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}}, wobei a=\frac{1}{2\sqrt{u}}
Dieser Startwert liegt nah genug, um mit dem Newton-Verfahren eine Lösung zu finden. In der Form
f\left ( a \right )=0
ist
f\left ( a \right )=2a \sinh \left ( \frac{h}{2a} \right )-c
und
{f}'\left ( a \right )=2 \sinh \left ( \frac{h}{2a} \right )-\frac{h}{a}\cosh \left ( \frac{h}{2a} \right )
Da
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 )}
Derzeit prüfe ich, ob die Folge bis auf 0,001 an ihren vorherigen Wert herankommt. Meist geschieht das nach 2 bis 4 Iterationen, manchmal dauert es aber über zehn, wenn die x-Koordinaten näher beieinander liegen. Householder-Verfahren höherer Ordnung wären vermutlich ohne viel zusätzlichen Rechenaufwand möglich, denn der Großteil der Rechenarbeit steckt in den Hyperbelfunktionen, und weitere Ableitungen der Funktion sollten nur weitere Vielfache davon ergeben. Da wir nun die Skalierung kennen, können wir die Verschiebung berechnen.

Die Verschiebung berechnen

Die Standardgleichung
y=a\cosh \left (\frac{x}{a} \right )
geht davon aus, dass die Kettenlinie symmetrisch zur y-Achse liegt und ihren tiefsten Punkt bei (0,a) hat, während die, die wir zeichnen, irgendwo auf dem Bildschirm liegen kann. Wir müssen also eine Verschiebung (p,q) bestimmen, die sie an die richtige Stelle bringt. Unsere Gleichung wird damit zu:
y-q=a\cosh \left (\frac{x-p}{a} \right )
(Zur Erinnerung: a ist an dieser Stelle bekannt, x und y sind Variablen, unbekannt sind also nur noch p und q.) Mein erster Gedanke war, in die Gleichung die bekannten Werte von x und y einzusetzen, also die beiden Endpunkte der Schnur. Theoretisch funktioniert das, aber die entstehenden Gleichungen lieferten trotz aller Vereinfachungsversuche im Lauf der Rechnung oft Zwischenwerte, mit denen Gleitkommazahlen doppelter Genauigkeit nicht zurechtkamen. Am Ende musste ich einen anderen Weg finden, der sich ohnehin als effizienter herausstellte. Erinnern wir uns: Ersetzt man in der Gleichung der Kettenlinie cosh durch sinh, erhält man die Länge der Kurve. Setzt man in diese Form die x-Koordinaten des linken und rechten Endpunkts ein, zusammen mit der bekannten gewünschten Länge der Schnur (s), ergibt sich die Gleichung:
a\sinh \left (\frac{x_{right}-p}{a} \right )-a\sinh \left (\frac{x_{left}-p}{a} \right )=s
Eine der Identitäten der Hyperbelfunktionen besagt:
\sinh x - \sinh y=2\cosh\left ( \frac{x+y}{2} \right ) \sinh \left ( \frac{x-y}{2} \right )
Daher gilt:
\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 )
Im Argument von cosh werden die rechte und die linke x-Koordinate addiert und durch 2 geteilt, das ergibt die Koordinate der Mitte, also können wir sie gleich dadurch ersetzen. Die Differenz im Argument von sinh ist derselbe Wert, den wir vorhin h genannt haben, sodass sich auch dieser sinh-Wert als bekannte Größe berechnen lässt. Die ursprüngliche Gleichung lässt sich dann umstellen zu:
\cosh\left ( \frac{x_{middle}-p}{a} \right ) = \frac{s}{2a\sinh \left ( \frac{h}{2a}\right )}
Jetzt können wir mit dem Areakosinus hyperbolicus nach p auflösen, allerdings ist Vorsicht geboten. Da cosh eine gerade Funktion ist, erhalten wir immer eine positive Lösung zurück, es gibt aber auch eine negative, also zwei mögliche Werte für p.
p =x_{middle}\pm a\cosh^{-1}\left ( \frac{s}{2a\sinh \left ( \frac{h}{2a}\right )} \right )
Welchen Wert wir brauchen, verraten die y-Koordinaten der Endpunkte. Ein Versuch mit einem Stück Schnur zeigt: Bei gleicher Schnurlänge und gleichen x-Koordinaten liegt der tiefste Punkt der Kurve (sofern er überhaupt zwischen den x-Koordinaten liegt) rechts der Mitte – gehört also zum positiven Areakosinus –, wenn die rechte y-Koordinate tiefer liegt als die linke, und umgekehrt. Ein weiteres Problem kann auftreten, wenn die y-Koordinaten gleich sind: Dann sollte der tiefste Punkt der Kurve genau in der Mitte liegen, das Argument von cosh sollte also null sein und der Wert, von dem man den Areakosinus nimmt, entsprechend eins. Wegen Gleitkomma-Ungenauigkeiten kommt oft etwas weniger als eins heraus, und der Areakosinus schlägt fehl. Es ist daher klug, vor dieser Rechnung zu prüfen, ob die y-Koordinaten gleich (oder fast gleich) sind – das spart Rechenzeit und vermeidet diesen Fehler –, und in diesem Fall einfach anzunehmen, dass p genau in der Mitte zwischen den Endpunkten liegt. Da nun a und p bekannt sind, lässt sich q berechnen, indem man die x- und y-Werte eines der Endpunkte in die Gleichung
y-q=a\cosh \left (\frac{x-p}{a} \right )
einsetzt und umstellt zu:
q=y-a\cosh \left (\frac{x-p}{a} \right )
Berechne daraus q, und du hast die endgültige Gleichung, mit der du deine Kurve zeichnen kannst.
y=a\cosh \left (\frac{x-p}{a} \right )+q
Die Kurve zu zeichnen ist jetzt einfach. Bei den meisten Schritten zur Berechnung der Parameter habe ich die Hyperbelfunktionen so genau wie möglich berechnet, weil schwer abzuschätzen war, wie sich Ungenauigkeiten fortpflanzen würden. In diesem letzten Schritt aber, der den Großteil der Rechnungen ausmacht, weißt du, wie groß das verwendete Argument ist und wie viel Ungenauigkeit du dir leisten kannst, und kannst vorhersagen, ob eine Taylorreihe bis zu einer bestimmten Potenz genau genug ist – und meiner Erfahrung nach lohnt sie sich meistens.