Интегрирование системы дифференциальных уравнений

Вдоль горизонтальной прямой движутся n материальных точек с массами \(m_k\) (\(k=1,\ldots,n\)), связанные пружинами с заданными жесткостями \(c_k\) и свободными длинами \(L_k\).

Дифференциальное уравнение движение первого тела: \[m_1 \ddot{x}_1 = - c_1 (x_{1}-L_1) + c_2 (x_{2}-x_{1}-L_{2})\]

Дифференциальное уравнение движение \(k\)-го тела (\(2 \leq k \leq n-1\)): \[m_k \ddot{x}_k = - c_k (x_{k}-x_{k-1}-L_k) + c_{k+1} (x_{k+1}-x_{k}-L_{k+1}), \quad k=2,\ldots,n-1\]

Дифференциальное уравнение движение \(n\)-го тела: \[m_n \ddot{x}_n = - c_n (x_{n}-x_{n-1}-L_n)\]

Файл-функция правых частей дифференциальных уравнений движения системы материальных точек

function dq = dqdt(t, q, p)
    % Матрица-столбец состояния механической системы имеет
    % размерность (2n x 1):
    % Количество тел (материальных точек)
    n = size(q,1)/2;
    % Первые n элементов столбца это координаты грузов:
    x  = q(1:n);
    % элементы с n+1 по 2n -- скорости тел:
    v  = q(n+1:2*n);
    % Столбец расстояний между соседними точками (для первой точки это сама координата x1)
    dx = [x(1);diff(x)];
    % Столбец сил растяжения пружин
    Fs = (dx-p.L).*p.c;
    % Столбец сил, действующих на точки
    F  = -Fs + [Fs(2:n); 0];
    % Столбец ускорений тел
    a    = F./p.m;
    dq = [v;a];
end

Файл-скрипт, запускающий процесс численного интегрирования

% Массы тел
p.m = [5;10;8];
% Свободные длины пружин
p.L = [1;1;1];
% Жесткости пружин
p.c = [1000;1000;500];
% Координатный столбец начальных условий
q0 = [1.1;2.5;3.2;0;0;0];
% Запускаем процесс интегрирования
[t, q] = ode45(@(t,q) dqdt(t,q,p),[0 5], q0);
% Графики изменений координат тел
plot(t,q(:,1:3));

Задание

1. Добавить в модель силу сухого трения, действующую на тела, полагая, что система движется по горизонтальной плоскости в поле силы тяжести (\(g = 9.81\) м/с\(^2\)).
2. Построить графики изменения расстояний между телами от времени.
3. Построить график изменения кинетической энергии системы от времени.
4. Построить график изменения потенциальной энергии системы от времени.

Использовать следующую модель зависимости коэффициента трения f от скорости (модель Штрикбека - Кулона):

Vbr = 0.002;
Vst = Vbr*sqrt(2);
Vc  = Vbr/10;
Fbr = 0.2;
Fc  = 0.15;
f   = @(v) sqrt(2)*exp(1)*(Fbr-Fc)*exp(-(v/Vst).^2).*v/Vst + Fc*tanh(v/Vc);