Законы ньютона во времени

[reqyz]

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

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

из второго закона ньютона получаем, что ускорение равно отношению скорости на время.
И, казалось бы что есть ответы на все вопросы, но..

верно ли я понимаю:
1) вектор направления двух притягивающихся объектов нужно нормализовывать?
2) вектор направления до другого тела нужно складывать с текущим вектором направления, или вектор направления не должен меняться, а только модифицироваться координаты.

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

procedure TWorld.NextTime(Time:Double);
var
  i1,i2:cardinal;
  U,R2,L:Double;
  V:TCoordinates;
begin
  if(Length(Bodies)=0)then
    exit;
  for i1:=0 to Length(Bodies)-2 do
    for i2:=i1+1 to Length(Bodies)-1 do
    begin
      R2:=GetDistance2(Bodies[i1].Coordinates,Bodies[i2].Coordinates);
      U:=(Bodies[i2].Mass/R2)*G*Time;
      with V do
      begin
        X:=Bodies[i2].Coordinates.X-Bodies[i1].Coordinates.X;
        Y:=Bodies[i2].Coordinates.Y-Bodies[i1].Coordinates.Y;
        Z:=Bodies[i2].Coordinates.Z-Bodies[i1].Coordinates.Z;
        L:=sqrt(X*X+Y*Y+Z*Z);
        X:=(X/L)*U;
        Y:=(Y/L)*U;
        Z:=(Z/L)*U;
      end;


      with Bodies[i1].Coordinates do
      begin
        X:=X+V.X;
        Y:=Y+V.Y;
        Z:=Z+V.Z;
        {L:=sqrt(X*X+Y*Y+Z*Z);
        X:=X/L;
        Y:=Y/L;
        Z:=Z/L};
      end;

      V.X:=-V.X;
      V.Y:=-V.Y;
      V.Z:=-V.Z;

      with Bodies[i2].Coordinates do
      begin
        X:=X+V.X;
        Y:=Y+V.Y;
        Z:=Z+V.Z;
      end;

    end; 
  for i1:=0 to Length(Bodies)-1 do
    with Bodies[i1] do
    begin
      Coordinates.X:=Coordinates.X+Route.X*Speed*Time;
      Coordinates.Y:=Coordinates.Y+Route.Y*Speed*Time;
      Coordinates.Z:=Coordinates.Z+Route.Z*Speed*Time;
    end;
end;

мне действительно нужна помощь, если не понятен мой вопрос, напишите, постараюсь переформулировать или пояснить

[gris]

А зачем отдельно задавать направления и следить за ними? Удобнее векторно задать координаты, скорость и ускорение каждого тела. Кстати, координаты точнее считать по фомуле равноускоренного, а не равномерного движения. Если массы не меняются и тел немного, то можно заранее просчитать массив произведений масс. Не вглядывался — а учитываете ли Вы третий закон Ньютона?

[reqyz]

gris в сообщении #925485 писал(а):

А зачем отдельно задавать направления и следить за ними? Удобнее векторно задать координаты, скорость и ускорение каждого тела. Кстати, координаты точнее считать по фомуле равноускоренного, а не равномерного движения. Если массы не меняются и тел немного, то можно заранее просчитать массив произведений масс. Не вглядывался — а учитываете ли Вы третий закон Ньютона?

координаты заданы вектором (x,y,z) направление тоже логично вектором (x,y,z)заданы так же скорость каждого тела, Зачем задавать ускорение? оно будет зависеть от времени через которое мы обновляем "мир", в коде я и использовал формулу равноускоренного движения

  TCoordinates=record
    X,Y,Z:double;
  end;

  TBody=record                //òåëî èëè ñðåäà
    Mass:double;              //ìàññà
    Coordinates:TCoordinates; //êîîðäèíàòû
    Route:TCoordinates;       //íàïðàâëåíèå
    Speed:double;             //ñêîðîñòü
    //Density:double;         //ïëîòíîñòü p=mass/V  V-îáüåì, ïîêà ÷òî ðàäèóñ
    Transparency:double;      //ïðîçðà÷íîñòü
    Color:cardinal;             //öâåò

    Radius:double;            //ðàäèóñ
  end;

U:=(Bodies[i2].Mass/R2)*G*Time;
  with V do
  begin
     X:=(X/L)*U;
     Y:=(Y/L)*U;
     Z:=(Z/L)*U;
  end;

  with Bodies[i1].Coordinates do
  begin
    X:=X+V.X;
    Y:=Y+V.Y;
    Z:=Z+V.Z;
  end;

хотя может где-то что-то напутал, отталкивание я тоже пока не добавлял

[arseniiv]

Почему вы скорость-то вектором не задаёте? Что такое «направление»?

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

[reqyz]

arseniiv в сообщении #925502 писал(а):

Почему вы скорость-то вектором не задаёте? Что такое «направление»?

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

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

[gris]

Я имел в виду, что направление у Вас фактически является скоростью, только Вы его нормализуете, а потом умножаете на модуль скорости, который непонятно, как получаете. Может быть, я просто плохо разглядел этот зелёный текст. Вы бы сначала написали математическую модель.
Про движение я вот что имел в виду. Вы для положения $X_0$ находите скорость и определяете $X_1=X_0+V_1\cdot \Delta t$, перед этим найдя ускорение в $X_0$. То есть у Вас получается $X_1=X_0+V_0\cdot \Delta t+A_0\cdot (\Delta t)^2$. А на двоечку поделить?

[reqyz]

gris в сообщении #925506 писал(а):

Я имел в виду, что направление у Вас фактически является скоростью, только Вы его нормализуете, а потом умножаете на модуль скорости, который непонятно, как получаете. Может быть, я просто плохо разглядел этот зелёный текст. Вы бы сначала написали математическую модель.
Про движение я вот что имел в виду. Вы для положения $X_0$ находите скорость и определяете $X_1=X_0+V_1\cdot \Delta t$, перед этим найдя ускорение в $X_0$. То есть у Вас получается $X_1=X_0+V_0\cdot \Delta t+A_0\cdot (\Delta t)^2$. А на двоечку поделить?

procedure TWorld.NextTime(Time:Double);
var
  i1,i2:cardinal;
  U,R2,L:Double;
  V:TCoordinates;
begin
  if(Length(Bodies)=0)then
    exit;
  for i1:=0 to Length(Bodies)-2 do        //прохожусь по всем телам
    for i2:=i1+1 to Length(Bodies)-1 do   //сравниваю со следующими телами(с предыдущими сравнивается в верхнем цикле)
    begin
      R2:=GetDistance2(Bodies[i1].Coordinates,Bodies[i2].Coordinates);    //нахожу квадрат расстояния мужду центрами тел
      U:=(Bodies[i2].Mass/R2)*G*Time;                                     //нахожу скорость через ускорение и время

      V:=VectorDifference(Bodies[i2].Coordinates,Bodies[i1].Coordinates); //нахожу направление от центра нашего тела до сравниваемого
      V:=NormalizeVector(V);                                              //нормализую вектор этого направления
      V:=VectorMultiplication(V,U);                                       //умножааю вектор на скорость, не знаю как по другому(
      Bodies[i1].Coordinates:=VectorSum(Bodies[i1].Coordinates,V);        //нахожу новые координаты нашего тела с учетом силы притяжения к сравниваемому телу
      V:=InverseVector(V);                                                //меняю направление вектора
      Bodies[i2].Coordinates:=VectorSum(Bodies[i1].Coordinates,V);        //нахожу новые координаты сравниваемого тела, относительно действия силы притяжения нашего тела

    end;
  for i1:=0 to Length(Bodies)-1 do
    with Bodies[i1] do
      VectorSum(Coordinates,VectorMultiplication(Speed,Time));            //нахожу новые координаты тел, учитывая их собственные скорости передвижения
end;

п.с. использовал второй и первый закон ньютона, почему $Delta t$должна делиться на двойку? может я не верно где то преобразовал или использую какую то не ту формулу, по сути у мен тут три формулы: $F = G \cdot {m_1 \cdot m_2\over R^2} ................\vec{a} = \frac{\vec{F}}{m}, .................\vec a={d\vec v \over dt}.$

[gris]

Ну да, я немножко поторопился. Произведения масс не нужны. Скорее всего, по аддитивности у Вас получается то, что нужно, но нет ли здесь лишних вычислений? Почему не посчитать для каждого тела векторную сумму ускорений по всем остальным телам, а потом уже не определить его новое положение. Только нужно следить за тем, чтобы ускорения считались по старым положениям и скоростям, а Вы их, по-моему, обновляете по мере продвижения по циклу.
Я имел в виду формулу $x=x_0+v_0t+\dfrac{a_0t^2}{2}$

[reqyz]

gris в сообщении #925522 писал(а):

Ну да, я немножко поторопился. Произведения масс не нужны. Скорее всего, по аддитивности у Вас получается то, что нужно, но нет ли здесь лишних вычислений? Почему не посчитать для каждого тела векторную сумму ускорений по всем остальным телам, а потом уже не определить его новое положение. Только нужно следить за тем, чтобы ускорения считались по старым положениям и скоростям, а Вы их, по-моему, обновляете по мере продвижения по циклу.
Я имел в виду формулу $x=x_0+v_0t+\dfrac{a_0t^2}{2}$

да, Вы правы я думал о том что бы считать векторную сумму для каждого тела, думаю действительно, можно пожертвовать ресурсами памяти для ускорения процессора и считать всё же суммарно. по остальным пунктам, можете подсказать, действительно ли мне надо менять координаты а не направление объекта? я сам не уверен как будет правильно..

[gris]

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

[Munin]

reqyz в сообщении #925465 писал(а):

верно ли я понимаю:
1) вектор направления двух притягивающихся объектов нужно нормализовывать?

Проще всего вектор направления вообще не использовать, а закон притяжения использовать в виде $\vec{F}=\dfrac{Gm_1m_2}{r^3}\vec{r}.$

reqyz в сообщении #925494 писал(а):

координаты заданы вектором (x,y,z) направление тоже логично вектором (x,y,z)

Для разных величин не стоит вводить одинаковые обозначения: запутаетесь.

reqyz в сообщении #925537 писал(а):

думаю действительно, можно пожертвовать ресурсами памяти для ускорения процессора и считать всё же суммарно.

Думать о ресурсах вам сейчас совершенно некритично.

Последовательность такая:
1. Добиться того, чтобы расчётная модель совпадала с физическими формулами, и выдавала "физически реалистичный, похожий на правду" результат.
2. Дальше можно заняться численными методами, и выбрать более удобный метод. Также сразу к этому шагу придётся перейти, если на шаге 1 возникнут численные проблемы: сходимость, устойчивость.
3. Только когда у вас будет хорошо отлаженная физическая модель и численный метод, можно заниматься уже оптимизациями в программном смысле. И то, возможно, вы поменяете метод, и все оптимизации станут не нужны.

-- 02.11.2014 21:58:33 --

Для листингов лучше использовать не тег code, а тег syntax lang="pascal".

[arseniiv]

Кстати ещё, паскали не запрещают писать каждый комментарий на отдельной строке, а не справа от кода (обычно при этом комментарий будет стоять уже перед описываемым кодом, а не после). Мнения могут быть разными, но, по-моему, так будет читаемее. А для комментария на одной строке с кодом нужна будет причина (он маленький, и кода тоже мало; закомментирован кусок строки и пр.).

[Munin]

Если комментарии даются к каждой строке - то удобнее справа.

Правда, в одном из листингов в комментариях слетела кодировка.