Визуализация траекторий движения небесных тел и релятивистских частиц в псевдотрехмерном пространстве
Орлов Иван, Зисман Михаил,
г. Челябинск, лицей №37, 11 класс
Введение
Проблема построения траекторий движения небесных тел долгое время интересует исследователей. Известно, что на основе второго закона Ньютона и уравнения гравитационного взаимодействия можно с успехом решить задачу двух тел. С возрастанием количества объектов решение этой задачи становится не под силу человеку. С появлением общей и частной теории относительности (ОТО и СТО соответственно) возрос интерес к исследованию траекторий движения частиц, перемещающихся с субсветовой скоростью в гравитационном и электромагнитном поле. Вначале все исследования такого движения были аналитическими. Компьютеры помогают решить задачу n тел в небесной механике и задачи ОТО и СТО.
Цель работы – моделирование движения релятивистских частиц и n тел в небесной механике, посредством компьютерных программ.
Трудности, связанные с разработкой подобных программ, заключались в отсутствии подходящих для этого языков программирования и необходимости создания собственного псевдотрехмерного пространства. В настоящей работе нам удалось создать программу, использующую собственное псевдотрехмерное пространство и не требующую высоких системных возможностей. В представленной работе нами рассмотрены наиболее «интересные» с точки зрения физики и небесной механики траектории движения тел и частиц.
Для достижения этой цели нами были проведены исследования возможностей компьютерного моделирования плоских и объемных траекторий, а также погрешностей, возникающих в процессе вычислений и их причин.
Математическая модель
Рассмотрим математическую модель задачи n тел в небесной механике. При ее построении учтем гравитационное взаимодействие (магнитное взаимодействие не учитывается, поскольку сильным магнитным полем обладают лишь весьма специфические небесные объекты, например, пульсары [1]). Приведенная ниже математическая модель лежит в основе главного исполнительного модуля программ. Для вычисления силы притяжения необходимо воспользоваться формулой гравитационного взаимодействия (F=GmM/r2), где m и M – массы взаимодействующих тел, r – расстояние между ними. Для задачи двух или более тел получение основных формул сводится ко второму закону Ньютона ma=F, где F – сила, действующая на тело, а – ускорение (жирный шрифт означает векторную величину).
Для получения необходимых формул распишем второй закон Ньютона в проекциях на соответствующие оси координат:
mAx=–GmM*x/r3 , mAy =–GmM*y/r3 , mAz =–GmM*z/r3, где A – ускорение, x, y, z – проекции радиус-вектора на соответствующие оси, r – расстояние между телами (см. рис.).
Так как ускорение есть производная скорости по времени, то заменим A на V¢, где V¢»(Vi+1–Vi)/h, h – шаг по шкале времени, Vi – значение скорости в i-тый момент времени. Шаг – это промежуток времени, на котором движение тела можно считать равноускоренным. Таким образом, чем меньше шаг, тем точней будет траектория. Применяемый метод расчета называется методом Эйлера [2].
Запишем основные уравнения, лежащие в основе численных расчетов, для задачи n тел.
Vi+1= Vi –hGM1(x-x1)/((x-x1)2+(y-y1)2+(z-z1)2)3/2 –hGM2(x-x2)/((x-x2)2+(y-y2)2+(z-z2)2)3/2–
–…–hGMn-1(x-xn-1)/ ((x-xn-1)2+(y-yn-1)2+(z-zn-1)2)3/2
Xi+1=Xi+hVi
Vi+1= Vi –hGM1(y-y1)/((x-x1)2+(y-y1)2+(z-z1)2)3/2 –hGM2(y-y2)/((x-x2)2+(y-y2)2+(z-z2)2)3/2–
–…–hGMn-1(y-yn-1)/ ((x-xn-1)2+(y-yn-1)2+(z-zn-1)2)3/2
Yi+1=Yi+hVi
Vi+1= Vi –hGM1(z-z1)/((x-x1)2+(y-y1)2+(z-z1)2)3/2 –hGM2(z-z2)/((x-x2)2+(y-y2)2+(z-z2)2)3/2–
–…–hGMn-1(z-zn-1)/ ((x-xn-1)2+(y-yn-1)2+(z-zn-1)2)3/2
Zi+1=Zi+hVi
Теперь выведем уравнения для расчета движения релятивистской частицы в гравитационном и электростатическом полях [3]. Известно, что производная импульса по времени равна силе, действующей на тело, значит dp/dt=F. Тогда, расписывая импульс и силу, получаем d(mV/(1-V2/c2)1/2)/dt=–(kQq+Gm0m)r)/r3(жирный шрифт означает векторную величину), где m и m0 – масса частицы и центрального тела, V – скорость частицы, Q и q – заряды центрального тела и частицы, r – расстояние от центрального тела до частицы, с – скорость света. Если разложить расстояние и скорость на составляющие, то мы получим следующее: r2=x2+y2, V2=Vx2+ Vy2, или V2=(x')2+(y')2, где x' и y' – первые производные от координат по времени. Распишем первоначальную систему для получения основных формул.
d(mx'/(1-(x'2+y'2)/c2)1/2)/dt=–(kQq+Gm0m)x/(x2+y2)3/2
d(my'/(1-(x'2+y'2)/c2)1/2)/dt=–(kQq+Gm0m)y/(x2+y2)3/2
Проведя необходимые преобразования, в левой части уравнений получим выражения cm(x''(c2-y'2)+x'y'y'')/(c2-x'2-y'2)3/2, и cm(y''(c2-x'2)+y'x'x'')/(c2-x'2-y'2)3/2, где x'' и y'' – вторые производные координат по времени. Выразив из левой части первого уравнения x'', подставив его во второе и проведя преобразования, получим окончательную систему уравнений для расчета скорости и координаты как функций времени:
Vxi+1=Vxi–h*(kQq+Gm0m)(c2-Vxi2-Vyi2)3/2(x(c2-Vyi2)-xyVyi)/cm(x2+y2)3/2((c2-Vxi2)(c2-Vyi2)–xVyi2Vxi)
Vyi+1=Vyi–h*(kQq+Gm0m)(c2-Vxi2-Vyi2)3/2(y(c2-Vyi2)-xyVxi)/cm(x2+y2)3/2((c2-Vxi2)(c2-Vyi2)–yVxi2Vyi)
Xi+1=Xi+h*Vxi
Yi+1=Yi+h*Vyi
Если центральное тело сферически симметричное, то его можно представить как материальную точку, в которой сосредоточена вся его масса. Если принять во внимание это свойство тела, то задачу для релятивистской частицы можно рассматривать в плоскости в силу закона сохранения момента импульса [3]. Поэтому представленная здесь система уравнений рассчитывает движение частицы в плоскости, а не в объеме.
На основе данной математической модели были составлены основные исполнительные модули программ. Алгоритм работы программы в общих чертах таков: программа просчитывает параметры заданной системы тел для каждого кадра и записывает положение всех тел в скрипт, который уже исполняется непосредственно в среде 3ds max (Скрипт – специальный код, исполняемый только в среде 3ds max). Либо, если установить опцию использования OpenGL, исполнение скрипта будет осуществляться в созданном нами псевдотрехмерном пространстве. Прорисовка тел и траекторий в 3ds max более наглядна, а недостаток заключается в больших требованиях к аппаратной части системы по сравнению с OpenGL. Детально принцип работы описан в приложении.
Рассмотрим некоторые примеры применения программы для задачи n тел.
1) Для начала рассмотрим систему четырех тел, расположенных в вершинах квадрата. Если массы всех тел будут равны, то, в случае нулевых начальных скоростей, результирующая сила для каждого тела будет направлена по диагонали квадрата, и их траектории движения сойдутся в центре, где расположен центр масс этой системы. Если массы будут неравными, то траектория тел сместится от центра квадрата. Теперь зададим начальные скорости всех тел, причем векторы скоростей направим под одним углом к сторонам квадрата (модули скоростей равны): тело №1 – (0;0;0) Vy=–10 м/с, m=1,5*1017кг; №2 – (105;0;0) Vx=–10 м/с, m=1,5*1017кг; №3 – (105;-105;0) Vy=10 м/с, m=1,5*1017кг; №4 – (0; –1e5;0) Vx=10 м/с, m=1,5*1017кг (параметры, стоящие в скобках, означают координаты тела в начальный момент времени; в дальнейшем все записи вида 1,5*1017 будут заменяться эквивалентной записью вида 1,5e17). При этих параметрах получатся следующие траектории (см. приложение, рис.1).
Казалось бы: мы задали одинаковые расстояния и скорости тел, а, в конце концов, система распадается. Данный эксперимент показал, что, несмотря на симметричность системы, она является нестабильной.
2) Теперь проверим на стабильность модель куба. Тела верхней грани куба пусть движутся «по часовой стрелке», а нижней – наоборот. Размеры квадратов и скорости тел пусть будут теми же, что и в предыдущей модели. На рисунке видно, что такая система тел распадается гораздо раньше, чем предыдущая, траектории движения тел более сложные (см. приложение, рис.2).
3) В ходе исследований мы убедились в существовании относительно стабильных систем тел (двойных звезд). Рассмотрим одну из таких систем. В центре масс мы расположили еще одно тело. Его предназначение – показать, что объект, расположенный в центре масс такой системы, тоже обладает относительной стабильностью, при симметричности модели. Параметры системы: №1 (центральное) – (0;0;0) Vxyz=0 m=1e10кг; №2 – (100;100;0) Vy=–5e-2 м/с, m=1e7кг; №3 – (-100;-100;0) Vy=5e-2 м/с, m=1e7кг (см. приложение, рис.3).
Для окончательной проверки всей программы решено было провести серию «контрольных» тестов на известных примерах, с целью выявления ошибок и погрешностей. Для тестов была выбрана наша Солнечная система. Рассмотрим некоторые из «контрольных» тестов:
4) В первом тесте нашей задачей было получить траектории движения нижних планет Солнечной системы. Все параметры планет задавались в соответствии с реальными, за исключением их размеров. При реальных радиусах планет, их невозможно будет увидеть, рассматривая всю систему целиком. Результат расчета траекторий получился близким к реальному. При шаге, равном 40000 секунд, относительная погрешность составила »2,24% (см. приложение, рис.4). На рисунке видно отклонение орбиты Меркурия от окружности (в действительности, эксцентриситет орбиты Меркурия равен 0,21).
5) Чтобы окончательно убедиться в правильности построения математической модели программы, мы проверили траекторию движения Луны вокруг Земли, относительно Солнца. Результат отличался от реального на »4,3% (см. приложение, рис.5). На рисунке видно, что траектория Луны представляет собой траекторию знакопеременной кривизны (в реальности это действительно так).
В приложении (рис.6) приведены траектории в задаче двух тел.
Теперь рассмотрим некоторые тесты программы для задач СТО.
1) В первом тесте мы решили проверить, будет ли отклоняться частица, летящая с субсветовой скоростью, в гравитационном и электростатическом поле. Параметры системы тел: частица – (1e7;1e5) m=1,6e-27 Vx=-2e8 Vy=0 q=1,67e-19 центральное тело – m0=2e30, Q=5 при шаге в 1e-5 секунд и времени просмотра в 0,1 сек. Получаем (см. приложение рис.7). В данном примере мы рассмотрели модель отклонения заряженной частицы вблизи звезды, подобной Солнцу, но саму частицу расположили в 1000 раз ближе, чем радиус Солнца, чтобы показать отклонение (если частица будет пролетать на расстоянии радиуса Солнца, то отклонение будет в 1,75 секунды). Отклонение же в нашем случае равно ≈2,860.
2) Во втором тесте скорость частицы много меньше скорости света. Параметры системы тел: частица – (-1,5e8;0) m=1,6e-27 Vx=0 Vy=3e6 q=3e-16 центральное тело – m0=2e30, Q=1 при шаге в 3e-1 секунд и времени просмотра в 8000 сек. получаем спираль (см. приложение рис.8).
3) В третьем тесте увеличим массу центрального тела: частица – (-1,5e8;0) m=1,6e-27 Vx=0 Vy=3e6 q=3e-16 центральное тело – m0=1,1e32, Q=1 при шаге в 3e-2 секунд (см. приложение рис.9).
4) В четвертом тесте увеличим заряд центрального тела и уменьшим заряд у частицы: частица – (-1,5e8;0) m=1,6e-27 Vx=0 Vy=3e6 q=2e-16 центральное тело – m0=2e30, Q=4 при шаге в 3e-1 секунд (см. приложение рис.10).
Рисунки для задачи n тел в небесной механике (см. приложение, рис1 – рис6) выполнены с помощью программы 3ds Max, поскольку в ней тела прорисовываются сплошными, и возможна одновременная прорисовка тел и их траекторий.
Погрешность, вносимая компьютером, составляет около 10-11%. Это происходит из-за невозможности учета всех знаков числа после запятой. По сравнению с ошибкой, вносимой шагом, ошибкой компьютера можно пренебречь. Ошибка, внесенная шагом (шаг равен 78840 секунд) в системе Солнце – Земля – Луна составляет »4,4%. Для исследования изменения ошибки, мы уменьшили шаг в 2 раза. При этом погрешность стала равна »2,2%. Она никогда не станет равной нулю, потому что для этого нам придется задавать шаг, стремящийся к нулю.
Основные выводы:
1. созданные программы позволяют исследовать траектории движения релятивистских частиц и n тел в небесной механике с использованием собственного псевдотрехмерного пространства, не требующего высоких системных возможностей;
2. представленная в работе модель программы применима как для индивидуального обучения, так и для научной работы.
Литература
1. Энциклопедия «Современное естествознание», том 4. «Астрофизика». М.: Магистр-пресс.2000, 250 с.
2. Энциклопедия «Математика». М.: Аванта+, 1998, 682 с.
3. , Лифшиц поля. М.: Наука, 1988. – 510с.
Приложение
Модель квадрата
Модель куба
Модель «двойной звезды»
Солнечная система (нижние планеты)
Траектория Земли и Луны
Траектории в системе двух тел
Отклонение частицы равно ≈2,860
Траектории движения релятивистской частицы в кулоновском поле (рис.8-рис.10)
Основные принципы построения программы
Ниже представлен код основного модуля программы, отвечающего за обработку исходных данных и получения координат тел в следующий момент времени.
Модуль разработан на языке Object Pascal в среде Delphi 6.
Рассмотрим основной алгоритм его работы:
1) Из оболочки вызывается процедура move с входными данными – массивом данных всех тел (тип а:game), количеством тел и шагом времени.
2) Выполняется процедура nc для текущего тела, определяющая положение и скорость данного тела в следующий момент времени по формулам, приведенным в основной части работы. Эта процедура выполняется для каждого тела отдельно.
3) Полученные данные записываются в выходной массив данных (тип b:game).
4) Из оболочки данные записываются в файл для дальнейшей обработки.
5) Входной массив данных a присваивается выходному массиву b. Текущее время увеличивается на время одного шага.
6) Повторяется с пункта 1. Если текущее время станет равным или превысит заданное время просчета, просчет данных завершается.
Примечания:
1) Массив данных для одного тела представляет собой набор чисел, описывающий его основные параметры (положение, скорость, масса и заряд).
2) Текущее время показывает состояние счетчика времени с начала просчета
unit Model_GK;
interface
type
vector=array [1..3] of extended;
coord=record
x, v:vector;
q, m:extended;
end;
game =array [1..400] of coord;
const k=9e9;
G=6.672e-11;
procedure move (a:game;var b:game; n:integer; h:extended); // предназначение – расчет
// положения тела в следующий момент времени. Исходные данные – входной массив
// данных, выходной массив данных, количество тел, шаг по шкале времени.
implementation
procedure move (a:game;var b:game; n:integer; h:extended);
var i:integer;
procedure nc(l:integer);
var j:integer;
e:extended;
begin
for j:=1 to n do
if (j<>l) then
begin
e:=sqr(a[l].x[1]-a[j].x[1])+sqr(a[l].x[2]-a[j].x[2])+sqr(a[l].x[3]-a[j].x[3]);
if (e<>0) and (a[l].m<>0) then b[l].v[1]:=b[l].v[1]-h*G*a[j].m*(a[l].x[1]-a[j].x[1])/sqrt(e*e*e)+h*k*a[j].q*a[l].q*(a[l].x[1]-a[j].x[1])/(sqrt(e*e*e)*a[l].m);
end;
b[l].x[1]:=b[l].x[1]+h*b[l].v[1];
for j:=1 to n do
if (j<>l) then
begin
e:=sqr(a[l].x[1]-a[j].x[1])+sqr(a[l].x[2]-a[j].x[2])+sqr(a[l].x[3]-a[j].x[3]);
if (e<>0) and (a[l].m<>0) then b[l].v[2]:=b[l].v[2]-h*G*a[j].m*(a[l].x[2]-a[j].x[2])/sqrt(e*e*e)+h*k*a[j].q*a[l].q*(a[l].x[2]-a[j].x[2])/(sqrt(e*e*e)*a[l].m);
end;
b[l].x[2]:=b[l].x[2]+h*b[l].v[2];
for j:=1 to n do
if (j<>l) then
begin
e:=sqr(a[l].x[1]-a[j].x[1])+sqr(a[l].x[2]-a[j].x[2])+sqr(a[l].x[3]-a[j].x[3]);
if (e<>0) and (a[l].m<>0) then b[l].v[3]:=b[l].v[3]-h*G*a[j].m*(a[l].x[3]-a[j].x[3])/sqrt(e*e*e)+h*k*a[j].q*a[l].q*(a[l].x[3]-a[j].x[3])/(sqrt(e*e*e)*a[l].m);
end;
b[l].x[3]:=b[l].x[3]+h*b[l].v[3];
end;
begin
b:=a;
for i:=1 to n do nc(i);
end;
begin
end.
Представленная программа разрабатывалась на компьютере AMD Duron 900 Mhz. Минимальные требования к системе, на которой будет осуществляться запуск программы: Pentium 100, Windows 95 OSR 2, 16 MB RAM, клавиатура, мышь, рекомендуется установленный 3ds Max v4-5.
