Майер роберт Валерьевич основы компьютерного моделирования: Учебное пособие. Глазов: ггпи, 2005. 25 с

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

---

Кроме нее на осциллятор действует сила упругости 
и сила вязкого 
трения 
, где -- скорость -ого осциллятора. По второму закону 
Ньютона, ускорение -ого осциллятора равно:
Определив ускорение , можно найти скорость и смещение -ой массы в 
следующий момент времени 
: 
3. Алгоритм.
Эта задача является частным случаем задачи о моделировании 
движения системы взаимодействующих между собой частиц. 
1. Задают параметры системы и начальные условия: число осцилляторов 
, 
их массы 
, жесткость упругих связей , жесткость пружины внутри 
осциллятора , шаг по времени 
начальные смещение 
и скорости 
, 
вынуждающие силы, действующие на различные материальные точки. 
2. Дают 
приращение 
по 
времени: 
переменной присваивают 
значение 
3. В цикле перебирают все 
осцилляторов и для каждого из них 
опеределяют равнодействующую силу, ускорение, скорость и смещение в 
следующий момент 

---

4. В цикле стирают предыдущее изображение каждого осциллятора и рисуют 
новое. 
5. Возвращение к операции 2. Если цикл по закончился, -- выход из цикла. 
4. Компьютерная 
программа.
Программа, 
представленная 
ниже, 
моделирует распространение импульса в случае, когда левый элемент среды 
совершает полколебания. 
program PROGRAMMA5; 
uses dos, crt, graph; 
Const m=0.5; r=0.1; k=0.01;
dt=0.001; q=100; N=50; 
Var teta, F, t, y : Real; 
i,xx,vv,aa,FF,tt,Gd,Gm : Integer; 
eta, ksi : array [0..N] of real; 
Procedure Graph_Init; 
begin 
Gd:= Detect; 
InitGraph(Gd, Gm, 'c:\bp\bgi'); 
if GraphResult <> grOk then Halt(1); 
end; 
Procedure Oscillator; 
begin 
F:=q*(ksi[i-1]-ksi[i])+q*(ksi[i+1]-ksi[i]); 
teta:=(F-r*eta[i]-k*ksi[i])/m; 
eta[i]:=eta[i]+teta*dt; ksi[i]:=ksi[i]+eta[i]*dt; 
end; 
BEGIN 
Graph_Init; 
Repeat begin t:=t+dt; 
For i:=1 to N do begin 
y:=ksi[i]; Oscillator; 
if 5*t<3.142 then ksi[1]:=20*sin(5*t) else ksi[1]:=0; 
ksi[N]:=0; {закреплен} 
{ksi[N]:=ksi[N-1]; незакреплен} 
setcolor(8); circle(10*i,240-round(y*10),2); 
setcolor(15); circle(10*i,240-round(ksi[i]*10),2); 
end; end; 
until KeyPressed; 
CloseGraph; 
END. 

---