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

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

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

d2x | Р dx _ q d2y і Р dy dt2 m dt ’ dt2 m dt

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

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

Таблица 3

Нестационарные дифференциальные уравнения

ode45

метод среднего порядка точности

ode23

метод низкого порядка точности

odel!3

метод переменного порядка точности

Стационарные дифференциальные уравнения

ode 15s

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

ode23s

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

ode23t

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

ode23tb

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

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

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

dX dVx -0

dY

~dt~Vy>

dVy -0 ^ = ^-уу~д-

— =VX, -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.
 
Посмотреть оригинал
< Пред   СОДЕРЖАНИЕ   ОРИГИНАЛ   След >