Een kettinglijn tekenen

Toen ik in Balloon Platform Defense touwtjes tussen ballonnen wilde tekenen, moest ik uitzoeken hoe je een touw tekent dat tussen twee punten hangt, de vorm die meer algemeen bekendstaat als kettinglijn. Alle voorbeelden die ik online vond, legden het probleem beperkingen op, zoals dat beide punten op dezelfde hoogte liggen, of gingen ervan uit dat je bepaalde kennis had, zoals de hoek waaronder het touw begint, en die was hier niet beschikbaar: alleen het begin- en eindpunt van het touw en de lengte ervan zijn bekend. Wikipedia geeft veel informatie over de vergelijkingen van de kettinglijn, en daar ben ik van uitgegaan. De vergelijking van de kettinglijn is:
y=a\cosh \left (\frac{x}{a} \right )=\frac{a\left ( e^{\frac{x}{a}}+e^{-\frac{x}{a}} \right )}{2}

Hierbij wordt aangenomen dat het laagste punt van het touw ligt waar het de y-as snijdt. In de praktijk komen er bij x en y constanten bij, die allebei nog bepaald moeten worden, omdat we niet weten waar het laagste punt zal liggen – een van de redenen waarom dit ingewikkelder is dan alle voorbeelden uit het lesboek. (Is het je opgevallen dat wiskundeproblemen uit het echte leven altijd veel ingewikkelder zijn dan de voorbeelden in het lesboek? Ik wacht nog steeds op het probleem uit het echte leven met een integraal die analytisch op te lossen is.) Maar dat moet wachten, want eerst moeten we de waarde van a bepalen, de constante in de vergelijking die in wezen bepaalt hoe smal of breed de kromme is. Let op: als de punten op dezelfde plaats op de x-as liggen, is a oneindig, dus dat geval moet je afgehandeld hebben voordat je zover komt – en ook het geval dat ze heel dicht bij elkaar liggen en a onberekenbaar groot is. (Bedenk dat ‘onberekenbaar groot’ hier elk getal is dat, als argument van de exponentiële functie, NaN of oneindig zou opleveren bij de zwevendekommaprecisie die je gebruikt.) En laat je code natuurlijk niet zover komen als de twee punten verder uit elkaar liggen dan de lengte van het touw – het resultaat wordt dan slecht.

De schaalfactor berekenen

Om a te berekenen, geeft Wikipedia de vergelijking (gebaseerd op een eigenschap van de kettinglijn: terwijl cosh de positie van de kromme geeft, geeft sinh de lengte ervan):
\sqrt{s^{2}-v^{2}}=2a\sinh \left ( \frac{h}{2a} \right )

waarbij s de lengte van het touw is, en h en v de (absolute) horizontale en verticale afstand tussen begin- en eindpunt. Dit zijn allemaal bekende waarden, zodat er nog maar één onbekende overblijft, maar die moet numeriek worden opgelost. Eerst probeerde ik dat zonder al te kieskeurig te zijn over de startwaarde, en of het convergeerde was een gok. Gelukkig blijkt het eenvoudig om een goede startwaarde te vinden. De Taylorreeks van de sinus hyperbolicus is:
\sinh x=x+\frac{x^{3}}{3!}+\frac{x^{5}}{5!}+\cdots

Met de substituties u = \frac{1}{4a^{2}} en c = \sqrt{s^{2}-v^{2}} wordt de op te lossen vergelijking:
c=\frac{1}{\sqrt{u}}\sinh \left ( h\sqrt{u} \right )

wat, met de eerste drie termen van de Taylorreeks voor sinh en na vereenvoudigen, geeft:
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!}

en herschreven:

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

Dit is een eenvoudige vierkantsvergelijking in u, en als je de coëfficiënten invult in de wortelformule van de middelbare school
x = \frac{-b\pm \sqrt{b^{2}-4ac}}{2a}
krijg je een goede startwaarde voor u, en daarmee voor 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}}, waarbij a=\frac{1}{2\sqrt{u}}
Deze startwaarde ligt dicht genoeg bij de oplossing om die met de methode van Newton-Raphson te vinden. In de vorm
f\left ( a \right )=0
is
f\left ( a \right )=2a \sinh \left ( \frac{h}{2a} \right )-c
en
{f}'\left ( a \right )=2 \sinh \left ( \frac{h}{2a} \right )-\frac{h}{a}\cosh \left ( \frac{h}{2a} \right )
Omdat
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 )}
Op dit moment controleer ik of de rij tot op 0,001 van de vorige waarde convergeert, wat meestal na 2 tot 4 iteraties gebeurt, al zijn er soms meer dan tien nodig als de x-coördinaten dichter bij elkaar liggen. Householder-methoden van hogere orde zouden waarschijnlijk mogelijk zijn zonder veel meer rekenwerk, omdat het meeste werk in het berekenen van de hyperbolische functies zit en verdere afgeleiden van de functie alleen meer veelvouden daarvan zouden moeten opleveren. Nu we de schaal kennen, kunnen we de verschuiving berekenen.

De verschuiving berekenen

De standaardvergelijking
y=a\cosh \left (\frac{x}{a} \right )
gaat ervan uit dat de kettinglijn symmetrisch is ten opzichte van de y-as en haar laagste punt in (0,a) heeft, terwijl de lijn die wij tekenen overal op het scherm kan liggen. We moeten dus een verschuiving (p,q) bepalen die haar op de juiste plek zet, zodat onze vergelijking wordt:
y-q=a\cosh \left (\frac{x-p}{a} \right )
(Ter herinnering: a is op dit punt bekend, x en y zijn variabelen, dus alleen p en q zijn nog onbekend.) Mijn eerste gedachte was om de bekende waarden van x en y in de vergelijking in te vullen – dus de twee eindpunten van het touw. In theorie werkt dat, maar de vergelijkingen die dat oplevert, gaven ondanks al mijn pogingen om ze te vereenvoudigen tijdens het rekenen vaak tussenwaarden waar zwevendekommagetallen met dubbele precisie niet mee overweg konden. Uiteindelijk moest ik een andere aanpak vinden, die bovendien efficiënter bleek. Bedenk dat je de lengte van de kromme krijgt als je in de vergelijking van de kettinglijn cosh door sinh vervangt. Vul je in die vorm de x-coördinaten van het linker- en rechtereindpunt in, samen met de bekende gewenste lengte van het touw (s), dan krijg je de vergelijking:
a\sinh \left (\frac{x_{right}-p}{a} \right )-a\sinh \left (\frac{x_{left}-p}{a} \right )=s
Een van de identiteiten van de hyperbolische functies zegt:
\sinh x - \sinh y=2\cosh\left ( \frac{x+y}{2} \right ) \sinh \left ( \frac{x-y}{2} \right )
Dus:
\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 )
In het argument van cosh worden de rechter- en linker-x-coördinaat opgeteld en door 2 gedeeld, wat de coördinaat van het midden geeft, dus die kunnen we er net zo goed voor in de plaats zetten. Het verschil in het argument van sinh is dezelfde waarde die we eerder h noemden, zodat ook die sinh als bekende grootheid te berekenen is, en de oorspronkelijke vergelijking wordt herschreven tot:
\cosh\left ( \frac{x_{middle}-p}{a} \right ) = \frac{s}{2a\sinh \left ( \frac{h}{2a}\right )}
Nu kunnen we p oplossen met de areaalcosinus hyperbolicus (arcosh), maar voorzichtig. Omdat cosh een even functie is, krijgen we altijd een positief antwoord terug, maar er is ook een negatief antwoord, wat twee mogelijke waarden voor p geeft.
p =x_{middle}\pm a\cosh^{-1}\left ( \frac{s}{2a\sinh \left ( \frac{h}{2a}\right )} \right )
Welke waarde we nodig hebben, zien we aan de y-coördinaten van de eindpunten. Een proefje met een stuk touw leert dat bij gelijke touwlengte en gelijke x-coördinaten het laagste punt van de kromme (als het al tussen de x-coördinaten ligt) rechts van het midden ligt – en dus bij de positieve arcosh hoort – als de rechter-y-coördinaat lager is dan de linker, en omgekeerd. Nog een probleem kan optreden als de y-coördinaten gelijk zijn: dan hoort het laagste punt van de kromme precies in het midden te liggen, dus het argument van cosh hoort nul te zijn en de waarde waarvan je de arcosh neemt dus één. Door onnauwkeurigheden van zwevendekommagetallen komt die vaak iets lager dan één uit, waardoor de arcosh een fout geeft. Het is verstandig om vóór deze berekening te controleren of de y-coördinaten gelijk (of bijna gelijk) zijn, wat rekenwerk scheelt en deze fout voorkomt, en in dat geval gewoon aan te nemen dat p precies halverwege de eindpunten ligt. Nu a en p bekend zijn, kun je q berekenen door de x- en y-waarde van een van de eindpunten in te vullen in de vergelijking
y-q=a\cosh \left (\frac{x-p}{a} \right )
en die te herschrijven tot:
q=y-a\cosh \left (\frac{x-p}{a} \right )
Bereken daarmee q, en je hebt de uiteindelijke vergelijking waarmee je je kromme kunt tekenen.
y=a\cosh \left (\frac{x-p}{a} \right )+q
De kromme tekenen is nu eenvoudig. Bij de meeste stappen voor het berekenen van de parameters heb ik de hyperbolische functies zo nauwkeurig mogelijk berekend, omdat moeilijk te zeggen was hoe een onnauwkeurigheid zich zou voortplanten. Maar bij deze laatste stap, waarin het meeste rekenwerk zit, weet je hoe groot het argument is dat je gebruikt en hoeveel onnauwkeurigheid je kunt verdragen, en kun je voorspellen of een Taylorreeks tot een bepaalde macht nauwkeurig genoeg is – en meestal is dat volgens mijn ervaring de moeite waard.