Система амортизации с двумя степенями свободы

Размещено на http://www.allbest.ru/

Оглавление

Введение

1. Постановка задачи

1.1 Техническая постановка

1.2 Обзор задач-прототипов

1.3 Математическая модель

2. Методы решения

2.1 Одношаговые методы

2.1.1 Методы Эйлера

2.1.2 Методы Рунге - Кутта

2.2 Многошаговые методы

2.3 Реализация методов решения систем нелинейных дифференциальных уравнений в пакете MATLAB

3. Анализ несимметричной нелинейной системы амортизации

3.1 Реализация решения задачи в пакете MATLAB

3.2 Исследование эффектов и поведения решения при различных значениях параметров и анализ полученных результатов

3.2.1 Исследование АЧХ в зависимости от значений расстояний от центра масс до точек прикрепления пружин

3.2.2 Исследование АЧХ в зависимости от значений жесткостей пружин

3.2.3 Исследование АЧХ в зависимости от значений диссипации

Заключение

Список использованной литературы



Введение

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

Исследование виброзащиты подразумевает анализ соответствующих математических моделей колебательных систем. Таким образом, необходимо также рассмотрение самой теории колебаний.

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

Любой колебательный процесс описывается с помощью математических уравнений или иначе математической моделью. Математической моделью колебаний являются обыкновенные дифференциальные уравнения. Различные типы колебаний и соответствующие им математические модели будут рассмотрены дальше.

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

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

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

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

Если говорить другими словами, виброизолятор способствует комфортному и плавному снижению амплитуды колебаний, а динамический гаситель поглощает энергию колебаний. Естественно, наиболее эффективным является рассмотрение объектов, которые будут включать в себя свойства гашения и изолирования. Так, например, у легкового автомобиля подвеска, а точнее ее стойка, состоит из амортизатора-гасителя и пружины-изолятора.

В транспортных задачах довольно-таки часто используются объекты «балочного» типа, а значит рассмотрение системы, состоящей из твердого тела с упругими опорами, играет немаловажную роль для инженерных работ.

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

Целью работы является построение математической модели и исследование несимметричной нелинейной системы амортизации.

Поставленные в данной работе задачи:

Построить математическую модель заданной системы;

Выбрать наиболее оптимальный метод решения модели;

Реализовать решение на программном уровне;

Исследовать поведение системы при различных значениях параметров;

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

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

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

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



1. Постановка задачи

1.1 Техническая постановка

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

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

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

Данная задача была исследована в другой постановке - симметричном относительно центра масс варианте (одинаковыми были расстояния от центра масс до мест приложения пружин d1=d2, коэффициенты упругости пружин c1=c2). Как известно, математически движение объектов с двумя степенями свободы описывается системой двух дифференциальных уравнений второго порядка с двумя неизвестными - обобщенными координатами. Остальные составляющие дифференциальных уравнений определяются динамикой самого объекта. Меняя физические величины системы, можно найти наиболее оптимальный вариант ее построения. Оптимальные величины мы будем выбирать такие, чтобы достичь поставленной цели, а именно, предоставить эффективную модель устройства виброзащиты.

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

Для исследования заданной системы необходимо построить ее математическую модель, с помощью математической модели построить амплитудно-частотную характеристику (далее АЧХ) системы, рассмотреть поведение АЧХ при различных значениях входных параметров, рассмотреть возможные во всех случаях эффекты, выяснить, какие параметры необходимо взять для снижения амплитуды колебаний.

1.2 Обзор задач-прототипов

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

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

По способам возникновения колебания бывают:

- свободные

- автоколебания

- вынужденные

- параметрические

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

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

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

Параметрические колебания возникают при периодическом изменении некоторого параметра системы. Частота колебаний определяется собственными параметрами системы, а соотношение собственной частоты системы и частоты изменения параметров системы определяет амплитуду колебаний. Всеми известный пример параметрических колебаний - качели.

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

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

В зависимости от вида уравнения модели могут быть:

- линейные (описываемые линейными дифференциальными уравнениями)

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

-существенно нелинейные (описываемые нелинейными дифференциальными уравнениями)

На самом простом примере рассмотрим все виды колебательных процессов, а также их математические модели.

Рис. 1. Математический маятник

Уравнение движения математического маятника выглядит следующим образом:



или,

где m - масса груза, l-длина нити, - частота системы, g-ускорение свободного падения, - угол отклонения нити от положения равновесия.

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

- линейной (пример линейного гармонического осциллятора)

-квазилинейной

- нелинейной

Это уравнения, описывающие свободные колебания или собственные (без учета диссипации).

Ограничимся линейным уравнением и рассмотрим оставшиеся способы возникновения колебательных процессов. Если к системе будет приложена внешняя сила F, которая не будет иметь колебательный характер, тогда будут автоколебания, уравнение будет выглядеть следующим образом:

Если к системе будет приложена внешняя сила, имеющая характер колебаний, то есть F=Wsinwt (где W-амплитуда колебаний приложенной силы, w-частота), колебания будут вынужденные, уравнение будет выглядеть:

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

Все уравнения рассмотрены без учета диссипации. Если мы для самого первого случая введем диссипацию (коэффициент диссипации - p), тогда колебания будут свободными и затухающими, уравнение буде выглядеть следующим образом:

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

Нашу систему можно разбить на более простые случаи:

-Гармонические колебания

- Затухающие линейные колебания



- Затухающие существенно нелинейные параметрические колебания

1.3 Математическая модель

Рис. 2. Иллюстрация нелинейной несимметричной системы амортизации с двумя степенями свободы

Использованные обозначения:

- с - центр масс;

- - вынуждающая сила;

- - амплитуда колебаний;

-w - частота колебаний;

- ц - перемещение вокруг центра масс, зависит от t;

- x- перемещение вдоль вертикальной оси, зависит от t;

- - расстояние от центра масс с до пружин i, ;

- - зависимость усилия от перемещения в пружине i,(s)=;

- - коэффициенты упругости, ;

- -коэффициенты при нелинейном члене;

- -коэффициенты диссипации.

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

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

- Движение центра масс вдоль вертикальной оси х.

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

- Вращение вокруг центра масс:

Во всех следующих формулах = ,=.

Получаем систему:



К системе учтем вид нелинейности F:

Раскроем скобки и приведем подобные :

Получим систему относительно x и :

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

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



2. Методы решения

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

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

или

Для второго уравнения решается задача Коши, то есть ищется частное решение уравнения, удовлетворяющее начальным условиям. Решением такого уравнения является пара x и y(x). При их подстановке уравнение должно обращаться в тождество.

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

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

Методы решения обыкновенных дифференциальных уравнений могут быть одношаговыми (одноступенчатыми) и многошаговыми (многоступенчатыми).

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

2.1 Одношаговые методы

2.1.1 Методы Эйлера

Метод Эйлера является методом Рунге-Кутта 1-го порядка. Этот метод является самым простым и наиболее распространенным в решении задач Коши, а также имеет различные модификации, поэтому его можно рассматривать отдельно. Дифференциальное уравнение первого порядка можно рассматривать как определение кривой через ее производную в плоскости (x,y). Любое дифференциальное уравнение задает наклон кривой как функцию от x и от y в любой точке. В начальный момент известна одна точка, через которую проходит кривая, а именно .Начиная с этой точки, вычисляем наклон кривой при и , продвигаемся на некоторое малое расстояние вдоль получившейся касательной. Шаг по х=h. Получаем и . Продолжая процедуру, получаем последовательность коротких отрезков прямой, которые являются достаточно хорошим приближением к искомой функции. Явный.

Если записывать явный метод Эйлера формулами, это будет выглядеть так:



, i=0,…,n-1

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

, i = 0, …, n-1

Метод Эйлера имеет ошибку пропорциональную h, неявный метод Эйлера имеет ошибку, пропорциональную

2.1.2 Методы Рунге-Кутта

Методы Рунге-Кутта согласуются с рядом Тейлора вплоть до членов порядка , где p-порядок метода. Требуют большого количества итераций, но не требуют вычисления производных, только вычисления самой функции.

Идея построения явных методов Рунге-Кутта p-го порядка заключается в получении приближения к значениям по формуле , где - некоторая функция, приближающая отрезок ряда Тейлора до р-го порядка и не содержащая частных производных функции . Если =, то получаем метод Эйлера. Для построения методов Рунге-Кутта высшего порядка функцию берут многопараметрическую и подбирают ее параметры сравнением с многочленом Тейлора для у(x), соответствующей желаемому порядку степени.

Метод Рунге-Кутта четвертого порядка выглядит следующим образом:

Ошибка метода пропорциональна .

2.2 Многошаговые методы

Одним из самых распространенных многоступенчатых методов является метод Адамса. Он заключается в том, что мы предполагаем уже найденными несколько значений решения на равномерной сетке . Формулы метода Адамса:

,

где, - интерполяционные полиномы Ньютона k-го порядка, интерполяция происходит из узла и узла соответственно.

Для получения экстраполяционного метода Адамса-Башфорта, необходимо сделать замену переменной , в итоге получим формулу:

,

где =, а - интерполяционный полином Ньютона, интерполяция происходит из узла .

Для получения интерполяционного метода Адамса-Моултона, необходимо аналогичным образом сделать другую замену переменной , в итоге получим формулу:

,

где =, а - интерполяционный полином Ньютона, интерполяция происходит из узла .

Точность многошаговых методов растет вместе с порядком порождаемых интерполяционных многочленов Ньютона, метод, порождаемый интерполяционным многочленом k-ой степени, является методом k+1 порядка точности.

Сравним метод Адамса с методом Рунге-Кутта того же порядка точности. Несмотря на то, что в методе Рунге-Кутта требуется вычисление 4 значений на каждом шаге, в то время, как в методе Адамса необходимо вычислить всего одно значение, гораздо удобнее все-таки использовать метод Рунге-Кутта, в связи с тем, что в методе Адамса расчет необходимо начинать с узла , а не . В практике же обычно известны только начальные условия, а вычисление необходимых для метода Адамса значений с помощью других методов значительно усложняет алгоритм решения. Еще один недостаток многошаговых методов заключается в том, что без изменения формул невозможно изменить изначально взятый шаг h.

2.3 Реализация методов решения систем нелинейных дифференциальных уравнений в пакете MATLAB

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

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

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

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

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

Самые распространенные:

ode23- явные одношаговые методы Рунге-Кутта второго и третьего порядка, достаточно высокая скорость при небольшой точности;

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

ode23t- метод трапеции с интерполяцией, применяется при решении задач, с почти гармоническими входными данными;

ode15s-многошаговый метод, использующий формулы численного дифференцирования, переменного порядка, применяется, если более точные методы не обеспечивают решения;

ode23s-одношаговый метод, использующий модифицированную формулу Розенброка, обеспечивает наиболее высокую скорость, но низкую точность вычислений;

ode133-многошаговый метод Адамса-Башфорта-Мултона переменного порядка, способен обеспечить высокую точность вычислений при достаточных знаниях начальных условий.

Решатели ode23t, ode15s, ode23s используются обычно для решения уравнений неявного вида. Мы выбираем решатель ode45, так как для нашей задачи он вполне подходит и является оптимальным.

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

functionf = F(t,y)

f=[(t,y)];

end

В М-файле f - обозначение функции, F(t,y)-название функции, в квадратных скобках (t,y)- исходная система дифференциальных уравнений, каждое уравнение системы пишется через точку с запятой без правой части.

Решатель может реализовываться в командном окне или может быть записан в скрипт.

[T,Y]=solver('F',[t0 tfinal],[y10 y20 y30 … yn0]);

В выше написанной формуле solver-выбранный решатель, 'F'-ссылка на решаемую систему, [t0 tfinal]-промежуток времени, за которое требуется решение, [y10 y20 y30 … yn0]-начальные условия.

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

Для начала берем самый простой случай:



Найдем его численное решение c помощью MATLAB.

Необходимо переписать систему в виде:

Здесь, , ,

Для начала создадим M-файл «s1»:

function f = s1(t,y)

global m I c1 c2 d1 d2 ;

f=[y(2);-((c1+c2)/m)*y(1);y(4);-((c1*d1^2+c2*d2^2)/I)*y(3)];

end

В другом файле «ss1» запишем программу для исследования системы:

warning off;

clear all

сlc

%Ввод значений переменных

global m I c1 c2 d1 d2 ;

m=input(' Enter m=' );

I=input(' Enter I=' );

d1=input(' Enter d1=' );

d2=input(' Enter d2=' );

c1=input(' Enter c1=' );

c2=input('Enter c2=' );

%Численное решение системы дифференциальных уравнений методом Рунге-Кутты 4 порядка (в данном случае ищем решение на интервале t=[0 10], с начальными условиями =0, =1, =0, =1):

[T,Y]=ode45('s1',[0 10],[0 1 0 1]);

%Построение графиков зависимостей от t(, ,, ):

subplot(2,2,1)

plot(T,Y(:,1))

title('y1(t)')

subplot(2,2,2)

plot(T,Y(:,2))

title('y2(t)')

subplot(2,2,3)

plot(T,Y(:,3))

title('y3(t)')

subplot(2,2,4)

plot(T,Y(:,4))

title('y4(t)')

%Построение фазовых портретов в отдельном окне(,):

figure

subplot(1,2,1)

plot(Y(:,1),Y(:,2))

title('y2(y1)')

subplot(1,2,2)

plot(Y(:,3),Y(:,4))

title('y4(y3)')

После запуска программы «ss1» предлагается ввод переменных. Будем вводить m=1, I=1, d1?d2, c1?c2.

Enter m=1

Enter I=1

Enter d1=1.5

Enter d2=2

Enter c1=1

Enter c2=2

Рис. 3. Слева графики зависимостей(, , , ), справа фазовые портреты (,) для случая гармонических колебаний

Теперь рассмотрим другой случай, введем диссипацию.

Переписываем систему в виде:



СоздаемМ-файл «s2»:

function f = s1(t,y)

global m I c1 c2 d1 d2 p1 p2;

f=[y(2);-((c1+c2)/m)*y(1)-((p1+p2)/m)*y(2);y(4);-((c1*d1^2+c2*d2^2)/I)*y(3)-

-((p1*d1^2+p2*d2^2)/I)*y(4)];

end

В другом файле «ss2» записываем программу для исследования системы:

clc

clear all

warning off;

global m I c1 c2 d1 d2 p1 p2;

m=input(' Enter m=' );

I=input(' Enter I=' );

d1=input(' Enter d1=' );

d2=input(' Enter d2=' );

c1=input(' Enter c1=' );

c2=input('Enter c2=' );

p1=input('Enter p1=' );

p2=input('Enter p2=' );

[T,Y]=ode45('s2',[0 10],[0 1 0 1]);

subplot(2,2,1)

plot(T,Y(:,1))

title('y1(t)')

subplot(2,2,2)

plot(T,Y(:,2))

title('y2(t)')

subplot(2,2,3)

plot(T,Y(:,3))

title('y3(t)')

subplot(2,2,4)

plot(T,Y(:,4))

title('y4(t)')

figure

subplot(1,2,1)

plot(Y(:,1),Y(:,2))

title('y2(y1)')

subplot(1,2,2)

plot(Y(:,3),Y(:,4))

title('y4(y3)')

Запускаем «ss2», вводим значения переменных.

Enter m=1

Enter I=1

Enter d1=1.7

Enter d2=1.2

Enter c1=1

Enter c2=1.4

Enter p1=0.3

Enter p2=0.1



Рис. 4. Слева графики зависимостей(, , , ), справа фазовые портреты (,) для случая затухающих колебаний

Теперь введем все оставшиеся слагаемые:

Переписываем систему:



Создаем соответствующий М-файл:

Function f = s3(t,y)

global m I c1 c2 d1 d2 p1 p2 h1 h2;

f=[y(2);-((c1+c2)/m)*y(1)-(3*(h1*d1^2+h2*d2^2)/m)*y(1)*y(3)^2-((h1+

+ h2)/m)*y(1)^3-(3*(h1*d1-h2*d2)/m)*y(3)*y(1)^2-((c1*d1-c2*d2)/m)*y(3)-

- ((h1*d1^3-h2*d2^3)/m)*y(3)^3-((p1+p2)/m)*y(2)-((p1*d1-p2*d2)/m)*y(4);y(4);-((c1*d1-c2*d2)/I)*y(1)-(3*(h1*d1^3-h2*d2^3)/I)*y(1)*y(3)^2-((h1*d1-

-h2*d2)/I)*y(1)^3-(3*(h1*d1^2+h2*d2^2)/I)*y(3)*y(1)^2-((c1*d1^2+

+ c2*d2^2)/I)*y(3)-((h1*d1^4+h2*d2^4)/I)*y(3)^3-((p1*d1-p2*d2)/I)*y(2) -

- ((p1*d1^2+p2*d2^2)/I)*y(4)];

end

Далее создаем файл для численного решения и построения графиков:

clc

clear all

warning off;

global m I c1 c2 d1 d2 p1 p2 h1 h2;

m=input(' Enter m=' );

I=input(' Enter I=' );

d1=input(' Enter d1=' );

d2=input(' Enter d2=' );

c1=input(' Enter c1=' );

c2=input('Enter c2=' );

p1=input('Enter p1=' );

p2=input('Enter p2=' );

h1=input('Enter h1=' );

h2=input('Enter h2=' );

[T,Y]=ode45('s3',[0 10],[0 1 0 1]);

subplot(2,2,1)

plot(T,Y(:,1))

title('y1(t)')

subplot(2,2,2)

plot(T,Y(:,2))

title('y2(t)')

subplot(2,2,3)

plot(T,Y(:,3))

title('y3(t)')

subplot(2,2,4)

plot(T,Y(:,4))

title('y4(t)')

figure

subplot(1,2,1)

plot(Y(:,1),Y(:,2))

title('y2(y1)')

subplot(1,2,2)

plot(Y(:,3),Y(:,4))

title('y4(y3)')

Предлагается ввести значения переменных:

Enter m=1

Enter I=1

Enter d1=1.4

Enter d2=1.3

Enter c1=1.6

Enter c2 = 1.2

Enter p1 = 0.1

Enter p2 = 0.3

Enter h1 = 3

Enter h2 = 2.7

Рис. 5. Слева графики зависимостей(, , , ), справа фазовые портреты (,) для нелинейного случая затухающих колебаний

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



3. Анализ несимметричной нелинейной системы амортизации

3.1 Реализация решения задачи в пакете MATLAB

Требуется переписать математическую систему, описывающую устройство амортизации, для правильной работы с ней в MATLAB. Все делается аналогичным образом, как уже было представлено в главе 2 п. 2.3 на простейших задачах. Система записывается в М-файл, который является функцией двух переменных: времени t и координаты y. Координата, в свою очередь, имеет 4 составляющих: y(1)=, y(2)=, y(3)=, y(4)=.

Создаем М-файл:

function f = s4(t,y)

global m I c1 c2 d1 d2 p1 p2 h1 h2 W w;

f=[y(2);-((c1+c2)/m)*y(1)-(3*(h1*d1^2+h2*d2^2)/m)*y(1)*y(3)^2-((h1+

+ h2)/m)*y(1)^3-(3*(h1*d1-h2*d2)/m)*y(3)*y(1)^2-((c1*d1-c2*d2)/m)*y(3) -

-((h1*d1^3-h2*d2^3)/m)*y(3)^3-((p1+p2)/m)*y(2)-((p1*d1-p2*d2)/m)*y(4)+ +W*sin(w*t); y(4);-((c1*d1-c2*d2)/I)*y(1)-(3*(h1*d1^3- h2*d2^3)/I)*y(1)*y(3)^2 -

- ((h1*d1-h2*d2)/I)*y(1)^3-(3*(h1*d1^2+h2*d2^2)/I)*y(3)*y(1)^2-((c1*d1^2++ c2*d2^2)/I)*y(3)-((h1*d1^4+h2*d2^4)/I)*y(3)^3-((p1*d1-p2*d2)/I)*y(2) - -((p1*d1^2+p2*d2^2)/I)*y(4)];

end

Далее создаем еще один файл, в котором пишем скрипт для решения системы методом Рунге-Кутта 4 порядка:

clc

clear all

warning off;

%Ввод значений переменных

globalmIc1 c2 d1 d2 p1 p2 h1 h2 Ww;

m = 1;

I = 1;

W = 1;

d1 = 10;

d2 = 1;

c1 = 1;

c2 = 1;

p1 = 0.1;

p2 = 0.1;

h1 = 1;

h2 = 5;

%Цикл по частоте внешнего воздействия для построения АЧХ

i = 1;

forw = 0:0.1:15

[T,Y]=ode45('s4',[0 10],[0 1 0 1]); %решение методом Рунге-Кутта

A(i) = max(Y(:,1)); %максимальное значение x(t)

A1(i) = max(Y(:,3)); %максимальное значение ц(t)

w_ot_i(i) = w;

i = i+1;

end

%Построение АЧХ для x(t) и для ц(t)

subplot(2,1,1)

plot(w_ot_i,A);

xlabel('Частота, w');

ylabel('Амплитуда, А');

title('АЧХ для x(t)');

grid on;

subplot(2,1,2)

plot(w_ot_i,A1);

xlabel('Частота, w');

ylabel('Амплитуда, А1');

title('AЧХ для Ф(t)');

grid on;

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

3.2 Исследование эффектов и поведения решения при различных значениях параметров и анализ полученных результатов

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

1) при одинаковых пружинах при большой диссипации возникают только вертикальные колебания, а крутильные не возникают, то есть угол поворота груза всегда равен нулю;

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

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

Ограничения на вводимые значения в связи с антисимметричной постановкой задачи: d1?d2 и c1?c2. Для начала возьмем жесткость пружин одинаковой, чтобы исследовать зависимость поведения АЧХ только от расстояний от центра масс до точек прикрепления пружин.

3.2.1 Исследование АЧХ в зависимости от значений расстояний от центра масс до точек прикрепления пружин

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



Рис. 6. Амплитудно-частотная характеристика при m = 1, I = 1, W = 1, d1 = 10, d2 = 1, c1 = 1, c2 = 1, p1 = 0.1, p2 = 0.1, h1 = 1, h2 = 5

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

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

Рис. 7. АЧХ при d1 = 9, d2 = 2

Рис. 9. АЧХ при d1 = 8, d2 = 3

Рис. 11. АЧХ при d1 = 7, d2 = 4



Рис. 8. АЧХ при d1 = 6, d2 = 5

Рис. 10. АЧХ при d1 = 5, d2 = 1



Рис. 12. АЧХ при d1 = 2, d2 = 1

Рис. 13. АЧХ при d1 = 1, d2 = 3



Рис. 15. АЧХ при d1 = 2, d2 = 5

Рис. 17. АЧХ при d1 = 3, d2 = 7



Рис. 14. АЧХ при d1 = 4, d2 = 9

Рис.16. АЧХ при d1 = 4, d2 = 6

Рис. 18. АЧХ при d1 = 1, d2 = 9

Для наглядности создадим таблицу.

Таблица 1 - Зависимость амплитуды найденного решения от значений расстояний от пружин до центра масс.

d1	d2	?d	Max(A(w))	w	Max(A1(w))	w
10	1	9	1.4	2.6	0.075	0
9	2	7	1.25	2.6	0.0539	0
8	3	5	1.25	2.5	0.0611	0
7	4	3	1.2	2.5	0.0679	0
6	5	1	1.2	2.6	0.113	0
5	1	4	1.25	2.4	0.0729	0
2	1	1	1.2	2.4	0.0893	0
1	3	2	1.7	2.1	0.349	2.1
2	5	3	1.6	2.1	0.165	2.1
3	7	4	1.6	2.1	0.15	2.1
4	9	5	1.6	2.1	0.081	2.1
4	6	2	1.48	2.1	0.087	2.1
1	8	7	1.9	2	0.17	2

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

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

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

При состоянии системы, близкому к симметричному, можно наблюдать, что в АЧХ крутильных колебаний в исследуемой области частот возникает только один резонанс. По мере приближения системы к симметричному состоянию область резонансных частот крутильных колебаний становится более узкой, значение резонансной амплитуды меньше, а при частотах, далеких от частоты резонанса, амплитуда остается постоянной. Если d2 брать больше d1 , то постоянное значение амплитуды крутильных колебаний приближается к нулевому значению.

Следует заметить, что при уменьшении разницы между расстояниями ?d, максимальные значения амплитуд колебаний функций тоже уменьшаются, то есть чем больше асимметрия системы, тем больше максимальные значения амплитуд.

Наибольшая амплитуда колебаний функций получилась в случае, когда d1=1 и d2=9, поэтому есть смысл проводить дальнейшее исследование с данными значениями расстояний от центра масс до пружин.

3.2.2 Исследование АЧХ в зависимости от значений жесткостей пружин

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



Рис. 19. АЧХ при c1 = 1, c2 = 5

Рис. 20. АЧХ при c1 = 1.5, c2 =4.5



Рис. 21. АЧХ при c1 = 2, c2 = 4

Рис.23. АЧХ при c1 = 2.5, c2 =3.5



Рис. 25. АЧХ при c1 = 4.5, c2 = 1

Рис. 22. АЧХ при c1 = 4, c2 = 1.5



Рис. 24. АЧХ при c1 = 3.5, c2 = 2

Рис. 26. АЧХ при c1 = 3, c2 = 2.5



Рис. 27. АЧХ при c1 = 1, c2 = 1

Рис. 29. АЧХ при c1 = 2, c2 = 2



Рис. 28. АЧХ при c1 = 3, c2 = 3

Рис. 30. АЧХ при c1 = 5, c2 = 5

Таблица 2 - Зависимость амплитуды найденного решения от значений жесткостей пружин

с1	с2	?с	Max(A(w))	w	Max(A1(w))	w
1	5	4	1.9	2	0.18	2
1.5	4.5	3	1.85	2.05	0.165	2.05
2	4	2	1.8	2.1	0.16	2.1
2.5	3.5	1	1.7	2.2	0.16	2.2
4.5	1	3.5	1.4	2.6	0.105	2.6
4	1.5	2.5	1.5	2.5	0.12	2.5
3.5	2	1.5	1.55	2.3	0.13	2.3
3	2.5	0.5	1.65	2.2	0.145	2.2
1	1	0	1.9	2	0.16	2
2	2	0	1.8	2.1	0.155	2.1
3	3	0	1.65	2.15	0.15	2.15
5	5	0	1.4	2.8	0.125	2.8

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

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

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

3.2.3 Исследование АЧХ в зависимости от значений диссипации

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



Рис. 7. АЧХ при p1 = 0.1, p2 = 0.1

Рис. 7. АЧХ при p1 = 0.2, p2 = 0.2



Рис. 7. АЧХ при p1 = 0.3, p2 = 0.3

Рис. 7. АЧХ при p1 = 0.4, p2 = 0.4



Рис. 7. АЧХ при p1 = 0.5, p2 = 0.5

Рис. 7. АЧХ при p1 = 0.6, p2 = 0.6



Рис. 7. АЧХ при p1 = 0.7, p2 = 0.7

Рис. 7. АЧХ при p1 = 0.8, p2 = 0.8



Рис. 7. АЧХ при p1 = 0.9, p2 = 0.9

Рис. 7. АЧХ при p1 = 1, p2 = 1



Рис. 7. АЧХ при p1 = 0.5, p2 = 0.6

Рис. 7. АЧХ при p1 = 0.7, p2 = 0.5



Рис. 7. АЧХ при p1 = 0.4, p2 = 0.7

Рис. 7. АЧХ при p1 = 0.6, p2 = 0.2



Рис. 7. АЧХ при p1 = 0.3, p2 = 0.8

Рис. 7. АЧХ при p1 = 1, p2 = 0.4



Рис. 7. АЧХ при p1 = 0.2, p2 = 0.9

Рис. 7. АЧХ при p1 = 0.9, p2 = 0.1

Таблица 3 - Зависимость амплитуды найденного решения от значений диссипаций

p1	p2	?p	Max(A(w))	w	Max(A1(w))	w
0.1	0.1	0	1.8	2	0.18	2
0.2	0.2	0	1.6	1.9	0.16	1.9
0.3	0.3	0	1.4	1.8	0.14	1.8
0.4	0.4	0	1.2	1.7	0.125	1.7
0.5	0.5	0	1	1.6	0.113	1.6
0.6	0.6	0	0.92	1.5	0.098	1.5
0.7	0.7	0	0.84	1.4	0.09	1.4
0.8	0.8	0	0.8	2	0.086	2
0.9	0.9	0	0.77	2	0.081	2
1	1	0	0.74	2	0.078	2
0.5	0.6	0.1	1	1.8	0.113	1.8
0.7	0.5	0.2	0.84	1.4	0.09	1.4
0.4	0.7	0.3	1.2	1.8	0.125	1.8
0.6	0.2	0.4	0.91	1.8	0.098	1.8
0.3	0.8	0.5	1.4	2	0.14	2
1	0.4	0.6	0.73	2.3	0.079	2.3
0.2	0.9	0.7	1.6	2	0.16	2
0.9	0.1	0.8	0.77	2	0.081	2
0.1	1	0.9	1.8	2	0.19	2

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

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

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

Заключение

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

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

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

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

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

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

Резонансы меньших порядков, по сравнению с главным, можно наблюдать практически на каждом графике АЧХ. При увеличении диссипации можно увидеть, как резонанс меньшего порядка «вытесняет» главный резонанс.

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

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

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

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



Список используемой литературы

1. Рабинович М.И., Трубецков Д.И. Введение в теорию колебаний и волн. -- Ижевск : НИЦ «Регулярная и хаотическая динамика», 2000. -- 560 с.

2. Бутенин Н.В., Неймарк Ю.И. Фуфаев Н.Л. Введение в теорию нелинейных колебаний. -- М.: Наука, 1987. -- 385 с.

3. Наместников С.М. Основы программирования в MatLab / Сборник лекций: УлГТУ, Ульяновск. - 2011. - 55 с.

4. Иглин С.П. Математические расчеты на базе MATLAB. СПб.: БХВ-Петербург, 2005. -640c.

5. Петров Л.Ф. Методы динамического анализа в экономике. М.: Инфа-М, 2010.-240с.

6. Водолазская И.В. Введение в систему MATLAB. Астрахань.:АГТУ, 2004.- 50с.

7. Ануфриев И.П. MATLAB 7.0. Наиболее полное руководство. СПб.: БХВ-Петербург, 2005-1104с.

8. Джон Г. Мэтьюз, Куртис Д. Финк. Численные методы. Использование MATLAB = Numerical Methods : Using MATLAB. - 3-е изд.-М.: Вильямс, 2001.-716с.Корненев Б.Г., Резников Л.М. Динамические гасители колебаний: теория и технические приложения. М.: Наука, 1988.-304с.

9. Бахвалов Н.С. Численные методы. М.: Наука, 2006.- 631с.

10. Вержбицкий В.М. Численные методы (математический анализ и обыкновенные дифференциальные уравнения): учебное пособие для вузов. М.: Высшая школа, 2001.-381с.

11. Корненев Б.Г., Резников Л.М. Динамические гасители колебаний: теория и технические приложения. М.: Наука, 1988.-304с.

Размещено на Allbest.ru

...