Задачи динамики

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

Задача 6.1

Частица, имеющая массу m и заряд е, влетает в однородное стационарное электрическое поле Е со скоростью v, перпендикулярной направлению поля. Определить траекторию движения частицы.

Для решения этой задачи следует в первую очередь записать уравнение Ньютона (второй закон). Уравнение, как известно, записывается для вектора, поэтому в действительности это система трех уравнений — для каждой из координатных осей свое уравнение. Систему координат следует задать оптимальным (наиболее удобным) образом. Поскольку в условии сказано, что векторы электрического поля и начальной скорости частицы перпендикулярны, можем выбрать систему координат так, чтобы ось X совпадала с направлением поля, ось У — с направлением начальной скорости, и, соответственно, ось Z будет перпендикулярной плоскости, в которой размещены вектор поля и начальной скорости. Начало отсчета выберем таким образом, чтобы в момент времени t=0 частица находилась в начале координат.

Теперь записываем уравнения движения. Сначала описываем движение вдоль оси X. Это единственное направление, вдоль которого на систему действует сила, равная произведению напряженности поля на величину заряда частицы.

Вдоль остальных двух осей силы не действуют (силы фавитации не учитываем).

Чтобы однозначно решить систему из приведенных выше трех уравнений, Цйеобходимо задать начальные условия. Поскольку уравнения имеют второй эрядок (порядок старшей производной), по каждой из координат необходимо задать два условия — начальное положение (соответствующая координата частицы) и проекция начальной скорости на координатную ось.

Начальное значение координат частицы для каждой из осей равно нулю выбрали начало системы координат!). Что касается проекций вектора Начальной скорости, то отличной от нуля будет только проекция на ось Y. Причем значение этой проекции равно модулю начальной скорости, т.е. v.

Удобства ради, начальные условия разобьем на три фуппы — в соответствии с числом координатных осей.

Таким образом, для осей имеем следующее.

В следующих двух переменных объединяем (в виде последовательности) уравнения (Eq_All) и начальные условия (InCon_All).

Для решения полученной системы с описанными выше начальными условиями используем процедуру dsolve(), а результат ее выполнения присваиваем в качестве значения переменной Res.

Сразу видим, что, поскольку z(t)=O, движение происходит в плоскости XY — в плоскости, где размещены векторы поля и начальной скорости.

Для дальнейшего анализа полезно записать временную зависимость координаты по каждой из осей х и Y в виде функций.

Исходя из этих зависимостей строим траекторию частицы. Для определенности будем исследовать движение частицы в течение первых 10 секунд. Но прежде чем строить траекторию, определим координаты частицы через каждую секунду после попадания ее в область действия поля. Эти координаты будем записывать в список List. С этой целью инициализируем соответствующую переменную, присвоив ей в качестве начального значения пустой список.

Дальше с помощью оператора цикла вычисляем координаты (х- и у-координату) через одну, две и так далее, до десятой секунды. Значения координат записываются в виде списка, и этот список добавляется в качестве элемента в список List, т.е. List является списком списков. Кроме того, следует отметить, что координаты точек нормированы соответствующим образом. Это сделано для того, чтобы точки в дальнейшем можно было отобразить на графике.

Можно просмотреть, как же в конечном счете выглядит список координат (при выполнении предыдущего оператора в области вывода ничего не отображается, поскольку соответствующая команда заканчивается двоеточием, а не точкой с запятой).

Теперь можно построить траекторию движения частицы. Квадратами на этой траектории будем отмечать положения частицы через каждую секунду. Шри этом первым параметром процедуры plot() указываем список с двумя тементами: первый элемент — это список для отображения заданной в параметрическом виде функции (это и есть траектория частицы), второй — переменная List, которая, как известно, содержит координаты частицы.

Однако графического отображения траектории мало — нужно определить уравнение этой траектории. Для этого из одной функциональной зависимости следует выразить время через координату и подставить полученное таким образом выражение в другое уравнение.

Чтобы отличить непосредственно координату от функциональной зависимости этой координаты от времени, введем переменные X и Y. Уравнение движения вдоль оси х может быть записано в новых переменных как X=x(t), a вдоль оси у — как Y=y(t).

Из последнего уравнения и выразим время через координату. Сделать это можно посредством команды isolate(Y=y(t),t). Далее выполняем подстановку в оставшемся уравнении с помощью процедуры subs(), указав это уравнение вторым параметром, а процедуру isolate)) — первым.

Таким образом, частица, попав в область действия поля, движется по параболе в плоскости, определяемой векторами поля и начальной скорости. Уравнение параболы представлено выше.

Задача 6.2

В некоторой области пространства одновременно имеются однородные и стационарные электрическое и магнитное поля, угол между которыми а. Частица с массой m и зарядом е, имеющая начальную скорость v, попадает в это пространство. Определить траекторию движения частицы.

Принципиально эта задача мало чем отличается от предыдущей. Ситуация, правда, несколько усложняется из-за наличия магнитного поля. Дело в том, что сила Лоренца (она определяет силу, действующую на частицу со стороны магнитного поля) выражается через векторное произведение. Трудность эта не принципиальна, но вычисления по использованной ранее схеме достаточно громоздки.

В данном случае решать будем векторное уравнение. Для этого понадобится определить две процедуры: для дифференцирования вектора и вычисления векторного произведения.

Процедура дифференцирования вектора (в наиболее простом варианте) выглядит следующим образом.

Данная процедура имеет два параметра — вектор (размерности 3), который и необходимо дифференцировать, и переменную, по которой следует вычислять производную. Процедура возвращает в качестве результата вектор, компоненты которого определяются через производные соответствующих компонентов вектора, указанного первым аргументом процедуры.

Следующая процедура vprod() нужна для вычисления векторного произведения двух векторов. В соответствии с правилом вычисления такого произведения, в теле процедуры компоненты результирующего вектора определены напрямую через компоненты векторов-параметров процедуры.

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

Ось Y направим вдоль вектора магнитного поля, а ось X — перпендикулярно плоскости векторов электрического и магнитных полей. Тогда проекция вектора электрического поля на ось X равна нулю, а на прочие оси определяется следующим образом.

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

Для определения координат частицы введем вектор г, который опишем следующим образом.

Другими словами, г является вектор-функцией, зависящей от времени.

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

Чтобы записать соответствующие уравнения, воспользуемся процедурой генерирования последовательностей seq().

Первым параметром процедуры seq() указано уравнение r(0)[i]=0. В его левой части содержится индексная переменная i, которая последовательно изменяется от 1 до 3 и определяет номер компонента функции-вектора г(). В качестве аргумента этой функции указан 0, поэтому в левой части уравнения соответствующие компоненты берутся в начальный момент. Результат (последовательность из трех уравнений) присваивается в качестве значения переменной IniCon_R.

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

Последовательность формируется все той же процедурой seq(). Для индексной переменной s указан список принимаемых значений (это — х, у и z).

Оператор D(s)(0) задает, в зависимости от принимаемого переменной s значения, первую производную от соответствующей координаты в начальный момент времени. Правая часть формируемого уравнения имеет вид v| |s. Так с помощью оператора объединения названий (11) символу V приписывается нужный суффикс.

С помощью следующей команды в переменной IniCon объединяются начальные условия для координат и скоростей.

Уравнение присваивается в качестве значения переменной VecEq. Левая часть уравнения представляет собой произведение массы на вторую производную от радиус-вектора (вектора координат), т.е. ускорения.

В правой части уравнения записан вектор действующей на частицу силы. Этот вектор равен сумме (умноженной на заряд) двух векторов: вектора электрического поля и вектора, являющегося векторным произведением скорости частицы на магнитную индукцию. Это векторное произведение записано с помощью описанной в самом начале процедуры vprod().

Далее записанное уравнение следует решить. Для этого формируем последовательность из трех элементов-уравнений. Каждое такое уравнение гописывает динамику частицы вдоль соответствующей координатной оси. «Чтобы составить эти уравнения, выбираем из правой и левой частей ис-рсодного векторного уравнения элементы списков с одинаковыми индексами и приравниваем их.

Доступ к левой части векторного уравнения осуществляется с помощью команды lhs(VecEq), а к правой — rhs(VecEq). На операнды эти части разби-гся посредством процедуры ор(). Первым операндом как в правой, так и левой части уравнения являются скалярные множители: для левой части — масса т, для правой — заряд е. На эти скаляры будут также множиться и довые уравнения. Вторые операнды правой и левой частей векторного уравнения — списки. Доступ к элементам списков реализуется путем указания вддекса этих элементов, например op(lhs(VecEq))[2][3] — третий элемент второго операнда левой части уравнения VecEq.

После последовательности из уравнений, формируемой процедурой seq(), ша переменная iniCon, в которой записаны все начальные условия. По-недовательность уравнений и начальных условий заключена в фигурные |скобки, после которых указано множество функций параметра t, относитель-шо которых следует решать систему.

Поступить, кроме прочего, можно следующим образом.

На Данном этапе значением переменной среды % является множество, элементами которого есть равенства, определяющие эволюцию частицы вдоль каждой оси. Уравнение для зависимости x(t) является третьим элементом множества. Это уравнение запишем в переменную Rx.

Внимание!
To, что уравнение для x(t) является именно третьим элементом множества решений, — факт достаточно случайный. От сеанса к сеансу это уравнение может оказаться и первым, и вторым элементом. Поэтому правильно указать индекс можно только после того, как решение системы уравнений отображено в области вывода. Кроме того, не рекомендуется присваивать решение системы в качестве значения переменной для дальнейшего вызова в формате переменная[индекс]. Дело в том, что при вызове переменной-множества порядок ее элементов может меняться. С точки зрения Maple переменная при этом не меняется, поскольку при изменении порядка элементов множества по определению полагают, что множество осталось неизменным. Здесь скрывается потенциальная опасность, и ее следует иметь в виду!

После выполнения предыдущей операции решение векторного уравнения записано уже в переменную %. Для оси Y уравнение будет следующим (переменная Ry).

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

В качестве частицы рассмотрим позитрон (то же, что электрон, но с положительным зарядом). Масса такой частицы (в килограммах) приведена ниже.

После этого можно строить траекторию частицы. Исследовать будем движение позитрона в течение 0.05 секунды (для позитрона и этого много). Поскольку строить предстоит параметрически заданную кривую в трехмерном пространстве, воспользуемся процедурой spacecurve() из пакета plots. Ее первым параметром является список с параметрическими зависимостями координат частицы вдоль каждой оси от времени, затем указывается интервал изменения параметра (в данном случае это время). Кроме того, указаны задающие ориентацию графика углы (orientation=[30,75]), тип отображаемых координатных осей (axes=FRAMED), число базовых точек (numpoints=200). Остальные опции (толщина и цвет линии, название графика, шрифт для названия) читателю должны быть уже знакомы.

Внимание!
В качестве отображаемых параметрических зависимостей указываются не сами уравнения (Rx.Ry.Rz), а их правые части (rhs(Rx),rhs(Ry),rhs(Rz)).

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

Задача 6.3

Груз массы М падает без начальной скорости с высоты Н на спиральную пружину. Под действием упавшего груза пружина сжимается на величину а. Вычислить время сжатия пружины (массой пружины и силами трения можно пренебречь).

В первую очередь задаем уравнения движения. Поскольку движение одномерное, уравнение будет всего одно. Координатную ось X, вдоль которой и будет двигаться шарик, направим вверх. Точку начала отсчета выберем так, чтобы она соответствовала верхнему свободному концу несжатой пружины (это будет точка 0). В общем виде, согласно второму закону Ньютона, записываем следующее (x(t) — координата шарика в момент времени t).

Однако для того, чтобы решить присвоенное в качестве значения переменной Eq уравнение, необходимо задать силу F(x), которая существенно зависит от того, долетел шарик до пружины или нет. Поскольку необходимо вычислить время сжатия пружины, движение шарика будем рассматривать начиная с того момента, когда шарик долетает до пружины. Другими словами, в момент времени t=0, по определению, координата шарика х(0)=0. Это, кстати, будет одним из двух начальных условий (второе — для скорости).

Поэтому равнодействующая двух сил, равная их сумме с учетом направленности, равна следующему.

Последнее выражение определяет действующую на шарик силу как функ-Щию координаты.

Далее, чтобы решить дифференциальное уравнение, необходимо определить еще одно начальное условие для скорости, т.е. необходимо определить .'скорость в нулевой момент времени — при падении шарика на пружину. [&ля этого воспользуемся законом сохранения энергии. Поскольку шарик начинал падать с высоты Н без начальной скорости, его полная энергия равнялась потенциальной и была равна М*д*Н. При падении шарика на пружину потенциальная энергия равна нулю (поскольку система координат выбрана так, что в момент столкновения с пружиной шарик находится на нулевом уровне). Однако кинетическая энергия отлична от нуля. Если скорость шарика в момент столкновения равна v, то кинетическая энергия, с одной стороны, равна M*v*2/2, а с другой — должна быть равна полной энергии М*д*Н. Из этого уравнения находим скорость шарика при столкновении с пружиной.

Для решения уравнения используем процедуру solve(), указав, что решать уравнение следует относительно переменной v.

В результате получаем два решения, которые отличаются только знаком. Выражение со знаком "минус" соответствует движению вниз (в направлении, Противоположном направлению координатной оси), а другое, соответственно, — движению вверх. В данном случае очевидно, что шарик движется вниз, поэтому именно отрицательное решение и выбираем (в списке полученных на предыдущем этапе решений оно является вторым).

После этого можно решать дифференциальное уравнение Eq. Если добавить к этому уравнению начальные условия, решение будет определено однозначно (уравнение с начальными условиями — это, как известно, задача Ко-ши). Решаемое уравнение и его начальные условия в процедуре dsolve() заключаются в фигурные скобки.

На заметку
Поскольку начальное значение для скорости — это производная в точке t=0, для записи этого условия (производной) используется оператор D- D(x)(0)=v.

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

Внимание!
Зависимость X справедлива до тех пор, пока шарик находится в контакте с пружиной Как только шарик от пружины отскочит, уравнение движения будет иным.

Зависимость скорости шарика от времени определяется через производную по времени от переменной X (которая, напоминаем, является зависимостью координаты шарика от времени).

Скорость нам понадобится вот для каких целей. Отсчет времени начинается с момента столкновения шарика с пружиной. Далее пружина сжимается на какую-то граничную величину (в условии задачи это а), после чего начинает разжиматься. Следовательно, время, в течение которого пружина сжимается, [равно по абсолютной величине (при данном выборе начала отсчета времени) шоменту, когда пружина прекращает сжатие и начинает разжиматься. Этот |момент характерен тем, что скорость шарика равна нулю!

Отсюда дальнейшие действия очевидны (почти!) — нужно определить «омент времени, когда производная от координаты (т.е. скорость) равна нулю. Однако здесь есть один нюанс. Дело в том, что полученное решение исходного внения и, соответственно, временная зависимость для скорости, содержат эигонометрические функции. При решении уравнения V=0 относительно t найдены не все решения (в данном случае только на интервале изменения рктангенса, и это решение не будет учитывать периодичность общего решения). Иногда в этом нет ничего страшного, но только не в этом случае. Другими Иовами, необходимо найти все решения уравнения V=0. С этой целью изменяем ачение переменной среды _EnvAllSolutions (переменная определяет, следует искать абсолютно все решения) на true.

Находим момент времени (если точнее, то моменты, но реальный интерес физический смысл имеет только один), в который скорость равна нулю. гшение присваиваем в качестве значения переменной Т.

Здесь переменная _Z1 обозначает любое целое число, поскольку данное сражение описывает общее решение, из которого следует выделить единст-енное. Определяется оно очень просто — это первое положительное реше-ше. Можно использовать для этого мощный "арсенал" Maple, однако в данном случае поступим проще.

Поскольку все переменные в аргументе арктангенса больше нуля, сам арктангенс возвращает значение в диапазоне от 0 до л/2, и этот арктангенс в выражение для Т входит со знаком "минус". Поэтому очевидно, что первым не отрицательным решением будет то, где _Z1=1. Это значение переменной _Z1 и присвоим. Используем процедуру subs(): первый параметр — равенство _Z1=1 (вместо переменной _Z1 следует подставить 1), второй параметр — Т (замену нужно делать в выражении для Т). Результат выполнения команды присваиваем переменной Т.

еперь переменная Т определяет тот самый момент — "единственный и неповторимый".

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

После этого выражение для времени сжатия пружины примет следующий вид.

Это выражение следует упростить. Чтобы упрощение возымело действие относительно самой переменной Т, а не просто было выведено упрощенное выражение, результат упрощения присвоим переменной Т. Кроме того, поскольку в выражении много радикалов, указываем опцию symbolic. В этом случае при вычислении корня квадратного из переменной в квадрате результатом будет сама эта переменная. Это удобно, когда используются положительные величины — нет необходимости вызывать процедуру assume ().

Таким образом, имеем следующее.

Это и есть ответ. Однако на достигнутом останавливаться не будем. Рассмотрим методы, с помощью которых для задач можно создавать иллюстрации (и не только!).

Для начала определим процедуру spring(), с помощью которой будет отоаться пружина. Данная процедура имеет три параметра: максимальное зачение координаты по вертикальной оси (Хтах), т.е. верхняя граница для гй картинки; закон движения верхнего свободного конца пружины L() — не выражение, а процедура (оператор); момент времени t, в который ото-эажается пружина. Другими словами, процедура создается "с перспектиэй" — чтобы можно было отображать состояние системы шарик-пружина ; произвольный момент времени.

На заметку
Наличие у процедуры spring() параметра Хтах объясняется следующим образом: в пьнейшем пружина будет отображаться "в компании" шарика. Поэтому размер рисунка i вертикали) зависит от высоты, на которой находится шарик. Наличие у процедуры змянутого параметра позволит затем связать его с положением шарика.

Далее, в процедуре используется одна локальная переменная (i) для запиоператора цикла и пять глобальных переменных. Подразумевается, что эти временные задаются вне процедуры до ее вызова. Назначение глобальных гменных таково: переменная хО определяет координату оси пружины (по Эризонтали), переменная Xmin определяет координату (по вертикальной оси) iero стационарного конца пружины, N задает число витков пружины, 1 — шу пружины в свободном состоянии и s — ее толщину в процентном отношении от длины 1.

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

Далее картинку будем отображать квадратной — ее размеры по вертикали и горизонтали будут совпадать (однако заранее они неизвестны). По горизонтали отсчет будет начинаться с 0. Пружину разместим посредине. Поскольку размер рисунка по горизонтали (и вертикали) равен Xmax-Xmin, ось пружины будет иметь координату (Xmax-Xmin)/2. Последней командой в теле процедуры является plot(). Первым параметром этой процедуры, определяющим отображаемую структуру, является команда формирования последовательности seq(), заключенная в квадратные скобки. Этой командой формируется последовательность точек, определяющая пружину. Горизонтальная координата каждой такой точки формируется поочередным добавлением/отниманием от координаты оси пружины хО половинной толщины пружины. Поскольку параметр s определяет толщину пружины в процентах к ее длине 1, непосредственно толщина пружины равна s*l/100 (половинная толщина отсюда равна s*l/200). Множитель (-l)*(i+l) нужен для чередования знаков "плюс" и "минус" при определении координаты.

Поскольку пружина состоит из N витков, в боковой проекции это соответствует N "зубцам", а значит, 2*N отрезкам, соединяющим угловые точки. Таким образом, на каждом шаге (при увеличении переменной i на единицу) вертикальная координата угловой точки должна увеличиваться на величину L(t)/(2*N), где L(t) — длина пружины в момент времени t. Начальной должна быть точка с координатой Xmin, конечной — точка с координатой Xmin+L(t). Из этих условий выбирается множитель (i-1) и диапазон изменения параметра (i=l. .2*N+1).

Параметрами процедуры plot() указаны диапазон отображения по гори-зонт&ли (0..Xmax-Xmin), вертикали (Xmin..Xmax), а также цвет пружины (синий) и толщина линии.

Внимание!
Выше в качестве отображаемой процедурой plot() структуры указан список точек. Поскольку по умолчанию значение опции style равняется LINE, отображаемые точки будут соединены ломаной линией строго в порядке их следования в списке. Таким образом, задаем только точки, а получаем — пружину. Чтобы отображались исключительно точки, следует указать style=POINT.

Теперь определим процедуру Sys_display(), которая будет отображать пружину вместе с шариком, и в этой процедуре будет использоваться описанная выше процедура отображения пружины spring ().

Параметрами процедуры Sys display() являются:

а) функциональная зависимость длины пружины от времени L; б) зависимость высоты, на которой находится шарик, от времени h (это оператор, как и L); в) зависимость скорости шарика от времени VI (зачем это нужно, объясняется ниже); г) момент времени t, в который отображается вся система.

Процедурой, помимо отображения пружины и шарика, будут отображаться | текстовые поля с указанием момента времени после начала падения, высоты 1 шарика в этот момент и его скорости. В связи с этим в процедуре ниже объяв-1ляются локальные переменные: Н — для значения высоты шарика в начальный |момент t=0, остальные переменные — для определения текстовых полей.

Первой командой в процедуре переменной Н присваивается значение: по-шльку высота шарика в момент t задается зависимостью h(t), начальная лсота равна h(0). После этого объявляются три текстовых поля, которые отображаются посредством процедуры textplot() из пакета plots.

На заметку
Чтобы воспользоваться процедурой, например, textplot() из пакета plots, можно постуодним из следующих способов: а) подключить пакет (with(plots)) и затем вызвать рйроцедуру (textplot()); б) ввести команду plots[textplot] (); в) воспользоваться командой 'plots/textplot1. Выше в описании процедуры Sys_display использована последняя форма вызова.

Так, переменная TextRegl является графическим объектом текстовое поле. Непосредственно отображаемый текст указан третьим элементом первого па-'раметра-списка процедуры textplot(). Выводимая текстовая строка формиру-?ется процедурой объединения cat(), которая имеет в данном случае три параметра: первый и третий — готовые строчные выражения (Ч = " и " сек" соответственно), а второй — также строка, но полученная преобразованием из численного выражения для времени t (команда convert(t,string)).

Первые два элемента списка — это координаты для вывода текста. Так, от левого края текст выводится на расстоянии трех четвертей от расстояния между начальной высотой шарика и длиной недеформированной пружины плюс четыре толщины последней (при желании можно придумать что-нибудь попроще). По высоте поле размещается на уровне начального положения шарика. Опция нужна для того, чтобы текст выравнивался вправо.

Два следующих поля принципиально мало чем отличаются от TextRegl. Каждое последующее поле размещается по отношению к предыдущему ниже на величину 0.2*Н (т.е. 20% от первоначальной высоты шарика). Поле TextReg2 содержит сведения о высоте шарика, которая определяется зависимостью h(t). Точно так же поле TextReg3 используется для вывода данных о скорости шарика (зависимость Vl(t)).

На заметку
Строго говоря, между зависимостями h(t) и vl(t) существует очевидная взаимосвязь: Vl(t)=D(h) (t). Поэтому, если задавать отдельно зависимость скорости шарика от времени, выполняется, казалось бы, ненужная работа. Однако это не совсем так. Дело в том, что, как будет показано ниже, зависимость координаты шарика от времени, равно как и длины пружины, имеет нетривиальный вид. Например, зависимость h(t) представляет собой "сшивку" двух зависимостей: одна — для свободного падения шарика, вторая — для процесса сжатия пружины вместе с шариком. Чтобы в теле процедуры можно было в дальнейшем в аналитическом виде вычислить производную от h(t), эту зависимость придется специальным образом описывать. Поэтому, не усложняя задачу, задаем зависимость скорости от времени отдельно.

Наконец, переменной TextReg в качестве значения присваиваем последовательность из трех описанных выше текстовых полей. Таким образом, весь выводимый текст запоминается в этой переменной.

Далее в процедуре используется процедура display() из пакета plots, которая отображает перечисленные в первом ее аргументе (этот аргумент может быть списком, множеством или массивом) графические структуры.

Совет
С процедурой display)) допускается использовать опцию insequence, которая может иметь значение true или false (последнее является значением по умолчанию). Если insequence=false, то отображаются сразу все графические объекты. Чтобы отображать эти объекты в той очередности, в которой они указаны (это относится в первую очередь к спискам), следует указать insequence=true. Это, кстати, один из методов создания анимации.

Первой в списке отображаемых структур является пружина (spring(H+4*s*l/100,L,t)). Два последних параметра этой процедуры особых комментариев не требуют. Первый, как известно, определяет границу по высоте. При ее выборе имели место следующие соображения. Во-первых, в задаче предполагается, что шарик точечный, т.е. его размеры при решении задачи во внимание не принимались. В целях наглядности, на рисунке размеры шарика сбудут существенными. Чтобы не вносить поправки на размер шарика в полученном уже решении, будем полагать, что шарик "сосредоточен" в самой |его нижней точке. Далее радиус шарика примем равным толщине пружины. |Следовательно, нижняя точка шарика первоначально находится на высоте Н, центр шарика — на высоте H+s*l/100, а верхняя точка — на высоте B+2*s*l/100. Сверху над шариком оставляем пространство толщиной с его |диаметр. Отсюда имеем для верхней границы значение H+4*s*l/100.

На заметку
В данном случае имя глобальной переменной (высота в условии задачи Н) совпадает : именем локальной, объявленной в процедуре Sys_display(), переменной. По большому /, поскольку и объявленная локальная, и такая же глобальная переменные обозначают ну и ту же величину, ничего страшного нет. Но даже если бы локальная и глобальная переменные обозначали разные величины, Maple корректно разграничивает область их использования: в теле процедуры используется локальная переменная, вне процеду-i — глобальная. Мало того, даже если не объявить переменную как локальную, Maple введет сообщение о том, что в процедуре неявно задана локальная переменная, и энно так ее и будет интерпретировать.

После пружины отображается шарик. Окружность в Maple может быть эрмирована процедурой circle () из пакета plottools. Первым параметром эцедуры указывается точка (центр окружности), а вторым — радиус. Однако в данном случае шарик неплохо было бы закрасить. Поэтому соз-ем последовательность окружностей: у них совпадают центры и последователь-изменяются радиусы. Последовательность формируется процедурой seq(), где звый параметр— команда создания окружности ('plottools/circle'O) с цен-эм в точке [x0,h(t)+s*l/100] (зависимость h(t) определяет динамику нижней «ки шарика, а значит, центр находится на высоте h(t)+s*l/100) и радиусом K05*s*l*i/100, а индексная переменная принимает значения в диапазоне i=l. .20. Гаким образом, радиус шарика (т.е. внешний, самый большой радиус) равен, как следовало ожидать, s*l/100. Кроме того, одна из опций процедуры circled за желто-коричневый цвет шарика (color=TAN).

Совет
М ожно поэкспериментировать с количеством внутренних окружностей, т.е. изменить диапазон изменения переменной i (при этом нужно изменить коэффициент 0.05 на 1/Nmax, Уде Nmax — верхняя граница диапазона изменения i). В этом случае можно наблюдать достаточно интересные визуальные эффекты. Еще один способ — задать для разных окружностей разный цвет, т.е. сделать цвет зависимым от индекса i.

Чтобы шарику было о что удариться, изобразим горизонтальную поверхность, f ограничивающую верхний конец пружины. Это будет обычная линия (процедура 'plottools/line'O). Параметрами процедуры указываются начальная ([хО-8*1/100,Xmin+L(t)]) и конечная ([x0+s*l/100,Xmin+L(t)]) точки. Обе точки имеют одинаковую координату по высоте, которая равна длине пружины в момент времени t, отсчитанной от координаты нижней точки пружины, т.е. Xmin+L(t). Горизонтальные координаты выбраны так, чтобы поверхность-подложка размещалась симметрично относительно оси пружины, а ширина подложки совпадала с диаметром шарика. Толщина отображаемой линии установлена опцией thickness=3, а сама линия изображена красным цветом (color=RED).

Наконец, последним отображаемым объектом является текстовое поле (на самом деле их три!) TextReg с данными о моменте времени, высоте шарика и его скорости.

Опция axes=FRAME устанавливает тип рамки рисунка (точка пересечения координатных осей в левом нижнем углу рисунка), а опция scaling=constrained задает масштаб, при котором единицы измерения по разным осям координат совпадают, иначе шарик будет похож на эллипс.

Внимание!
Как отмечалось выше, отображаемые процедурой display!) объекты могут быть организованы в виде списка (массива) или множества. В данном случае это список (то есть последовательность отображаемых объектов заключена в квадратные скобки). При этом объекты отображаются именно в той последовательности, в какой они указаны. Это важно, поскольку процедурой spring () определяется ряд глобальных переменных! Эта процедура должна выполняться первой, т.е. первой должна отображаться пружина. Оформлять графические объекты в виде множества настоятельно не рекомендуется, поскольку элементы множества перебираются, в принципе, в произвольном порядке. Это может стать причиной серьезных недоразумений.

На этом описание процедуры Sys_display() заканчивается. Далее следует определиться с законом движения шарика и, разумеется, пружины.

Выше была вычислена зависимость высоты шарика от времени после падения его на пружину (переменная X). Ниже эта зависимость представлена в виде функции Y(t).

Внимание!
В приведенной выше зависимости время отсчитывается от момента падения шарика на пружину!

Как уже отмечалось ранее, зависимость высоты шарика от времени является "сшивкой" двух зависимостей: первая — это зависимость для свободно падающего шарика, вторая — зависимость для шарика, упавшего на пружину и движущегося далее вместе со сжимающейся пружиной.

Для дальнейшего анализа необходимо определить время свободного паде-1ния шарика с высоты Н (как долго шарик падает, прежде чем столкнуться с шружиной). Время это, как известно, равно следующему.

Для большей конкретности зададим такие значения для параметров задачи котя это можно сделать и позже) — длина пружины будет равна 5 метрам.

Толщина пружины будет составлять 20% от длины недеформированной ружины, т.е. 1 метр.

Теперь задаем функциональную зависимость высоты шарика от времени. Зависимость запишем так, чтобы можно было определить положение шарика в произвольный момент времени. Полезными будут следующие рассуждения.
Понятно, что движение шарика (при отсутствии трения) периодично — шарик будет периодически подпрыгивать, падать, ударяться о пружину, сжимать ее, двигаться вверх при разжатии пружины и снова подпрыгивать. Время от прыжка до прыжка (т.е. период), очевидно, равно удвоенному времени свободного падения шарика (to) и времени сжатия пружины (Т), т.е. период равен 2*(T+tO). Чтобы восстановить динамику системы, достаточно знать ее динамику на интервале времени от 0 до (T+tO). В силу симметрии уравнений механики и периодичности движения, в произвольный момент времени t положение и скорость шарика могут быть вычислены согласно следующим правилам. Во-первых, от интервала t можно отнять целое число периодов 2*(T+tO), при этом положение и скорость шарика не изменятся. Во-вторых, если t>(T+tO), то положение и модуль скорости шарика будут такими же, как в момент времени 2*(T+tO)-t, но только скорость имеет противоположный знак.
Зависимость, определяющая положение шарика в произвольный момент времени t, описана ниже как процедура h().

В теле процедуры используются две локальные переменные ТТ и j. В качестве значения переменной j присваивается остаток от целочисленного деления (irem()) результата выполнения операции trunc(t/(T+tO)) на 2. Функция trunc() вычисляет целую часть выражения, указанного в качестве ее аргумента. В данном случае с помощью функции trunc() устанавливается, сколько целых полупериодов (T+tO) укладывается в интервале времени t. Если это число четное, то целочисленный остаток от его деления на 2 (значение переменной j) равен 0. В противном случае значение j равно 1.
Если j=0, то локальной переменной ТТ присваивается значение (T+tO)*frac(t/(T+tO)), т.е. переходим к локальному времени, равному остатку от вычитания из параметра процедуры t целого числа периодов (первое правило вычисления высоты и скорости шарика).

На заметку
Результатом выполнения функция frac() является дробная часть ее аргумента. Таким образом, результат команды frac(t/(T+tO)) — остаток от вычитания из t целого числа полупериодов, но только в отношении к интервалу времени (T+tO). Чтобы получить время в абсолютных единицах, данное число следует умножить на (T+tO).

Если же в интервале времени t вмещается нечетное число полупериодов (j=l), используем второе правило вычисления высоты и скорости шарика и [присваиваем локальному параметру ТТ значение (T+tO) * (1-fгас (t/ (T+tO))).

После этого достаточно вычислить положение и скорость шарика в момент ТТ, который заведомо не превышает интервал (T+tO). I Далее, если ТТ не превышает времени свободного падения шарика (to), рследует использовать формулу для свободного падения тела с высоты Н без начальной скорости. В противном случае зависимость высоты шарика от времени дается зависимостью Y(). Однако поскольку в послед-ией зависимости предполагается, что время отсчитывается от момента парения шарика на пружину, аргументом следует указать параметр (TT-tO), ито и сделано.
Принципиально процедура L(), позволяющая вычислить длину пружины, вт описанной только что процедуры h() отличается мало. Локальные переменные и их значения абсолютно такие же, как и в предыдущем случае. От-иичие только в последнем условном операторе.

Если шарик еще не долетел до пружины, последняя находится в недефор-рированном состоянии и ее длина равна 1. После того как шарик столкнулся пружиной, она начинает сжиматься. Высота шарика определяется зависимостью Y(TT-tO), и эта высота меньше нуля (поскольку шарик находится нисе уровня верхнего конца недеформированной пружины). Поэтому длина Пружины изменяется на такую же величину и равна l+Y(TT,-tO). Осталось задать зависимость скорости шарика от времени. Эта зависимость определяется процедурой VI ().

Локальные переменные ТТ и j определяются так же, как и в предыдущих случаях. Если шарик находится "в свободном полете" (TT<tO), модуль скорости шарика равен д*ТТ, а знак определяется так: если в интервале t вмещается четное число периодов (j=0), вектор скорости направлен вниз и скорость отрицательна; в противном случае (j=l) она положительна, т.е. шарик движется вверх. Отсюда появляется множитель (-l)*(j+l): он равен 1 при j=l и -1 при j»0. В противном случае скорость вычисляется как производная (в виде оператора D(Y) ()) от функции У() в момент времени (TT-tO). Наличие множителя j обусловлено теми же причинами, что и в предыдущем случае.

На заметку
В описанных выше процедурах при сравнении параметра ТТ с to использована команда evalf (). Сделано это "на всякий случай" — разница (TT-t0) преобразуется в формат числа с плавающей точкой, поскольку это упрощает процедуру сравнения данного выражения с нулем.

Системы с колебаниями

Часто приходится рассматривать механические системы, которые, будучи выведены из состояния равновесия, возвращаются к нему. Если при этом силы, возвращающие систему в состояние равновесия, прямо пропорциональны "смещению" ее элементов относительно положения равновесия, то в таких системах имеют место гармонические колебания, которые при наличии диссипативных сил являются затухающими. Рассмотрим примеры.

Задача 6.4

Маятник состоит из жесткого стержня длины 1 и массы m на конце. К стержню прикреплены две пружины с жесткостью к на расстоянии а от точки подвеса. Найти частоты малых колебаний маятника. Массой стержня пренебречь.

В первую очередь создадим фафическую иллюстрацию для данной задачи. Поскольку Maple является командной средой, то и рисунок следует создавать с помощью команд — этим Maple отличается от фафических редакторов, где рисунки создаются путем непосредственного рисования.

Для начала опишем процедуру wall (), с помощью которой затем отобразим верхние и боковые стенки. Стенку представим в виде прямой линии со штрихами с внутренней стороны.

Параметры процедуры таковы: а) левая крайняя точка Р, с помощью которой будет осуществляться "привязка" отображаемой стенки к рисунку; б) длина стенки L; в) угол phi поворота стенки относительно базовой точки Р.

Далее в процедуре объявляются локальные переменные S, N, i и Res. Переменная N определяет число штрихов (в данном случае — 50). Переменная | S является массивом размерности N+1, а нумерация его элементов начинается |, с 0. Сами элементы массива S — это фафические объекты линия. Нулевым I элементом S[0] является основная линия, определяющая стенку. Линия эта создается процедурой line() из пакета plottools. Как известно, аргументами этой процедуры указываются начальная и конечная точки, которые соединяет прямая. Начальной является точка Р. Поскольку в первоначальном варианте прямая отображается горизонтально, по вертикальной оси координата конечной точки, по сравнению с начальной, измениться не должна, а вот горизонтальная координата должна увеличиться на длину линии, т.е. на L. Задается конечная точка следующим образом. Любая точка в Maple является списком из двух элементов: первый элемент — это горизонтальная координата, второй — вертикальная. Поэтому горизонтальная координата начальной точки определяется командой ор(Р)[1], вертикальная— командой ор(Р)[2]. Соответственно, конечная точка— это список [op(P)[l]+L,op(P)[2]]. После этого выполняется заполнение прочих элементов массива S. Эти элементы определяют "штрихи" — набор равноудаленных линий, направленных вовнутрь стенки под углом 45 фадусов к основной линии. Если число таких штрихов N, а длина стенки L, то горизонтальная координата начальной точки каждого "штриха", по сравнению с предыдущим, увеличивается на величину L/N и равна горизонтальной координате конечной точки следующего штриха. Вертикальная координата начальной точки штриха совпадает с координатой основной линии (первоначально она размещается по горизонтали), а вертикальная координата конечной точки, по сравнению с начальной точкой, больше на L/N. В соответствии с вышесказанным и формируются штрихи.

В процедуре предусмотрена возможность поворота всей описанной выше конструкции на один и тот же угол. Полученная в результате такого поворота структура будет присвоена в качестве значения переменной Res, которая инициализируется как пустой список. После этого, при последовательном переборе элементов массива S, эти элементы, повернутые на угол phi вокруг точки Р, добавляются в список Res. Поворот осуществляется процедурой rotate() из пакета plottools. В качестве первого, "разворачиваемого", аргумента указывается элемент S[i] (типа plot), где индексная переменная i пробегает значения в диапазоне от 0 до N. После этого указывается угол (phi), на который выполняется поворот (против часовой стрелки), и точка, вокруг которой следует поворот выполнять (Р).
После описания процедуры wall() инициализируем переменную Wall (последовательность), элементы которой — три стенки (две вертикальные и одна горизонтальная).

Все события будут происходить в квадрате размером 1x1. Поэтому левая вертикальная стенка имеет длину 1 и получается из горизонтальной поворотом вокруг точки [0,0] на угол 90 градусов (Pi/2) против часовой стрелки. "Потолок" получается, если горизонтальную линию перенести по вертикали на расстояние 1, т.е. базовой будет точка [0,1]. Последняя, правая стенка, получается поворотом на угол 90 градусов по часовой стрелке (-Pi/2) вокруг точки [1,1].

Для отображения шарика на стержне описывается процедура С(). Процедура имеет два параметра: длина стержня (L) и угол отклонения стержня (влево) от вертикали (alpha). Кроме того, в процедуре объявляется две локальные переменные С1 и С2. Первая переменная нужна для записи в нее объекта стержня (в качестве значения этой переменной присваивается соответствующий графический объект), а вторая — для шарика.

Стержень создается процедурой 'plottools/line'O в виде линии приемлемой толщины (опция thickness=3), где первым параметром (начальная точка) указана точка [0.5,1] (стержень подвешен по центру картинки), а вторым параметром (конечная точка)— точка [0.5-L*sin(alpha),l-L*cos(alpha)] — таковы координаты конца стержня длиной L при отклонении его влево на угол alpha. Эта же конечная точка — центр подвешенного на стержне шарика. Шарик создается с помощью команды seq('plottools/circle'([0.5-sin(alpha) ,1-L*cos(alpha)],0.004*i),i=l. .10), которая формирует последовательность из 10 окружностей с центром в упомянутой выше точке
' (0.5-L*sin(alpha),l-L*cos(alpha)] и радиусами, дискретно увеличивающимися (шаг 0.004) от 0.004 до 0.04. В результате выполнения процедуры формируется последовательность из двух элементов С1,С2, т.е. стержня и шарика (графические объекты).
Кроме шарика на стержне, нужно отобразить две пружины. Начнем с процедуры для выведения на экран одной пружины.
Пружина все время будет ориентирована по горизонтали, так что для корректного ее отображения следует знать два параметра: точку Р (левую) фиксации пружины и ее длину L.
Локальные переменные процедуры определяют число витков (N), половинную толщину (thick) пружины в абсолютных единицах (задавать толщину в процентах от длины пружины смысла не имеет, поскольку при сжатии (растягивании) пружины ее толщина будет уменьшаться (увеличиваться)), а также набор базовых точек для отображения пружины (S).

После инициализации переменной S в виде пустого списка, в этот список вносятся точки, формирующие пружину. Чтобы по базовым точкам сформировать пружину, используется процедура plot().
Теперь процедуру spring () используем для отображения двух соединенных друг с другом пружин. С этой целью определяем процедуру S(), параметрами которой являются высота фиксации пружин Н и угол alpha отклонения стержня, на котором подвешен шарик (этот угол нужен, чтобы определить длину каждой из пружин). Сформированные в процедуре графические объекты (левая и правая пружины) присваиваются в качестве значения локальным переменным S1 и S2 соответственно.
Левая пружина зафиксирована на левой стенке (горизонтальная координата равна 0) на высоте Н — это вертикальная координата. Если пружина расположена на высоте Н, расстояние между точкой подвеса стержня и точкой крепления пружины к стержню равно (1-Н) (в терминах задачи это параметр а). При отклонении стержня влево на угол alpha горизонтальная координата точки крепления пружины к стержню сместится в том же направлении на величину (1-H)*sin(alpha), и длина пружины будет равна 0.5-(l-H)*sin(alpha).

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

Что касается правой пружины, то здесь нужно учесть следующее. Согласно определению процедуры spring(), ее параметром указывается левая точка пружины. Поэтому для правой пружины точкой фиксации является точка крепления ее к стержню — и эта точка совпадает с точкой крепления со стержнем левой пружины, т.е. это точка [0.5-(l-H)*sin(alpha),H]. Длина же правой пружины увеличивается на величину, на которую уменьшается длина левой пружины. Эта величина равна 0.5+(l-H)*sin(alpha).

Результат выполнения процедуры S() формируется в виде последовательности.
Наконец, создаем процедуру Picture() для формирования всей картинки. У этой процедуры будет всего один аргумент — угол alpha отклонения стержня (влево). Высота крепления пружин и длина стержня определяются локальными переменными Н и L соответственно.

Картинка формируется процедурой 'plots/display'() с перечислением всех отображаемых объектов: трех стенок (переменная Wall), стержня и шарика (формируются командой С(L,alpha)), а также двух пружин (команда S(H,alpha)). Кроме того, чтобы масштабы по горизонтали и вертикали совпадали, указано значение опции scaling=constrained.
Теперь можно посмотреть, как система выглядит в состоянии равновесия (когда угол отклонения стержня равен 0).

После этого приступаем непосредственно к решению задачи. Заметим, что, поскольку в системе не действуют диссипативные силы, полная энергия системы не меняется. Полная энергия системы состоит из:

а) кинетической энергии шарика; б) потенциальной энергии деформации пружин; в) потенциальной энергии шарика в гравитационном поле.

Начнем с кинетической энергии, которая определяется через модуль скорости шарика. Его скорость, очевидно, направлена по касательной к описываемой шариком траектории (дуга окружности) и определяется через производную по времени от угла отклонения.

Энергию системы определим как оператор. В частности, оператор кинетической энергии шарика, зависящий от функции alpha (эта функция задает эволюцию системы во времени), действует на параметр t следующим образом.

Внимание!
В приведенном выше описании конструкция Ekin:=alpha-> означает, что оператору Ekin в зависимости от параметра alpha ставится в соответствие то, что находится после стрелки (->). После стрелки находится код t->m*(l*D(alpha)(t))A2/2, определяющий действие Ekin на аргумент t. Таким образом показано, что оператору Ekin в соответствие ставится действие (поэтому Ekin и является оператором).

Практически так же определяется потенциальная энергия сжатой и разжатой пружин.

На заметку
Если пружина сжимается или разжимается на величину д: , то потенциальная энергия деформации пружины равна Е = кх 2/2, где к — жесткость пружины. Если расстояние от точки подвеса стержня до точки крепления пружины равно а , а стержень отклонился на угол а , то при малых углах отклонения можем полагать, что сжатие одной пружины и, разумеется, растяжение другой одинаковы и равны х = аа{(). Поскольку пружин две и каждая дает свой вклад в энергию деформации, полная энергия деформации пружин равна Е = ка2а2. Что касается потенциальной энергии шарика в поле силы тяжести, определяемой ниже, то удобно нулевой уровень для такой энергии выбрать на высоте, где щврик и стержень находятся в равновесии. Тогда при отклонении на угол а стержня шарик поднимется, по сравнению с первоначальным положением, на величину h*l(\-cos(a)) и его потенциальная энергия будет равна Е = mgl(\-cos(a(/))).

Прирост потенциальной энергии шарика в гравитационном поле при отклонении стержня от вертикали дается следующей зависимостью.

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

Следует заметить, что полная энергия Etot определена как переменная, а не как оператор.

Как уже отмечалось, полная энергия системы (переменная Etot) при эволюции последней не меняется. Это значит, что если продифференцировать выражение для полной энергии по времени, то эта производная с неизбежностью должна равняться нулю. Таким образом, получаем уравнение движения.

Учитывая, что в задаче ищется нетривиальное решение, можно сократить полученное уравнение на производную.

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

Замена осуществляется с помощью процедуры subs(). В частности, указано, что в уравнении Eq_l выражение sin(alpha(t)) следует заменить на результат преобразования в полином (процедура convert () с опцией polynom) разложения в ряд Тейлора выражения sin (alpha) в окрестности нуля (опция alpha=0) с остатком ряда второго порядка по аргументу. При этом сразу после процедуры convert () указан заключенный в скобки параметр t. Дело в том, что угол alpha, по которому выполняется разложение в ряд, является функцией времени. Использованная конструкция реализуется по следующей схеме: вычисляется разложение в ряд по alpha, затем ряд преобразуется в полином, который действует как оператор на параметр t. Другими словами, в данной записи alpha интерпретируется вычислительным ядром Maple как оператор.

После этого полученное уравнение Eq_2 упрощаем, разделив правую и левую части на (ш*1Л2) и сгруппировав слагаемые при alpha(t).

Решить это дифференциальное уравнение особого труда не представляет.

Переменные среды _С1 и _С2 определяются из начальных условий. Эти начальные условия можно было сразу указать в процедуре dsolve() еще при решении уравнения. Однако в данном случае поступим иначе. Так, переменной ехрг в качестве значения присвоим полученную при решении дифференциального уравнения Eq 2 временную зависимость угла отклонения стержня.

Соответственно, угловая скорость стержня определяется через производную по времени от переменной ехрг.

Теперь временной переменной t присвоим начальное значение, т.е. 0. t:=0;

Если в начальный момент угол равен А, а скорость (угловая) равна v, то переменные С1 и _С2 можно найти как решение системы уравнений.

После этого в переменной ехрг заменяем переменные среды _С1 и _С2 на полученные выше значения для них.

Внимание!
Первым параметром процедуры subs() следует указать равенства, согласно которым выполняется замена в выражении, указанном вторым параметром. В данном случае первым параметром является переменная среды 1%, значение которой равнорезультату предпоследней операции, т.е. это множество, полученное в результате выполнения команды, помеченной звездочкой (*). В соответствии с равенствами, указанными в качестве элементов этого множества, и выполняется замена.

В следующей задаче вопрос о гармонических колебаниях разрешается несколько иначе.

Задача 6.5

Тело массы М, соединенное с пружиной жесткости к, другой конец которой закреплен неподвижно, может двигаться без трения по горизонтальной плоскости. К телу прикреплен математический маятник массы m и длины 1. Найти функцию Лагранжа системы и частоты малых колебаний.

Сразу подключаем нужные пакеты.

Внимание!
Появление сообщений с предупреждениями после подключения пакетов объясняется следующим образом. После подключения пакета, как известно, для вызова любой проце дуры из этого пакета необходимости указывать принадлежность последней к пакету нет — и это удобно. Однако иногда подключается сразу несколько пакетов, в которых имеются процедуры с одинаковыми названиями. В таких случаях по умолчанию при соот ветствующем вызове используется процедура из пакета, подключенного последним. Например, стандартная процедура Maple (доступная без подключения каких бы то ни было пакетов) changecoords)) выполняет преобразование указанного в качестве первого ее аргумента выражения, записанного в декартовых координатах, в переменные других координатных систем (необходимая координатная система также указывается как аргумент процедуры changecoords ()). Процедура plots[changecoords ] () (процедура с таким же названием, но из пакета plots) предназначена для отображения графических структур в системах координат, альтернативных декартовой. При подключении пакета plots исходная процедура changecoords () переопределяется, о чем и выводится сообщение. Похожая ситуация имеет место и с процедурой arrow)). Что процедура arrow() из пакета plots, что из пакета plottools — обе предназначены для отображения стрелок Однако, в отличие процедуры plots [arrow] (), которая сразу выводит картинку со стрелкой, процедура plottools [arrow] () только формирует соответствующий объект, и для его отображения необходимо использовать процедуру display)). Ниже показано, как можно использовать указанные процедуры. Например, можно сначала подключить пакет plots.

В появляющемся в данном случае предупреждении упоминается только процедура changecoords)), как и должно быть. Теперь можно отобразить три стрелки единичной длиНы (параметр {[0,0,1], [0,1,0], [1,0,0]}), размещенные вдоль координатных ос центром в начале координат (параметр [0,0,0]).



На следующем этапе подключаем пакет plottools.

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

В Maple 9 сообщение будет таким: Error, (in arrow) expecting at least 5 arguments, but got 4 — (Ошибка, (в arrow) ожидается не меньше 5 аргументов, а их Л). Однако от этого мало что меняется. Правильным теперь будет, например, такой вызов.



В пакете plottools в процедуре arrow() в качестве аргументов указываются, кроме прочего, начальная ([0,0]) и конечная ([1,1]) точки, толщина линии (.05), толщина "наконечника" (. 15) и его длина (.2). Как уже отмечалось, чтобы стрелка отображалась, нужно использовать процедуру display{). Вызов только процедуры arrow() сформирует графический объект, но не отобразит его.



После того как подключен пакет plottools, процедуру arrow () из пакета plots вызвать все же можно, но для этого следует явно указать пакет.

Результат выполнения этой команды будет таким же, как и в самом первом случае.

Определим процедуру body() для отображения тела (бруска) на горизонтальной плоскости.

Параметр процедуры х задает горизонтальную координату центра симметрии бруска относительно граничной левой стенки. В этой процедуре используются две глобальные переменные Н и L, определяющие высоту и ширину бруска соответственно.

Следующей командой формируется изображение пружины, которое затем присваивается переменной s.

Основная линия вертикальной стенки начинается в точке [0,0] и заканчивается в точке [0,Sft], т.е. имеет длину Sft — и этот объект в качестве значения присваивается локальной переменной А.
Штрихи формируются как массив S, элементами которого являются тонкие (thickness=0) линии, направленные под углом 135 градусов к вертикали, т.е. внутрь стенки и вниз. Вертикальная координата окончания каждого штриха (лежащего на основной линии) на величину Sf t/N больше, чем у предыдущего.
Вертикальная стенка создается процедурой display)), первый аргумент которой — базовая линия стенки А, а второй — последовательность, формируемая из элементов массива S. Фрагмент горизонтальной плоскости (переменная hWall) получается из вертикальной стенки путем отражения относительно прямой, проходящей через точки [0,0] и [1,1]. Осуществляется отражение посредством процедуры reflect() из пакета plottools. Параметрами этой процедуры указывают графическую структуру, которую следует отражать (hWall), и точку или линию, относительно которой выполняется отражение. Прямая линия задается двумя точками.
Полученная таким образом часть горизонтальной плоскости имеет длину Sft, поэтому ее нужно "удлинить". Делается это посредством параллельного переноса созданного фрагмента вдоль горизонтали. Такая операция реализуется в рамках условного оператора while.

На заметку
Синтаксис вызова условного оператора выглядит следующим образом: while условие do опрератор_1 ... оператор_Я end do. Если выполняется условие, то последовательно выполняются операторы, размещенные между ключевыми фразами do и end do, после чего снова проверяется условие, и если оно верно, то снова выполняются операторы, и т.д.

Перед вызовом этого условного оператора инициализируется переменная-счетчик R, с помощью которой можно корректно определить, когда следует закончить выполнение условного оператора. Переменная эта будет "помнить" длину фрагмента горизонтальной плоскости. Первоначальное значение этой переменной, как несложно догадаться, должно равняться Sft. Картинку будем создавать такой, чтобы размер по горизонтали (и, в принципе, по вертикали) был не менее 1, но не более 2. Поэтому процесс "дополнения" горизонтальной плоскости продолжается до тех пор, пока R<1. Именно это условие и указано в условном операторе while. В теле оператора размещено всего две команды. Первая "удлиняет" горизонтальную плоскость. Делается это следующим образом: переменной hWall в качестве значения присваивается последовательность, первый элемент которой — сама переменная hWall (т.е. ее прежнее значение), а второй элемент — результат выполнения процедуры translate(). У этой процедуры три параметра-аргумента: первый задает графическую структуру, которую следует "перенести", а два других определяют направление переноса. Так, в данном случае вторым параметром указана переменная Sft, а третьим — 0. Это значит, что графическую структуру следует перенести на расстояние Sft вдоль горизонтальной оси и на расстояние 0 по вертикали. Первый аргумент процедуры translate () — результат отражения вертикальной стенки относительно биссектрисы левого нижнего угла.

Вторая команда в цикле while нужна для изменения значения переменной-счетчика R — после описанной выше манипуляции длина горизонтальной плоскости становится больше на величину Sft.

На следующем этапе формируется полная картинка. Выполняется это с помощью процедуры display(), и результат присваивается в качестве значения переменной Wall.

Наконец, вся картинка переносится вверх на расстояние (1-Sft) и по горизонтали на расстояние длины штриха, т.е. на величину Sft/N. Такой перенос выполняется с очень простой целью — чтобы координаты по вертикали и горизонтали были положительными. Если картинку не смещать, принципи-[, ально ничего не изменится — внешний вид системы будет таким же, разница будет только в выборе начала отсчета.

На заметку
В данном примере координатная система не отображается. Поэтому, строго говоря, картинку можно было и не смещать.

Используя последнюю описанную процедуру, можно отобразить внешний вид рассматриваемой системы. Однако для этого следует инициализировать ряд параметров. Так, высоту бруска примем равной

Теперь можно отобразить всю систему. Координату центральной точки бруска примем равной 0.5, а угол отклонения стержня — я/8.

Внимание!
Процедурой wall_and_all() картинка формируется, но не отображается! Для ее отображения следует использовать процедуру display!) из пакета plots.

На заметку
Если выделить рисунок и из раскрывающегося меню выбрать подменю Axes (Оси), а затем одну из команд выбора системы координат (любая, кроме None (Отсутствуют): Boxed (В рамке), Frame (Точка пересечения в левом углу) или Normal (Обычные)), на рисунке будут отображены координатные оси. В этом случае несложно заметить, что центральная точка бруска расположена несколько правее метки 0.5. Дело в том, что координатацентральной точки бруска равна 0.5 до переноса картинки на (1-Sft) вверх и Sft/N вправо. После переноса вправо координата центральной точки равна 0.5+Sft/N.

Далее приступаем непосредственно к решению задачи. Но прежде следуе "сбросить" значение, присвоенное переменной 1. Эта переменная использу ется непосредственно в решении, а решение необходимо получить в сим вольном виде. Поэтому выполняем следующую команду.

Для определения функции Лагранжа нужно задать кинетическую и потен циальную энергии системы. Кинетической энергией обладают движущийс: брусок и совершающий колебания шарик (стержень и пружину считаем неве сомыми, поэтому кинетической энергией они не обладают).
Кинетическая энергия Т1 () бруска равна половинному произведенш массы бруска М на его скорость v. Определяем данную зависимость ка функцию от скорости.

Кинетическая энергия шарика определяется не так просто, поскольку з; кон его движения есть суперпозиция движения центра бруска и непосредс венно шарика относительно бруска — ведь колебания шарик совершает отн< сительно бруска, а относительно неподвижной системы координат движеш шарика намного сложнее. В общем случае кинетическая энергия шарика Т2 зависит от трех параметров: скорости бруска v, угловой скорости шарш omega в системе координат, связанной с бруском, и угла phi отклонен стержня относительно вертикали.

На заметку
Относительная линейная скорость колебательных движений шарика равна произведению угловой скорости шарика на длину стержня и направлена по касательной к описываемой шари ком траектории, т.е. перпендикулярно стержню. Из элементарных геометрических соображений очевидно, что проекция этой скорости на горизонтальную ось дается умножением на коа нус угла отклонения стержня от вертикали, а на вертикальную — умножением на синус. Чтоб1 найти абсолютную скорость шарика, нужно сложить соответствующие проекции скорости бру ска и относительной скорости шарика на горизонталь и вертикаль. Поскольку брусок движете вдоль горизонтали, проекция его скорости на горизонталь совпадает с этой скоростью, а вертикальная проекция равна нулю. В выражении для кинетической энергии используется квадрат абсолютной скорости. Как известно, квадрат вектора равен сумме квадратов его проекци! на перпендикулярные оси. Этим свойством и воспользуемся ниже.

Таким образом, кинетическая энергия шарика определяется как функцу трех параметров.

Кинетическая энергия системы Т() равна сумме кинетических энергий бруска и шарика и также является функцией трех параметров: скорости бруска v, угловой скорости шарика omega и угла phi отклонения стержня от вертикали.

Далее определяем потенциальную энергию системы, которая имеет две составляющие: энергию сжатия пружины VI () и энергию шарика в поле тяжести V2(). Энергия деформации пружины определяется смещением пружины х относительно положения равновесия.

При отклонении стержня на угол phi приращение потенциальной энергии шарика, как известно, определяется следующим образом.

Суммарная потенциальная энергия зависит как от смещения бруска х относительно положения равновесия, так и от угла отклонения стержня phi относительно вертикали.

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

Для дальнейшего анализа удобно присвоить выражение для функции Ла-гранжа системы в качестве значения переменной Lg. При этом обобщенными координатами являются смещение х бруска относительно положения равновесия и угол отклонения стержня phi от вертикали. Соответственно, скорость бруска v является производной по времени от х, а угловая скорость omega — производной по времени от phi.

Чтобы по функции Лагранжа определить уравнения движения системы, необходимо вычислить частные производные от функции Лагранжа по каждой из ее переменных — всего четыре выражения. Каждое такое выражение присвоим в качестве значения элементам таблицы А, индексы которых будут обозначать переменную, по которой вычислялась производная.
Дальнейшая процедура подразумевает вычисление производных по времени от тех выражений (элементов таблицы А), которые являются производными по скоростям от функции Лагранжа (элементы Av и Ат). В самих же
уравнениях следует учесть, что первые два параметра в функции Лагранжа являются производными по времени от двух последних, и, разумеется, все они — функции времени. Поэтому после вычисления частных производных сразу же в полученных выражениях делаем ряд замен: указываем явно, что х и phi зависят от времени, a v и omega есть производные от х и phi соответственно и также зависят от времени.
Реализуется такая замена ниже с помощью оператора цикла for, где переменная s пробегает значения из набора v, omega, x и phi. Элементу с индексом s присваивается в качестве значения результат дифференцирования функции Лагранжа по переменной s (второй параметр процедуры subs () — команда diff (Lg,s)). При этом сразу указано, что координаты следует считать функцией времени, а скорости выражены через производные от координат. Замена осуществляется процедурой subs ().

Последнее уравнение можно существенно упростить.

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

Упрощенную таким образом систему уравнений можно решить в общем виде. Ниже на этот случай приведена соответствующая команда. Результат ее выполнения, из которого можно определить частоты колебаний, не приводится в силу его исключительной громоздкости.

Однако для определения частот систему дифференциальных уравнений можно и не решать.
Для начала сделаем в уравнениях замену: вместо функции x(t) будем пользовать y(t), которая связана со "старыми" координатами х и phi со ношением y(t)=x(t)+l*phi(t). Фактически y(t) — это координата (гориэ тальная) шарика в неподвижной системе координат.
Таким образом, получаем следующее.

Это уравнение неплохо было бы упростить.

Из последнего уравнения, как нетрудно заметить, можно выразить у) отклонения. Сделаем это (команда solve(Eq2new,phi(t))), и полученное выражение для функции phi(t) подставим в первое уравнение.

Упрощаем результат.

Составим характеристическое уравнение. Поэтому в уравнении Eqn выполняем соответствующую замену, после чего уравнение сокращаем на экспоненту (еще умножаем на д) и упрощаем.

Полученное уравнение определяет частоты собственных колебаний системы. Эти частоты (точнее, частоты в квадрате) можно найти.

Чтобы нагляднее представить, как выглядит система "в действии", создадим анимационную картинку. Но прежде присвоим параметрам, используемым в задаче, конкретные значения.

Для отображения динамики системы во времени необходимо знать закон движения. Для этого решим исходную систему дифференциальных уравнений (Eql, Eq2). Решение будем искать в численном виде, для чего результат выполнения процедуры dsolve() указывается аргументом процедуры преобразования в формат числа с плавающей точкой evalf ().
Система уравнений решается для случая, когда в начальный момент времени брусок отклоняется от положения равновесия на расстояние 0.2, а стержень — на угол я/16. Начальные скорости шарика и бруска равны нулю.

Внимание!
Указывать индексы в последних двух командах следует крайне осторожно! От сеанса к сеансу решения системы дифференциальных уравнений выводятся в разной последовательности, поэтому приходится соответствующим образом менять и индексы.

После этого формируем кадры анимации. Динамику системы будем от слеживать на протяжении 2,5 секунд с интервалом 0,1 секунды. Каждый кад{ будем записывать в массив о. Ниже приведен соответствующий код.

Момент времени последовательно, с шагом 0.1, изменяется от 0 до 2. Положение равновесия бруска взято равным 0.5, а чтобы было лучше ориеь тироваться в ситуации, отображается, кроме картинки, еще и момент времни (процедура textplot()). Отображаемый текст формируется процедуре объединения cat().
Наконец, в процедуре display!), помимо первого аргумента — списка (эт важно!) кадров, указывается еще и опция insequence=true.

В данной главе для отображения анимации использовалась процедура dis-play(). В отличие от процедуры animate(), рассмотренной в предыдущей главе, процедура display!) более удобна в том смысле, что позволяет отображать последовательность, по большому счету, произвольных картинок. Например, если при отображении динамики системы с помощью процедуры animate () интервал времени разбивается на равные промежутки, то при использовании процедуры display() имеется возможность отображать состояние системы через промежутки времени разной протяженности. Это иногда бывает полезно.

Заключительные замечания

Рассмотренные в данной главе темы служат не столько иллюстрацией методов решения физических задач, сколько являются примером того, как в Maple можно графически интерпретировать получаемые результаты. С другой стороны, не следует воспринимать то, о чем шла речь, как некий путеводитель по графическим утилитам системы Maple. Это иллюстрация, всего лишь один из возможных способов применения системы в решении прикладных задач. Возможности Maple несоизмеримо шире. Желающие могут в этом убедиться сами.

Еще одной темой, важность которой не следует недооценивать, являются возможности Maple в сфере численно-аналитических методов. Как будет показано в следующей главе, Maple и здесь не уступит любому доступному на сегодня инженерному пакету.

Контрольные вопросы

1. Какая структура определяется командой <а|Ь|с>: а) команда некорректна; б) вектор с компонентами а, Ь и с; в) последовательность с компонентами a, b и с; г) сразу три переменные: a, b и с. 2. Процедура spacecurve () используется для: а) построения окружностей; б) построения линий; в) построения параметрически заданных кривых в трехмерном пространстве; г) построения прямоугольников. 3. Процедура textplot() используется для: а) создания заголовков графиков; б) отображения текстовых полей; в) отображения легенд графиков; г) определения названий координатных осей. 4. Какие из приведенных ниже высказываний относительно процедуры display! ) верны: а) процедура нужна для отображения графических структур; б) процедура нужна для отображения результатов выполнения команд пользователя; в) это эквивалент процедуры plot(); г) это эквивалент процедуры animate (). 5. "Повернуть" фафический объект можно с помощью процедуры: а) reflect(); б) changecoords(); в) translate(); г) rotate().