\documentclass[11pt,a4paper,final,onecolumn]{article}
\usepackage[cp1251]{inputenc}
\usepackage[russian]{babel}
\usepackage{graphics}
\usepackage{amsmath}
\begin{document}
\section{Постановка задачи}
%---------------------------------------------------
Предположим, что тело (например спутник) находится в гравитационном поле
 массивного тела (планеты или звезды) на расстоянии $r_0$ от
его центра. В некоторый момент времени ему сообщили скорость $v_0$ под углом
 $\alpha$ к линии, соединяющей центр
массивного тела и вращающийся объект. Предположив, что сила тяготения
 определяется соотношением
$$F=\frac{k}{r^{n-1}},$$ необходимо определить траекторию движения.
\section{Решение}
%---------------------------------------------------
%---------------------------------------------------
\subsection{Допущения}
%---------------------------------------------------
Прежде всего необходимо отметить, что мы будем рассматривать скорости, которые
 значительно меньше скорости
света в вакууме. Кроме того, размеры вращающегося тела пренебрежимо малы, по
 сравнению с массивным телом. Иногда
вместо словосочетания `массивное тело` мы будем употреблять слово `планета`, а
 вместо `вращающееся тело` слово `спутник`.
\subsection{Вывод уравнения}
%---------------------------------------------------
Поскольку действует центральная сила, то, взяв начало координат в центре
 планеты, будем иметь
$$ \vec{r} \times \vec{F} =0,
$$
\center {или}
$$\frac{d}{dt}(\vec{r} \times m\vec{v})=\vec{r}\times \vec{F}.$$
\flushleft{Это уравнение дает решение}
\begin{equation}
\label{e1}
\vec{r}\times \vec{v}=\vec{\zeta}=\displaystyle{\overrightarrow{\textup{const.}}}
\end{equation}
%Рисунок с площадью
%--------------------------------------------
\begin{figure}[h!]
  \label{square}

  %\centerline{\includegraphics{/home/geckelberryfinn/field.ps}}
    \centerline{\includegraphics{/home/geckelberryfinn/figure1.ps}}

  \caption{}
\end{figure}
\flushleft
Рассмотрим теперь более подробно движение вокруг планеты.
Как видно из рисунка изменение площади, заметываемой радиус-вектором спутника,
 можно представить так:
$$\Delta (\vec{S})=\overrightarrow{OMM`}=\frac{1}{2}(\vec{r}\times\Delta\vec{r}).$$
Но тогда, деля обе части равенства на $dt$ и устремляя к пределу, получим
$$\frac{d}{dt}(\vec{S})=\lim_{dt \to 0}(\vec{r} \times \frac{d}{dt}\vec{r}) ,$$
\center{или}
$$2\vec{V_{S}}=\vec{r}\times\vec{v}.$$
\flushleft
Тогда из выражения \ref{e1}  и последнего равенства следует
\begin{equation}
\label{e2}
|r\times v|=2\frac{d}{dt}S=\zeta.
\end{equation}
%Рисунок с обоснованием перехода к полярным координатам
%--------------------------------------------
\begin{figure}[h!]
  \label{sqpolar}
 %\centerline{ \includegraphics{/home/geckelberryfinn/sqpolar.ps}}
 \centerline{ \includegraphics{/home/geckelberryfinn/figure2.ps}}
 \caption{}
\end{figure}
\flushleft
Перейдем к полярным координатам, исходя из геометрических соображений. Изменение
 площади можно приближенно
выразить как площадь прямоугольного треугольника $OM_1M_2$ (рис. 2), которая
 равна:
$$
dS=\frac{1}{2}r^2d\phi,
$$
\flushleft
откуда получаем
\begin{equation}\label{polar}
\frac{d}{dt}S=\frac{1}{2}r^2\frac{d}{dt}\phi.
\end{equation}
Из \ref{e2} и \ref{polar} видно, что
\begin{equation}
\label{sqr}
r^2\frac{d\phi}{dt}=\zeta.
\end{equation}
Далее представим квадрат скорости спутника следующим образом:
\begin{equation}\label{speeds}
v^2=(\frac{d}{dr}r)^2+r^2(\frac{d}{dt}\phi)^2=v_{r}^{2}+v_{p}^{2}.
\end{equation}
Преобразуем выражения для $v_r$ и $v_p$, используя равенство \ref{sqr}:
\begin{equation}
\label{vr}
v_r=\frac{dr}{dt}=\frac{dr}{d\phi}\frac{d\phi}{dt}=\frac{\zeta}{r^2}\frac{dr}{d\phi},
\end{equation}
\begin{equation}
v_p=r\frac{d\phi}{dt}=\frac{\zeta}{r}.
\end{equation}
Сделаем замену переменной
$$u=\frac{1}{r}.$$
Тогда выражение \ref{speeds} для скорости примет вид
\begin{equation} \label{speednew}
  v^2=\zeta^2\left[\left(\frac{du}{d\phi}\right	)^2+u^2\right].
\end{equation}
По теореме об изменении кинетической энергии, учитывая, что движение происходит
 в поле центральной силы, будем иметь
$$d\left(m\frac{v^2}{2}\right)=F_rdr.$$
Деля обе части этого равенства на $d\phi$, получим
$$\frac{m}{2}\frac{d(v^2)}{d\phi}=F_r\frac{dr}{d\phi}.$$
Используя формулу \ref{speednew} , запишем
$$\frac{m\zeta^2}{2}\frac{d}{d\phi}\left[\left(\frac{du}{d\phi}\right)^2+u^2\right]=F_r\frac{dr}{d\phi}.$$
Учитывая введенную замену, имеем
\begin{equation}\label{afterreplace}
\frac{m\zeta^2}{2}\left(2\frac{du}{d\phi}\frac{d^2u}{d\phi^2}+2u\frac{du}{d\phi}
\right)=-F_r\frac{1}{u^2}\frac{du}{d\phi}.
\end{equation}
Окончательно имеем (если $\frac{du}{d\phi}\ne 0$ - этот случай тривиален и будет
 рассмотрен отдельно)
\begin{equation}\label{maineq}
m\zeta^2u^2\left(\frac{d^2u}{d\phi^2}+u\right)=-F_r.
\end{equation}
Константу $\zeta$ найдем из начальных условий, используя выражение \ref{e1}:
$$
\zeta=|r \times v|=r_0 \cdot v_0 \cdot \sin{\alpha},
$$
где $r_0$ - расстояние от спутника до центра планеты,  $v_0$ - начальная
 скорость (напомним, что мы рассматриваем только
случаи, когда $v_0 \ll c, c$  - скорость света в вакууме,  а $\alpha$  - угол
между вектором скорости и радиус-вектором, соединяющим центр планеты и спутника.
 \\
Поскольку в нашем случае сила $F_r$ имеет вид
$$F_r=\frac{-km}{r^{n-1}},$$
то после разрешения уравнения \ref{maineq} относительно производной получим
\begin{equation}
\label{final}
\frac{d^2}{d\phi^2}u(\phi)=\frac{ku(\phi)^{n-3}}{(r_0v_0 \sin{\alpha}
 )^2}-u(\phi).
\end{equation}
\subsection{Решение уравнения}
%---------------------------------------------------
Введем новую переменную
$$
z=\frac{du}{d\phi}.
$$
Тогда уравнение \ref{final} преобразуется к виду
\begin{equation}\label{solve}
\frac{zdz}{du}=\frac{ku^{n-3}(\phi)}{(r_0v_0\sin{\alpha})^2}-u(\phi).
\end{equation}
Из этого уравнения видно, что при $n \ne 2$ решение \ref{solve} будет иметь
 вид:
$$
z^2=\frac{2ku^{n-2}(\phi)}{n-2}-u^2(\phi)+\sigma,
$$
где $\sigma$- постоянная интегрирования, которую можно найти из начальных условий. В
 самом деле, возвращаясь к исходным
обозначениям, получаем
\begin{equation}
\label{int2}
\frac{du}{d\phi}=\pm \sqrt{\frac{2ku^{n-2}(\phi)}{n-2}-u^2(\phi)+\sigma}.
\end{equation}
Но так как
$$
\frac{du}{d\phi}=-\frac{1}{r^2}\frac{dr}{d\phi},
$$
кроме того, из уравнения \ref{vr} следует, что
$$
\frac{dr}{d\phi}=\frac{r^2}{\zeta}\frac{dr}{dt}.
$$
В начальный момент времени последнее выражение имеет вид (учитывая полученное
 нами значения для $\zeta$)
$$
\frac{dr}{d\phi}=\frac{v_0\cos{\alpha}\cdot r_0^2}{v_0
 \sin{\alpha}}=r_0^2\ctg{\alpha}.
$$
Теперь можно записать выражение для $\sigma$
\begin{equation}\label{sigma}
\sigma=\frac{1}{r_0^2}-\ctg^2{\alpha}-\frac{2k}{r_0^n(n-2)(v_0\sin{\alpha})^2}.
\end{equation}
Обратимся теперь к выражению \ref{int2}. Его структура такова, что дальнейшее
 интегрирование в квадратурах возможно лишь в
случаях $n=3$ и $n=4$. Выпишем для каждого случая эти решения.
\subsubsection{Тривиальный случай $\frac{du}{d\phi}=0$}
Заметим, для начала, что в процессе вывода уравнения \ref{final} мы
 воспользовались неэквивалентным преобразованием, поделив уравнение
\ref{afterreplace} на $\frac{du}{d\phi}$. При этом мы могли потерять решение
 $\frac{du}{d\phi}=0$. Это возможно в случае, когда $r=\textup{const},$
т.е. траектория представляет собой окружность. Это тривиальный случай, нас же
 будет интересовать вариант, когда $\frac{du}{d\phi}\ne0$.
\subsubsection{Случай $n=3$}
%-------------------n=3-------------------------------
Перепишем выражения \ref{int2}  и \ref{sigma} следующим образом:
\begin{gather}\label{n3}
\frac{du}{d\phi}=\pm \sqrt{\frac{2ku(\phi)}{\zeta^2}-u^2(\phi)+\sigma_{n=3}},\\
\label{sigma3}
\sigma_{n=3}=\frac{1}{r_0^2}-\frac{2k}{r_0^3(v_0\sin{\alpha})^2}-\ctg^2{\alpha}.
\end{gather}
Уравнение \ref{n3} дает решение
$$
\phi=\arcsin{\frac{u(\phi)-\frac{k}{\zeta}}{\sqrt{\sigma_{n=3}+(\frac{k}{\zeta})^2}}}+\Phi,
$$
где
$$
\Phi=-\arcsin \frac{\frac{1}{r_0}-\frac{k}{\zeta}}{\sqrt{\sigma_{n=3}+(\frac{k}{\zeta})^2}}
$$
\center{или}
$$
r(\phi)=\frac{1}{\sqrt{\sigma_{n=3}+\left(\frac{k}{\zeta}\right)^2}\cdot\sin (\phi-\Phi)+\frac{k}{\zeta}}. $$
%--------------------n=4-----------------------------------
\flushleft
\subsubsection{Случай n=4}
При таком значении параметра $n$ выражения \ref{int2} и \ref{sigma} принимают
 вид
\begin{gather}\label{n4}
\frac{du}{d\phi}=\pm \sqrt{u^2(\phi)(k-1)+\sigma_{n=4}}, \\
\label{sigma4}
\sigma_{n=4}=\frac{1}{r_0^2}-\ctg^2{\alpha}-\frac{k}{r_0^4(v_0\sin{\alpha})^2}.
\end{gather}
Решение уравнения \ref{n4} выглядит следующим образом:
\begin{equation}
\phi=\frac{1}{\sqrt{k-1}}\ln{\Phi'\cdot\left(u(\phi)\sqrt{k-1}+\sqrt{u^2(\phi)(k-1)+\sigma_{n=4}}\right)},
\end{equation}
где
$$
\Phi'=\frac{1}{\left(\frac{1}{r_0}\sqrt{k-1}+\sqrt{\frac{1}{r_0^2}(k-1)+\sigma_{n=4}}\right)}
$$
%----------------------n=2---------------------------------
\subsubsection{Случай n=2}
Рассмотрим теперь двумерный случай. При этом все выражения, котрые мы получили
для $\sigma$ и для $\frac{du}{d\phi}$ не имеют смысла, поэтому вернемся к
уравнению \ref{solve}. При $n=2$ оно выглядит следующим образом:
\begin{equation} \label{n=2}
z\frac{dz}{du}=\frac{k}{u(\phi)(r_0v_0\sin \alpha)^2}-u.
\end{equation}
После интегрирования получим
\begin{equation}
\frac{du}{d\phi}=\pm \sqrt{\frac{2k\ln u(\phi)}{(r_0v_0\sin\alpha)^2}-u^2(\phi)+\sigma_{n=2}},
\end{equation}
где
$$
\sigma_{n=2}=\ctg^2 \alpha+\frac{2k\ln r_0}{(r_0 v_0 \sin \alpha)^2}+\frac{1}{r_0^2}.
$$
Дальнейшее интегрирование в квадратурах невозможно, поэтому мы проведем численное решение этого
и других неинтегрируемых в квадаратурах случаев.
\subsubsection{О методах, используемых при численном интегрировании}
Для численного решения дифференциального уравнения \ref{final} мы пользовались математическим пакетом MathCad 2001 Pro, при этом
воспользовавшись встроенной функцией \textup{odesolve}.
%----------------------------------------------------------
%----------------------Рисунки------------------------------
\begin{figure}[p]
\begin{tabular}[t]{|p{100mm}|p{40mm}|}
\hline
 вид решения &начальные условия \\ \hline
 \includegraphics{fig1_s.ps} & \vbox{$\alpha=\frac{\pi}{4}$,k=3, $r_0=2$, $v_0=0.31$, n=2 \vspace{30mm}}\\ \hline
 \includegraphics{fig2_s.ps} & \vbox{$\alpha=\frac{\pi}{2}$,k=3, $r_0=2$, $v_0=2 $, n=2\vspace{30mm}} \\ \hline
 \includegraphics{fig3_s.ps} & \vbox{$\alpha=\frac{\pi}{3}$,k=3, $r_0=2$, $v_0=2 $, n=3 \vspace{30mm}}\\ \hline
\end{tabular}
\end{figure}

\begin{figure}[p]
\begin{tabular}[t]{|p{100mm}|p{40mm}|}
\hline
 \includegraphics{fig4_s.ps} & \vbox{$\alpha=\frac{\pi}{2}$,k=3, $r_0=2$, $v_0=2 $, n=3 \vspace{30mm}} \\ \hline
 \includegraphics{fig5_s.ps} & \vbox{$\alpha=\frac{\pi}{2}$,k=3, $r_0=2$, $v_0=0 $, n=3\vspace{30mm}} \\ \hline
 \includegraphics{fig6_s.ps} & \vbox{$\alpha=\frac{\pi}{4}$,k=3, $r_0=2$, $v_0=0.31 $, n=4 \vspace{30mm}}\\ \hline
\end{tabular}
\end{figure}

\begin{figure}[p]
\begin{tabular}[t]{|p{100mm}|p{40mm}|}
\hline
 \includegraphics{fig7_s.ps} &  \vbox{$\alpha=\frac{\pi}{4}$,k=3, $r_0=2$, $v_0=1 $, n=4 \vspace{30mm} }  \\ \hline
 \includegraphics{fig_out.ps} & \vbox{$\alpha=\frac{\pi}{3}$, k=3, $r_0=2$, $v_0=0.971 $, n=4 \vspace{50mm}} \\ \hline
\end{tabular}
\end{figure}

\end{document}
