Составление программы на основе обработки символьных выражений

       В нашем задании уравнения движения точки содержат «элементарные» функции, встроенные в набор возможностей программы МАТЛАБ. Поэтому будем решать задачу с помощью символьных преобразований. Объявляем 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) в конце строки. Выводимую переменную сделаем вектором-строкой, содержащей ряд чисел, разделенных пробелом или запятой. Компоненты  вектора заключаем в квадратные скобки. Ниже приведена вся наша процедура.

 

Листинг простейшей программы

%task_k1: evaluate the path, veloc. and accel.

%for given x(t),y(t),t1

clear all

syms x y t;

novar = input('Task K.1   Введи номер варианта = ');

x = input('Введи x = ');

y = input('Введи y = ');

t1 = input('Введи момент t1 = ');

mv = input('Введи масштаб изображения скорости mv = ');

ma = input('Введи масштаб изображения ускорения ma = ');

tedge = input('Введи правую границу промежутка tedge = ');

if t1>tedge

tedge = 2*t1;

end%if t1

ezplot(x,y,[0,tedge]), axis equal, hold on;

t = t1;

x1 = subs(x); y1 = subs(y);

plot(x1,y1, 'ko'), hold on;

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(x1, y1, vx1*mv, vy1*mv, 'g'); hold on;

quiver(x1, y1, ax1*ma, ay1*ma, 'r'); hold off;

if v1>0.000001

at1 = (vx1*ax1+vy1*ay1)/v1;

an1 = sqrt((a1+at1)*abs(a1-at1));

  if an1<0.000001

     ro1=Inf;

  else

     ro1 = v1*v1/an1;

  end%if an1

else

at1 = Nan;

an1 = Nan;

ro1 = Nan;

end%if v1

no_var = novar

vx1_vy1_v1 = [vx1 vy1 v1]

ax1_ay1_a1 = [ax1 ay1 a1]

at1_an1_ro1 = [at1 an1 ro1]

    Сохраняем программу (save as) под именем, скажем, task_k1. Набираем и запускаем команду <task_k1> в командном окне. В этом же окне отобразятся результаты счета, а

8

в графическом окне – соответствующая фигура. Используя инструменты этого окна, можно вставлять надписи, менять толщину и цвет линий и проч.

Замечания к составлению листинга.

1) На одной строке могут находиться несколько команд.

2) Если команда не умещается в одну строку, ставим многоточие «…» и переносим команду на другую строку.

3) Некоторые (не вычисляемые) команды можно заканчивать запятой или пробелом.