Движение при наличии сопротивления среды

Рассмотрим движение тела в среде с силой сопротивления, пропорциональной скорости движения.

В этом случае уравнения движения имеют вид

d²x | Р dx _ q d²y і Р dy dt² m dt ’ dt² m dt

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

Для численного решения обыкновенных уравнений существует много методов, наиболее известны из них метод Эйлера, метод Адамса, семейство методов Рунге - Кутты [3]. Отличие между методами заключается в разных способах конечно-разностного представления производной и соответственно разной степени точности получаемых расчетов. В MATLAB численное решение дифференциальных уравнений реализовано следующим набором функций (полный список функций можно посмотреть в руководстве):

Таблица 3

В простейшем случае функции имеют следующий шаблон вызова: [переменная, функция] = ode4 5 (символьная производная функции, диапазон изменения переменной, начальные условия).

При использовании MATLAB нужно помнить, что система способна решать ОДУ только первого порядка. Поэтому уравнения более высокого порядка должны быть сначала приведены к системе уравнений первого порядка. В нашем случае

dX dV_x -0

dY

~dt~^Vy>

dV_y -0 ^ = ^-^уу~^д-

— =V_X, -7^ = — • ^ = 0, dt ^x dtm

Рассмотрим алгоритм решения.

Пример: численное решение системы ОДУ второго порядка: % Сброс всех переменных clear all;

% Начальные условия

V0=10;

Alpha0=60;

BetaO=l;

Mass=l;

X0=0;

Y0=0;

g=9.8;

% Формируем начальные данные для системы ОДУ Condition=[0; V0*cosd(Alpha0);0;V0*sind(Alpha0)];

% Задаем промежуток времени tspan=[0 2];

% Решаем систему

[t,R] = ode45(@MyDiffEquation, tspan, Condition);

%Построение графиков

figure();

hold on;

plot(t,R, 'Line Width2);

legend ('X(t)'Vx(t)'Y(t)'Vy(t)');

figure(); hold on; plot(R(:, 1),R(:,3), 'Line Width ',2); legend ('Y(X));

% Формируем систему уравнений

%R(1) — соответствует x(t),

% R(2) — соответствует Vx=dx/dt.

% R(3) - соответствует y(t),

% R(4) — соответствует Vy=dy/dt.

function dRdt=MyDiffEquation(t,R) ?0=10;

Alpha0=60;

Beta=l;

Mass=l;

X0=0;

Y0=0;

g=9.8;

dRdt=[R(2); -R(2)*Beta/Mass; R(4); -R(4)*Beta/Mass-g]; end

Функция MyDiffEquation(t,R) получает в качестве аргументов названия переменной и системы дифференциальных уравнений и возвращает символьные выражения для производных функций системы. В соответствии с правилами MATLAB функции, созданные пользователем, размещаются в конце программного кода.

Сформированные выражения передаются для решения функции ode45. Символ стоящий перед первым аргументом (функцией MyDiffEquation), говорит о том, что будут переданы не численные значения, а аналитическая формула или алгоритм вычисления величины. Одновременно функции ode45 передаются диапазон изменения переменной и начальные значения. В результате применения функции ode45 формируется массив данных, столбцы которого содержат зависимость от времени для искомых функций.

В результате выполнения программы будут построены графики зависимости кинематических величин от времени и график траектории движения (рис. 8 и 9).

Рис. 8. Зависимость кинематических величин от времени

Рис. 9. Траектория движения

В GNU Octave для решения ОДУ используется функция Isode ( ). Шаблон вызова функции имеет вид:

[х, istate, msg]=lsode (символьная производная функции, начальные условия, диапазон изменения переменной).

Массив данных с результатами численного решения будет передан переменной х, переменная istate содержит код правильности выполнения

(по умолчанию 2), переменная msg содержит дополнительное сообщение системы о процессе решения. При использовании функции Isode GNU Octave самостоятельно подберет алгоритм решения, при желании пользователь может сделать это самостоятельно с помощью функции Isode options (). Можно определять метод решения, параметры сходимости, шаг и др. При использовании специализированного пакета odepkg появляется дополнительный набор специализированных функций.

Задания к лабораторному практикуму

• 1. По заданию преподавателя запустите MATLAB или GNU Octave.
• 2. Напишите программу для построения траектории движения тела, брошенного под углом к горизонту в среде без сопротивления.
• 3. Решите численно систему ОДУ, описывающих одномерное движение тела в среде с сопротивлением. По заданию преподавателя проведите расчет для нескольких вариантов силы сопротивления.
• 4. Постройте графические зависимости кинематических характеристик от времени, экспортируйте полученные графики в формат .pdf.