Как нарисовать цепную линию
Когда мне понадобилось рисовать нити между воздушными шарами в Balloon Platform Defense, пришлось разобраться, как нарисовать нить, висящую между двумя точками, — форму, известную как цепная линия. Все примеры, которые я нашёл в интернете, либо упрощали задачу, например помещая обе точки на одну высоту, либо предполагали знание того, что здесь неизвестно, например угла, под которым начинается нить. Известны только начальная и конечная точки нити и её длина. В Википедии много сведений об уравнениях цепной линии, и с них я начал. Уравнение цепной линии:y=a\cosh \left (\frac{x}{a} \right )=\frac{a\left ( e^{\frac{x}{a}}+e^{-\frac{x}{a}} \right )}{2}
Здесь предполагается, что самая низкая точка нити лежит там, где она пересекает ось y. На практике к x и y добавляются константы, которые тоже нужно найти, ведь мы не знаем, где окажется самая низкая точка, — это одна из причин, почему задача сложнее любых примеров из учебника. (Вы замечали, что математические задачи из реальной жизни всегда гораздо сложнее примеров из учебника? Я всё ещё жду реальной задачи с интегралом, который берётся аналитически.) Но это подождёт, потому что сначала нужно найти значение a — константы в уравнении, которая, по сути, определяет, насколько узкой или широкой будет кривая. Осторожно: если точки находятся в одном и том же месте оси x, то a бесконечно. Этот случай нужно обработать до того, как вы дойдёте до этого места, — как и случай, когда точки так близки, что a становится невычислимо большим. (Помните, что «невычислимо большое» здесь — это любое число, которое в качестве аргумента показательной функции при используемой точности чисел с плавающей запятой дало бы NaN или бесконечность.) И конечно, ваш код вообще не должен доходить до этого места, если точки находятся дальше друг от друга, чем длина нити, — результат был бы бессмысленным.
Вычисление масштабного множителя
Для вычисления a Википедия даёт следующее уравнение (оно основано на свойстве цепной линии: если cosh задаёт положение кривой, то sinh задаёт её длину):\sqrt{s^{2}-v^{2}}=2a\sinh \left ( \frac{h}{2a} \right )
Здесь s — длина нити, а h и v — горизонтальное и вертикальное расстояния (по модулю) между начальной и конечной точками. Всё это известные величины, так что неизвестная остаётся одна, но найти её можно только численно. Сначала я пробовал делать это, не заботясь о начальном приближении, и сходился ли метод, было делом случая. К счастью, хорошее начальное приближение найти несложно. Ряд Тейлора для гиперболического синуса:
\sinh x=x+\frac{x^{3}}{3!}+\frac{x^{5}}{5!}+\cdots
После подстановок u = \frac{1}{4a^{2}} и c = \sqrt{s^{2}-v^{2}} решаемое уравнение принимает вид:
c=\frac{1}{\sqrt{u}}\sinh \left ( h\sqrt{u} \right )
Если взять первые три члена ряда Тейлора для sinh и упростить, получим:
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!}
и после перестановки:
\frac{h^{5}}{120}u^{2}+\frac{h^{3}}{6}u+\left ( h-c \right )=0
Это простое квадратное уравнение относительно u, и если подставить его коэффициенты в знакомую со школы формулу корней квадратного уравнения x = \frac{-b\pm \sqrt{b^{2}-4ac}}{2a}
, получится хорошее начальное приближение для u, а значит, и для 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}}, где a=\frac{1}{2\sqrt{u}}
Это приближение достаточно близко, чтобы найти решение методом Ньютона. В форме f\left ( a \right )=0
имеем f\left ( a \right )=2a \sinh \left ( \frac{h}{2a} \right )-c
и {f}'\left ( a \right )=2 \sinh \left ( \frac{h}{2a} \right )-\frac{h}{a}\cosh \left ( \frac{h}{2a} \right )
Поскольку 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 )}
Сейчас я проверяю, отличается ли очередной член последовательности от предыдущего меньше чем на 0,001. Обычно это происходит через 2–4 итерации, но иногда, когда x-координаты ближе друг к другу, требуется больше десяти. Методы Хаусхолдера более высокого порядка, вероятно, можно было бы применить без особых дополнительных затрат, потому что основная часть вычислений приходится на гиперболические функции, а дальнейшие производные функции должны давать лишь новые кратные этих же величин. Теперь, когда масштаб известен, можно вычислить смещение.
Вычисление смещения
Стандартное уравнениеy=a\cosh \left (\frac{x}{a} \right )
предполагает, что цепная линия симметрична относительно оси y и её самая низкая точка находится в (0,a), тогда как та, которую мы рисуем, может оказаться где угодно на экране. Значит, нужно найти смещение (p,q), которое поставит её на нужное место. Наше уравнение становится таким:
y-q=a\cosh \left (\frac{x-p}{a} \right )
(Напомню: a к этому моменту известно, x и y — переменные, так что неизвестными остаются только p и q.) Сначала я хотел подставить в уравнение известные значения x и y, то есть две конечные точки нити. Теоретически это работает, но получающиеся уравнения, несмотря на все попытки их упростить по ходу дела, часто давали промежуточные значения, с которыми не справлялись числа двойной точности. В итоге пришлось найти другой путь, который к тому же оказался эффективнее.
Вспомним: если в уравнении цепной линии заменить cosh на sinh, получится длина кривой. Подставив в эту форму x-координаты левой и правой конечных точек вместе с известной нужной длиной нити (s), получаем уравнение:
a\sinh \left (\frac{x_{right}-p}{a} \right )-a\sinh \left (\frac{x_{left}-p}{a} \right )=s
Одно из тождеств гиперболических функций гласит:
\sinh x - \sinh y=2\cosh\left ( \frac{x+y}{2} \right ) \sinh \left ( \frac{x-y}{2} \right )
Следовательно:
\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 )
В аргументе cosh правая и левая x-координаты складываются и делятся на 2, что даёт координату середины, так что их можно сразу ею заменить. Разность в аргументе sinh — та самая величина, которую мы раньше назвали h, поэтому и это значение sinh вычисляется как известная величина. Тогда исходное уравнение можно переписать так:
\cosh\left ( \frac{x_{middle}-p}{a} \right ) = \frac{s}{2a\sinh \left ( \frac{h}{2a}\right )}
Теперь можно найти p с помощью обратного гиперболического косинуса, но здесь нужна осторожность. Поскольку cosh — чётная функция, мы всегда получаем положительное решение, но есть и отрицательное, то есть у p два возможных значения.
p =x_{middle}\pm a\cosh^{-1}\left ( \frac{s}{2a\sinh \left ( \frac{h}{2a}\right )} \right )
Какое из них нужно, подскажут y-координаты конечных точек. Опыт с куском верёвки показывает: при той же длине нити и тех же x-координатах самая низкая точка кривой (если она вообще лежит между x-координатами) находится правее середины — то есть соответствует положительному значению обратного гиперболического косинуса, — когда правая y-координата ниже левой, и наоборот.
Ещё одна проблема может возникнуть, когда y-координаты равны: тогда самая низкая точка кривой должна быть ровно посередине, аргумент cosh должен быть равен нулю, а значение, от которого берётся обратный гиперболический косинус, — единице. Из-за неточностей вычислений с плавающей запятой часто получается чуть меньше единицы, и обратный гиперболический косинус не вычисляется. Поэтому разумно перед этим вычислением проверить, равны ли (или почти равны) y-координаты, — это экономит время и позволяет избежать ошибки, — и в таком случае просто считать, что p находится ровно посередине между конечными точками.
Теперь, когда a и p известны, можно вычислить q, подставив значения x и y одной из конечных точек в уравнение
y-q=a\cosh \left (\frac{x-p}{a} \right )
и переписав его так:
q=y-a\cosh \left (\frac{x-p}{a} \right )
Вычислите отсюда q — и у вас есть окончательное уравнение, по которому можно нарисовать кривую.
y=a\cosh \left (\frac{x-p}{a} \right )+q
Нарисовать кривую теперь просто. На большинстве шагов вычисления параметров я вычислял гиперболические функции как можно точнее, потому что трудно было предсказать, как будут распространяться погрешности. Но на этом последнем шаге, на который приходится большая часть вычислений, вы знаете, насколько велик используемый аргумент и какую погрешность можете себе позволить, и можете предсказать, будет ли ряд Тейлора до определённой степени достаточно точным, — и, по моему опыту, обычно он того стоит.