Задача Кеплера

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

Разделы:

1. Вводим систему масштабов, которая делает задачу безразмерной. Устанавливаются три шкалы:

AU (среднее расстояние Солнце–Земля) – мера расстояния;
масса Солнца M_s – мера массы;
год (время, необходимое для того, чтобы Земля обошла Солнце) – мера времени.

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

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

3. Кратко обсудим проблему рассеяния.

Задача о Солнечной системе

Параметры Солнечной системы:

солнечная масса – 2 · 1030 кг;
гравитационная постоянная – 6.67 · 10 СИ (Н · м2/кг2);
масса Земли – 6 · 1024 кг;
расстояние Солнце–Земля – 1.5 · 1011 м;
время обращения Земли – 3.15 · 107 с.

restart;
G:=6.67*10^(-11);

M_s:=1.99*10^(30);

M_e:=5.97*10^(24);

G*M_s*M_e;

AU:=1.5*10^(11);

G*M_s*M_e/AU^2;

Притяжение Земли к Солнцу довольно велико, если выразить его в ньютонах, но ускорение в м/с2 не очень велико:

G*M_s/AU^2;

Обратите внимание, что в центре Земли гравитационное ускорение уравновешивается центробежной силой из-за вращения вокруг Солнца. На поверхности Земли на объект действует только приливная сила.

Период задается так:

T_e:=evalf(365.26*24*60*60);

v_e:=evalf(2*Pi*AU)/T_e;

Земля движется вокруг Солнца со скоростью около 30 км/с. Введем систему единиц, в которой радиус орбиты Земли равен единице, масса измеряется кратно массе Солнца, а время измеряется в единицах периода орбиты. Выбираем три размерные величины и устанавливаем шкалы:

M = M ' · M_s;

T = T ' · T_e;

X = X ' · AU.

Перепишем уравнения с учетом этой замены. Аналогичный прием используется для масштабирования орбиты атома в квантовой механике, разве что в качестве масштаба массы используется масса электрона. В задаче Кеплера этого не требуется, потому что масса Солнца выпадает из уравнения движения. В квантовой задаче это не так, потому что сила притяжения электрическая, а не гравитационная, т. е. она не пропорциональна массе более легкой частицы (в кантовой задаче – электрона).

Найдем единицу скорости в задаче движения Земли вокруг Солнца, учитывая введенный выше масштаб:

unitLength:=AU; # introduced unit #1

unitTime:=T_e; # introduced unit #2

unitVel:=unitLength/unitTime; # a derived unit!

Скорости V оказываются кратны 4,7 км/с, что означает, что крейсерская скорость движения Земли вокруг Солнца равна 2π!

v_e/unitVel;

Чему равна единица силы? Естественной единицей должна быть сила, испытываемая Землей под действием силы гравитации. Единица ускорения – это, очевидно, единица скорости, деленная на единицу времени:

unitAcc:=unitVel/unitTime; # a derived unit!

unitMass:=M_s; # introduced unit #3 !

unitForce:=unitMass*unitAcc; # a derived unit!

Зная единицу для ускорения, запишем ускорение, которое Земля испытывает от гравитационного притяжения к Солнцу, в этом же масштабе:

(G*M_s/AU^2)/unitAcc;

Сила, действующая на Землю, в этом масштабе:

(G*M_s*M_e/AU^2)/unitForce;

Гравитационная постоянная запишется так:

unitG:=unitForce*(unitLength)^2/unitMass^2;

Запишем центробежную силу в системе СИ: в точке, где скорость и координата направлены под прямым углом друг к другу, величина углового момента равна произведению координаты и скорости на массу:

L0:=(M_e*AU*v_e);

L0/(unitMass*unitVel*unitLength); # angular momentum in our units:

subs(r=AU,-diff(L0^2/(2*M_e*r^2),r));

Гравитационное притяжение:

-G*M_s*M_e/AU^2;

Равновесие гравитационного притяжения и центробежной силы не идеально, а соответствует числу значащих цифр в расчете. Это означает, что орбита Земли – не идеальный круг, а величина AU усреднена по эллипсу с эксцентриситетом, близким к нулю. Эксцентриситет Земли равен 0.017, Марса – 0.092, Меркурия – 0.2.

Запишем условие равновесия. Центробежную силу вычислим дифференцированием центробежного потенциала, записанного как L2/(2·m·r2):

subs(r=AU,diff(L0^2/(2*M_e*r^2),r))*(unitMass*unitLength^3/
(unitMass*unitLength*unitVel)^2);

(-G*M_s*M_e/AU^2)*(unitLength^2/(unitMass^2*unitG));

G/unitG;

Это показывает масштаб величин G · M_s.

Стартуем заново и смотрим (в декартовой системе координат), как Земля (масса m) движется вокруг неподвижного гравитационного центра притяжения:

restart; GM:=39.17: m:=6*10^(24)/(2*10^(30)):

Зададим для потенциала показатель в степенной зависимости от радиуса (при желании здесь можно поиграть):

n:=10/10;

V:=-GM*m/r^n;

Из закона сохранения ориентации вектора углового момента следует, что движение происходит в плоскости, т. е. задача является двумерной, в координатах x и y:

V:=subs(r=sqrt(x^2+y^2),V);

Сила рассчитывается так:

Fx:=-diff(V,x);

Fy:=-diff(V,y);

Видно, что обе составляющие силы зависят как от x, так и от y. Сделаем это более очевидным, применив два выражения к функциям:

FX:=unapply(Fx,x,y): FY:=unapply(Fy,x,y):

Теперь можно оценить вектор силы в некоторой точке. Покажем, что сила является центральной, т. е. направлена из начала системы координат:

F0:=[FX(1,2),FY(1,2)];

В точке r = [1, 2] y-компонента вектора силы вдвое больше х-компоненты, чего достаточно для того, чтобы вектор силы был пропорционален вектору координаты. Уравнение Ньютона дает два связанных уравнения второго порядка:

NE1:=m*diff(x(t),t$2)=FX(x(t),y(t));

NE2:=m*diff(y(t),t$2)=FY(x(t),y(t));

Зададим начальное условие, выбрав точку на положительной части оси x, обладающую ненулевой компонентой скорости в направлении y, но не имеющую компоненты скорости в направлении x, т. е. векторы начального положения и скорости находятся под прямым углом:

x0:=1: y0:=0: vx0:=0: vy0:=6.28: # это условия для орбиты Земли
x0:=1: y0:=0: vx0:=0: vy0:=4.: # эти условия дают заметную эллиптичность
IC:=x(0)=x0,D(x)(0)=vx0,y(0)=y0,D(y)(0)=vy0;

sol:=dsolve({NE1,NE2,IC},{x(t),y(t)},numeric, output=listprocedure):
X:=subs(sol,x(t)): Y:=subs(sol,y(t)):
[X(1),Y(1)];

Можно рассчитать траекторию и создать график. Запишем X(t) and Y(t) в кавычках, потому что они вычисляются численно только для числовых значений t. Результат представлен на рис. 4.11.1.

plot(['X(t)','Y(t)'],t=0..3,title="x(t) and
y(t)",color=[blue,green],thickness=3);

Рис. 4.11.1

Параметрический график:

plot(['X(t)','Y(t)',t=0..3],title="parametric plot",
color=red,thickness=3,scaling=constrained);

Задание 4.11.1

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

Задание 4.11.2

Измените в потенциале зависимость от расстояния, немного отличающуюся от 1/r. Посмотрите, что произойдет с движением.

Чтобы лучше представить, что происходит, изобразим векторные величины вдоль траектории (рис. 4.11.2, 4.11.3):

Np:=30; T:=0.5;

Lr:=[seq([X(T/Np*i),Y(T/Np*i)],i=1..Np)]:
with(plots):
PLT:=plot(Lr,style=point):
display(PLT,scaling=constrained);

Рис. 4.11.2

VX:=subs(sol,diff(x(t),t)): VY:=subs(sol,diff(y(t),t)):
with(plottools):
scale_v:=0.1: scale_F:=300: # подгонка - методом проб и ошибок
for i from 1 to Np do: t_i:=i*T/Np:
PLv[i] := arrow([X(t_i),Y(t_i)],
evalm(scale_v*vector([VX(t_i),VY(t_i)])), .05, .2, .1,
color=green):
od:
for i from 1 to Np do: t_i:=i*T/Np:
PLF[i] := arrow([X(t_i),Y(t_i)],
evalm(scale_F*vector([FX(X(t_i),Y(t_i)),FY(X(t_i),Y(t_i))])),
.05, .2, .1, color=red):
od:
for i from 1 to Np do: t_i:=i*T/Np:
PLr[i] := arrow([0,0],vector([X(t_i),
Y(t_i)]), .05, .2, .1, color=blue): od:
display(seq(display(PLT,PLv[i],PLF[i],PLr[i]),i=1..Np),
insequence=true,scaling=constrained,view=[-1..1,-1..1]);

Рис. 4.11.3

Покажем, что происходит с потенциальной и кинетической энергиями как функциями времени (рис. 4.11.4):

T_kin:=t->m/2*(VX(t)^2+VY(t)^2);

T_kin(0),T_kin(0.25),T_kin(0.5);

V_pot:=t->subs(x=X(t),y=Y(t),V);

E_tot:=t->T_kin(t)+V_pot(t);

E_tot(0),E_tot(0.25),E_tot(0.5);

plot(['T_kin(t)','V_pot(t)','E_tot(t)'],t=0..1,
color=[red,blue,green],thickness=3,axes=boxed);

Рис. 4.11.4

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

Period:=0.5;

T_avg:=evalf(Int(T_kin,0..Period,5,_NCrule))/Period;

E_avg:=evalf(Int(E_tot,0..Period,5,_NCrule))/Period;

V_avg:=evalf(Int(V_pot,0..Period,5,_NCrule))/Period;

V_avg/E_avg;

Задание 4.11.3

Выполните такие же расчеты с другими начальными условиями и проверьте, что соотношения между средними значениями есть свойство ньютоновского закона притяжения (обратный квадрат расстояния). По графику проверьте правильность выбора периода. Попробуйте посмотреть, что случится с немного другими степенными зависимостями, например: n = 12/10, n = 8/10 и т. п.

Задача о радиальном движении

Радиальная задача для r(t) следует из сохранения углового момента импульса. Она позволяет определять траекторию из одного ДУ.

Сначала рассчитаем угловой момент импульса:

with(linalg):
L0:=evalm(m*crossprod([x0,y0,0],[vx0,vy0,0]));

t0:=2;

evalm(m*crossprod([X(t0),Y(t0),0],[VX(t0),VY(t0),0]));

Действительно, вычисленная траектория приблизительно сохраняет угловой момент импульса. Это – проверка точности численного алгоритма. Вектор углового момента перпендикулярен плоскости (xy). Напомним, что евклидовой нормой, или нормой Фробениуса называется квадратный корень из суммы квадратов элементов. Она есть в пакете linalg.

L_mag:=norm(L0,frobenius);

Одномерная задача для зависимости расстояния от центра притяжения от времени:

Vr:=-GM*m/r^n;

Veff:=Vr+L_mag^2/(2*m*r^2);

Feff:=unapply(-diff(Veff,r),r);

DE:=m*diff(r(t),t$2)=Feff(r(t));

r0:=sqrt(x0^2+y0^2);

vr0:=0; # это записано для конкретного выбора вектора направления нормали и векторов скорости.

IC_r:=r(0)=r0,D(r)(0)=vr0:
#sol_r:=dsolve({DE,IC_r},r(t)); # плохо работает
sol_r:=dsolve(DE,r(t));

Получившийся результат выглядит не очень пригодным (эта громоздкая выдача не показана – увидите сами). График показан на рис. 4.11.5.

sol_r:=dsolve({DE,IC_r},r(t),numeric,output=listprocedure):
R:=subs(sol_r,r(t)):
plot('R(t)',t=0..2,thickness=3);

Рис. 4.11.5

В случае орбиты Земли центробежные силы почти уравновешены гравитационными. При другом начальном условии результат дает существенное изменение радиального расстояния.

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

Упростим ситуацию, взяв потенциал притяжения вида –k/rn, и переместим в этом силовом поле единичную массу (график представлен на рис. 4.11.6).

restart;
n:=3/2; # сделайте это при n=3 и посмотрите что получается; затем при n=2.0001; затем при n=1.999...

Необходимое для равномерного кругового движения условие устойчивости: центростремительное ускорение обеспечивается силой притяжения.

STC:=0=Lsq/(m*r0^3)-n*k/r0^(n+1);

Lsq:=solve(STC,Lsq);

V_eff:=Lsq/(2*m*r^2)-k/r^n;

F_eff:=unapply(diff(-V_eff,r),r);

NE:=m*diff(r(t),t$2)=F_eff(r(t));

m:=1; k:=1;

r0:=1;

IC:=r(0)=r0,D(r)(0)=-1/100;

sol:=dsolve({NE,IC},r(t),numeric,output=listprocedure):
R:=subs(sol,r(t)):
plot('R(t)',t=0..100,thickness=3);

Рис. 4.11.6

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

P1:=plot(1/(2*r^2)-1/r,r=0..10,color=red):
P2:=plot(1/(2*r^2)-1/r^3,r=0..10,color=blue):
plots[display](P1,P2,view=[0..5,-2..1],
thickness=3,title="Effective Potential");

Синим на рис. 4.11.7 показано нестабильное решение.

Рис. 4.11.7

Нам также нужен график радиальной силы:

P3:=plot(-diff(1/(2*r^2)-1/r,r),
r=0..10,color=red,numpoints=300):
P4:=plot(-diff(1/(2*r^2)-1/r^3,r),
r=0..10,color=blue,numpoints=300):
plots[display](P3,P4,view=[0..5,-1..0.5],
thickness=3,title="Effective Force");

Красная кривая на рис. 4.11.7 показывает график возвращающей силы, где точка равновесия находится в точке r = 1.

Еще одна ситуация, в которой круговые орбиты устойчивы, – это случай удерживающих (линейных или более мощных) потенциалов в трех измерениях (рис. 4.11.8, 4.1.9):

P1:=plot(1/(2*r^2)+r,r=0..10,color=red):
P2:=plot(1/(2*r^2)+r^2,r=0..10,color=blue):
plots[display](P1,P2,view=[0..5,0..10],
thickness=3, title="Effective Potential");

Рис. 4.11.8

P3:=plot(-diff(1/(2*r^2)+r,r),r=0..10, color=red,numpoints=300):
P4:=plot(-diff(1/(2*r^2)+r^2,r),r=0..10, color=blue,numpoints=300):
plots[display](P3,P4,view=[0..5,-10..10], thickness=3,title="Effective Force");

Рис. 4.11.9

Задача о рассеянии

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

Предположим, что частица находится очень далеко от центра рассеяния, поэтому притяжение можно игнорировать. «Свободная» частица характеризуется:

вектором асимптотической скорости и связанной с ним кинетической энергией;
вектором координаты, который вместе с вектором скорости и массой частицы определяет вектор углового момента.

Применим Maple для получения численных решений:

restart; with(linalg):
m:=1: r0:=[-100,5,0]: v0:=[1,0,0]:
L0:=evalm(m*crossprod(r0,v0));

Lsq:=dotprod(L0,L0);

n:=1;

V:=-1/r^n;

V:=subs(r=sqrt(x^2+y^2),V);

T0:=m/2*dotprod(v0,v0);

E0:=evalf(T0+subs(x=r0[1],y=r0[2],V));

Видно, что на расстоянии 100 единиц частица не является полностью «свободной». Фактически закон силы (квадрат обратного расстояния) не спадает достаточно быстро на больших расстояниях, чтобы по-настоящему позволить частице быть совсем свободной при конечном расстоянии (рис. 4.11.10, 4.11.11):

F_x:=unapply(diff(-V,x),x,y);

F_y:=unapply(diff(-V,y),x,y):
NE1:=m*diff(x(t),t$2)=F_x(x(t),y(t));

NE2:=m*diff(y(t),t$2)=F_y(x(t),y(t)):
IC:=x(0)=r0[1],D(x)(0)=v0[1],y(0)=r0[2],D(y)(0)=v0[2];

sol:=dsolve({NE1,NE2,IC},{x(t),y(t)},numeric, output=listprocedure):
X:=subs(sol,x(t)): Y:=subs(sol,y(t)):
PL0:=plot([[0,0]],style=point,symbol=cross, color=red,symbolsize=20):
PL1:=plot(['X(t)','Y(t)',t=0..200],axes=boxed, color=green,thickness=3): plots[display](PL1,PL0);