Maple 9.5/10 в математике, физике и образовании - Владимир Дьяконов

Рис. 11.12. Пример вычисления временной зависимости свободных колебаний и построения их графика

Нередко о характере колебаний удобно судить по фазовому портрету колебаний. Он задается графиком в параметрической форме, при которой по одной оси откладывается зависимость у(t), а по другой — ее производная. Это показано на рис. 11.13. Фазовый портрет в данном случае представляет собой сворачивающуюся спираль.

Рис. 11.13. Фазовый портрет затухающих свободных колебаний

11.1.6. Улучшенное моделирование колебаний при синусоидальном воздействии

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

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

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

Рис. 11.15. Результаты моделирования цепи второго порядка при синусоидальном воздействии

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

11.1.7. Улучшенное моделирование колебаний при пилообразном воздействии

Рассмотрим методику улучшенного моделирования еще на одном примере — вычислении реакции системы при пилообразном воздействии. На рис. 11.16 показано задание такого воздействия с помощью функции floor. Для упрощения расчетных выражений амплитуда и период воздействия взяты равными я. Поскольку в данном случае аналитическое решение получить невозможно (функция floor не позволяет этого), то заменим воздействие рядом Фурье. Его коэффициенты также представлены на рис. 11.16.

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

На рис. 11.17 представлены графики воздействия в идеальном случае и при его представлении рядом Фурье с пятью гармоники. Показано также аналитическое решение для временной зависимости y(t) при таком воздействии.

Рис. 11.17. Воздействие и временная зависимость реакции системы при пилообразной форме воздействия

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

Рис. 11.18. Реакция системы на пилообразное воздействие и фазовый портрет колебании при таком воздействии

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

11.1.8. Анализ и моделирование линейных систем операторным методом

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

Для начала зададим инициализацию применяемых пакетов расширения

> restart:with(plots): readlib(spline): with(inttrans):

Warning, the name changecoords has been redefined

Далее зададим операторные выражения для коэффициента передачи системы G и входного сигнала R (в виде единичного перепада) и вычислим с упрощением их произведение:

> G := K/(M*s^2+C*s+1); R := 1/s;

> X := simplify(R*G);

Теперь, используя обратное преобразование Лапласа, найдем временную зависимость реакции системы в аналитическом (что наиболее ценно) виде:

> h := simplify(invlaplace(X, s, t));

Теперь мы можем построить график этой зависимости для конкретных значений М, С и K:

> h1 : = subs(M=1,C=0.75,K=1,h);

h1 := 0.5393598900 I (-1.854019622 I + 1.854049622 I e(-0.3750000000 t) cosh(0.9270248110 I t) + 0.75 e(-0.3750000000 t) sinh(0.9270248110 I t))

> linresp := plot(h1, t=0..20, axes=boxed, color=black): display(linresp);

Вид этой зависимости представлен на рис. 11.19. Он соответствует реакции системы второго порядка для случая затухающих колебаний.

Рис. 11.19. Одна из временных зависимостей реакции системы второго порядка

А теперь зададимся целью наглядно проиллюстрировать изменение временной зависимости реакции системы при изменении параметра С от 0 до 2 при М=1 и K=1. Для этого выполним следующие вполне очевидные команды:

> х := subs(M=1, K=1, h);

> plot3d(x, С=0..2, t=0..20, axes=boxed);

Соответствующий график показан на рис. 11.20. Он прекрасно иллюстрирует переход от апериодического режима при С=2 к колебательному при С= 0 при изменении времени от 0 до 20.

Рис. 11.20. Динамика развития колебаний в системе при изменении параметра С

Аналогичным образом можно построить трехмерный образ временной зависимости реакции системы для М=1, С=0.25 и изменении параметра K от 0 до 3. Для этого надо выполнить команды:

> x1 := subs(М=1, С=0.25, xt);

x1 := 0.5039526307 I K (-1.984313483 I + 1.984313483 I e(-0.1250000000 t) cosh(0.9921567415 I t) + 0.25 е(-0.1250000000 t) sinh(0.9921567415 I t))

> plot3d(x1, K=0..3, t=0..20, axes=boxed);

Диаграмма временных зависимостей представлена на рис. 11.21.

Рис. 11.21. Динамика развития колебаний в системе при изменении параметра K

Представленные на рис. 11.20 и 11.21 диаграммы дают весьма наглядное представление о динамике поведения рассмотренной системы. Но еще важнее то, что просто изменением операторной записи G и R по описанной методике можно анализировать и наглядно представлять работу множества линейных систем.

11.2. Моделирование динамических задач и систем

11.2.1. Расчет траектории камня с учетом сопротивления воздуха

Вы хотите метнуть камень в огород вашего вредного соседа? Разумеется во время его отсутствия. Давайте промоделируем эту ситуацию, предположив два актуальных случая: дело происходит на Земле в условиях, когда наша планета лишилась воздуха и когда, слава богу, он все же есть. В первом случае сопротивления воздуха нет, а в другом сопротивление воздуха есть и его надо учитывать. Иначе камень упадет в ваш огород, а не в огород соседа!

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

Итак, пусть подвернувшиеся под руку камни с массой 500 и 100 грамм брошены под углом 45 к горизонту со скоростью VO=20 м/с. Найдем их баллистические траектории, если сила сопротивления воздуха Fmp=A*V, где А=0,1 Н∙с/м. Сравним их с траекториями, получающейся без учета сопротивления воздуха. Документ с решением этой задачи, описанным ниже, представлен в файле balist.

Начнем с подключения пакета plots, нужного для визуализации данной задачи:

> restart; with(plots):

Warning, the name changecoords has been redefined

Составим параметрические уравнения для проекций скорости на оси координат:

> Vox:=Vc*cos(alpha);Voy:=Vo*sin(alpha);

Vox:= Vo cos(α) Voy:= Vo sin(α)

Мы рассматриваем два случая: камень массой 500 г и камень массой 100 г. Поскольку для каждого случая мы предусматриваем расчет в двух вариантах (с учетом сопротивления воздуха и без такого учета), то мы должны составить 4 системы дифференциальных уравнений (ДУ). Каждая система состоит из двух ДУ второго порядка и вид этих систем известен из курса физики. Ниже представлено задание этих систем ДУ (для первой системы дан вывод ее вида):