Брахистохрона под землей

Мотивация

В сборниках задач по физике есть известная задача о свободном падении тела сквозь прямой тоннель, проходящий через центр Земли и соединяющий диаметрально противоположные точки. Для решения нужно знать, что сила притяжения к центру Земли в таком тоннеле уменьшается линейно до нуля, и тело, упавшее в такой тоннель, будет совершать гармонические колебания. Время движения составляет половину периода и для модели однородной Земли составляет 42 минуты. Более того, это время не меняется для любого прямого тоннеля, идущего по хорде (не проходящего через центр), если трением можно пренебречь.

Также в популярных книгах и других источниках есть задача о брахистохроне — кривой скорейшего спуска. В ней требуется найти форму желоба, проложенного между двумя фиксированными точками A и B, чтобы время скатывания тела из A в B было наименьшим. В однородном гравитационном поле брахистохроной является дуга циклоиды.

Вполне естественным кажется желание скрестить эти две задачи и выяснить, какой формы нужно проложить тоннель под землей, чтобы в идеальной ситуации без трения между точками A и B на поверхности (скажем, между двумя городами) поезд проезжал за наименьшее время. Эта новая задача имеет аналитическое решение. Её ответом является гипоциклоида — кривая, которую описывают точки окружности радиуса $$a$$, катящейся без проскальзывания внутри другой окружности радиуса $$b$$. Вот как записывается параметрическое уравнение гипоциклоиды и выглядит ее чертеж:

$$\left\{\begin{array}{lr} x=(b-a)\cos \theta+a\cos(b-a)\theta,\\ y=(b-a)\sin\theta-a\sin(b-a)\theta. \end{array}\right.$$

$$\begin{tikzpicture} \def\clr{black} \def\a{1} \def\b{5} \newcommand{\xt}[1]{(\b-\a)*cos(#1)+\a*cos((\b-\a)*#1/\a} \newcommand{\yt}[1]{(\b-\a)*sin(#1)-\a*sin((\b-\a)*#1/\a} \draw[cyan,very thin] (-\b-1,-\b-1) grid (\b+1,\b+1); \draw[->] (-\b-1,0) -- (\b+1,0); \draw[->] (0,-\b-1) -- (0,\b+1); \draw[\clr,thick] (0,0) circle (\b); \draw[\clr,dashed,very thin] (0,0) circle (\b-\a); \draw[line width=2pt,orange!80!red] plot[samples=60,domain=0:\a*360,smooth,variable=\t] ({\xt{\t}},{\yt{\t}}); \def\t0{40} \draw[\clr!50!black] (\t0:\b-\a) circle (\a); \draw[purple,fill] ({\xt{\t0}},{\yt{\t0}}) circle (2pt) -- (\t0:\b-\a) circle (1pt) node [midway, sloped, above] {\scriptsize $a$} -- (0,0) circle (1pt) node [midway, sloped, above] {\scriptsize $b-a$} node [below right] {$O$}; \node[\clr,fill=white,text width=1cm,align=center] at (-\b,\b) {$a=\a$ $b=\b$}; \end{tikzpicture}$$

Я задумался об этой задаче в старших классах или на первом курсе института, но тогда, разумеется, у меня не хватало знаний для решения в аналитическом виде. Я вернулся к задаче на четвертом курсе и смог довести вычисления до конца, не обращаясь к источникам, но и не без помощи Maple. Эти вычисления были одними из самых сложных среди всех задач с явным аналитическим ответом, которые я решал.

Физическая сторона задачи

$$\begin{tikzpicture} \def\a{1} \def\b{5} \draw ([shift=(-10:\b)]0) arc (-10:90:\b); \draw[->] (\b,0) node[right]{$A$} -- node[pos=0.92,above,inner sep=2] {$\varphi$} (0,0) -- (\b*0.35,\b*0.6) node[below,pos=0.9]{$\vec{r}$} ; \node[above] at (\b*0.32,\b*0.95) {$B$}; \draw[line width=1pt,orange!80!red] plot[samples=12,domain=0:72,smooth,variable=\t] ({(\b-\a)*cos(\t)+\a*cos((\b-\a)*\t/\a},{(\b-\a)*sin(\t)-\a*sin((\b-\a)*\t/\a}); \end{tikzpicture}$$

Движение поезда по подземному тоннелю будем рассматривать в полярных координатах. Считаем потенциал внутри Земли осцилляторным: $$U\sim r^2$$. Запишем закон сохранения энергии с учетом нулевой скорости поезда на поверхности:

$${mv^2\over 2}+{m\omega^2r^2\over 2}={m\omega^2r_0^2\over 2},$$(1)

где $$r_0$$ — радиус Земли, $$\omega=2\pi/T$$ — круговая частота колебаний, соответствующая периоду в 84 минуты. Время движения через тоннель $$\inline t=\int{dl\over v}$$. Выражаем скорость $$\inline v=\omega\sqrt{r_0^2-r^2}$$ из (1) и для времени движения поезда через тоннель получаем следующее выражение:

$$t=\int{dl\over v}={1\over\omega}\int{\sqrt{dr^2+r^2d\varphi^2}\over\sqrt{r_0^2-r^2}}={1\over\omega}\int{\sqrt{r'^2+r^2}\over\sqrt{r_0^2-r^2}}d\varphi,$$

Здесь мы выразили элемент длины пути через полярные координаты. Таким образом, задача свелась к поиску функции $$r(\varphi)$$ с фиксированным значением $$r(\varphi_A)=r(\varphi_B)=r_0$$ на концах отрезка, для которой функционал

$$\omega t=\int\limits_{\varphi_A}^{\varphi_B}\sqrt{{r'^2+r^2}\over{r_0^2-r^2}}d\varphi$$(2)

принимает наименьшее значение.

Математическая сторона задачи

Задача по минимизации интеграла (2) на множестве достаточно хороших с точки зрения физика функций (скажем, дифференцируемых сколько нужно раз) — типовая задача вариационного исчисления. Заметим, что подынтегральное выражение в (2) можно переписать в виде функции:

$${\cal L}(r,r')=\sqrt{{r' ^2+r^2}\over{r_0^2-r^2}}$$(3)

Из курса дифференциальных уравнений известно, что экстремальные значения (2) нужно искать на решениях дифференциального уравнения Эйлера — Лагранжа:

$${\partial{\cal L}\over\partial r}-{d\over d\varphi}{\partial{\cal L}\over\partial r'}=0.$$

Если вы вооружитесь бумагой и ручкой и попробуете подставить сюда (3), то быстро заметите сложность получившегося дифференциального уравнения со второй производной и нагромождением корней и осознаете тупиковость этого подхода. В системах компьютерной алгебры есть встроенные пакеты для решения вариационных задач. Вот пример для Maple, где вы видите трехэтажное дифференциальное уравнение:

Кроме того, Maple заботливо подсказывает, что у получившегося дифференциального уравнения есть первый интеграл. И тут вы вспоминаете, что на курсах теоретической механики и дифференциальных уравнений всё-таки рассказывают некоторые полезные вещи. В частности, у уравнений Эйлера — Лагранжа есть интегралы движения. В нашем случае интеграл движения имеет вид

$${\cal L}-r'{\partial{\cal L}\over\partial r'}=const.$$

Если вам проще воспринимать физику, а не математику, вы можете представить некоторую систему, описываемую функцией Лагранжа (3). Движение во «времени» $$\varphi$$ от «момента» $$\varphi_A$$ к $$\varphi_B$$ будет происходить по такой траектории $$r(\varphi)$$, на которой «действие» (2) примет экстремальное значение. Так как функция Лагранжа $$\cal L$$ не зависит от «времени» $$\varphi$$, то у системы есть сохраняющаяся величина — «энергия». Первый интеграл выше как раз и дает величину этой «энергии».

Выполним дифференцирование:

$$k={\cal L}-r'{\partial{\cal L}\over\partial r'}=\sqrt{{r' ^2+r^2}\over{r_0^2-r^2}} -r'{r'\over\sqrt{r_0^2-r^2}\sqrt{r'^2+r^2}}={r^2\over\sqrt{r_0^2-r^2}\sqrt{r'^2+r^2}}.$$

Мы пришли к дифференциальному уравнению первого порядка $$\inline r^2=k\sqrt{r'^2+r^2}\sqrt{r_0^2-r^2}$$, в котором разделяются переменные и получается следующий интеграл:

$$\varphi=\int{dr\over\sqrt{\cfrac{r^4}{k^2(r_0^2-r^2)}-r^2}}.$$

Его можно вычислить аналитически. Например, Maple выдает такой результат:

Чтобы доказать, что это выражение описывает гипоциклоиду, нужно взять её параметрическое уравнение и проделать достаточно громоздкие преобразования обратных тригонометрических функций, которые я не буду сейчас приводить.

В интернете опять кто-то неправ, или учимся избегать ошибок

Когда я занимался этой задачей, мне попадалась статья с нестандартными вычислениями: вместо полярных координат использовались декартовы. Тогда мне показалось, что автору очень повезло: уравнение гипоциклоиды в параметрической форме написано в таком виде, что параметр оказался пропорциональным физическому времени. При подготовке этой заметки я опять наткнулся на ту же статью. Сейчас эти нестандартные вычисления кажутся мне ошибочными.

Я предположил, что это чья-то студенческая работа с не очень высокими требованиями. Решил подробнее узнать об авторе, Аманде Максхэм, и нашел заметку в ее блоге как раз об этой статье. Она пишет, что эта статья выдается на первой странице в гугле по запросу «gravity train». Раз так, то тем более нужно разобраться, скрывается ли за нестандартными вычислениями какой-то хитрый прием, или это просто ошибка. Причем моя цель не в том, чтобы показать, что в интернете опять кто-то неправ, а в том, чтобы научиться отличать правильный способ рассуждения от неправильного.

Вот ключевой момент статьи, в котором делается попытка доказать, что гипоциклоида — искомое решение:

Вопросы начинаются с формулы (19), в которой пропущен дифференциал, и поэтому не совсем понятно, по какому параметру дальше берутся производные. Эти же производные сохраняются вплоть до (22), куда подставляются $$x(t)$$ и $$y(t)$$ из (24) и (25). Но параметр $$t$$ в (24) и (25) не может быть временем $$t_\text{физ}$$, так как он является аргументом тригонометрических функций и должен быть безразмерным. Получается, это некоторый параметр, задающий положение точки на траектории, но не являющийся временем.

(Попутно возникает мини-вопрос о том, откуда взялась разность $$x'^2-y'^2$$ вместо соответствующей суммы в правых частях (20) и (21), но эта опечатка в четырех местах не распространяется на дальнейшие рассуждения, и мы к ней придираться не будем.)

Еще возникает вопрос о степенях свободы. В уравнениях (19) и (22) неизвестные величины — $$(x, y, t)$$. Таким образом, эти уравнения содержат лишнюю степень свободы по сравнению с неизвестными $$(r, \varphi)$$ из нашего рассмотрения (и из предыдущего раздела работы Аманды). Эта избыточность мешает корректному применению формализма Лагранжа — Эйлера. Какие признаки этого мы видим?

1. (22) — не единственное возможное следствие из (20) и (21). Если поделить (20) на (21), получим, что $$(y'/x')^2=(dy/dx)^2=const$$, а это есть уравнение, описывающее не гипоциклоиду, а прямую траекторию. Почему это следствие из (20) и (21) менее правильное, чем (22)?

2. (22) с точностью до константы $$C$$ по форме совпадает с законом сохранения энергии, с которого начиналось решение. Действительно, выражение в левой части напоминает квадрат скорости, а в правой — разницу квадратичных потенциалов, то есть изменение потенциальной энергии. Присутствие константы $$C$$ объясняется тем, что производные в (22) берутся не по физическому времени $$t_\text{физ}$$, а по другому параметру $$t$$, который в $$C$$ раз меньше: $$t_\text{физ}=Ct$$. Если бы связь $$t_\text{физ}=Ct$$ нарушилась, то (22) противоречило бы закону сохранения энергии. Таким образом, (22) можно получить без применения лагранжева формализма, просто напрямую из закона сохранения энергии с помощью замены $$t_\text{физ}=Ct$$. Следовательно, (22) не может иметь большего физического содержания, чем закон сохранения энергии, и не несет информацию о траектории, минимизирующей время.

3. (24) и (25) формально являются решением (22), так как зануляют его при выполнении (26). Но они не являются единственным решением. Фактически, любое физически возможное движение из-за сохранения энергии будет удовлетворять уравнению (22). Важно лишь выбрать параметризацию этого движения так, чтобы физическое время было пропорционально параметру. Возьмем для примера движение по хорде $$x(t)=A$$, $$y(t)=B\cos t$$, где $$\inline A^2+B^2=R^2$$, подставим в (22) и при должном выборе константы $$C$$ (то есть при выборе правильного движения во времени вдоль хорды) получим тождество:

$$RB^2\sin^2 t=gC^2(R^2-A^2-B^2\cos^2 t)=gC^2B^2\sin^2 t.$$

Таким образом, осцилляторное движение по тоннелю в виде прямой хорды также является решением (22). Но оно отнюдь не минимизирует функционал (19) для общего времени движения между двумя точками.

Я так и не смог понять, можно ли провести рассуждения в декартовых координатах с избыточной степенью свободы и получить правильный ответ. Понятно желание использовать декартовы координаты, а не полярные, так как в них уравнение гипоциклоиды выглядит проще. Но интуитивно кажется, что лишняя степень свободы, появившаяся в функционале времени и соответствующая выбору закона движения тела вдоль траектории, не дает при применении лагранжева формализма достаточных ограничений на траекторию.

Я оставил комментарий в блоге Аманды с вопросами. Но никакой реакции не последовало.

Напишите свой комментарий, если видите ошибку в моих рассуждениях, или если знаете, как сделать корректное рассмотрение в декартовых координатах.

Задача о взвешенном выборе и случайной величине: единственность решения

В прошлый раз мы решали задачу о взвешенной сортировке и показали, что существует функция $$f_w(x)$$, такая что максимальное значение этой функции $$\inline\max_{i=1,2,...n}\left\{f_{w_i}(x_i)\right\}$$ будет достигнуто на k-той паре $$(w_k, x_k)$$ с вероятностью, пропорциональной $$w_k$$, где $$x_i$$ — значение случайной величины, равномерно распределенной на единичном отрезке $$[0, 1]$$.

Мы потребовали непрерывности и монотонности функции $$f_w(x)$$ по $$x$$, а также зафиксировали значения на концах отрезка: $$f_w(0) = a, f_w(1) = b$$. В этих предположениях показали, что функция должна удовлетворять условию

$$ \int\limits_0^1dx_1\,f^{-1}_{w_2}\left(f_{w_1}(x_1)\right)\cdot f^{-1}_{w_3}\left(f_{w_1}(x_1)\right)\cdot\ldots\cdot f^{-1}_{w_n}\left(f_{w_1}(x_1)\right) ={w_1\over w_1+w_2+w_3+\ldots+w_n} $$(1)

для любых наборов значений $$w_i>0$$. Мы нашли целый класс подходящих под это условие функций:

$$f_w(x)=h(x^{1/w}),$$(2)

где $$h(y)$$ — любая строго монотонная функция.

Теперь докажем, что других решений у этой задачи не существует.

Теорема о единственности. Пусть для каждого значения $$w>0$$ функция $$f_w(x)$$ непрерывна и строго монотонна по $$x$$ на отрезке $$[0,1]$$, принимает на концах отрезка фиксированные значения $$f_w(0) = a, f_w(1) = b$$ и удовлетворяет (1). Тогда найдется такая строго монотонная функция $$h(y)$$, для которой выполняется тождество (2).

Доказательство будет конструктивным: мы построим конкретный пример $$h(y)$$ на основе вида функции $$f_w(x)$$.

Лемма 1. Обратная функция $$f^{-1}_w(s)$$ при каждом фиксированном значении аргумента $$s\in[a,b]$$ непрерывна по параметру $$w$$ на всем допустимом множестве значений этого параметра $$w>0$$.

Доказательство от противного. Пусть существует точка $$w_2$$, в которой $$f^{-1}_w(s)$$ не является непрерывной. Это значит, что

$$\exists\varepsilon>0\ \forall\delta_1>0\ \exists s_1\in[a,b]\ \exists w_3:|w_2-w_3|<\delta_1,\left|f^{-1}_{w_2}(s_1)-f^{-1}_{w_3}(s_1)\right|>\varepsilon.$$(3)

Так как функция $$f^{-1}_w(s)$$ непрерывна по $$s$$ на отрезке $$[a,b]$$, она равномерно непрерывна на нем в силу теоремы Кантора — Гейне. Из равномерной непрерывности следует, что $$|f^{-1}_{w_2}(s)-f^{-1}_{w_3}(s)|>\varepsilon/2$$ в некоторой $$\delta_2$$-окрестности точки $$s_1$$, причем $$\delta_2$$ зависит только от $$\varepsilon$$, не от $$s_1$$. (Иными словами, разрыв по $$w$$ будет наблюдаться не только в одной точке $$s_1$$, но и в ее окрестности.)

Аналогично функция $$f_w(x)$$ непрерывна и равномерно непрерывна на отрезке $$[0,1]$$. По $$\delta_2$$-окрестности точки $$s_1$$ можно найти $$\delta_3$$-окрестность точки $$x_1=f^{-1}_{w_1}(s_1)$$, для каждого $$x$$ из которой $$\left|f^{-1}_{w_2}(f_{w_1}(x))-f^{-1}_{w_3}(f_{w_1}(x))\right|>\varepsilon/2$$ (значение $$w_1$$ выбирается произвольно и фиксируется для следующих рассуждений).

Это значит, что следующий интеграл от квадрата разности функций в соседних точках $$w_2$$ и $$w_3$$ ограничен снизу ненулевой величиной, не зависящей от $$\delta_1$$:

$$\int\limits_a^bdx\left[f^{-1}_{w_2}(f_{w_1}(x))- f^{-1}_{w_3}(f_{w_1}(x))\right]^2>{\delta_3\varepsilon^2\over 4}$$

Теперь непосредственно вычислим интеграл, раскрыв скобки и воспользовавшись условием (1):

$$ \begin{align*} &\int\limits_a^bdx\left[f^{-1}_{w_2}(f_{w_1}(x))- f^{-1}_{w_3}(f_{w_1}(x))\right]^2=\\ =&\!\int\limits_a^b\!dx\left[f^{-1}_{w_2}(f_{w_1}(x))\right]^2-2\!\int\limits_a^b\!dx\,f^{-1}_{w_2}(f_{w_1}(x))\, f^{-1}_{w_3}(f_{w_1}(x))+\!\int\limits_a^b\!dx\left[ f^{-1}_{w_3}(f_{w_1}(x))\right]^2=\\ =&{w_1\over w_1+2w_2}-2{w_1\over w_1+w_2+w_3}+{w_1\over w_1+2w_3}=\\ =&{2w_1\over(w_1+2w_2)(w_1+2w_3)(w_1+w_2+w_3)}(w_2-w_3)^2=\\ =&{2w_1\over(w_1+2w_2)(w_1+2w_2-2(w_2-w_3))(w_1+2w_2-(w_2-w_3))}(w_2-w_3)^2. \end{align*} $$

Оценка этого выражения снизу дает неравенство:

$${2w_1(w_2-w_3)^2\over(w_1+2w_2)(w_1+2w_2-2(w_2-w_3))(w_1+2w_2-(w_2-w_3))}>{\delta_3\varepsilon^2\over 4}.$$

Но согласно отрицанию определения непрерывности (3) $$\delta_1$$ можно выбрать сколь угодно малым, при этом $$\varepsilon$$, $$\delta_3$$ и $$w_1$$ зависят только от выбора $$w_2$$ и остаются фиксированными. Так как $$|w_2-w_3|<\delta_1$$, выбором $$\delta_1$$ также можно сделать левую часть неравенства сколь угодно малой при фиксированной правой. Противоречие. $$\square$$

Лемма 2. Для любых значений $$k,m\in\mathbb{N}$$, $$w>0$$, $$s\in[a,b]$$ выполняется равенство

$$\left[f^{-1}_{w/k}(s)\right]^k=\left[f^{-1}_{w/m}(s)\right]^m.$$

Доказательство. Введем в (1) новую переменную $$s=f_{w_1}(x_1)$$:

$$ \int\limits_a^bds\ {\partial f^{-1}_{w_1}(s)\over\partial s}\ f^{-1}_{w_2}(s)\cdot f^{-1}_{w_3}(s)\cdot\ldots\cdot f^{-1}_{w_n}(s) ={w_1\over w_1+w_2+w_3+\ldots+w_n}. $$

Разделим набор весов $$w_2$$, $$w_3$$, … $$w_n$$ на две группы количеством $$k$$ и $$m$$, в каждой из которых веса совпадают и равны $$w/k$$ и $$w/m$$ соответственно. Тогда:

$$ \int\limits_a^bds\ {\partial f^{-1}_{w_1}(s)\over\partial s}\ \left[f^{-1}_{w/k}(s)\right]^k \left[f^{-1}_{w/m}(s)\right]^m ={w_1\over w_1+2w}\quad\forall k,m\in\mathbb{N}. $$

Последний интеграл можно интерпретировать как скалярное произведение функций $$[f^{-1}_{w/k}(s)]^k$$ и $$[f^{-1}_{w/m}(s)]^m$$ в пространстве с положительной в силу монотонности $$f^{-1}_{w_1}$$ метрикой $$\partial f^{-1}_{w_1}(s)/\partial s$$. По неравенству Коши — Буняковского скалярное произведение не превосходит произведения норм, и равенство достигается при коллинеарности сомножителей. Из последнего уравнения видно, что скалярное произведение совпадает с квадратом нормы каждого элемента. Поэтому в нашем случае имеет место равенство. Действительно, вычислим норму от квадрата разности функций:

$$ \begin{align*} &\int\limits_a^bds\ {\partial f^{-1}_{w_1}(s)\over\partial s}\ \left\{\left[f^{-1}_{w/k}(s)\right]^k- \left[f^{-1}_{w/m}(s)\right]^m\right\}^2=\\ =&\int\limits_a^bds\ {\partial f^{-1}_{w_1}(s)\over\partial s}\ \left\{\left[f^{-1}_{w/k}(s)\right]^{2k} -2\left[f^{-1}_{w/k}(s)\right]^k\left[f^{-1}_{w/m}(s)\right]^m +\left[f^{-1}_{w/m}(s)\right]^{2m}\right\}=0. \end{align*} $$

Если предполагать противное — непрерывные функции не совпадают хотя бы в одной точке — интеграл от квадрата их разности не может быть нулевым. Противоречие.

Несмотря на присутствие в выражениях производной $$\partial f^{-1}_{w_1}(s)/\partial s$$, дифференцируемость по $$s$$ в этом доказательстве не требуется. Замена переменной $$x_1$$ на $$s$$ была проведена для удобства и наглядности. Аналогичные рассуждения имеют место и без замены переменных. $$\square$$

Следствие 1. Для любых значений $$r\in\mathbb{Q}_+$$, $$w>0$$, $$s\in[a,b]$$ выполняется равенство

$$f^{-1}_w(s)=\left[f^{-1}_{w/r}(s)\right]^r.$$

Доказательство. Представим положительное рациональное число $$r=m/k$$ как отношение натуральных чисел. В утверждении леммы 2 выполним замену $$w\to wk$$ и возведем обе части в степень $$1/k$$:

$$f^{-1}_w(s)=\left[f^{-1}_{wk/m}(s)\right]^{m/k}=\left[f^{-1}_{w/r}(s)\right]^r.\ \square$$

Следствие 2. Для любых значений $$p\in\mathbb{R}_+$$, $$w>0$$, $$s\in[a,b]$$ выполняется равенство

$$f^{-1}_w(s)=\left[f^{-1}_{w/p}(s)\right]^p.$$

Доказательство. Рассмотрим выражение $$[f^{-1}_{w/p}(s)]^p$$ как функцию от $$p$$. По следствию 1 в рациональных точках $$p$$ ее значение равно $$f^{-1}_w(s)$$. По лемме 1 функция непрерывна по параметру. Поэтому ее значение в иррациональных точках $$p$$ также равно $$f^{-1}_w(s)$$. $$\square$$

Доказательство теоремы. Из следствия 2 для любых $$p>0$$ выполняется равенство $$f_w(x)=f_{w/p}(x^{1/p})$$. Пусть $$p=w$$. Тогда

$$f_w(x)=f_1(x^{1/w})=h(x^{1/w}),$$

где $$h(y)\equiv f_1(y)$$ — строго монотонная функция. $$\square$$

Задача о взвешенном выборе и случайной величине

Условие

Пусть заданы n положительных чисел $$w_1$$, $$w_2$$, … $$w_n$$. Для каждого из них выберем значение $$x_i$$ случайной величины, равномерно распределенной на единичном интервале (0, 1). Существует ли функция $$f_w(x)$$, такая что максимальное значение этой функции $$\inline\max_{i=1,2,...n}\left\{f_{w_i}(x_i)\right\}$$ достигается на k-той выбранной паре $$(w_k, x_k)$$ с вероятностью, пропорциональной $$w_k$$?

Мотивация

Задача возникла из следующей программистской проблемы. Пусть в базе данных есть таблица с колонкой весов w. Нужно выбрать из этой таблицы случайную строку с вероятностью, пропорциональной значению веса w.

Можно было бы провести одно испытание, выбрав значение случайной величины $$x$$, и взять среди сумм $$w_1_$$, $$w_1 + w_2$$, $$w_1 + w_2 + w_3$$, … первую, превосходящую $$x(w_1 + w_2 + ... + w_n)$$. Но такой перебор сделать сложнее в синтаксисе SQL. Проще отсортировать таблицу по некоторой функции от веса и случайного числа:

SELECT * FROM table ORDER BY f(w, random()) DESC LIMIT 1

Какую функцию нужно применить к весу и случайному числу, чтобы значение этой функции на некоторой строке оказалось наибольшим с вероятностью, пропорциональной весу?

Решение

Допустим, что функция $$f_w(x)$$ существует. Пусть в первом выборе случайная величина приняла значение $$x_1$$. Вычислим, с какой вероятностью значение $$f_{w_1}(x_1)$$ будет наибольшим. Нам нужно взять меру множества значений $$x_i$$, на котором выполняется система неравенств:

$$ \begin{equation*} \begin{cases} f_{w_1}(x_1) > f_{w_2}(x_2),\\ f_{w_1}(x_1) > f_{w_3}(x_3),\\ \cdots \end{cases} \end{equation*} $$(1)

Предполагаем, что $$f_w(x)$$ непрерывна и монотонно возрастает по $$x$$. Тогда существует обратная функция $$f^{-1}_w(y)$$. Чтобы упростить рассмотрение на границах, предположим, что $$f_w(x)$$ принимает одинаковые значения на концах отрезка для всех значений параметра:

$$ f_w(0) = a,\quad f_w(1) = b. $$(2)

Тогда решением системы (1) будет множество

$$0<x_i<f^{-1}_{w_i}\left(f_{w_1}(x_1)\right),\quad i>1.$$

Меру этого множества нужно проинтегрировать по всем значениям $$x_1$$, чтобы вычислить вероятность получения максимального значения в первом выборе. С другой стороны, эта вероятность нормирована на 1 и по условию пропорциональна $$w_1$$. Выводим следующее ограничение на функцию $$f_w(x)$$ для произвольных значений $$w_i$$:

$$ \int\limits_0^1dx_1\,f^{-1}_{w_2}\left(f_{w_1}(x_1)\right)\cdot f^{-1}_{w_3}\left(f_{w_1}(x_1)\right)\cdot\ldots\cdot f^{-1}_{w_n}\left(f_{w_1}(x_1)\right) ={w_1\over w_1+w_2+w_3+\ldots+w_n}. $$(3)

Воспользуемся некоторыми эвристическими соображениями, которые позволяют угадать вид нашей функции. В интеграле можно перейти к новой переменной $$s=f_{w_1}(x_1)$$. Тогда подынтегральное выражение есть произведение величин $$f^{-1}_{w_2}(s)\cdot f^{-1}_{w_3}(s)\cdot\ldots\cdot f^{-1}_{w_n}(s)$$, а соответствующие веса в правой части входят только в виде суммы $$w_2+w_3+\ldots+w_n$$. Логично предположить, что обратная функция обладает показательным свойством по параметру:

$$f^{-1}_{w_2}(s)\cdot f^{-1}_{w_3}(s)\cdot\ldots\cdot f^{-1}_{w_n}(s) =f^{-1}_{w_2+w_3+\ldots+w_n}(s).$$(4)

Если обратная функция $$f^{-1}_w(y)$$ еще и непрерывна по этому параметру $$w$$, она имеет следующий вид: $$x=f^{-1}_w(y)=(g(y))^w$$, где $$g(y)$$ — некоторая новая монотонно возрастающая непрерывная функция. Отсюда следует, что прямая функция $$y=f_w(x)=g^{-1}(x^{1/w})$$.

Подставим полученный вид функций в интеграл из (3):

$$ \begin{align*} P_1&=\int\limits_0^1dx_1\,\left(g\left(g^{-1}(x_1^{1/w_1})\right)\right)^{w_2}\cdot \left(g\left(g^{-1}(x_1^{1/w_1})\right)\right)^{w_3}\cdot\ldots=\\ &=\int\limits_0^1dx_1\,\left(x_1^{1/w_1}\right)^{w_2+w_3+\ldots}= {1\over 1 +{w_2+w_3+\ldots\over w_1}}x_1^{1 +{w_2+w_3+\ldots\over w_1}}\Biggr|_0^1= {w_1\over w_1+w_2+\ldots+w_n}. \end{align*} $$

Таким образом, подстановка $$f_w(x)=g^{-1}(x^{1/w})$$ обращает уравнение (3) в тождество.

Вывод

Функция $$f_w(x)=x^{1/w}$$, а также ее композиция с любой монотонной функцией подходят для выбора из вариантов с заданными весами с помощью описанного в условии способа.

Обсудим единственность этого решения. Помимо естественного требования непрерывности и монотонности, мы наложили дополнительные ограничения (2) и (4). В этих предположениях решение единственно. В следующий раз мы рассмотрим вопрос единственности подробнее, а именно, появляются ли другие решения, если не требовать выполнения (2) и (4), а также можно ли вывести (2) и (4) непосредственно из условия или из более сильных ограничений на искомую функцию вроде дифференцируемости.

Хоккейная задача — 2

Продолжим решать «хоккейную задачу». В прошлом посте мы рассмотрели вариант задачи с искусственной поддержкой постоянной скорости вращения кольца. Перейдем теперь к варианту, когда за счет трения замедляется не только поступательное движение, но и вращение.

Напомним, что динамика кольца описывается уравнениями

$$\begin{align*} {du\over dt}&=-\int\limits_0^{2\pi}{d\alpha\over 2\pi}\,{u-v\sin\alpha\over\sqrt{u^2+v^2-2uv\sin\alpha}},\\ {dv\over dt}&=-\int\limits_0^{2\pi}{d\alpha\over 2\pi}\,{v-u\sin\alpha\over\sqrt{u^2+v^2-2uv\sin\alpha}}, \end{align*}$$(1)

где $$v$$ — скорость поступательного движения, $$u=\omega R$$ — скорость вращательного движения кольца.

Система уравнений (1) симметрична относительно замены $$u$$ на $$v$$. Если в начальный момент $$u=v$$, то по соображениям симметрии это соотношение будет сохраняться между скоростями вплоть до остановки. Другой случай, который мы рассмотрим — это $$u>v$$. Противоположный случай $$u<v$$ будет вытекать из него простой заменой $$u$$ на $$v$$ из тех же соображений симметрии.

Вырожденный случай u = v

Оба уравнения системы принимают один и тот же вид

$${dv\over dt}=-\int\limits_0^{2\pi}{d\alpha\over 2\pi\sqrt2}\,\sqrt{1-\sin\alpha}={2\over\pi}.$$

Скорости совместного движения падают линейно со временем, но в $$\pi/2=1,\!57$$ раз медленнее отдельного вращения или отдельного поступательного движения.

Случай u > v

К сожалению, интегралы в системе (1) сводятся к эллиптическим, что не оставляет надежды решить систему аналитически. Остается численное решение.

Прямая попытка решить уравнение в Maple проваливается. Определенные интегралы заменяются нагромождением эллиптических интегралов, приводить которое здесь нет смысла. При построении графика мы видим сообщение об ошибке

Warning, cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

Чтобы упростить выражения, накладываем ограничение $$u(t)\ge v(t)\ge 0$$ на каждый из интегралов.

sys_ode := {
    diff(u(t), t) = -`assuming`(
        [int(
            (u(t)-v(t)*sin(alpha))/sqrt(u(t)*u(t)+v(t)*v(t)-2*u(t)*v(t)*sin(alpha)),
            alpha = 0 .. 2*Pi
        )],
        [u(t) >= v(t) and v(t) >= 0, u(t) >= 0]
    )/(2*Pi),
    diff(v(t), t) = -`assuming`(
        [int(
            (v(t)-u(t)*sin(alpha))/sqrt(u(t)*u(t)+v(t)*v(t)-2*u(t)*v(t)*sin(alpha)),
            alpha = 0 .. 2*Pi
        )],
        [u(t) >= v(t) and v(t) >= 0, u(t) >= 0]
    )/(2*Pi), u(0) = 10, v(0) = 1
};

Ограничения устраняют сообщение об ошибке. Maple надолго задумывается, но выводит пустой график. Чтобы понять причину, выведем решение уравнения в какой-то точке.

p := dsolve(sys_ode, type = numeric):
p(1);

$$[t=1.,\\u(t)=10.9976102709142-2.84005015829721 10^{-8}{\rm I},\\v(t)=1.04887416587408+2.92542849266058 10^{-8}{\rm I}]$$

Накопление ошибки округления приводит к появлению ненулевой мнимой части в искомых функциях, из-за чего график пустой. Чтобы избежать появления мнимой части, возьмем в явном виде действительную часть правых частей уравнений:

sys_ode := {
    diff(u(t), t) = -Re(`assuming`(
        [int(
            (u(t)-v(t)*sin(alpha))/sqrt(u(t)*u(t)+v(t)*v(t)-2*u(t)*v(t)*sin(alpha)),
            alpha = 0 .. 2*Pi
        )],
        [u(t) >= v(t) and v(t) >= 0, u(t) >= 0]
    ))/(2*Pi),
    diff(v(t), t) = -Re(`assuming`(
        [Re(int(
            (v(t)-u(t)*sin(alpha))/sqrt(u(t)*u(t)+v(t)*v(t)-2*u(t)*v(t)*sin(alpha)),
            alpha = 0 .. 2*Pi
        )],
        [u(t) >= v(t) and v(t) >= 0, u(t) >= 0]
    ))/(2*Pi), u(0) = 10, v(0) = 1
};

Наконец, Maple, долго думая, всё же рисует график.

plots[odeplot](p, [
    [t, (10-t), color=gray],
    [t, sqrt(1-(1/10)*t), color=gray],
    [t, u(t), color = blue],
    [t, v(t), color = red]
], 0 .. 11);

При этом мы видим предупреждение о точке остановки вычислений, которая, очевидно, совпадает с моментом окончания движения:

Warning, cannot evaluate the solution further right of 10.155665, probably a singularity

Для численного решения мы выбрали начальную поступательную скорость 1 и скорость вращения 10. Без вращения кольцо перемещается в течение 1 единицы времени на расстояние 0,5. Наличие вращения привело к тому, что поступательное движение сохранялось более 10 единиц времени.

Линеаризованное решение при u >> v

График $$u(t)$$ — синяя линия — практически совпадает с прямой. График $$v(t)$$ — красная линия — напоминает параболу. Такое поведение предсказывается линеаризованным решением в пределе $$u\gg v$$. Пренебрегая $$v$$ в первом уравнении системы (1) и раскладывая второе по степеням $$p=v/u$$, получаем

$$\begin{align*} {du\over dt}&=-1,\\ {dv\over dt}&=-{v\over 2u}. \end{align*}$$(2)

Решением этой системы являются функции $$u(t)=u_0-t$$, $$v(t)=v_0\sqrt{1-t/u_0}$$. Их графики изображены серыми линиями.

Как видим из численного решения и графика, приближение верно описывает характер движения, пока $$u$$ и $$v$$ не становятся сравнимыми. Вот график для случая, когда одна скорость в два раза больше другой:

Здесь отличие от линеаризованного решения заметно практически сразу.

Путь до остановки: линеаризованное приближение

Пройденный кольцом путь дается площадью криволинейной трапеции под красной линией. В линеаризованном приближении криволинейная трапеция ограничена параболой и, как легко видеть, имеет площадь $$l=2u_0v_0/3$$. Если бы кольцо не вращалось, оно прошло бы расстояние $$l_0=v_0^2/2$$. Таким образом, чтобы увеличить проходимый путь в 5 раз, нужно закрутить кольцо с угловой скоростью

$$\omega_0={u_0\over R}={15\over 4}{v_0\over R}.$$

Напомним, что для прохождения кольцом того же пути при постоянном поддержании вращения его скорость должна быть в три раза меньше: $$\omega_0=(5/4)\,{v_0/R}$$.

Путь до остановки: численное решение

Теперь вычислим проходимый путь из первоначальных уравнений, без линеаризации.

Maple почему-то не нашел точку остановки («сингулярности») для начального значения 15/4 и некорректно продолжил решение. Приближенное положение этой точки находим подбором.

dsol2p := dsolve(sys_ode, type = numeric, output = listprocedure):
v := eval(v(t), dsol2p):
v(4.04);

$$0.00664275658541608$$

v(4.05);

$$0.000276558860981139$$

v(4.06);

$$0.00946581059191108$$

int(v(t), t = 0 .. 4.05);

$$2.5535588257745343$$

Последний результат Maple вычислял больше 10 минут. Как видим, для начального значения 15/4 погрешность линеаризованного приближения составляет 2%.

Вывод

Мы решили «хоккейную задачу» в исходной упрощенной формулировке и без упрощений в линеаризованном приближении. С помощью системы компьютерной алгебры Maple убедились, что для условия этой задачи линеаризованное приближение дает ошибку на несколько процентов.

Хоккейная задача

В 2007 году на физтеховской олимпиаде по физике была такая задача:

Тонкое кольцо лежит на шероховатой горизонтальной поверхности. После толчка в направлении центра кольца оно перемещается на некоторое расстояние. Когда это кольцо раскрутили до некоторой угловой скорости (поддерживаемой постоянной за всё время движения), то при той же начальной скорости кольцо прошло в $$k$$ раз большее расстояние. Как было раскручено это кольцо? (попытайтесь в ответе найти нелинейную поправку). (В.С. Булыгин)

На примере этой задачи я хочу показать, как использовать системы компьютерной алгебры, в частности, Maple.

В условии поддержка постоянной угловой скорости вращения выглядит искусственно. Это требование упрощает задачу, чтобы ее можно было решить на олимпиаде. В этом посте рассмотрим такую формулировку, а в следующем — потерю через трение не только поступательной скорости, но и вращательной.

Физическая сторона задачи

Решение задачи разобрано на Элементах (там ее назвали «хоккейной задачей»). Физическая часть решения не требует выхода за рамки школьных знаний, на ней останавливаться не будем. А вот на математике остановимся подробнее. Начнем с системы уравнений, выведенной по ссылке выше:

$$\begin{align*} {du\over dt}&=-\mu g\int\limits_0^{2\pi}{d\alpha\over 2\pi}\,{u-v\sin\alpha\over\sqrt{u^2+v^2-2uv\sin\alpha}},\\ {dv\over dt}&=-\mu g\int\limits_0^{2\pi}{d\alpha\over 2\pi}\,{v-u\sin\alpha\over\sqrt{u^2+v^2-2uv\sin\alpha}}. \end{align*}$$(1)

Здесь $$v$$ — скорость поступательного движения, $$u=\omega R$$ — скорость вращательного движения кольца.

Уравнения таковы, что соответствующим выбором единицы измерения времени мы можем избавиться от несущественного множителя $$\mu g$$, поэтому в рассуждениях ниже мы его опустим.

Предельный режим

В предположении $$u=const$$ уравнение на $$v$$ имеет простой предельный режим, когда $$v\to 0$$. Тогда раскладывая подынтегральное выражение в ряд по $$v$$ с точностью до второго порядка малости и интегрируя, видим, что убывание скорости пропорционально самой скорости. При малых скоростях кольцо останавливается по экспоненте, как будто трение не сухое, а жидкостное.

Случай постоянной скорости вращения

По условию угловая скорость вращения поддерживается постоянной, а линейная скорость падает до 0. Логично предположить, что мы имеем дело с режимом $$0<v<u=\text{const}$$. Пусть $$p=v/u$$. Тогда

$$u{dp\over dt}&=-\int\limits_0^{2\pi}{d\alpha\over 2\pi}\,{p-\sin\alpha\over\sqrt{1+p^2-2p\sin\alpha}}.$$(2)

До остановки кольцо проходит расстояние

$$l=\int\limits_0^{t_\text{ост}} v\,dt=u\int\limits_0^{t_\text{ост}} p\,dt.$$(3)

Подставим $$dt$$ из (2) в (3):

$$u^2\int\limits_{p_0}^0{p\,dp\over-\int\limits_0^{2\pi}{d\alpha\over 2\pi}\,{p-\sin\alpha\over\sqrt{1+p^2-2p\sin\alpha}}}=u\int\limits_0^{t_\text{ост}} p\,dt=l.$$(4)

Когда вращения нет, кольцо движется равнозамедленно и проходит расстояние $$l_0=v_0^2/2$$. Вращающееся кольцо проходит в $$k=5$$ раз большее расстояние $$l$$:

$$l=kl_0=k\,{v_0^2\over 2}=ku^2\,{p_0^2\over 2}.$$

Подставим $$l$$ в (4) и поделим на $$p_0^2$$:

$${1\over p_0^2}\int\limits^{p_0}_0{p\,dp\over\int\limits_0^{2\pi}{d\alpha\over 2\pi}\,{p-\sin\alpha\over\sqrt{1+p^2-2p\sin\alpha}}}={k\over 2}.$$(5)

Это уравнение дает нам соотношение между начальной скоростью вращения кольца $$\omega_0=v_0/(p_0R)$$ и ростом проходимого до остановки расстояния $$k$$.

Численное решение в Maple

Внутренний интеграл выражается через эллиптические интегралы. Однако Maple ничего не знает про знак $$p$$ и выдает слишком громоздкое выражение, содержащее фрагменты вроде $$\sqrt{-{2p/(p-1)^2}}$$, csgn(p-1). Подскажем очевидное ограничение на $$p$$:

`assuming`([
    int(
        (p-sin(alpha))/(2*Pi*sqrt(p^2 + 1 - 2*p*sin(alpha))),
        alpha = 0 .. 2*Pi
    )
], [p > 0]):
factor(%);

$${\left(p-1\right)\mathrm{EllipticK}\left(\dfrac{2\sqrt{p}}{p+1}\right)+\left(p+1\right)\mathrm{EllipticE}\left(\dfrac{2\sqrt{p}}{p+1}\right)\over p\pi}$$

После этого легко подобрать ответ $$p_0$$ с нужной точностью:

p0 := 0.78:
(int(
    p^2 * Pi / (
        (p-1) * EllipticK(2*sqrt(p)/(p+1)) +
        (p+1) * EllipticE(2*sqrt(p)/(p+1))
    ),
    p = 0 .. p0
)) / p0^2;

$$2.491446599$$

В начальный момент скорости поступательного движения и вращения были связаны соотношением $$v_0\approx0,\!78\,\omega_0R$$.

Приближенный ответ можно получить не только численно, но и из разложения внутреннего интеграла в ряд

expand(series(`assuming`([int(...)], [p >= 0]), p = 0, 8));

$${1\over 2}p+{1\over 16}p^3+{3\over 128}p^5+O(p^7)$$

Вычислим левую часть (5) при $$p_0=0,\!78$$, последовательно уточняя разложение внутреннего интеграла:

$$\begin{align*} &{1\over 0,\!78^2}\int\limits_0^{0,78}{p\over{\frac{1}{2}p}}\,dp=2,\!564102564,\\ &{1\over 0,\!78^2}\int\limits_0^{0,78}{p\over{\frac{1}{2}p+{1\over 16}p^3}}\,dp=2,\!501916372,\\ &{1\over 0,\!78^2}\int\limits_0^{0,78}{p\over{\frac{1}{2}p+{1\over 16}p^3+{3\over 128}p^5}}\,dp=2,\!493976851. \end{align*}$$

Хотя начальное значение параметра $$p_0=v_0/u_0=0,\!78$$ нельзя назвать малым, линеаризация по этому параметру дает погрешность менее 3%. Для линейного приближения ответ составил бы $$p_0=4/k=0,\!8\,$$. А учет кубической поправки уже дает погрешность меньше процента.

туда →