- •Применение информационных технологий при решении задач теоретической механики кинематика точки
- •191028, С.-Петербург, ул. Моховая, 26
- •Введение
- •Начало работы над программой
- •Проект программы для обработки массивов значений координат и времени
- •Составление программы на основе обработки символьных выражений
- •Листинг простейшей программы
- •Движение точки по пространственной кривой
- •Анимация
- •Дополнительные задания
- •Заключение
- •Список литературы
Проект программы для обработки массивов значений координат и времени
Рассмотрим два способа решения нашей задачи: с помощью обработки числовых массивов и с помощью символьных преобразований математических выражений. Идея первого способа состоит в том, что траекторию точки можно приближенно построить, задав последовательность из n точек плоскости и соединив их отрезками прямых. Для этого следует создать два вектора: X = [X(1) X(2) … X(n)] – вектор абсцисс и Y = [Y(1) Y(2) … Y(n)] – вектор ординат узловых точек траектории и вызвать встроенную функцию построения графика путем линейной интерполяции: plot(X,Y). Дополнительными входными параметрами функции plot могут быть установки вида, цвета, толщины линии и маркеров, помечающих узловые точки и др. Чтобы узнать об этом подробнее, можно вызвать позицию меню «Help» или значок «?» в инструментарии и поискать описание команды plot .
Чтобы найти и сохранить значения скорости точки на разных этапах движения, создаем вектор узловых значений времени T = [T(1) T(2) … T(n)]. Применив далее команду Dx = diff(X), получаем вектор Dx конечных разностей [X(2) – X(1), X(3) – X(2), … X(n) – X(n-1)]. Создав аналогично вектор Dy, а затем вектор разностей значений времени Dt, находим средние значения проекций скорости точки на отдельных прямолинейных участках: Vx(k) = Dx(k)/Dt(k), Vy(k) = Dy(k)/Dt(k), из которых создаем векторы Vx и Vy длиной n-1. Построив далее векторы конечных разностей следующего, второго порядка DVx, DVy и DDt, создаем векторы ( длиной n-2) средних ускорений Ax, Ay. Далее по известным формулам рассчтываем вектор касательных ускорений, вектор нормальных ускорений, вектор значений радиуса кривизны траектории.
Чтобы найти какой-то кинематический параметр, например, vx в некоторый момент времени, заключенный между двумя узловыми точками, можно посчитать его как среднюю скорость и приравнять соответствующему элементу вектора Vx. Более точный счет способен дать метод интерполирования (см. «Методы вычислений»).
Составление программы на основе обработки символьных выражений
В нашем задании уравнения движения точки содержат «элементарные» функции, встроенные в набор возможностей программы МАТЛАБ. Поэтому будем решать задачу с помощью символьных преобразований. Объявляем x, y, t переменными символьного типа (иначе по умолчанию они считались бы вещественными переменными типа double):
syms x y t;
Символ « ; » (semicolon), закрывающей строку, означает, что результат операции не выводится на экран в окне «Command Window». Набираем команду интерактивного ввода номера варианта:
novar = input('task_K1 Введи номер варианта = ');
Впоследствии, при выполнении этой инструкции, в командном окне появится
5
текст, содержащийся между апострофами. Программа будет ожидать ввода. Там, куда указывает курсор, надо будет набрать номер варианта и нажать «Enter». Аналогично пишем команды ввода выражений и числа . При выполнении программы по указанию курсора будем вводить соответствующие математические выражения. Число вводится с десятичной точкой (а не запятой!). Аргумент функции указывается в круглых скобках. Например, если , будем набирать выражение 2.5*sin(pi * t / 3). Число - встроенная в МАТЛАБ константа, обозначаемая «pi».
X = input('Введи x = ');
y = input('Введи y = ');
t1 = input('Введи t1 = ');
mv = input(‘Введи scale mv = ‘);
ma = input(‘Введи scale ma = ‘);
Масштабы mv, ma для построения векторов скорости и ускорения будем подбирать, глядя на результат построения и повторяя выполнение программы. Сначала можно назначит их равными единице.
Зададимся далее промежутком времени , на котором будем исследовать движение точки:
tedge = input('Введи правую границу промежутка tedge = ');
Обеспечим выполнение условия . Отдавая приоритет величине , отодвигаем границу в случае невыполнения данного условия:
if t1>tedge
tedge = 2*t1;
end%if t1
Здесь использован условный оператор, начинающийся ключевым словом if и заканчивающийся словом end.
Если траектория – замкнутая кривая типа фигуры Лиссажу, то желательно изобразить всю фигуру целиком, и назначаем , где - общий период функций x(t),y(t). Так, например, функция имеет период , функция имеет период .
Вызываем функцию ezplot построения графика на основе символьных выражений. В скобках на строго определенных местах указываем входные параметры функции: имена функций x, y и границы промежутка . Границы промежутка представлены вектором из двух элементов. Элементы вектора заключены в квадратных скобках.
ezplot(x,y,[0,tedge]), axis equal, hold on;
При этом создается графический файл figure. Команда axis equal задает одинаковый масштаб на осях Ox, Oy. В заголовок рисунка автоматически выводятся формулы . Команда hold on означает сохранение графика при добавлении следующего графического объекта (вектора скорости). Если не задать входной параметр
6
[0,tedge], то по умолчанию параметр t примет значения от 0 до 2π.
Начинаем счет. Рассчитываем числовые значения символьных переменных t, x, y с
помощью команды subs (substitute) и строим расчетную точку на траектории с помощью функции plot. Выражение 'ko' означает, что точка будет обозначена кружком «о» черного цвета «k».
T = t1;
x1 = subs(x); y1 = subs(y);
plot(x1,y1, ‘ko’); hold on;
Используем оператор diff для нахождения скорости и ускорения точки. Будучи применен к символьной функции x(t), он создает символьную производную dx/dt. Находим проекции скорости и ускорения, а также модули этих векторов:
vx = diff(x); vy = diff(y);
ax = diff(vx); ay = diff(vy);
vx1 = subs(vx); vy1 = subs(vy);
ax1 = subs(ax); ay1 = subs(ay);
v1 = sqrt(vx1*vx1+vy1*vy1);
a1 = sqrt(ax1*ax1+ay1*ay1);
Для рисования стрелок, изображающих векторы скорости и ускорения, вызываем команду-функцию quiver («колчан стрел»):
quiver(x1, y1, vx1*mv, vy1*mv, ‘g’), hold on;
quiver(x1, y1, ax1*ma, ay1*ma, ‘r’);
Первые два аргумента этой функции – координаты начала вектора-стрелки, вторые два – проекции вектора. В апострофах указан шифр цвета: 'g' (green)– зеленый, 'r'- красный. Дальше по известным из теории формулам высчитываем касательное ускорение, нормальное и радиус кривизны траектории в расчетной точке. Особым назовем случай, когда скорость точки близка к нулю. Из формулы для касательного ускорения следует, что при условии касательное ускорение является неопределенностью типа 0/0, которую надо раскрывать. Выведем из рассмотрения особый случай с помощью условного оператора «if» (тело которого непременно должно заканчиваться ключевым словом «end») и организуем две ветви расчета : стандартную и особую. В особом случае естественным составляющим ускорения и радиусу кривизны присвоим значение «NaN» («Not a Number»). В стандартной версии еще один «подводный камень» - случай, когда при расчете радиуса кривизны знаменатель близок к нулю. Потому организуем две стандартные ветви расчета: расчет в случае, когда радиус кривизны конечен, и в случае, когда он бесконечно велик. В последнем случае придаем ему значение «Inf» - «бесконечность».
If v1>1e-12
at1 = (vx1*ax1+vy1*ay1)/v1;
an1 = sqrt((a1+at1)*abs(a1-at1));
if an1<1e-12
ro1=Inf;
else
ro1 = v1*v1/an1;
end%if an1
else
7
at1 = NaN;
an1 = NaN;
ro1 = NaN;
end%if v1
Для удобства чтения записываем строки со сдвигом, помечая, какой «end» к какому «if» относится. Число 1e-12 равно . В конце программы организуем вывод результатов счета в командное окно. Именуем выводимую переменную одним (пусть и длинным) словом и выводим, не проставляя символ « ; » (semicolon) в конце строки. Выводимую переменную сделаем вектором-строкой, содержащей ряд чисел, разделенных пробелом или запятой. Компоненты вектора заключаем в квадратные скобки. Ниже приведена вся наша процедура.