Randomize timer

Задача 5. Проверка равномерного распределения энергии по термодинамическим степеням свободы.

'программа моделирования равномерного распределения по степеням свободы

'одномерная модель: на прямой маятник и частица

'левая граница случайно перемещается

pi = ATN(1) * 4

m1 = 1: k = 8 'масса частицы и коэффициент жесткости пружины маятника

INPUT "m2/m1= ", m2 'отношение масс

L = 1 'длина всего отрезка

x1 = 0: x2 =.8 'начальные координаты

v1 = SQR((2 - k * (L - x2) ^ 2) / m1) 'начальная скорость при E=1

v2 = 0: a2 = k * (L - x2) / m2 'ускорение частицы маятника

v0 = m1 * v1 / (m1 + m2) 'скорость центра масс

t = 2 * pi * SQR(m2 / k) 'период колебаний маятника, опорное время

nstep = 1000: dt = t / nstep 'число шагов и время шага

ni = nstep * 2000 'полное число шагов

E0 = 1 'начальная энергия

gr = 0 'левая граница

FOR i = 1 TO ni 'основной цикл движения

'сдвиг свободной частицы

x1 = x1 + v1 * dt

'флуктуация левой стенки и вкл. ключа после ее сдвига

IF x1 >.1 * L AND kl = 0 THEN gr =.05 * L * (2 * RND - 1): kl = 1

'отражение от левой границы

IF x1 < gr THEN x1 = 2 * gr - x1: v1 = ABS(v1): kl = 0

' предварительный сдвиг частицы маятника

dx2 = v2 * dt + a2 * dt ^ 2 / 2

'предсказание ускорения после сдвига

a22 = k * (L - x2 - dx2) / m2

'основной сдвиг частицы маятника

'во втором приближении с поправкой ускорения

x2 = x2 + v2 * dt + (a2 + a22) / 2 * dt ^ 2 / 2

'проверка удара, изменение скоростей и счетчик ударов

IF x1 > x2 THEN x1 = x2: v1 = 2 * v0 - v1: v2 = 2 * v0 - v2: n = n + 1

'половинка изменение скорости, определение ускорения, еще половинка

v2 = v2 + a2 * dt / 2: a2 = k * (L - x2) / m2: v2 = v2 + a2 * dt / 2

'вычисление энергии

E = m1 * v1 ^ 2 / 2 + m2 * v2 ^ 2 / 2 + k * (L - x2) ^ 2 / 2

'коэффициент изменения энергии из-за неточности алгоритма

kE = SQR(E0 / E)

'правка скоростей

v1 = v1 * kE: v2 = v2 * kE

'скорость центра масс

1 v0 = (m1 * v1 + m2 * v2) / (m1 + m2)

вычисление средних величин

x1cp = x1cp + x1 / ni: x2cp = x2cp + x2 / ni

x12cp = x12cp + x1 ^ 2 / ni: x22cp = x22cp + x2 ^ 2 / ni

E1cp = E1cp + m1 * v1 ^ 2 / 2 / ni

E2cp = E2cp + (m2 * v2 ^ 2 / 2 + k * (L - x2) ^ 2 / 2) / ni