МАЙЕР Р.В.

ИСПОЛЬЗОВАНИЕ КОМПЬЮТЕРНЫХ МОДЕЛЕЙ ПРИ ИЗУЧЕНИИ МЕХАНИКИ

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

1. ОСНОВЫ КИНЕМАТИКИ

Задача 1.1. Первые t1=8 с точка двигалась с ускорением a1(t)=1*t. До момента t2=30 с ускорение точки a2=0,2 м/с, после чего движение происходит с ускорением a3=-4 м/с. Постройте графики зависимостей координаты, скорости и ускорения от времени, если начальные координата и скорость точки равны -200 м и -15 м/с.

 ====== Программа ПР -- 1.1.
SCREEN 12: dt = .002: x = -200: v = -15           'QBASIC
WHILE INKEY$ = ""
IF t < 8 THEN a = 1 * t ELSE IF t < 30 THEN a = .2 ELSE a = -4
v = v + a * dt: x = x + v * dt: t = t + dt
CIRCLE (10*t+10, 240-.5*x), 2 : CIRCLE (10*t+10, 240-10*v), 1
CIRCLE (10 * t + 10, 240 - 10 * a), 1
WEND

Задача 1.2. Точка движется по закону: x=3 cos(t) , y=2 cos(3.7t+2). Вычислите координаты x, y, проекции и модули скорости vx, vy, v и ускорения ax, ay, a, нормальное an и тангенциальное aτ ускорения в последовательные моменты времени t=i Δt, i = 1, 2, 3, ...

 ===== Программа ПР -- 1.2. 
dt = .2: FOR t = 0 TO 10 STEP dt                  'QBASIC
x0 = x: y0 = y: vx0 = vx: vy0 = vy
x = 3 * COS(t): y = 2 * COS(3.7 * t + 2)
vx = (x - x0) / dt: vy = (y - y0) / dt
ax = (vx - vx0) / dt: ay = (vy - vy0) / dt
vv = SQR(vx * vx + vy * vy): aa = SQR(ax * ax + ay * ay)
cosa = (ax * vx + ay * vy) / (aa * vv): at = aa * cosa:
IF at < aa THEN an = SQR(aa * aa - (at * at))
PRINT t; x; y; vv; aa; at; an
NEXT

2. ДИНАМИКА СИСТЕМЫ ЧАСТИЦ

Задача 2.1. Камень бросили под углом α к горизонту со скоростью v0. Сила сопротивления воздуха пропорциональна скорости: F=rv. Определите координаты и скорость камня в последующие моменты времени t. Постройте положения камня через равные промежутки времени.

 ===== Программа ПР -- 2.1. 
SCREEN 12                                         'QBASIC
m = 1: r = .5: g = 10: dt = .002: x = 0: y = 1: vx = 30: vy = 30
LINE (0, 400)-(640, 400): LINE (20, 0)-(20, 480)
WHILE y > 0
t = t + dt: k = k + 1
ax = -r * vx / m: ay = (-m * g - r * vy) / m
vx = vx + ax * dt: x = x + vx * dt
vy = vy + ay * dt: y = y + vy * dt
PSET (20 + 10 * x, 400 - 10 * y), 15
IF k / 50 = INT(k / 50) THEN CIRCLE (20 + 10 * x, 400 - 10 * y), 4
WEND 

Задача 2.2. Тело брошено под углом к горизонту. Постройте графики зависимости координаты y тела и проекций его скорости от времени. Убедитесь в том, что время подъема меньше времени спуска. Сила сопротивления воздуха пропорциональна скорости.

 ===== Программа ПР -- 2.2. 
SCREEN 12                                         'QBASIC
m = 1: r = 2: g = 10: dt = .002: x = 0: y = .01
v = 10: A = 3.14 / 5
vx = v * COS(A): vy = v * SIN(A)
LINE (0, 300)-(640, 300): LINE (20, 0)-(20, 480)
WHILE y > 0
t = t + dt: ax = -r * vx / m: ay = (-m * g - r * vy) / m
vx = vx + ax * dt: x = x + vx * dt
vy = vy + ay * dt: y = y + vy * dt
CIRCLE (20 + 500 * t, 300 - 25 * vx), 1
CIRCLE (20 + 500 * t, 300 - 25 * vy), 1
CIRCLE (20 + 500 * t, 300 - 150 * y), 1
WEND 

3. ДИНАМИКА ТВЕРДОГО ТЕЛА

Задача 3.1. Имеется неоднородная пластина, ограниченная функцией y=x1/2, осью абсцисс, прямой x=2, плотность которой равна ρ(x,y)=2+0,4y+0,2x2. Найдите момент инерции пластины относительно оси Ox.

 ===== Программа ПР -- 3.1. 
dx = .05: dy = .05: h = .02                       'QBASIC
FOR x = 0 TO 2 STEP dy
WHILE y < SQR(x)
y = y + dy: rho = 2 + .4 * y + .2 * x * x
dV = dx * dy * h: I = I + rho * dV * y * y
WEND: NEXT: PRINT I  

Задача 3.2. К шкиву массой m радиуса R, установленному на горизонтальной оси, с помощью невесомого стержня длиной l прикреплен груз массы m1. На шкив намотана нить, к концу которой привязано тело массой m2. При t=0 система покоилась. Рассчитайте угловые ускорение, скорость и перемещение шкива в последующие моменты времени, если тормозящий момент пропорционален скорости.

 ===== Программа ПР -- 3.2. 
SCREEN 12                                         'QBASIC
m1 = .1: m2 = .3: m = .1: R = .05: l = .1: dt = .001
FOR t = 0 TO 8 STEP dt
eps = (9.8 * m2 * R - 9.8 * m1 * l * SIN(fi) - .01 * w) / 
                          (m1 * l * l + (m2 + m / 2) * R * R)
fi = fi + w * dt: w = w + eps * dt
CIRCLE (10 + 100 * t, 400 - 10 * w), 1
CIRCLE (10 + 100 * t, 400 - 5 * fi), 1
CIRCLE (10 + 100 * t, 400 - 1 * eps), 1
NEXT  

4. ЗАКОНЫ СОХРАНЕНИЯ

Задача 4.1. Рассчитайте движение ракеты, удаляющейся по прямой от Земли, если масса горючего уменьшается со скоростью b=const. Масса оболочки равна m_0, масса заправленной ракеты m>m0. Относительная скорость вытекания газов равна vотн. Постройте графики x(t), vx(t), m(t).

 ===== Программа ПР -- 4.1. 
SCREEN 12:  dt = .1: x = 50: m1 = 10              'QBASIC
b = 1: dt = .01: vot = 4: m0 = 1
LINE (0, 420)-(640, 420)
FOR t = 0 TO 100 STEP dt
m1 = m1 - b * dt: m2 = m1 - b * dt
IF m2 < m0 THEN b = 0
Fr = vot * (m1 - m2) / dt: F = 3000 / (x * x)
a = (Fr - F) / m1: v = v + a * dt: x = x + v * dt
CIRCLE (5 + 20 * t, 420 - 2.5 * x), 1
CIRCLE (5 + 20 * t, 420 - 18 * v), 1
CIRCLE (5 + 20 * t, 420 - 20 * m1), 1
NEXT   

Задача 4.2. Промоделируйте движение двух МТ, взаимодействующих друг с другом: а) с силами отталкивания; б) с силами притяжения.

 ===== Программа ПР -- 4.2. 
SCREEN 12:  dt = .1: m1 = 1: m2 = 2               'QBASIC
x1 = 0: y1 = 200: v1x = 2: x2 = 640: y2 = 220: v2x = -2:
FOR t = 1 TO 1000 STEP dt: k = k + 1
ll = (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)
cosA = (x2 - x1) / SQR(ll): sinA = (y2 - y1) / SQR(ll)
F = 500 / ll: a1x = -F * cosA / m1: a1y = -F * sinA / m1
v1x = v1x + a1x * dt: x1 = x1 + v1x * dt
v1y = v1y + a1y * dt: y1 = y1 + v1y * dt
a2x = F * cosA / m2: a2y = F * sinA / m2
v2x = v2x + a2x * dt: x2 = x2 + v2x * dt
v2y = v2y + a2y * dt: y2 = y2 + v2y * dt
PSET (x1, y1): PSET (x2, y2)
IF k / 100 = INT(k / 100) 
    THEN CIRCLE (x1, y1), 1: CIRCLE (x2, y2), 1
IF k / 100 = INT(k / 100) 
    THEN CIRCLE (x1, y1), 2: CIRCLE (x2, y2), 2
IF k / 100 = INT(k / 100) 
    THEN CIRCLE (x1, y1), 3: CIRCLE (x2, y2), 3
IF k / 100 = INT(k / 100) 
    THEN CIRCLE (x1, y1), 4: CIRCLE (x2, y2), 4
NEXT  
5. СИЛЫ В МЕХАНИКЕ

Задача 5.1. Определите силу гравитационного притяжения, действующую со стороны шарообразного тела радиуса R на материальную точку массой m, находящуюся на расстоянии z' от центра. Плотность шара ρ(r)=100/r, где r --- расстояние от его центра.

 ===== Программа ПР -- 5.1. 
CLS: dr = .1: dtheta = .1: dfi = .2               'QBASIC
m = 1: rho1 = 100: z1 = 20
FOR r = 0 TO 3 STEP dr
FOR theta = 0 TO 3.14 STEP dtheta
FOR fi = 0 TO 6.28 STEP dfi
IF r < 1 THEN rho = 100 ELSE rho = 100 / r
z = r * COS(theta): y = r * SIN(theta) * SIN(fi)
x = r * SIN(theta) * COS(fi)
dm = rho * r * r * dtheta * dfi * dr
rast = SQR(x * x + y * y + (z1 - z) * (z1 - z))
l = z1 - r * COS(theta) 
F = F + m * dm * l / (rast * rast * rast)
NEXT: NEXT: PRINT F, r: NEXT   

Задача 5.2. Определите силу гравитационного притяжения, действующую между неоднородным диском радиусом R, толщиной h и однородным стержнем массой m1, концы которого имеют координаты a и b>R. Зависимость плотности диска от координаты: ρ(r)=100/r, где r --- расстояние от его центра O.

 ===== Программа ПР -- 5.2. 
CLS: dr = .1: dfi = .05: rho1 = 100               'QBASIC
dx1 = .1: dm1 = .2 * dx1: h = .1
FOR x1 = 5 TO 6 STEP dx1
FOR r = 0 TO 3 STEP dr
FOR fi = 0 TO 6.28 STEP dfi
IF r < 1 THEN rho = 100 ELSE rho = 100 / r
x = r * COS(fi): y = r * SIN(fi): dm = rho * h * r * dr * dfi
rast = SQR((x - x1) * (x - x1) + y * y): l = x1 - r * COS(fi)
F = F + dm1 * dm * l / (rast * rast * rast)
NEXT: NEXT: PRINT F, r, x1: NEXT 

Задача 5.3. Три МТ массами m1, m2, m3 имеют координаты (x1, y1, z1), (x2, y2, z2), (x3, y3, z3). Рассчитайте напряженность поля во всех точках плоскости, содержащей эти МТ.

 ===== Программа ПР -- 5.3. 
uses crt, graph;                                  'PASCAL
var EC,D,DrV,MV,i,j,m1,m2,m3,x1,y1,x2,y2,x3,y3 : integer;
cosa1,cosa2,cosa3,sina1,sina2,sina3,
        a,a1,a2,r1,r2,r3,x,y,Gx,Gy,pi: real; Label metka;
BEGIN  
DrV:=Detect; InitGraph(DrV,MV,'c:\bp\bgi');
EC:=GraphResult; if EC<>grOK then Halt(1);
pi:=arctan(1)*4; x1:=150; y1:=240; m1:=15; 
x2:=300; y2:=100; m2:=20; x3:=450; y3:=400; m3:=10;
circle(x1,y1,2);circle(x2,y2,2);circle(x3,y3,2);
circle(x1,y1,4);circle(x2,y2,4);circle(x3,y3,4);
circle(x1,y1,6);circle(x2,y2,6);circle(x3,y3,6);
for i:=1 to 100 do for j:=1 to 100 do
begin x:=35*i; y:=35*j;
r1:=sqrt((x-x1)*(x-x1)+(y-y1)*(y-y1));
r2:=sqrt((x-x2)*(x-x2)+(y-y2)*(y-y2));
r3:=sqrt((x-x3)*(x-x3)+(y-y3)*(y-y3));
if (r1=0)or(r2=0)or(r3=0) then goto metka;
cosa1:=(x-x1)/r1; sina1:=(y-y1)/r1; 
cosa2:=(x-x2)/r2; sina2:=(y-y2)/r2;
cosa3:=(x-x3)/r3; sina3:=(y-y3)/r3;
Gx:=m1*cosa1/(r1*r1)+m2*cosa2/(r2*r2)+m3*cosa3/(r3*r3);
Gy:=m1*sina1/(r1*r1)+m2*sina2/(r2*r2)+m3*sina3/(r3*r3);
if Gx<>0 then a:=arctan(Gy/Gx) else a:=pi/2;
line(35*i-round(12*cos(a)),35*j-round(12*sin(a)),
                 35*i+round(12*cos(a)),35*j+round(12*sin(a)));
line(35*i-round(12*cos(a))+1,35*j-round(12*sin(a)),
               35*i+round(12*cos(a))+1,35*j+round(12*sin(a)));
line(35*i-round(12*cos(a)),35*j-round(12*sin(a))+1,
               35*i+round(12*cos(a)),35*j+round(12*sin(a))+1);
line(35*i-round(12*cos(a)),35*j-round(12*sin(a))-1,
               35*i+round(12*cos(a)),35*j+round(12*sin(a))-1);
metka: end; Repeat until Keypressed; CloseGraph;
END.    

Задача 5.4. Рассчитайте силовые линии гравитационного поля, созданного несколькими материальными точками в безграничной среде. Используйте метод покоординатного спуска.

 ===== Программа ПР -- 5.4. 
uses crt, graph;
var EC,D,DrV,MV,i,j,m1,m2,m3,x1,y1,x2,y2,x3,y3,k : integer;
x0,y0,amax,max,fi1,fi,a,r1,r2,r3,x,y,pi: real;
Label metka;
Procedure Raschet;
begin
Repeat r1:=sqrt((x-x1)*(x-x1)+(y-y1)*(y-y1));
r2:=sqrt((x-x2)*(x-x2)+(y-y2)*(y-y2)); 
r3:=sqrt((x-x3)*(x-x3)+(y-y3)*(y-y3));
fi:=-(m1/r1)-(m2/r2)-(m3/r3); max:=0; x0:=x; y0:=y;
for i:=1 to 90 do begin a:=2*pi*i/90; x:=x0+cos(a); y:=y0+sin(a);
r1:=sqrt((x-x1)*(x-x1)+(y-y1)*(y-y1)); 
r2:=sqrt((x-x2)*(x-x2)+(y-y2)*(y-y2));
r3:=sqrt((x-x3)*(x-x3)+(y-y3)*(y-y3)); 
fi1:=-(m1/r1)-(m2/r2)-(m3/r3);
if fi1-fi>max then begin max:=fi1-fi; amax:=a; end;
end; k:=k+1; x:=x0+cos(amax); y:=y0+sin(amax);
circle(round(x),round(y),1);
until k>500;
end;
BEGIN  DrV:=Detect; InitGraph(DrV,MV,'c:\bp\bgi');
       EC:=GraphResult; if EC<>grOK then Halt(1);
pi:=arctan(1)*4; x1:=250; y1:=240; m1:=10;
x2:=300; y2:=100; m2:=20; x3:=400; y3:=300; m3:=10;
circle(x1,y1,2);circle(x2,y2,2);circle(x3,y3,2);
for j:=1 to 20 do begin
x:=250+1*cos(0.314*j); y:=240+1*sin(0.314*j); Raschet; k:=0;
end; Repeat until Keypressed; CloseGraph;
END.    

Задача 5.5. Рассчитайте распределение потенциала гравитационного поля нескольких МТ на плоскости, содержащей эти МТ.

 ===== Программа ПР -- 5.5. 
uses crt, graph;                                  'PASCAL
var r1,r2,r3,x1,x2,x3,y1,y2,y3,fi,m1,m2,m3 :real; 
Gd,Gm,i,j :integer;
BEGIN
Gd:= Detect; InitGraph(Gd, Gm, 'c:\bp\bgi\');
if GraphResult <> grOk then Halt(1);
m1:=15; m2:=23; m3:=8; { m1:=8; m2:=46;}
x1:=180.5; y1:=180.5; x2:=400.5; 
y2:=220.5; x3:=240.5; y3:=340.5;
For i:=0 to 640 do For j:=0 to 480 do
begin 
r1:=sqrt(sqr(x1-i)+sqr(y1-j)); r2:=sqrt(sqr(x2-i)+sqr(y2-j));
r3:=sqrt(sqr(x3-i)+sqr(y3-j)); fi:=m1/r1+m2/r2+m3/r3;
{putpixel(i,j,round((fi+1)*20));} 
if round(fi/0.05)-fi/0.05<0.002 then putpixel(i,j,15);
end; 
Repeat until Keypressed; CloseGraph;
END.    

7. АНАЛИТИЧЕСКАЯ ДИНАМИКА

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

 ===== Программа ПР -- 7.1. 
uses crt, graph;                                  'PASCAL
Const m= 1; k = 4; r = 0.5; dt = 0.002; pi = 3.1415926; 
var i,j,Gd,Gm: integer; 
x,y,z,v,a,f,zz,zzz,t,dx,dy,dz,vx,vy,vz: real;
BEGIN Gd := Detect;  InitGraph(Gd, Gm, 'c:\bp\bgi');
      if GraphResult <> grOk then Halt(1);
line(180,0,180,480);  line(0,280,640,280);
For i:=1 to 40 do For j:=1 to 40 do begin
x:=1+0.04*i; v:=0.04*j; t:=0;
Repeat t:= t + dt; 
a:=(f-k*sin(x)-r*v)/m; v:=v+a*dt; x:=x+v*dt;
until t>5;
putpixel(round(50*x)+180,280-round(50*v),15);
circle(round(50*x)+180,280-round(50*v),1);
end; Repeat until KeyPressed; CloseGraph;
END.     

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

Расползание фазового пятна в случае незатухающих колебаний математического маятника:

Сжатие фазового пятна в случае затухающих колебаний:

Задача 7.2. Камень брошен вертикально вверх со скоростью v=10 м/c так, что в момент t1=0 его координата y1=0, а в момент t2=1 c она равна y2=5 м. Численным методом найдите действие S для действительного пути. Варьируя коэффициенты в уравнении движения y=C1t2/2+C2t+C3, определите действие S' при движении по окольным путям, и убедитесь в том, что оно всегда больше.

 ===== Программа ПР -- 7.2. 
uses crt, graph;                                  'PASCAL
const m=1; dt=0.00005; g=10;
var y,vy,s,ds,ds1,t,C1,c2,c3,c4:real; Gd,Gm,i : integer;
BEGIN  
Gd:= Detect; InitGraph(Gd, Gm, 'c:\bp\bgi');
if GraphResult <> grOk then Halt(1);
line(0,0,640,0); line(0,0,0,480);
For i:=1 to 70 do
begin S:=0; t:=0; C1:=-13+i*0.2; C2:=5-C1/2;
Repeat y:=C2*t+C1*t*t/2; vy:=C2+C1*t;
dS:=m*(vy*vy/2-g*y)*dt; s:=s+ds; t:=t+dt;
until t>1; circle(round(C1*40)+500,round(-S*20),1);
end; Repeat until keypressed;
END.     

Задача 7.3. Электрон движется из точки A области с потенциалом φ1 в точку B другой области с потенциалом φ2. При переходе плоской границы раздела областей происходит "преломление" траектории в точке C. Задавшись скоростью v1 в первой области, рассчитать действие S при различных положениях точки C, если время движения из A в B фиксировано.

 ===== Программа ПР -- 7.3. 
uses crt, graph;                                  'PASCAL
var h2,h1,pi,t,alpha,sinusA,sinusB,x,v1,v2,
                 k,z,q,s,ds,t1,t2,dx: real;
    i: integer; Gd, Gm : integer; x1,x2,y1,y2: integer;
const m=0.1; h=0.1; g=10; a=2;
BEGIN 
Gd:= Detect;  InitGraph(Gd, Gm, 'c:\bp\bgi\');
if GraphResult <> grOk then Halt(1);
pi:=3.1415; h1:=1; v1:=1; alpha:=0.5;
v2:=sqrt(v1*v1+2*g*h); sinusB:=v1*sin(alpha)/v2;
x:=h1*sin(alpha)/cos(alpha);
h2:=(a-x)*sqrt(1-sinusB*sinusB)/sinusB;
t1:=sqrt(h1*h1+sqr(x))/v1; t2:=sqrt(h2*h2+sqr(a-x))/v2;
t:=t1+t2; line(round(x*230+35),150,round(x*230+35),380);
line(10,450,600,450); line(35,50,35,470); x:=0; dx:= 0.01;
For i:=0 to 200 do
begin x:=x+dx; t1:=sqrt(h1*h1+x*x)/v1; t2:=t-t1;
v2:=sqrt(h2*h2+(a-x)*(a-x))/t2;
s:=(m*v1*v1/2-m*g*h)*t1+m*v2*v2/2*t2;
circle(round(x*230+35),round(-S*1500+800),2);
sinusA:=x/sqrt(h1*h1+x*x);
sinusB:=(a-x)/sqrt(h2*h2+(a-x)*(a-x));
circle(round(x*230+35),round(-250*sinusA*v1+450),1);
circle(round(x*230+35),round(-250*sinusB*v2+450),1);
end; Repeat until keypressed; CloseGraph;
END.   

8. КОЛЕБАТЕЛЬНОЕ ДВИЖЕНИЕ

Задача 8.1. Колебательная система состоит из тела массой m, находящегося на пружине с жесткостью k в вязкой среде с коэффициентом сопротивления r. Систему вывели из состояния равновесия и предоставили самой себе. Изучите возникающие колебания.

 ===== Программа ПР -- 8.1. 
SCREEN 12                                         'QBASIC
m = .1: k = 10: r = .1: w = 9.5: dt = .001: x = 15      
WHILE INKEY$ = ""
t = t + dt : F = 0: REM F = 10 * SIN(w * t)
a = (F - k * x - r * v) / m : v = v + a * dt 
x = x + v * dt
CIRCLE (100 * t + 10, 240 - 5 * x), 1 
CIRCLE (100 * t + 10, 240 - 1 * v), 2
CIRCLE (100 * t + 10, 100 - 1 * F), 1 
'CIRCLE (320 + 10 * x, 240 - 1 * v), 1
WEND  

Задача 8.2. На пружинный маятник, состоящий из тела на пружине, действует внешняя гармоническая сила. Получите графики зависимостей x(t), vx(t), ax(t), постройте фазовый портрет системы.

Задача 8.3. Исследуйте автоколебательную систему, которая посредством положительной обратной связи сама регулирует поступление энергии от источника. Пусть вблизи положения равновесия, когда -b < x < b и тело движется в направлении оси OX, на него действует постоянная сила Fm. Постройте графики x(t), vx(t), ax(t), и фазовый портрет системы.

 ===== Программа ПР -- 8.3. 
SCREEN 12                                        'QBASIC
m = 1: k = 1: r = .03: dt = .001: x = 0: v = .01  
LINE (0, 240)-(640, 240): LINE (0, 241)-(640, 241)
WHILE INKEY$ = ""
IF (.3 - ABS(x) > 0) AND (v > 0) THEN F = 3 ELSE F = 0
t = t + dt : a = (F - k * x - r * v) / m 
v = v + a * dt : x = x + v * dt
CIRCLE (15 * t + 10, 240 - 15 * x), 2 
CIRCLE (15 * t + 10, 240 - 40 * v), 1
'CIRCLE (320 + 50 * x, 240 - 40 * v), 1
WEND       

Задача 8.4. Промоделируйте фрикционный маятник Фроуда, состоящий из физического маятника, висящего на вращающемся валу. Сила трения между валом и маятником убывает по экспоненциальному закону Fтр=F0e-β |v - x't|. Постройте график φ(t), ω(t) и фазовую кривую.

 ===== Программа ПР -- 8.4. 
SCREEN 12                                         'QBASIC
m = 1: k = 1: r = .05: dt = .005: vv = 2: x = 0: v = .01
LINE (10, 0)-(10, 480): LINE (0, 240)-(640, 240)
WHILE INKEY$ = ""
IF ((vv - v) > 0) THEN z = 1 ELSE z = -1
F = 1.6 * z * EXP(-.1 * ABS(v - vv))
t = t + dt: a = (F - k * x - r * v) / m
v = v + a * dt: x = x + v * dt
CIRCLE (10 * t + 10, 240 - 20 * x), 2
CIRCLE (10 * t + 10, 240 - 100 * v), 1
CIRCLE (10 * t + 10, 240 - 100 * vv), 1
CIRCLE (180 + 100 * x, 240 - 100 * v), 2
WEND: END 

9. ВОЛНОВОЕ ДВИЖЕНИЕ

Задача 9.1. Промоделируйте распространение одиночного импульса в одномерной упругой среде, его отражение и прохождение границы раздела двух сред.

 ===== Программа ПР -- 9.1. 
uses crt, graph; const n=200; h=1; dt=0.05;       'PASCAL
var  i, j, DriverVar,ModeVar, ErrorCode : integer;
eta,xi,xxii : array[1..N] of real; t, vv : real;
Procedure Graph_Init;
begin {- Инициализация графики -}
DriverVar:=Detect;   
InitGraph(DriverVar,ModeVar,'c:\bp\bgi');
ErrorCode:=GraphResult; if ErrorCode<>grOK then Halt(1);
end;
Procedure Raschet; {Расчет смещения}
  begin for i:=2 to N-1 do begin xxii[i]:=xi[i];
  if i<N/2 then vv:=4 else vv:=0.5;
  eta[i]:=eta[i]+vv*(xi[i+1]-2*xi[i]+xi[i-1])/(h*h)*dt;  
  end;
  for i:=2 to N-1 do xi[i]:=xi[i]+eta[i]*dt;
  xi[N]:=0;      { Конец закреплен}
{ xi[N]:=xi[N-1];}{ незакрепленный}  end;
Procedure Draw; 
  begin        {- Вывод на экран -}
  setcolor(black); line(i*3-3,240-round(xxii[i-1]*50),
                           i*3,240-round(xxii[i]*50));
    setcolor(white); line(i*3-3,240-round(xi[i-1]*50),
                       i*3,240-round(xi[i]*50));  end;
BEGIN      {- Основная программа -}
  Graph_Init;  Repeat t:=t+dt;
    if t<6.28 then xi[1]:=2*sin(t) else xi[1]:=0;
    Raschet; For i:=1 to N do Draw;
  until KeyPressed; CloseGraph;
END.  

Задача 9.2. Имеется прямоугольная пластина, некоторые элементы которой совершает колебания, а другие --- не способны колебаться. Края пластины закреплены. Изучите распространение, отражение, интерференцию и дифракцию упругой волны в двумерной среде.

Задача 9.3. Рассчитайте интерференционную картину от двух источников, удаленных друг от друга на расстоянии 2d, в плоскости, содержащей источники. Интенсивности волн равны I, длина волны λ.

 ===== Программа ПР -- 9.2. 
SCREEN 12: pi = 3.1415: Lambda = 60: d = 60       'QBASIC
x1 = 320 - d: x2 = 320 + d: y1 = 340
FOR x = 1 TO 640 STEP 2
FOR y = 1 TO 480 STEP 2
l1 = SQR((x - x1) * (x - x1) + (y - y1) * (y - y1))
l2 = SQR((x - x2) * (x - x2) + (y - y1) * (y - y1))
I1 = 5000 / (l1 * l1): I2 = 5000 / (l2 * l2)
I = I1 + I2 + 2 * SQR(I1 * I2) * COS(2 * pi * (l1 - l2) / Lambda)
IF I > .2 THEN PSET (x, y)
IF I > 2 THEN CIRCLE (x, y), 1
NEXT
NEXT 


10. МЕХАНИКА ЖИДКОСТИ И ГАЗА

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

 ===== Программа ПР -- 10.1. 
uses crt, graph; const n=140; m=50;               'PASCAL
var i,ii,j,jj,k,DriverVar, ModeVar, ErrorCode : integer;
psi: array[1..N, 1..M] of real;
procedure Raschet; {---- Расчет потенциала ----}
begin 
psi[i,j]:=(psi[i+1,j]+psi[i-1,j]+psi[i,j+1]+psi[i,j-1])/4;
end;
procedure Gran_usl; {---- Граничные условия ----}
begin for i:=1 to N do
begin psi[i,2]:=-200; psi[i,M-1]:=200; 
      psi[i,1]:=-200; psi[i,M]:=200;
end;
for j:=1 to M do begin psi[N-1,j]:=-204+8*j; 
                         psi[N,j]:=-204+8*j; end;
for j:=1 to 25 do begin
psi[2,j]:=-204+16*j; psi[N-1,j]:=-204+8*j;
psi[1,j]:=-204+16*j; psi[N,j]:=-204+8*j; end;
for i:=1 to N do for j:=1 to M do
if (j>25)and((i-60)*(i-60)+(j-25)*(j-25)<200) 
                                then psi[i,j]:=0;
for i:=1 to N do for j:=1 to M do
if (abs(j-25)<15)and(abs(i-110)<5) then psi[i,j]:=0;
for i:=1 to N do for j:=1 to M do
if (j>25)and(i<30) then psi[i,j]:=200;
end;
procedure Draw; {---- Вывод на экран ----}
begin
if psi[i,j]>-200 then setcolor(14);
if psi[i,j]>-150 then setcolor(0);
if psi[i,j]>-100 then setcolor(2);  
if psi[i,j]>-50 then setcolor(4);
if psi[i,j]>0  then setcolor(6);    
if psi[i,j]>50  then setcolor(8);
if psi[i,j]>100  then setcolor(10); 
if psi[i,j]>150 then setcolor(12);
{setcolor(round((psi[i,j]+200)/35));}
{if round(psi[i,j]/75)25)and((i-60)*(i-60)+(j-25)*(j-25)<200) 
                then setcolor(0);
if (j>25)and(i<30) then setcolor(0);
if (abs(j-25)<15)and(abs(i-110)<5) then setcolor(0);
rectangle(i*4+50,j*4,i*4+53,j*4+3);
rectangle(i*4+51,j*4+1,i*4+52,j*4+2);
end;
BEGIN             {---- Основная программа ----}
begin DriverVar:=Detect; 
InitGraph(DriverVar,ModeVar,'c:\bp\bgi');
ErrorCode:=GraphResult; 
if ErrorCode <> grOK then Halt(1); end;
Repeat inc(k);
Gran_usl; For i:=2 to N-1 do For j:=2 to M-1 do Raschet;
Gran_usl; For j:=2 to M-1 do For i:=2 to N-1 do Raschet;
Gran_usl; For i:=2 to N-1 do For jj:=2 to M-1 do
                          begin j:=M+1-jj; Raschet; end;
Gran_usl; For j:=2 to M-1 do For ii:=2 to N-1 do
                          begin i:=N+1-ii; Raschet; end;
if k/50=round(k/50) then begin cleardevice;
For i:=2 to N-1 do For j:=2 to M-1 do Draw; end;
until KeyPressed; CloseGraph;
END.

Задача 10.2. Исследуйте установившееся безвихревое течение вязкой жидкости в бесконечно длинной трубе (канале) постоянного сечения, внутри которой движется бесконечно длинное тело.

 ===== Программа ПР -- 10.2. 
uses crt, graph;                                  'PASCAL
const n= 90; m=58; h=1; dt=0.02;                  
var  i,j,ii,jj,kk,DV,MV,EC : integer;
vv, v: array[1..N, 1..M] of real; 
procedure Gr_usl; {---- Граничные условия ----}
begin
for i:=1 to N do for j:=1 to M do v[i,j]:=vv[i,j];
for i:=2 to N-1 do for j:=2 to M-1 do begin
if (j<5)and(abs(5-i)<5) then v[i,j]:=0;
if (j<40)and(abs(i-50)<1+20*sqr(cos(j/25))) then v[i,j]:=300;
end; for i:=1 to N do v[i,M]:=0;
for j:=1 to M do begin v[1,j]:=0; v[N,j]:=0; end;
end;                                    
procedure Raschet; {---- Расчет скорости ----}
begin vv[i,j]:=v[i,j]+(v[i,j+1]-2*v[i,j]+v[i,j-1])*dt/(h*h)
+(v[i+1,j]-2*v[i,j]+v[i-1,j])*dt/(h*h); end;
procedure Draw;{---- Вывод на экран ----}
begin
setcolor(round(v[i,j]/40));
if v[i,j]<1 then setcolor(9);
{if round(v[i,j]/60) grOK then Halt(1);
Repeat Gr_usl; 
       For i:=1 to N-1 do For j:=2 to M-1 do Raschet;
       Gr_usl; For jj:=2 to M-1 do For ii:=1 to N-1 do
               begin i:=N+1-ii; j:=M+1-jj; Raschet; end;
       For i:=1 to N do vv[i,1]:=vv[i,2];
       kk:=kk+1; if kk/30=round(kk/30) then
       For i:=2 to N-1 do For j:=2 to M-1 do Draw;
until KeyPressed; CloseGraph;
END.          

Задача 10.3. Рассчитайте течение вязкой жидкости в бесконечно длинной трубе произвольного сечения при наличии разности давлений на концах трубы.


11. ОСНОВЫ НЕРЕГУЛЯРНОЙ ДИНАМИКИ

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

 ===== Программа ПР -- 11.1. 
uses dos, crt, graph;                             'PASCAL
const dt=0.03;
var m,Fx,Fy,x,y,vx,vy,xx,yy,x1,y1 : real;
    Gd,Gm,i,j :integer; ax,ay,F,l : real;
BEGIN
Gd:=Detect; InitGraph(Gd, Gm, 'c:\bp\bgi');
if GraphResult <> grOk then Halt(1);
m:=0.001; x1:=150; y1:=150; x:=100; 
y:=100; vx:=30; vy:=10;
Repeat
 Fx:=0; Fy:=0; l:=sqrt(sqr(x1-x)+sqr(y1-y)); 
 if l<10 then F:=-0.03;
 if x<80 then Fx:=0.05; if x>220 then Fx:=-0.05;
 if y<80 then Fy:=0.05; if y>220 then Fy:=-0.05;
 Fx:=Fx+F*(x1-x)/l; Fy:=Fy+F*(y1-y)/l; ax:=Fx/m; ay:=Fy/m;
 xx:=x; yy:=y; vx:=vx+ax*dt; vy:=vy+ay*dt;  
 x:=x+vx*dt; y:=y+vy*dt; 
 delay(50); circle(50+round(x),240-round(vx),1);
{setcolor(8);  circle(round(xx),round(yy),2);}
 setcolor(15); circle(round(x),round(y),2);
until KeyPressed; 
Repeat until keypressed; CloseGraph;
END.                

Задача 11.2. Логистическое отображение задается рекурсивной процедурой: xi+1=axi(1-xi). Изменяя бифуркационный параметр a, исследуйте процесс перехода в хаотический режим. Как ведет себя показатель устойчивости Ляпунова?

 ===== Программа ПР -- 11.2. 
uses crt, graph;                                  'PASCAL
const N=500; dt=0.001;
var ErrorCode,Detect,DriverVar,ModeVar, i, j : integer;
l,a,x,e,y,z,delta: real;
BEGIN
  DriverVar:=Detect; InitGraph(DriverVar,ModeVar,'c:\bp\bgi');
  ErrorCode:=GraphResult; if ErrorCode<>grOK then Halt(1);
  line(0,460,640,460); line(0,50,640,50); a:=2; e:=0.0001;
Repeat x:=0.05; y:=x+e; z:=x-e; a:=a+0.001; delta:=0;
for i:=1 to 600 do if abs(a*x*(1-x))<10 then begin
x:=a*x*(1-x); y:=a*y*(1-y); z:=a*z*(1-z);
if i>300 then delta:=delta+abs(y-x+z-x);
if i>100 then putpixel(round(300*a)-600,460-round(300*x),15); 
end;
l:=ln(0.0000000000000000000001+(delta/2/e/300));
circle(round(300*a)-600,50-round(2*l),1);
until (a>4) or (Keypressed);
Repeat until Keypressed; CloseGraph;
END.   

Задача 11.3. Промоделируйте колебания маятника в потенциальной яме, у которой на месте одного постепенно появляется два углубления, разделенные плавно увеличивающимся потенциальным барьером. Изучите бифуркацию типа "вил".

 ===== Программа ПР -- 11.3. 
uses crt, graph;                                  'PASCAL
const pi = 3.1415926; dt = 0.001; m= 1; 
      k = 4; r = 0.9; w = 2.3;
var i,j,Gd,Gm: integer; b,f,fi,t,z,zz,x,v,a: real;
BEGIN Randomize;
Gd := Detect; InitGraph(Gd, Gm, 'c:\bp\bgi');
if GraphResult <> grOk then Halt(1);
for j:=0 to 1000 do begin
b:=-10+0.025*j; x:= 0; v:= 0; t:=0; fi:=random(10);
for i:=0 to 10000 do begin t:= t + dt; f:=2*sin(w*t+fi);
a:=(f-k*(x*x*x-b*x)-r*v)/m; v:=v+a*dt; x:=x+v*dt; z:=SIN(w*t);
IF (z>0)and(zz<0)and(t>1) THEN 
                 circle(round(20*b)+200, 240-round(50*x),1);
IF (z>0)and(zz<0)and(t>1) THEN 
                 circle(round(20*b)+200, 240-round(50*x),2);
zz:= z;
end; end; Repeat until KeyPressed; CloseGraph;
END.   

Задача 11.4. Промоделируйте перемешивание фазового объема в случае свободных незатухающих колебаний маятника Дафинга.

 ===== Программа ПР -- 11.4. 
uses crt, graph;                                  'PASCAL
Const m= 1; k = 4; dt = 0.002; pi = 3.1415926;
var i,j, Gd, Gm: integer; 
v,a,f,zz,zzz,t, x,y,z,dx,dy,dz,vx,vy,vz: real;
BEGIN Gd := Detect;  InitGraph(Gd, Gm, 'c:\bp\bgi');
if GraphResult <> grOk then Halt(1);
line(200,0,200,480); line(0,280,640,280);
For i:=1 to 40 do For j:=1 to 40 do begin
x:=0.04*i; v:=0.04*j; t:=0;
Repeat t:=t+dt; a:=-(x*x*x-x)/m; v:=v+a*dt; x:=x+v*dt;
until t>2; circle(round(50*x)+200,280-round(50*v),1);
end; Repeat until KeyPressed; CloseGraph;
END.         

Задача 11.5. Получите сечения Пуанкаре для маятника Дафинга.

 ===== Программа ПР -- 11.5. 
uses crt, graph;                                  'PASCAL
const m = 1; k = 4; r = 0.5; w = 2.3;
      dt = 0.002; pi = 3.1415926;
var  Gd, Gm: integer; f,t,x,v,a,z,zz: real;
BEGIN x:= 0; v:= 0;
  Gd := Detect; InitGraph(Gd, Gm, 'c:\bp\bgi');
  if GraphResult <> grOk then Halt(1);
  line(200,0,200,480); line(0,280,640,280);
Repeat t:= t + dt; f:=1*sin(w*t);
a:=(f-k*(x*x*x-x)-r*v)/m; v:=v+a*dt; 
x:=x+v*dt; z:=SIN(w*t+3*pi/4);
IF (z>0)and(zz < 0) 
  THEN circle(round(100*x) + 200, 280 - round(100*v),1);
zz:= z; until (t>10000);
Repeat until KeyPressed; CloseGraph;
END.    

Задача 11.6. Исследуйте странный аттрактор Реслера, задаваемый системой уравнений:

x' = - (y + z), y' = x + y/5, z' = 1/5 + z( x - μ).

 ===== Программа ПР -- 11.6. 
uses crt, Graph;                                 'PASCAL
Const a=15; b=30; c=2.2; dt=0.001; m=4.9;
var Gd,Gm: integer; t,x,y,z,dx,dy,dz,vx,vy,vz: real;
BEGIN Gd := Detect; InitGraph(Gd, Gm, 'c:\bp\bgi');
if GraphResult <> grOk then Halt(1); x:=2; y:=2; z:=1;
line(300,0,300,480); line(0,220,640,220);
Repeat t:=t+dt; dx:=-(y+z)*dt; x:=x+dx; vx:=dx/dt;
dy:=(x+y/5)*dt; y:=y+dy; vy:=dy/dt;
dz:=(0.2+z*(x-m))*dt; z:=z+dz; vz:=dz/dt;
if t>20 then putpixel(300+round(14*x),220-round(13*y),15);
until KeyPressed;
CloseGraph;
END.  

Задача 11.7. Промоделируйте странный аттрактор Лоренца, уравнения которого имеют вид:

x' = a ( - x + y), y' = bx - y - xz, z' = - cz + xy.

 ===== Программа ПР -- 11.7. 
uses crt, graph;                                  'PASCAL
const a=15; b=30; c=2.2; dt=0.0001;                    
var  Gd, Gm: integer;
zz,zzz,t, x,y,z,dx,dy,dz,vx,vy,vz: real;
BEGIN  Gd:=Detect; InitGraph(Gd,Gm,'c:\bp\bgi');
       if GraphResult <> grOk then Halt(1);
  x:=20; y:=20; z:=5;
  Repeat t:=t+dt; dx:=a*(-x+y)*dt; x:=x+dx; vx:=dx/dt;
  dy:=(b*x-y-x*z)*dt; y:=y+dy; vy:=dy/dt;
  dz:=(-c*z+x*y)*dt; z:=z+dz;  vz:=dz/dt;
  putpixel(220+round(10*x),220-round(1*vx),15);
{ putpixel(220+round(10*y),220-round(1*vy),15);}
{ putpixel(100+round(10*z),220-round(1*vz),15);}
until KeyPressed; CloseGraph;
END.  

12. НЕКОТОРЫЕ ЗАДАЧИ МЕХАНИКИ

Задача 12.1. Рассчитайте траекторию тела в гравитационном поле Земли при различных начальных скоростях движения.

 ===== Программа ПР -- 12.1. 
SCREEN 12                                         'QBASIC
m1 = 10: m2 = 1000: dt = .0005:
LINE (0, 240)-(640, 240): LINE (320, 0)-(320, 480)
CIRCLE (320, 240), 100: CIRCLE (320, 240), 102
FOR i = 0 TO 5
x = 0: y = 20: vy = 0: vx = 2 * i + 1
FOR t = 0 TO 10 STEP dt
r = SQR(x * x + y * y): F = m1 * m2 / (r * r)
Fx = -F * x / r: Fy = -F * y / r: ax = Fx / m1: ay = Fy / m1
vx = vx + ax * dt: x = x + vx * dt
vy = vy + ay * dt: y = y + vy * dt
CIRCLE (320 + 10 * x, 240 - 10 * y), 1
NEXT: NEXT  

Задача 12.2. Промоделируйте движение тела в центрально симметричном поле, напряженность которого обратно пропорциональна r2/3. Убедитесь в том, что траекторией является незамкнутая кривая.

 ===== Программа ПР -- 12.2. 
SCREEN 12                                         'QBASIC
m1 = 10: m2 = 1000: dt = .0002 
x = 0: y = 10: vx = 10: vy = 0 
LINE (0, 240)-(640, 240): LINE (320, 0)-(320, 480)
WHILE t < 50: t = t + dt
r = SQR(x * x + y * y): F = m1 * m2 / SQR(r * r * r)
Fx = -F * x / r: Fy = -F * y / r: ax = Fx / m1: ay = Fy / m1
vx = vx + ax * dt: x = x + vx * dt
vy = vy + ay * dt: y = y + vy * dt
CIRCLE (320 + 15 * x, 240 - 15 * y), 1
WEND       

Задача 12.3. Планета вращается вокруг Солнца по эллиптической орбите. Постройте графики зависимости расстояния, линейной скорости и секториальной скорости планеты от времени. Подтвердите второй закон Кеплера, согласно которому секториальная скорость планеты остается постоянной.

 ===== Программа ПР -- 12.3. 
SCREEN 12: m1 = 10: m2 = 1000: dt = .002          'QBASIC
LINE (0, 400)-(640, 400): LINE (20, 0)-(20, 480)
x = 0: y = 20: vy = 0: vx = 4: r1 = SQR(x * x + y * y)
FOR t = 0 TO 30 STEP dt
r = SQR(x * x + y * y): F = m1 * m2 / (r * r)
Fx = -F * x / r: Fy = -F * y / r
ax = Fx / m1: ay = Fy / m1
vx = vx + ax * dt: x = x + vx * dt
vy = vy + ay * dt: y = y + vy * dt
dl = SQR((x - x1) * (x - x1) + (y - y1) * 
                      (y - y1) - (r - r1) * (r - r1))
dS = dl * r / 2: w = dS / (10 * dt)
x1 = x: y1 = y: r1 = r: v = SQR(vx * vx + vy * vy)
CIRCLE (20 + 20 * t, 400 - 60 * w), 1
CIRCLE (20 + 20 * t, 400 - 10 * r), 2
CIRCLE (20 + 20 * t, 400 - 10 * v), 1
NEXT   

Задача 12.4. По направлению к массивному положительно заряженному ядру движутся альфа--частицы (опыт Резерфорда). Постройте траектории движения частиц при различных значениях прицельного параметра ρ.

 ===== Программа ПР -- 12.4. 
SCREEN 12: m1 = 1: m2 = 200: dt = .001            'QBASIC
LINE (0, 400)-(640, 400): LINE (320, 0)-(320, 480)
FOR i = 0 TO 30 STEP 3
x = -80: y = i: vx = 5: vy = 0
FOR t = 0 TO 30 STEP dt
r = SQR(x * x + y * y): F = -m1 * m2 / (r * r)
Fx = -F * x / r: Fy = -F * y / r: ax = Fx / m1: ay = Fy / m1
vx = vx + ax * dt: x = x + vx * dt
vy = vy + ay * dt: y = y + vy * dt
CIRCLE (320 + 5 * x, 400 - 5 * y), 1
NEXT: NEXT   

Задача 12.5. Изучите затухающие колебания двух одинаковых маятников, каждый из которых состоит из тела массой m и пружины жесткостью k, связанных между собой упругой связью жесткостью q.

 ===== Программа ПР -- 12.5. 
SCREEN 12                                         'QBASIC
dt = .001: r = .2: k = 100: q = 8
t = 0: v = 15: x1 = .4: m = 1
DO
F = -q * (x1 - x2): a1 = (-F - r * v1 - k * x1) / m
x1 = x1 + v1 * dt
v1 = v1 + a1 * dt: a2 = (F - r * v2 - k * x2) / m
x2 = x2 + v2 * dt: 
v2 = v2 + a2 * dt: t = t + dt
CIRCLE (10 + 20 * t, 150 - 100 * x1), 1
CIRCLE (10 + 20 * t, 350 - 100 * x2), 1
LOOP UNTIL t > 40
END 

13. ТЕОРИЯ ОТНОСИТЕЛЬНОСТИ

Задача 13.1. В точке с координатами (x1, 0, 0) в момент времени t (инерциальная система отсчета ИСО K1) произошло событие A. Определите пространственную и временную координаты события A в ИСО K2, движущейся со скоростью v=0,6c в направлении оси Ox. Вычислите и сравните пространственно--временные интервалы между событиями O и A в ИСО K и K'.

 ===== Программа ПР -- 13.1. 
CLS                                               'QBASIC
x1 = 2: x0 = 1: beta = .6: gamma = 1 / SQR(1 - beta * beta)
x11 = (x1 + beta * x0) * gamma: x01 = (x0 + beta * x1) * gamma
ss = SQR(x1 * x1 - x0 * x0): ss1 = SQR(x11 * x11 - x01 * x01)
PRINT x1, x0, x11, x01: PRINT "Интервал", ss, ss1
END  


Майер Р.В. Математические начала современной теории механического движения: Учебное пособие. --- Глазов: ГГПИ, 2007. --- 164 с.


ВВЕРХ