Рассмотрим два тела c массами m₁ и m₂, взаимодействующих по закону Ньютона

(3.5.1)

(3.5.2)

Эти уравнения описывают компоненты движения:

(a) Примените эти уравнения движения, учитывая, что и т. д., получите 12 ДУ первого порядка этой системы и запишите их.

(b) Решите их в Maple, используя G = 1, m₁ = 1, m₂ = 2 и начальные условия:

x1(0) = 1	x2(0) = –1
y1(0) = 0.5	y2(0) = –0.3
z1(0) = –0.3	z2(0) = 0.6
vx1(0) = 0.65	vx2(0) = –0.45
vy1(0) = 0.2	vy2(0) = 0.3
vz1(0) = 0.1	vz2(0) = –0.3

Запустите эти решения для интервала от t = 0 до t = 50. Получив массивы решений, интерполируйте их на новые массивы, равномерно расположенные во времени (xle, yle, zle, x2e, y2e, ... c N = 5 * длина (t)). Анимируйте движение двух масс. Для этого используйте интерполированные массивы таким образом, чтобы увидеть, что при сближении друг с другом тела ускоряются. Рисуйте орбиты кусками по пяти точкам данных. По одной точке можно построить орбиту как последовательность точек, а по большему числу точек орбита будет проявляться рывками.

Ваша картинка будет похожа на рис. 3.5.1, на котором тела совершают сложный танец.

Рис. 3.3.3. Траектории двух тел, притягивающихся по закону обратных квадратов

Для вычислений с трехмерными векторами удобнее перевести отдельные x-, y-, z-массивы в матрицы. Например, вектор r₁ будет матрицей с тремя столбцами и с числом строк, равным числу шагов по времени в выводе из процедуры решения ДУ (ode45). Трехмерные векторы задаются так (операторы придется переписать на языке Maple):

r1 = [x1, y1, z1] ;

r2 = [x2, y2, z2] ;

v1 = [vx1, vy1, vz1] ;

v2 = [vx2, vy2, vz2] ;

Вот версии этих векторов для данных, распределенных с равномерным шагом по времени:

r1e = [x1e, y1e, z1e] ;

r2e = [x2e, y2e, z2e] ;

v1e = [vx1e, vy1e, vz1e] ;

v2e = [vx2e, vy2e, vz2e] ;

Совет: посмотрите возможности пакета Physics в [2].

В физике обычно ищется простейшее описание движения, поэтому r1 и r2 заменяются на относительное положение центра масс:

(3.5.3)

(a) Для рассмотренного выше случая создайте 3D-графики R и и покажите, что это очень простое движение.

С помощью определенных выше матриц определите векторы R, r, V и v-подобный (операторы надо переписать на языке Maple):

R=(ml*rl+m2*r2)/(ml+m2);

V=(ml*vl+m2*v2)/(ml+m2);

r=rl–r2;

v=vl–v2;

(равномерное распределение во времени)

Re=(ml*rle+m2*r2e)/(ml+m2);

Ve=(ml*vle+m2*v2e)/(ml+m2);

re=rle–r2e;

ve=vle–v2e;

Постройте трехмерные графики.

(b) Создайте 3D-график вектора разности r и посмотрите на изменения этого вектора (в какой плоскости находится, как выглядит траектория... – эллипс?).

(c) Вычислите угловой момент центра масс (это объяснит принадлежность разностного вектора плоскости)

(3.5.4)

И численно покажите, что этот вектор постоянен во времени. Вектор L также вычисляется как матрица (операторы надо переписать на языке Maple; команда cross означает векторное произведение):

L=ml*cross(rl-R,vl-V)+m2*cross(r2-R,v2-V); (строки – время, столбцы – x-, y-, z-компоненты)

Постройте каждую компоненту в зависимости от времени: plot(t,L(:,1),t,L(:,2),t,L(:,3))

(d) Покажите графически, что L перпендикулярно r1 – R и r2 – R (и, следовательно, r1 – r2). Для вычисления скалярных произведений примените Maple (посмотрите возможности пакета Physics в [2]). Осторожно применяйте скалярное произведение к матрицам (у вас строки – это время, столбцы – это x-, y-, z-компоненты). Поэтому надо сказать Maple: нужны скалярные произведения координатных компонент; надо работать с индексами столбцов, например: dot1=dot(r1–R,L,2); dot2=dot(r2–R,L,2)).

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

Применим Maple к законам Кеплера.

Первый закон Кеплера. Разностный вектор r описывает эллипс с фокусом в положении второго тела (по определению r – это вектор, соединяющий точки нахождения масс: от m₂ к m₁, поэтому r = 0 в точке m₂).

Второй закон Кеплера. Воображаемая линия, проведенная от m₂ к m₁ (вектор r), описывает равные площади в равные интервалы времени.

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

(3.5.5)

где μ – приведенная масса:

(3.5.6)

Е – полная энергия центра масс:

(3.5.7)

(напомним: r = r1 – r2 и v = v1 – v2).

Вначале смоделируем Солнечную систему по Кеплеру. Зададим m₂ = 1 и m₁ = 1e^–8. Кроме того, можно избежать многих приключений, заменив начальные условия так, чтобы движение происходило только в плоскости xy и эллипс был бы ориентирован нужным образом. Для этого задайте все z-компоненты координат и скоростей равными 0. Затем сделайте Солнце неподвижным, задав для m₂ все координаты и скорости равными 0. В конце задайте для m₁ начальные условия такими, чтобы удобно сориентировать эллиптическую орбиту в плоскости xy (рис. 3.5.2).

Рис. 3.5.2. Прецессия орбит двух тел, притягивающихся по квадратичному закону

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

(a) Проверка первого закона Кеплера.

Напомним: расстояние от центра эллипса до его фокусов , где a – большая полуось, b – малая полуось. Для доказательства того, что это эллипс, вначале найдите правильные начальные условия (орбита должна быть в плоскости xy и, соответственно, ориентирована по осям координат). Параметрическое уравнение эллипса x = a cos s; y = b sin s; s ∈ [0, 2π].