Управление процессом численного интегрирования (триггеры/события)
Рассмотрим движение материальной точки, положение которой определяется двумя координатами \(x\) и \(y\) в вертикальной плоскости в ограниченном пространстве: \[x_{min} \leq x \leq x_{max}, \quad y_{min} \leq y \leq y_{max}.\]
Предположим, что точка движется под действием силы тяжести и ее уравнения движения имеют вид: \[\left\{ \begin{aligned} \dot x & = & V_x \\ \dot y & = & V_y \\ \dot{V}_x & = & 0 \\ \dot{V}_y & = & -g \end{aligned} \right.\]
Пусть при достижении правой или левой границы происходит отскок материальной точки от стенки (границы) с потерей энергии, так что: \[V_x^+ = - k V_x^-, \quad V_y^+ = \cdot V_y^-,\]
где \(V_x^-\), \(V_x^+\) – проекции скорости точки на ось x до и после удара соответственно, \(V_y^-\), \(V_y^+\) – проекции скорости точки на ось y до и после удара. \(0 < k \leq 1\) – коэффициент восстановления.
Также, при ударе точки о верхнюю \(y = y_{max}, \, V_y>0\) или нижнюю границы \(y = y_{min}, \, V_y<0\): \[V_x^+ = V_x^-, \quad V_y^+ = -k \cdot V_y^-,\]
При моделировании движения такой механической системы в MATLAB необходимо останавливать процесс численного интегрирования при достижении любой из границ и запускать следующий процесс численного интегрирования с новыми начальными условиями (скоростями точки).
Создадим файл-скрипт и объявим в нем параметры системы при помощи структуры p, функцию правых частей дифференциальных уравнений pool_dqdt и начальные условия движения материальной точки:
clc; clear all;
p.xmin = -2;
p.xmax = +5;
p.ymin = -2;
p.ymax = +2;
p.g = 9.807;
p.vmin = 0.1;
% Функция правых частей
pool_dqdt = @(t,q,p) [q(3:4);0;-p.g];
% Начальные условия
x0 = 0;
y0 = 0;
vx0 = 10*cos(45*pi/180);
vy0 = 10*sin(45*pi/180);
q0 = [x0;y0;vx0;vy0];
Создадим еще один файл – файл-функцию pool_event_function, аргументами которой является время t, соответствующий этому времени координатный столбец вектора состояния системы \([x,y,V_x,V_y]\) и структура с параметрами системы.
В функции определяется вектор функция value, состоящая из 4 элементов. При пересечении границ поля точкой один элементов этого вектора становится равным нулю. Переменная direction определяет направление изменения компонентов вектор-функции, объявленной в переменной value, которые контролируется. Например, если значение x-p.xmin (первый элемент массива value) пересечёт ноль при убывании x direction(1)=-1, то процесс интегрирования необходимо остановить isterminal(1)=1. Вторые элементы массивов value, direction, isterminal отвечают за контроль условия \[x(t) - x_{max} = 0, \quad \dot x(t) > 0\]
Последний пятый элемент “следит” за уменьшением скорости движения материальной точки. Если эта скорость падает (direction(5)=-1) до значений меньше (по модулю) p.vmin, то процесс интегрирования также останавливается.
function [value,isterminal,direction] = pool_event_function(t, q, p)
x = q(1);
y = q(2);
% Скорость точки (модуль)
v = norm(q(3:4));
% Контрольные функции
value = [x-p.xmin x-p.xmax y-p.ymin y-p.ymax v-p.vmin];
% Направление изменений функций при пересечении нуля
direction = [-1 +1 -1 +1 -1];
% Остановить процесс интегрирования [1 - да, 0 - нет],
isterminal = [1 1 1 1 1];
end
Вернемся в главный файл-скрипт. При помощи функции odeset создадим структуру параметров интегрирования opt, указав, что в процессе интегрирования необходимо контролировать события, определенные в функции pool_event_function:
opt = odeset('Events', @(t,q) pool_event_function(t,q,p));
Далее зададим интервал интегрирования и объявим пока пустые массивы для хранения результатов интегрирования:
t0 = 0 ; dt = 20;
t = []; q = [];
Интегрирование производится интервалами по dt = 20 секунд (от t0 до t0+dt). Первый запуск процесса численного интегрирования выполняется для интервала от 0 до 20 секунд с начальными условиями q0. При “срабатывании” любого из пяти событий, определенных в функции события pool_event_function, процесс интегрирования останавливается, при этом в переменную te записывается время наступления события, в переменную qe вектор состояния системы (строка: \([ x,y,V_x,V_y ]\) ) в этом момент, в переменную ie – номер сработавшего события или номера событий, если “сработало” несколько условий сразу. При помощи конструкции switch, case определяется номер события и в соответствии с ним изменяются начальные условия – меняются значения скоростей с учетом коэффициента восстановления. Для пятого события определяется, не находится ли точка, скорость которой уменьшилась до уровня менее p.vmin, ниже уровня p.ymin+0.01. Если да, то в переменную stop записывается значение true. На выходе из блока switch Задается новое значение начала интервала интегрирования и результаты интегрирования присоединяются к массивам t и q.
stop = false;
for i=1:100
[ti,qi,te,qe,ie] = ode45(@(t,q) pool_dqdt(t,q,p), [t0:0.1:t0+dt], q0, opt);
switch ie
case 1
q0 = qe'.*[1; 1; -0.9; 1];
case 2
q0 = qe'.*[1; 1; -0.9; 1];
case 3
q0 = qe'.*[1; 1; 1.0 ;-0.9];
case 4
q0 = qe'.*[1; 1; 1.0 ;-0.9];
case 5
if qe(2) < p.ymin+0.01
stop = true;
end
q0 = qe';
end
t0 = te;
q = [q;qi];
t = [t;ti];
if stop
break;
end
end
В следующем блоке кода главного файл-скрипта из массива решения удаляются строки с одинаковыми значениями времени при помощи функции unique. Это необходимо сделать для того, чтобы можно было выполнять интерполирование по таблице полученного решения. Функция интерполирования требует монотонного роста таблицы аргумента (времени) интерполируемой табличной функции.
[t, iu] = unique(t);
q = q(iu,:);
Далее создается видеофайл, в который записываются кадры движения материальной точки. Весь интервал времени от 0 до max(t) разбивается на 800 точек равномерно при помощи функции linspace. При помощи функции интерполяции interp1 определяется положение точки в каждый момент времени ti (от 0 до max(t)).
v = VideoWriter('pool.avi');
open(v);
for ti = linspace(0,max(t),800)
r = interp1(t,q(:,1:2),[(ti-0.5)*(ti>0.5):0.02:ti]);
plot(r(:,1),r(:,2),'*-','LineWidth',1);
daspect([1,1,1]);
xlim([p.xmin,p.xmax]);
ylim([p.ymin,p.ymax]);
frame = getframe(gcf);
writeVideo(v,frame);
end
close(v);
