Модель движения системы материальных точек

1. Задача. Имеется система N материальных точек с массами m_i , i = 1, 2, ..., N, взаимодействующих друг с другом с внутренними силами, на каждую из которых действует внешняя сила. Исходя из начальных координат x_i, y_i, и скоростей v_xi, v_yi, определите координаты и скорости материальных точек в последующие моменты времени.

2. Теория. Рассмотрим механическую систему из N материальных точек. Основной закон динамики:

где F'_ij -- внутренняя сила, действующая на i-ую материальную точку со стороны j -ой материальной точки, F_i -- равнодействующая внешних сил, действующих на i - ую материальную точку со стороны тел не входящих в систему.

Дифференциальное уравнение второго порядка может быть представлено двумя дифференциальными уравнениями первого порядка. Имеем:

Зная внешние и внутренние силы, действующие на каждую материальную точку, можно определить их ускорения. Исходя из координат и скоростей точки в момент времени t, можно расчитать координаты и скорости точки в следующий момент времени t + Δ t.

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

3. Алгоритм.

1. Задают число материальных точек N, их массы m_i , координаты x_i , y_i и проекции начальных скоростей v_ix , v_iy , силовое поле F_x = F_x (x,y), F_y = F_y (x,y), а также шаг по времени Δ t.

2. Начало цикла по t. Дают приращение по времени: переменной t присваивают значение t + Δ t.

3. Определяют проекции F_xi , F_yi равнодействующей всех внешних и внутренних сил, действующих на каждую i - ую материальную точку в момент t + Δt, и записывают их в массивы.

4. В цикле переобозначают координаты всех материальных точек, записывая их в массивы xx[i], yy[i].

5. В цикле перебирают все материальные точки и определяют проекции ускорения, скорости и координаты для каждой из них в момент t + Δt:

a_xi (t + Δt) = F_xi (t + Δt)/m_i ,

v_ix (t + Δ t) = v_ix (t) + a_ix (t + Δt)Δt,

x_i (t + Δ t) = x_i (t) + v_ix (t + Δt)Δt.

По аналогичным формулам вычисляют проекции на ось OY. Результаты записывают в массивы x[i], y[i], vx[i], vy[i].

6. Стирают изображения материальных точек в предыдущий момент времени t, координаты которых сохранены в массивах xx[i], yy[i].

7. На экране строят точки в следующий момент t + Δt, либо рисуют графики или выводят результат в числовом виде.

8. Возвращение к операции 2. Если цикл по t закончился, -- выход из цикла.

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

program PROGRAMMA4;uses dos, crt, graph;const N=50; dt=0.01;var m,Fx,Fy,x,y,vx,vy,xx,yy : array[1..N] of real; Gd, Gm, i, j : integer; ax, ay, F, l : real;label Metka, metka1;Procedure Init_Graph; {Инициализация графики}begin Gd:=Detect; InitGraph(Gd, Gm, 'c:\bp\bgi'); if GraphResult <> grOk then Halt(1);end;Procedure Sila; {Вычисление действующих сил}label Metka;begin For i:=1 to N do begin Fx[i]:=0; Fy[i]:=0; end; For i:=1 to N do for j:=1 to N do begin if j=i then goto Metka; l:=sqrt(sqr(x[i]-x[j])+sqr(y[i]-y[j])); if l<2 then l:=2; F:=-50000*m[i]*m[j]/sqr(l)+500000*m[i]*m[j]/sqr(l*l); Fx[i]:=Fx[i]+F*(x[i]-x[j])/l; Fy[i]:=Fy[i]+F*(y[i]-y[j])/l+m[i]*10; Metka: end;end;Procedure Nach_uslov;begin Randomize; {Задание случайных координат и скоростей} for i:=1 to N do begin m[i]:=2; x[i]:=random(280)+60; y[i]:=random(280)+60; vy[i]:=random(30)-15; vx[i]:=random(30)-15;end; end;BEGINInit_Graph; Nach_uslov;Repeat Sila;for i:=1 to N do begin xx[i]:=x[i]; yy[i]:=y[i]; {Запись предыдущих координат} ax:=Fx[i]/m[i]; ay:=Fy[i]/m[i]; {Вычисление ускорений,} vx[i]:=vx[i]+ax*dt; vy[i]:=vy[i]+ay*dt; {скоростей,} x[i]:=x[i]+vx[i]*dt; y[i]:=y[i]+vy[i]*dt; {координат} if (x[i]<50)or(x[i]>350) then vx[i]:=-vx[i];{отражение} if (y[i]<50)or(y[i]>350) then vy[i]:=-vy[i];{от стенок} end; delay(500); setcolor(8); for i:=1 to N do circle(round(xx[i]),round(yy[i]),2);setcolor(15); for i:=1 to N do circle(round(x[i]),round(y[i]),2);until KeyPressed; Repeat until keypressed; CloseGraph;END.

Тема 3. Математические системы. Реализация алгоритма для математических систем