Математический анализ в Maple 9

мельник салат


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

При решении задач динамики, как правило, приходится интегрировать уравнения движения материальной точки или системы тел. Эти уравнения являются дифференциальными уравнениями второго порядка, и для их решения широко используются методы, рассмотренные в главе 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)).

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

Два следующих поля принципиально мало чем отличаются от 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) преобразуется в формат числа с плавающей точкой, поскольку это упрощает процедуру сравнения данного выражения с нулем.