«Mоделирование процессов и технических систем в MATLAB Учебный курс Киев – 2004 2 УДК 681.3.06(075.8) ББК 32.973.26-018.2 Я73 Л17 Лазарев Юрий Федорович Л17 MatLAB 6.5. Mатематическое моделирование физических процессов и ...»

Юрий ЛАЗАРЕВ

_

Mоделирование процессов

и технических систем

в MATLAB

Учебный курс

Киев – 2004

2

УДК 681.3.06(075.8)

ББК 32.973.26-018.2 Я73

Л17

Лазарев Юрий Федорович

Л17 MatLAB 6.5. Mатематическое моделирование физических

процессов и технических систем: Учебный курс. - К.: 2004. - 474 с.

Изложены основные особенности проведения вычислений в среде MatLAB как в режиме калькулятора, так и в программном режиме. Ознакомление с системой рассчитано на начинающего. Приведены сведения об основных командах, операторах, функциях и процедурах MatLAB. Изложение ведется таким образом, чтобы пользователь мог сразу применить полученные знания для проведения вычислений. Пособие содержит много примеров, которые поясняют и иллюстрируют работу по использованию процедур. Рассмотрена работа с некоторыми наиболее важными для инженеров пакетами прикладных программ MatLAB (Signal Toolbox, Control и SimuLink).

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

Илл. 487. Библиогр. 21 назв.

Содержание Предисловие Введение Урок 1. MatLAB как научный калькулятор 1.1. Командное окно 1.2. Операции с числами 1.3. Простейшие операции с векторами и матрицами 1.4. Функции прикладной численной математики 1.5. Построение простейших графиков 1.6. Операторы управления вычислительным процессом 1.7. Вопросы для самопроверки Урок 2. Программирование в среде MatLAB 2.1. Функции функций 2.2. Создание М-файлов 2.3. Создание простейших файлов-функций (процедур) 2.4. Создание Script-файлов 2.5. Графическое оформление результатов 2.6. Создание функций от функций 2.7. Программа моделирования движения маятника 2.8. Вопросы для самопроверки Урок 3. Интерфейс системы MatLAB 3.1. Команды связи с операционной средой 3.2. Использование MatLAB при оформлении текстового документа 3.3. Использование в MatLAB файлов данных 3.4. Вопросы для самопроверки Урок 4. Классы вычислительных объектов 4.1. Основные классы объектов 4.2. Производные классы MatLAB 4.3. Пример создания нового класса polynom 4.4. Создание методов нового класса 4.5. Вопросы для самопроверки Урок 5. Цифровая обработка сигналов (пакет Signal Processing Toolbox) 5.1. Формирование типовых процессов 5.2. Общие средства фильтрации. Формирование случайных процессов 5.3. Процедуры спектрального (частотного) и статистического анализа процессов 5.4. Проектирование фильтров 5.5. Графические и интерактивные средства 5.6. Вопросы для самопроверки Урок 6. Исследование линейных стационарных систем(пакет Control) 6.1. Общая характеристика процедур пакета Control 6.2. Ввод и преобразование моделей 6.3. Получение информации о модели 6.4. Анализ системы 6.5. Интерактивный "обозреватель" ltiview 6.6. Синтез системы Урок 7. Основы визуального моделирование динамических систем Урок 9. Моделирование динамики аэрокосмических объектов 9.2. Моделирование свободного углового движения космического аппарата 9.3. Моделирование управляемого углового движения космического аппарата Урок 10. Моделирование електродинамических систем Урок 11. Моделирование машин и механизмов Предисловие В последние годы в университетских и инженерно-технических кругах мира получила широкое распространение новая компьютерная система осуществления математических расчетов - система MatLAB. Более того, в настоящее время система MatLAB принята в качестве официального вычислительного средства при подготовке и оформлении инженерной документации и научных публикаций. В чем причина такой популярности этой системы?

Главные преимущества "языка технических вычислений" MatLAB, выгодно выделяющие его среди других существующих ныне математических систем и пакетов, состоят в следующем:

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

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

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

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

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

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

Книга состоит из одиннадцати глав (уроков). Она задумана как учебное пособие для студентов высших технических учебных заведений и естественнонаучных специальностей университетов. Первые четыре урока могут быть рекомендованы в качестве пособия по учебным дисциплинам "Моделирование процессов и систем" и "Математическое моделирование на ЭВМ". Остальная часть может быть использована как пособие по курсовому и дипломному проектированию. В целом книгу можно рассматривать и как введение в MatLAB для инженеров и научных работников в области проектирования технических систем.

В Уроке 1 читатель знакомится с возможностями системы в режиме научного калькулятора. Здесь помещены сведения об основных операторах, командах, функциях и процедурах системы. В Уроке2 описаны правила и примеры составления программ на языке MatLAB. Кроме того, в нем представлены некоторые дополнительные процедуры, которые помогают рационально организовать вычислительный процесс. Урок 3 содержит перечень некоторых процедур и команд общего назначения, которые связывают систему MatLAB с операционной системой компьютера, а также со средствами, позволяющими использовать возможности MatLAB при оформлении документов в текстовом редакторе WORD. Здесь же описан механизм формирования и чтения файлов данных в среде MatLAB.

Важной частью MatLAB, которая позволяет приспосабливать систему к задачам пользователя, является возможность образования новых классов вычислительных объектов. С понятием классов вычислительных объектов в MatLAB и правилами создания новых классов пользователь ознакомится в Уроке 4.

В Уроке 5 сосредоточены сведения об особенностях использования процедур цифровой обработки сигналов пакета SIGNAL. Содержание Урока 6 - начальное ознакомление с особенностями работы с процедурами анализа и синтеза линейных стационарных систем автоматического управления пакета CONTROL.

Урок 7 знакомит с ядром пакета SimuLink интерактивного (визуального) моделирования динамических систем во временной области. Содержанием Урока 8 является более глубокое ознакомление с важнейшими средствами SimuLink, позволяющими обеспечить эффективное взаимодействие SimuLink со средой MatLAB.

Наконец Уроки 9, 10 и 11 посвящены ознакомлению читателя с основами использования трех дополнительных библиотек пакета SimuLink: моделирования динамики аэрокосмических объектов (Aerospace Blockset), моделирования електротехнических систем (SimPowerSystems) и моделирования динамики механизмов (SimMechanics).

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

В пособии использованы материалы из изданий, указанных в списке литературы.

Учебный курс, в основном, ориентирован на русифицированную версию MatLAB 6.5.0.180913a (R13) Система MatLAB создана фирмой MathWork Inc. (США, г. Нейтик, штат Массачузетс). Хотя впервые эта система начала использоваться в конце 70-х годов, расцвет ее применения начался в конце 80-х, в особенности после появления на рынке версии 4.0. Последние версии MatLAB, - это чрезвычайно развитые системы, которые содержат огромную совокупность процедур и функций, необходимых инженеру и научному работнику для осуществления сложных численных расчетов, моделирования поведения технических и физических систем, оформления результатов этих расчетов в наглядном виде.

MatLAB (сокращение от MATrix LABoratory - матричная лаборатория) представляет собой интерактивную компьютерную систему для выполнения инженерных и научных расчетов, ориентированную на работу с массивами данных. Система предполагает возможность обращения к программам, которые написаны на языках FORTRAN, C и C++.

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

Основной объект системы MatLAB - прямоугольный числовой массив (матрица), который допускает комплексные элементы. Использование матриц не требует явного указания их размеров. MatLAB позволяет решать многие вычислительные задачи за значительно меньшее время, чем то, которое необходимо для написания соответствующих программ на языках FORTRAN, BASIC и C.

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

В базовый набор слов системы входят: спецзнаки; знаки арифметических и логических операций;

арифметические, тригонометрические и некоторые специальные математические функции; функции быстрого преобразования Фурье и фильтрации; векторные и матричные функции; средства для работы с комплексными числами; операторы построения графиков в декартовой и полярной системах координат, трехмерных поверхностей и т.п. То есть MatLAB предоставляет пользователю большой набор готовых средств (более половины из них - внешние расширения в виде m-файлов).

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

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

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

Работа в среде MatLab может осуществляться в двух режимах:

- в режиме калькулятора, когда вычисления осуществляются сразу после набора очередного оператора или команды MatLab; при этом значение результатов вычисления могут присваиваться некоторым переменным, или результаты получаются непосредственно, без присваивания (как в обычных калькуляторах);

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

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

MatLAB имеет черты разных известных языков программирования высокого уровня.

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

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

Принципы сохранения значений переменных в MatLAB более всего приближаются к тем, которые присущи языку FORTRAN, а именно: все переменные являются локальными, то есть действуют лишь в границах той программной единицы (процедуры, функции или главной, управляющей программы), где им присвоены некоторые конкретные значения. При переходе к выполнению другой программной единицы, значения переменных предыдущей программной единицы или совсем теряются (в случае, если выполненная программная единица является процедурой или функцией), или становится недосягаемыми (если выполненная программа - управляющая). В отличие от языков BASIC и Pascal, в языке MatLAB нет глобальных переменных, действие которых распространялась бы на все программные единицы. Но при этом язык MatLAB обладает особенностью, которая отличает его от других языков. В отличие от интерпретатора BASIC, интерпретатор MatLAB позволяет в одном и том же сеансе работы выполнять несколько самостоятельных программ, причем все переменные, используемые в этих программах, являются общими для этих программ и образуют общее рабочее пространство (Work Space). Это позволяет более рационально организовывать сложные (громоздкие) вычисления по принципу, который напоминает оверлейные структуры.

Указанные особенности MatLAB делают ее весьма гибкой и удобной в пользовании вычислительной системой.

Урок 1. MatLab как научный калькулятор Простейшие операции с векторами и матрицами Функции прикладной численной математики Построение простейших графиков Операторы управления вычислительным процессом Работа в среде MatLAB может осуществляться в двух режимах:

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

путем запуска из командного окна специально написанной на языке MatLAB (M-языке) программы;

такой режим можно назвать программным.

В командном режиме пользователю доступны практически все возможности и функции MatLAB, за небольшим исключением. Систему MatLAB обоснованно относят к одному из наиболее мощных научных калькуляторов, которому доступны практически все численные средства решения научных и инженерных задач, разработанные в научных организациях мира на настоящий момент. Использование этих средств в командном режиме является, в большинстве случаев, очень простым. Результат получают сразу в самом командном окне в наглядной простой форме или в графическом виде в дополнительном графическом окне. Поэтому знакомство с MatLAB и освоение приемов работы в ее среде целесообразно начать с изучения возможностей системы именно в командном режиме 1.1. Командное окно После вызова MatLAB 6.5 из среды Windows на экране появляется окно MATЛAБ, представленное на рис 1.1. В нем могут отображаться несколько окон. Главным из них является Окно команд, или так называемое командное окно среды MatLAB.

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

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

В верхней части окна (под заголовком) расположена строка меню,. Чтобы открыть какое-либо меню, следует установить на нем курсор мыши и нажать ее левую кнопку. Здесь отметим лишь, что для выхода из среды MatLAB достаточно открыть меню Файл и выбрать в нем команду Выход из MATLAB, или просто закрыть командное окно, щелкнув мышью на кнопке закрытия окна с изображением крестика.

1.2. Операции с числами Основные объекты системы MatLAB – числа. Операции с ними лежат в основе всей работы с этой системой.

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

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

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

Если, например, ввести в командном окне MatLAB строку 1. 20357651е -17, то после нажатия клавиши Enter в этом окне появится запись, показанная на рис. 1.3.

Видно, что выведенное на экран число не совпадает с введенным. Это обусловлено тем, что результат вычислений в MatLAB выводится в виде (формате), который определяется предварительно установленным форматом представления чисел. Этот формат может быть установлен с помощью команды Файл Предпочтения. После ее вызова на экране появится одноименное окно (рис. 1.4).

Один из участков этого окна имеет название Числовой формат. Он предназначен для установки и изменения формата представления чисел, которые выводятся в командное окно в процессе расчетов. Предусмотрены такие форматы (рис 1.5):

Short краткая запись (применяется по умолчанию) с фиксированной (default) запятой;

Long длинная запись с фиксированной запятой;

Short E краткая запись в формате с плавающей запятой;

Long E длинная запись в формате с плавающей запятой;

Short G вторая форма краткой записи в формате с плавающей запятой;

Long G вторая форма длинной записи в формате с плавающей запятой;

Hex запись в виде шестнадцатеричного числа;

Rational запись в виде рациональной дроби.

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

Как видно из рис. 1.3, число, которое выведено на экран, не совпадает с введенным числом. Это обусловлено тем, что установленный по умолчанию формат представления чисел (Short Е) не позволяет вывести больше значащих цифр. На самом деле введенное число сохраняется внутри MatLAB со всеми введенными его цифрами. Например, если выбрать формат Long E, то, повторяя те же действия, получим представление числа, где все цифры отображены верно (рис. 1.6).

Следует помнить:

- введенное число и результаты всех вычислений в системе MatLAB сохраняются в памяти ПК с относительной погрешностью около 2.10 -16 (т. е. с точными значениями в 15 десятичных разрядах);

- диапазон представления модуля действительных чисел лежит в диапазоне между 10308 и 10 +308.

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

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

(4.5)2*7.23 - 3.14*10. Если после ввода с клавиатуры этой последовательности нажать клавишу Enter, в командном окне возникнет результат выполнения в виде, представленном на рис. 1.7, т. е. на экран под именем системной переменной ans выводится результат действия последнего выполненного оператора.

Вообще вывод промежуточной информации в командное окно подчиняется таким правилам:

если запись оператора не заканчивается символом «;», результат действия этого оператора сразу же выводится в командное окно;

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

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

в общем случае форма представления результата в командном окне имеет вид:

Пример. Пусть нужно вычислить выражение (25+17)*7. Это можно сделать таким образом. Сначала набираем последовательность 25+17 и нажимаем Enter. Получаем на экране результат в виде ans = 42. Теперь записываем последовательность ans*7 и нажимаем Enter. Получаем ans = 294 (рис. 1.8).

Чтобы предотвратить выведение промежуточного результата действия 25+17, достаточно после записи этой последовательности добавить символ « ; ». Тогда получим результаты в виде, представленном на рис. 1.9.

Используя MatLAB как калькулятор, можно использовать имена переменных для записи промежуточных результатов в память ПК. Для этого служит операция присваивания, которая вводится знаком равенства « = » в соответствия со схемой:

Имя переменной может содержать до 30 символов и не должно совпадать с именами функций, процедур системы и системных переменных. При этом система различает большие и малые буквы в переменных. Так, имена amenu, Аmenu, aMenu в MatLAB обозначают разные переменные.

Выражение справа от знака присваивания может быть просто числом, арифметическим выражением, строкой символов (тогда эти символы нужно заключить в апострофы) или символьным выражением. Если выражение не заканчивается символом « ; », после нажатия клавиши Enter в командном окне возникнет результат выполнения в виде:

Например, если ввести в командное окно строку x = 25 + 17, на экране появится запись (рис. 1.10):

Система MatLAB имеет несколько имен переменных, которые используются самой системой и входят в состав зарезервированных (эти переменные можно использовать в математических выражениях):

i, j мнимая единица (корень квадратный из -1);

inf обозначение машинной бесконечности;

NaN обозначение неопределенного результата (например, типа 0/0 или inf/inf);

eps погрешность операций над числами с плавающей запятой;

ans результат последней операции без знака присваивания;

realmax и максимально и минимально возможные величины числа, которые могут быть исrealmin пользованы.

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

Для обозначения мнимой единицы в языке MatLAB зарезервированы два имени i и j. Ввод с клавиатуры значения комплексного числа осуществляется путем записи в командное окно строки вида:

= + i[j]*, где ДЧ — действительная часть комплексного числа, МЧ — мнимая часть. Из примера, приведенного на рис.

1.11, видно, в каком виде система выводит комплексные числа на экран (и на печать).

1.2.4. Элементарные математические функции Общая форма использования функции в MatLAB такова:

В языке MatLAB предусмотрены следующие элементарные математические функции.

sinh(Z) гиперболический синус;

asin(Z) asinh(Z) обратный гиперболический синус;

cosh(Z) гиперболический косинус;

acos(Z) acosh(Z) обратный гиперболический косинус;

tanh(Z) гиперболический тангенс;

atan(Z) четырехквадрантный арктангенс (угол в диапазоне (-, + ] между горизонтальным atan2(X,Y) правым лучом и лучом, который проходит через точку с координатами X и Y);

atanh(Z) обратный гиперболический тангенс;

sech(Z) гиперболический секанс;

asech(Z) обратный гиперболический секанс;

csch(Z) гиперболический косеканс;

acsch(Z) обратный гиперболический косеканс;

coth(Z) гиперболический котангенс;

acoth(Z) обратный гиперболический котангенс.

log(Z) натуральный логарифм;

log10(Z) десятичный логарифм;

sqrt(Z) квадратный корень из числа Z;

fix(Z) округление к ближайшему целому в сторону нуля;

floor(Z) округление к ближайшему целому в сторону отрицательной бесконечности;

ceil(Z) округление к ближайшему целому в сторону положительной бесконечности;

round(Z) обычное округление числа Z к ближайшему целому;

mod(X,Y) целочисленное деление X на Y;

rem(X,Y) вычисление остатка от деления X на Y;

sign(Z) вычисление сигнум-функции числа Z (0 при Z=0, -1 при Z0) 1.2.5. Специальные математические функции Кроме элементарных, в языке MatLAB предусмотрен целый ряд специальных математических функций. Ниже приведен перечень и краткое содержание этих функций. Правила обращения к ним и использования пользователь может отыскать в описаниях этих функций, которые выводятся на экран, если набрать команду help и указать в той же строке имя функции.

cart2sph преобразование декартовых координат в сферические;

cart2pol преобразование декартовых координат в полярные;

pol2cart преобразование полярных координат в декартовые;

sph2cart преобразование сферических координат в декартовые.

besselj функция Бесселя первого рода;

bessely функция Бесселя второго рода;

besseli модифицированная функция Бесселя первого рода;

besselk модифицированная функция Бесселя второго рода.

betainc усеченная Бета- функция;

betaln логарифм Бета- функции gammainc усеченная Гамма-функция;

gammaln логарифм Гамма- функции ellipj эллиптические функции Якоби;

ellipke полный эллиптический интеграл;

expint функция экспоненциального интеграла.

erfc дополнительная функция ошибок;

erfcx масштабированная дополнительная функция ошибок;

erfinv обратная функция ошибок.

legendre обобщенная функция Лежандра;

rat представление числа в виде рациональной дроби;

rats представление чисел в виде рациональной дроби.

1.2.6. Элементарные действия с комплексными числами Простейшие действия с комплексными числами - сложение, вычитание, умножение, деление и возведение в степень - осуществляются при помощи обычных арифметических знаков +, -, *, /, \ и ^ соответственно. Примеры использования приведены на рис. 1.12.

ПРИМЕЧАНИЕ В приведенном фрагменте использована функция disp (от слова 'дисплей'), которая тоже выводит в командное окно результаты вычислений или некоторый текст. При этом численный результат, как видно, выводится уже без указания имени переменной или ans.

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

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

В MatLAB есть несколько дополнительных функций, рассчитанных только на комплексный аргумент:

real(Z) выделяет действительную часть комплексного аргумента Z;

imag(Z) выделяет мнимую часть комплексного аргумента;

вычисляет значение аргумента комплексного числа Z (в радианах в диапазоне от - до angle(Z) conj(Z) выдает число, комплексно сопряженное относительно Z.

Примеры приведены на рис. 1.13.

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

v = [ -1, -1+2i,-5,4,5i,-1-2i,-5i] disp(cplxpair(v)) Приспособленность большинства функций MatLAB к оперированию с комплексными числами позволяет значительно проще строить вычисления с действительными числами, результат которых является комплексным, например, находить комплексные корни квадратных уравнений.

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

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

Например, запись строки V = [ 1.2 -0.3 1.2e-5] задает вектор V, который содержит три элемента со значениями 1.2, -0.3 и 1.2е-5 (рис. 1.14):

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

Длинный вектор можно вводить частями, которые потом объединять с помощью операции объединения векторов в строку : v = [ v1 v2 ]. Например (см. рис.1.15):

Язык MatLAB дает пользователю возможность сокращенного введения вектора, значения элементов которого составляют арифметическую прогрессию. Если обозначить nz - начальное значение этой прогрессии (значение первого элемента вектора), kz - конечное значение прогрессии (значение последнего элемента вектора), а h разность прогрессии (шаг), то вектор можно ввести с помощью короткой записи V = nz : h : kz. Например, введение строки V = - 0.1 : 0.3 : 1.4 приведет к такому результату (рис. 1.16):

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

Например, команда приводит к формированию такого вектора Так вводятся векторы-строки. Вектор-столбец вводится аналогично, но значения элементов отделяются знаком « ; ».

Ввод значений элементов матрицы осуществляется в MatLAB в квадратных скобках, по строкам. При этом элементы строки матрицы один от другого отделяются пробелом или запятой, а строки одна от другой отделяются знаком « ; » (рис. 1.17).

1.3.2. Формирование векторов и матриц MatLAB имеет несколько функций, которые позволяют формировать векторы и матрицы некоторого определенного вида. К таким функциям относятся:

zeros(М,N) — создает матрицу размером (МN) с нулевыми элементами, например:

» zeros(3,5) ones(М,N) — создает матрицу размером (МN) с единичными элементами, например:

» ones(3,5) eye(М,N) — создает единичную матрицу размером (МN), т. е. с единицами по главной диагонали и остальными нулевыми элементами, например:

» eye(3,5) rand(М,N) - создает матрицу размером (МN) из случайных чисел, равномерно распределенных в диапазоне от 0 до 1, например:

» rand(3,5) - создает матрицу размером (М N) из случайных чисел, распределенных по нормальному (гаrandn(М,N) уссовому) закону с нулевым математическим ожиданием и стандартным (среднеквадратичным) отклонением, равным единице, например:

» randn(3,5) - создает матрицу Адамара размером (N N), например:

hadamard(N) » hadamard(4) - создает матрицу Гильберта размером (N N), например:

hilb(N) » hilb(4) invhilb(N) - создает обратную матрицу Гильберта размером (N N), например:

» invhilb(4) - создает матрицу Паскаля размером (N N), например:

pascal(N) » pascal(5) В языке MatLAB предусмотрено несколько функций, которые позволяют формировать матрицу на основе другой (заданной) или используя некоторый заданный вектор. К таким функциям принадлежат:

- формирует матрицу, переставляя столбцы известной матрицы А относительно вертикальной fliplr(A) оси, т. е. меняя местами левую и правую стороны матрицы, например:

» fliplr(A) flipud(A) - переставляет строки заданной матрицы А относительно горизонтальной оси, т. е. меняя местами верхнюю и нижнюю стороны матрицы, например:

» flipud(A) - формирует матрицу путем "поворота" заданной матрицы А на 90 градусов против часовой rot90(A) стрелки:

» rot90(A) - образует матрицу размером (m n) путем выборки элементов заданной матрицы А по reshape(A,m,n) столбцам и последующего распределения этих элементов по n столбцам, каждый из которых содержит m элементов; при этом число элементов матрицы А должно равняться m*n, например:

» reshape(A,2,9) - образует нижнюю треугольную матрицу на основе матрицы А путем обнуления ее элементов tril(A) выше главной диагонали:

» tril(A) triu(A) - образует верхнюю треугольную матрицу на основе матрицы А путем обнуления ее элементов ниже главной диагонали:

» triu(A) hankel(V) - образует квадратную матрицу Ганкеля, первый столбец которой совпадает с заданным вектором V, например:

>> V = [-5 6 7 4] » hankel(V) Процедура diag(х) - формирует или извлекает диагональ матрицы.

Если х - вектор, то функция diag(х) создает квадратную матрицу с вектором х на главной диагонали:

» diag(V) Чтобы установить заданный вектор на другую диагональ, при обращении к функции необходимо указать еще один параметр (целое число) - номер диагонали (при этом диагонали отсчитываются от главной вверх), например:

Если х - матрица, то функция diag(х) создает вектор-столбец, который состоит из элементов главной диагонали заданной матрицы х, например, для матрицы А, указанной перед примером применения процедуры fliplr:

» diag(A) Если при этом указать дополнительно номер диагонали, то можно получить вектор-столбец из элементов любой диагонали матрицы х, например:

» diag(A,3) Функция zeros(1,N) формирует (создает) вектор-строку из N нулевых элементов. Аналогично zeros(N,1) создает вектор-столбец из N нулей.

Векторы, значения элементов которых являются случайными равномерно распределенными, формируются таким образом: rand(1,n) - для вектора-строки и rand(m,1) - для вектора-столбца.

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

Пусть имеем некоторую матрицу А:

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

>> A(2,3) Если нужно, наоборот, установить на это место некоторое число, например, Иногда нужно создать меньшую матрицу из большей, формируя ее путем извлечения из последней матрицы элементов ее нескольких строк и столбцов. Или, наоборот, вставить меньшую матрицу таким образом, чтобы она стала определенной частью матрицы большего размера. Это в MatLAB делается с помощью знака двоеточия (« : »).

Рассмотрим эти операции на примерах.

Пусть нужно создать вектор V1, состоящий из элементов третьего столбца последней матрицы А. Для этого произведем такие действия:

>> V1 = A(:, 3) Чтобы создать вектор V2, состоящий из элементов второй строки матрицы А, поступают так:

>> V2 = A(2, : ) Допустим, что необходимо из матрицы А образовать матрицу В размером (2 2), которая состоит из элементов левого нижнего угла матрицы А. Тогда делают так:

>> B = A(2:3, 1:2) Аналогично можно вставить матрицу В в верхнюю середину матрицы А:

>> A(1:2,2:3)=B Как видно, для этого вместо указания номеров элементов матрицы можно указывать диапазон изменения этих номеров путем указания нижней и верхней границ, разделяя их двоеточием.

Примечание. Если верхней границей изменения номеров элементов матрицы является ее размер в Например:

>> A(2:end,2:end) Эти операции очень удобны для формирования матриц, большинство элементов которых одинаковы, в частности, так называемых разреженных матриц, которые состоят, в основном, из нулей, за исключением отдельных элементов. Для примера рассмотрим формирование разреженной матрицы размером (5 7) с единичными элементами в ее центре:

>> A = zeros(5,7);

>> B = ones(3,3);

>> A(2:4,3:5)=B "Растянуть" матрицу (А) в единый вектор (V) можно с помощью обычной записи "V = A(:)". При этом создается вектор-столбец с количеством элементов (m n), в котором столбцы заданной матрицы размещены сверху вниз в порядке самих столбцов:

» v = A(:) Наконец, "расширить" матрицу, составляя ее из отдельных заданных матриц ("блоков") можно тоже довольно просто. Если заданы несколько матриц-блоков А1, А2,... АN с одинаковым количеством строк, то из них можно "слепить" единую матрицу А, объединяя блоки в одну "строку" таким образом:

A = [A1, A2,..., AN] Эту операцию называют горизонтальной конкатенацией (сцеплением) матриц. Вертикальная конкатенация матриц реализуется (при условии, что все составные блоки-матрицы имеют одинаковое количество столбцов) аналогично, путем применения для отделения блоков вместо запятой точки с запятой:

A = [A1; A2;... ; AN].

Приведем примеры. Пример горизонтальной конкатенации :

>> A2 = [10;11;12];

>> A3 = [14 15; 16 17; 18 19];

>> A = [A1, A2, A3] Пример вертикальной конкатенации:

>> B3 = [17 18 19 20 21];

>> B = [ B1; B2; B3] 1.3.4. Действия над векторами Будем различать две группы действий над векторами:

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

б) действия по преобразованию элементов - это действия, которые преобразуют элементы вектора, но не являются операциями, разрешенными математикой.

Сложение векторов. Как известно, суммироваться могут только векторы одинакового типа (т. е. такие, которые оба являются или векторами-строками, или векторами-столбцами), имеющие одинаковую длину (т. е. одинаковое количество элементов). Если X и Y - именно такие векторы, то их сумму Z можно получить, введя команду Z = X + Y, например:

»v=x+y Аналогично с помощью арифметического знака « - » осуществляется вычитание векторов, имеющих одинаковую структуру (Z = X - Y).

Например:

»v=x-y Транспонирование вектора осуществляется применением знака апострофа, который записывается сразу после имени транспонируемого вектора. Например:

Умножение вектора на число осуществляется в MatLAB с помощью знака арифметического умножения « * »

таким образом: Z = X*r или Z = r*X, где r - некоторое действительное число.

Пример:

» v = 2*x Умножение двух векторов определено в математике только для векторов одинакового размера (длины) и лишь тогда, когда один из векторов-множителей - строка, а второй - столбец. Иначе говоря, если векторы X и Y являются строками, то математическое смысл имеют лишь две формы умножения этих векторов: U = X' * Y и V = X * Y'. При этом в первом случае результатом будет квадратная матрица, а во втором - число. В МatLAB умножение векторов осуществляется применением обычного знака умножения ' * ', который записывается между множителями-векторами.

Пример:

Для трехкомпонентных векторов в MatLAB предусмотрена функция cross, которая позволяет найти векторное произведение двух векторов. При этом, если заданы два трехкомпонентных вектора v1 и v2, достаточно ввести оператор cross(v1, v2).

Пример:

» cross(v1,v2) На этом перечень допустимых математических операций с векторами исчерпывается.

В языке MatLAB предусмотрен ряд операций, которые преобразуют заданный вектор в другой того же размера и типа, хотя не являются операциями с вектором как c математическим объектом. К таким операциям относятся, например, все элементарные математические функции, приведенные в разделе 1.2.4 и которые зависят от одного аргумента. В языке MatLAB запись, например, вида Y = sin(X), где X - некоторый известный вектор, приводит к формированию нового вектора Y, имеющего тот же тип и размер, элементы которого равны синусам соответствующих элементов вектора-аргумента X. Например:

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

Добавление (отнимание) числа к (из) каждому элементу вектора. Осуществляется с помощью знака + (–).

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

Поэлементное деление векторов. Осуществляется с помощью совокупности знаков «./ ». Результат - вектор, каждый элемент которого является частным от деления соответствующего элемента первого вектора на соответствующий элемент второго вектора.

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

Поэлементное возведение в степень. Осуществляется с помощью совокупности знаков «.^». Результат - вектор, каждый элемент которого является соответствующим элементом первого вектора, возведенным в степень, величина которой равна значению соответствующего элемента второго вектора.

Примеры:

» x = [1,2,3,4,5]; y = [-2,1,4,0,5];

» disp(x + 2) » disp(y - 3) » disp(x. *y) » disp(x. /y) Warning: Divide by zero » disp(x. \y) » disp(x. ^y) Вышеуказанные операции позволяют очень просто вычислять (а затем – строить их графики) сложные математические функции, не используя при этом операторов цикла, т. е. осуществлять построение графиков в режиме калькулятора. Для этого достаточно задать значение аргумента как арифметическую прогрессию так, как это было показано в п. 1.3.1, а потом записать нужную функцию, используя знаки поэлементного преобразования векторов.

Например, пусть нужно вычислить значения функции:

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

» y = a * exp(-h*x). * sin(x) Columns 1 through Columns 8 through 1.3.5. Поэлементное преобразование матриц Для поэлементного преобразования матрицы пригодны все указанные ранее в п. 1.2.4 алгебраические функции.

Каждая такая функция формирует матрицу того же размера, что и заданная, каждый элемент которой вычисляется как указанная функция от соответствующего элемента заданной матрицы. Кроме этого, в MatLAB определены операции поэлементного умножения матриц одинакового размера (совокупностью знаков «.* », записываемой между именами перемножаемых матриц), поэлементного деления (совокупности «./ » и «.\ »), поэлементного возведения в степень (совокупность «.^ » ), когда каждый элемент первой матрицы возводится в степень, равную значению соответствующего элемента второй матрицы.

Приведем несколько примеров:

» A = [1,2,3,4,5; -2, 3, 1, 4, 0] » B = [-1,3,5,-2,1; 1,8,-3,-1,2] » sin(A) »A.*B »A./B »A.\B Warning: Divide by zero »A.^B Оригинальной в языке MatLAB является операция прибавления к матрице числа. Она записывается следующим образом: A + x, или х + A (А - матрица, а x - число). Такой операции нет в математике. В MatLAB она эквивалентна совокупности операций где Е - обозначение матрицы, все элементы которой равны единице, тех же размеров, что и матрица А. Например:

»A+ 1.3.6. Матричные действия над матрицами К матричным действиям над матрицами относят такие операции, которые используются в матричном исчислении в математике и не противоречат ему.

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

- при сложении или вычитании матрицы должны иметь одинаковые размеры;

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

Пример сложения и вычитания:

»A+B »A-B Пример умножения на число:

» 5*A » A* Пример транспонирования матрицы:

Пример умножения матрицы на матрицу:

» A' * B Функция обращения матрицы - inv(A) - вычисляет матрицу, обратную заданной матрице А. Исходная матрица А должна быть квадратной, а ее определитель не должен равняться нулю.

Приведем пример:

» inv(C) -2. 6000e-001 1. 0000e- -8. 1739e-002 3. 4783e- Проверим правильность выполнения операции обращения, применяя ее еще раз к полученному результату:

» inv(ans) -4. 0000e+001 1. 1500e+ -9. 4000e+001 2. 9900e+ Как видим, мы получили исходную матрицу С, что является признаком правильности выполнения обращения матрицы.

Возведение матрицы в целую степень осуществляется в MatLAB с помощью знака «^»: А^n. При этом матрица должна быть квадратной, а n - целым (положительным или отрицательным) числом. Это матричное действие эквивалентно умножению матрицы А на себя n раз (если n - положительно) или умножению обратной матрицы на себя (при n отрицательном).

Приведем пример:

» A^ » A^(-2) 1.5385e-001 -7.6923e-002 3.0769e- 7. 6923e-002 3. 0769e-001 -4. 6154e- 2. 1328e-018 -1. 5385e-001 3. 8462e- Оригинальными в языке MatLAB являются две новые, не определяемые в математике функции деления матриц. При этом вводятся понятие деления матриц слева напрво и деления матриц справа налево. Первая операция записывается с помощью знака « / », а вторая - « \ ».

Операция В / A эквивалентна последовательности действий B * inv(A), где функция inv осуществляет обращение матрицы. Ее удобно использовать для решения матричного уравнения:

Х А = В.

Аналогично операция A\B равносильна совокупности операций inv(A)*B, которая представляет собой решение матричного уравнения:

Для примера рассмотрим задачу отыскания корней системы линейных алгебраических уравнений:

В среде MatLAB это можно сделать таким образом:

» B = [ 14;-15;-4] »x=A\B 1.3.7. Матричные функции Вычисление матричной экспоненты (eА) осуществляется с помощью функций expm, expm1, expm2, expm3.

Эти функции следует отличать от прежде рассмотренной функции exp(A), которая формирует матрицу, каждый элемент которой равняется е в степени, равной соответствующему элементу матрицы А.

Функция expm является встроенной функцией MatLAB. Функция expm1(A) реализована как М-файл, который вычисляет матричную экспоненту путем использования разложения Паде матрицы А. Функция еxpm2(A) вычисляет матричную экспоненту, используя разложение Тейлора матрицы А. Функция expm3(A) вычисляет матричную экспоненту на основе использования спектрального разложения А.

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

» A = [1,2,3; 0, -1,5;7, -4,1] » expm(A) » expm1(A) » expm2(A) » expm3(A) Функция logm(А) осуществляет обратную операцию - логарифмирование матрицы по натуральному основанию, например:

» B = expm3(A) » logm(B) Функция sqrtm(А) вычисляет такую матрицу Y, что Y*Y = A:

» Y = sqrtm(A) »Y*Y 1.4. Функции прикладной численной математики К преимуществам системы MatLAB относится то, что она содержит в своем составе большое число функций и процедур, выполняющих стандартные математические операции, используемые в прикладной (инженерной) математике. Сюда можно отнести операции с полиномами, обработку данных измерений, функции линейной алгебры, аппроксимацию и интерполяцию данных, векторную фильтрацию и спектральный анализ сигналов.

Далее ознакомимся с важнейшими из них.

1.4.1. Операции с полиномами В системе MatLAB предусмотрены некоторые дополнительные возможности математического оперирования с полиномами.

Полином (многочлен) как функция определяется выражением:

В среде MatLAB полином задается и сохраняется в виде вектора, элементами которого являются коэффициенты полинома от an до a0 в указанном порядке:

P = [an... a2 a1 a0 ].

Введение полинома в MatLAB осуществляется так же, как и ввод вектора длиной n+1, где n - порядок полинома.

Умножение полиномов. Произведением двух полиномов степеней n и m соответственно, как известно, называют полином степени n+m, коэффициенты которого определяют простым перемножением этих двух полиномов. Фактически операция умножения двух полиномов сводится к построению расширенного вектора коэффициентов по заданным векторам коэффициентов полиномов-сомножителей. Эту операцию в математике называют сверткой векторов (а сам вектор, получаемый в результате такой процедуры - вектором-сверткой двух векторов). В MatLAB ее осуществляет функция conv(P1, P2).

Аналогично, функция deconv(P1, P2) осуществляет деление полинома P1 на полином P2, т. е. обратную свертку векторов P1 и P2. Она определяет коэффициенты полинома, который является частным от деления P на P2.

Пример :

» p = conv(p1,p2) » deconv(p,p1) В общем случае деление двух полиномов приводит к получению двух полиномов - полинома-результата (частного) и полинома-остатка. Чтобы получить оба этого полинома, следует оформить обращение к функции таким образом:

[Q,R] = deconv(B,A).

Тогда результат будет выдан в виде вектора Q с остатком в виде вектора R таким образом, что будет выполнено соотношение B = conv(A,Q) + R.

В системе MatLAB предусмотрена функция roots(P), которая вычисляет вектор, элементы которого являются корнями заданного полинома Р.

Пусть нужно найти корни полинома:

Ниже показано, как просто это сделать:

» p = [1,8,31,80,94,20];

» disp(roots(p)) Обратная операция - построение вектора р коэффициентов полинома по заданному вектору его корней - осуществляется функцией poly:

p = poly(r).

Здесь r - заданный вектор значений корней, p - вычисленный вектор коэффициентов полинома. Приведем пример:

» p = [1,8,31,80,94,20] » r = roots(p) » p1 = poly(r) Заметим, что получаемый вектор не содержит старшего коэффициента, который по умолчанию полагается равным единице.

Эта же функция в случае, если аргументом ее является некоторая квадратная матрица А размером (n n), строит вектор характеристического полинома этой матрицы. Обращение p = poly(A) формирует вектор p коэффициентов характеристического полинома p(s) = det(s*E - A) = p1sn +... + pn*s + pn+1, где Е - обозначение единичной матрицы размером (n n).

Рассмотрим пример:

» p = poly(A) Для вычисления значения полинома по заданному значению его аргумента в MatLAB предусмотрена функция polyval. Обращение к ней осуществляется по схеме:

y = polyval(p,x), где p - заданный вектор коэффициентов полинома, а x - заданное значение аргумента. Пример:

» y = polyval(p,2) Если в качестве аргумента полинома указана матрица Х, то функция polyval(p,X) вычисляет матрицу Y, каждый элемент которой является значением указанного полинома при значении аргумента, равном соответствующему элементу матрицы Х, например:

» disp(polyval(p,X)) В этом случае функция вычисляет значение полинома для каждого элемента матрицы Х, и поэтому размеры исходной и конечной матриц одинаковы size(Y) = size(X).

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

dp = polyder(p) по заданному полиному р вычисляет вектор dp, элементы которого являются коэффициентами полинома-производной от заданного:

» dp = polyder(p) dp = polyder(p1,p2) вычисляет вектор dp, элементы которого являются коэффициентами полиномапроизводной от произведения двух полиномов р1 и р2:

» p1 = [1,8,31,80,94,20];

» p2 = [1,2,16];

» p = conv(p1,p2) » dp = polyder(p) » dp1 = polyder(p1,p2) [q,p] = polyder(p1,p2) вычисляет производную от отношения (p1/p2) двух полиномов р1 и р2 и выдает результат в виде отношения (q/p) полиномов q и p:

» p1 = [1,8,31,80,94,20];

» p2 = [1,2,16];

» [q,p] = polyder(p1,p2) » z = deconv (q,p) » y = deconv(p1,p2) » z1 = polyder(y) 1.4.2. Обработка данных измерений Система MatLAB дает пользователю дополнительные возможности для обработки данных, которые заданы в векторной или матричной форме.

Допустим, что есть некоторая зависимость y(x), заданная рядом точек Ее можно задать в командном окне MatLAB как матрицу xydata, содержащую две строки - значения x и значения y:

>> xydata =[2 4 6 8 10; 5.5 6.3 6.8 8 8.6] На примере этой зависимости рассмотрим основные средства для обработки данных.

Функция size(xydata) предназначена для определения числа строк и столбцов матрицы xydata. Она формирует вектор [n, p], содержащий эти величины:

>> size(xydata) Обращение к ней вида >> [n, p] = size(xydata);

позволяет сохранить в памяти машины и использовать потом при дальнейших вычислениях данные о числе строк n и столбцов p этой матрицы:

>> n, p С помощью этой функции можно установить длину и тип (строка или столбец) вектора:

» v = xydata(:) » n = size(v) » v1 = v' » size(v') Функция max(V), где V - некоторый вектор, выдает значение максимального элемента этого вектора. Аналогично, функция min(V) извлекает минимальный элемент вектора V. Функции mean(V) и std(V) определяют, соответственно, среднее значение и среднеквадратичное (стандартное) отклонение от него значений элементов вектора V.

Функция сортировки sort(V) формирует вектор, элементы которого расположены в порядке возрастания их значений.

Функция sum(V) вычисляет сумму элементов вектора V.

Функция prod(V) выдает произведение всех элементов вектора V.

Функция cumsum(V) формирует вектор того же типа и размера, любой элемент которого является суммой всех предшествующих элементов вектора V (вектор кумулятивной суммы).

Функция cumprod(V) создает вектор, элементы которого являются произведением всех предшествующих элементов вектора V.

Функция diff(V) создает вектор, который имеет размер на единицу меньший размера вектора V, а элементы его являются разностью между соседними элементами вектора V.

Применение описанных функций проиллюстрировано ниже.

» v = [ 1, 0.1, 0.5, 0.1, 0.1,0.4 ];

» disp(size(v)) » disp(max(v)) » disp(min(v)) » disp(mean(v)) » disp(std(v)) » disp(sort(v)) » disp(sum(v)) » disp(prod(v)) » disp(cumsum(v)) » disp(cumprod(v)) » disp(diff(v)) Если указать второй выходной параметр, то можно получить дополнительную информацию об индексе первого элемента, значение которого является максимальным или минимальным:

>> [M,n]=max(v) >> [N,m]=min(v) Интегрирование методом трапеций осуществляет процедура trapz. Обращение к ней вида trapz(x,y) приводит к вычислению площади под графиком функции y(x), в котором соседние точки, заданные векторами х и у, соединены отрезками прямых. Если первый вектор х не указан в обращении, по умолчанию допускается, что шаг интегрирования равняется единице (т. е. вектор х представляет собой вектор из номеров элементов вектора у).

Пример. Вычислим интеграл от функции y = sin(x) в диапазоне от 0 до. Его точное значение равно 2.

Возьмем равномерную сетку аргумента из 100 элементов. Тогда вычисления сведутся к совокупности операций:

» x = (0 : 0.01:1)*pi;

» y = sin(x);

» disp(trapz(x,y)) Те же функции size, max, min, mean, std, sort, sum, prod, cumsum, cumprod, diff могут быть применены и к матрицам. Основным отличием использования в качестве аргументов этих функций именно матриц является то, что соответствующие описанные выше операции ведутся не по отношению к строкам матриц, а к каждому из столбцов заданной матрицы. Т. е. каждый столбец матрицы А рассматривается как переменная, а каждая строка - как отдельное наблюдение. Так, в результате применения функций max, min, mean, std получаются векторы-строки с количеством элементов, которое равняется количеству столбцов заданной матрицы. Каждый элемент содержит, соответственно, максимальные, минимальное, среднее или среднеквадратичное значения элементов соответствующего столбца заданной матрицы.

Приведем примеры. Пусть имеем 3 величины y1, y2 и y3, измеренные при некоторых пяти значениях аргумента (они не указаны). Тогда данные измерений образуют 3 вектора по 5 элементов:

>> y1 = [ 5.5 6.3 6.8 8 8.6];

>> y2 = [-1. 2 0.5 -0. 6 1 0.1];

>> y3 = [ 3.4 5.6 0 8.4 10.3] ;.

Сформируем из них матрицу измерений так, чтобы векторы y1, y2 и y3 образовывали столбцы этой матрицы:

» A = [ y1', y2', y3'] Применим к этой матрице измерений описанные функции. Получим » size(A) » max(A) » min(A) » mean(A) » std(A) Если при обращении к функциям max и min указать второй выходной параметр, то он даст информацию о номерах строк, где находятся в соответствующем столбце первые элементы с максимальным (или минимальным) значением. Например:

>> [M,n]=max(A) M= 8.6000 1.0000 10. >> [N,m]=min(A) Функция sort сортирует элементы любого из столбцов матрицы. Результатом является матрица того же размера.

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

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

Наконец, функция diff создает из заданной матрицы размером (m n) матрицу размером ((m-1) n), элементы которой являются разностью между элементами соседних строк начальной матрицы.

Применяя эти процедуры к принятой матрице измерений, получим:

» sort(A) » sum(A) » prod(A) » cumsum(A) » cumprod(A) » diff(A) Рассмотрим некоторые другие функции, предоставляемые пользователю системой MatLAB.

Функция cov(A) вычисляет матрицу ковариаций измерений. При этом получают симметричную квадратную матрицу с количеством строк и столбцов, равным количеству измеренных величин, т. е. количеству столбцов матрицы измерений. Например, при применении к принятой матрице измерений она дает такой результат:

» cov(A) На диагонали матрицы ковариаций размещены дисперсии измеренных величин, а вне ее - взаимные корреляционные моменты этих величин.

Функция corrcoeff(A) вычисляет матрицу коэффициентов корреляции при тех же условиях. Элементы матрицы S = corrcoef(A) связаны с элементами матрицы ковариаций C=cov(A) таким соотношением:

Пример:

» corrcoef(A) 1.4.3. Функции линейной алгебры Традиционно к линейной алгебре относят такие задачи, как обращение и псевдообращение матрицы, спектральное и сингулярное разложения матриц, вычисление собственных значений и векторов, сингулярных чисел матриц, вычисление функций от матриц. Ознакомимся с некоторыми основными функциями MatLAB в этой области.

Функция k = cond(A) вычисляет и выдает число обусловленности матрицы относительно операции обращения, которое равняется отношению максимального сингулярного числа матрицы к минимальному.

Функция k = norm(v,p) вычисляет р-норму вектора v по формуле:

k = sum(abs(v). p)(1/p), где р - целое положительное число. Если аргумент р при обращении к функции не указан, вычисляется 2норма (р=2).

Функция k = norm(А,p) вычисляет р-норму матрицы, где р = 1,2, 'fro' или inf. Если аргумент р не указан, вычисляется 2-норма. При этом справедливы такие соотношения:

norm(A,1) = max(sum(abs(A)));

norm(A,inf) = max(sum(abs(A')));

norm(A,’fro') = sqrt(sum(diag(A'*A)));

norm(A) = norm(A,2) = max(A).

Функция rd = rcond(А) вычисляет величину, обратную значению числа обусловленности матрицы А относительно 1-нормы. Если матрица А хорошо обусловлена, значение rd близко к единице. Если же она плохо обусловленная, rd приближается к нулю.

Функция r = rank(A) вычисляет ранг матрицы, который определяется как количество сингулярных чисел матрицы, превышающие порог max(size(A))*norm(A)*eps.

Приведем примеры применения этих функций:

» disp(cond(A)) » disp(norm(A,1)) » disp(norm(A)) » disp(rcond(A)) » disp(rank(A)) Процедура d = det(A) вычисляет определитель квадратной матрицы на основе треугольного разложения методом исключения Гаусса.

Функция t = trace(A) вычисляет след матрицы А, равный сумме ее диагональных элементов.

Q = null(A) вычисляет ортонормированный базис нуль-пространства матрицы А.

Q = orth(A) выдает ортонормированный базис матрицы А.

Процедура R = rref(A) осуществляет приведение матрицы к треугольному виду на основе метода исключения Гаусса с частичным выбором ведущего элемента.

Примеры:

» disp(det(A)) » disp(trace(A)) » disp(null(A)) » disp(orth(A)) » disp(rref(A)) Функция R=chol(A) осуществляет разложение Холецького для симметричных действительных и комплексных эрмитовых матриц. Например:

» disp(chol(A)) Функция lu(A) осуществляет LU-разложение матрицы А в виде произведения нижней треугольной матрицы L (возможно, с перестановками) и верхней треугольной матрицы U, так что A = L * U.

Обращение к этой функции вида [ L, U, P ] = lu(A) позволяет получить три составляющие этого разложения - нижнюю треугольную матрицу L, верхнюю треугольную U и матрицу перестановок P такие, что P * A =L * U.

Приведем пример:

» disp(lu(A)) » [ L, U, P] = lu(A);

Из него вытекает, что в первом, упрощенном варианте обращения функция выдает комбинацию из матриц L и Обращение матрицы осуществляется с помощью функции inv(А):

» disp(inv(A)) Процедура pinv(А) находит матрицу, псевдообратную матрице А, которая имеет размеры матрицы А’ и удовлетворяет условиям Например:

» P = pinv(A) » A*P*A, » P*A*P Для квадратных матриц эта операция равнозначна обычному обращению.

Процедура [ Q, R, P ] = qr(A) осуществляет разложение матрицы А на три - унитарную матрицу Q, верхнюю треугольную R с диагональными элементами, уменьшающимися по модулю, и матрицу перестановок P - такие что Например:

» [Q,R,P] = qr(A) Определение характеристического полинома матрицы A можно осуществить с помощью функции poly(A).

Обращение к ней вида p=poly(A) дает возможность найти вектор-строку p коэффициентов характеристического полинома p(s) = det(s*E - A) = p1*sn +... + pn*s + pn+1, где Е - обозначение единичной матрицы размером (n n). Например :

» p = poly(A) Вычисление собственных значений и собственных векторов матрицы осуществляет процедура eig(А).

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

[ R, D ] = eig(A), то в результате получают диагональную матрицу D собственных значений и матрицу R правых собственных векторов, которые удовлетворяют условию A*R = R * D.

Эти векторы являются нормированными так, что норма любого из них равна единице. Приведем пример:

» disp(eig(A)) » [ R,D] = eig(A) Сингулярное разложение матрицы осуществляет процедура svd(А). Упрощенное обращение к ней позволяет получить сингулярные числа матрицы А. Более сложное обращение вида:

[ U, S, V ] = svd(A) позволяет получить три матрицы - U, которая состоит из ортонормированных собственных векторов, отвечающих наибольшим собственным значениям матрицы А*АT; V - из ортонормированных собственных векторов матрицы АT*А и S - диагональную матрицу, которая содержит неотрицательные значения квадратных корней из собственных значений матрицы АT*А (их называют сингулярными числами). Эти матрицы удовлетворяют соотношению:

Рассмотрим пример:

» disp(svd(A)) » [U,S,V] = svd(A) Приведение матрицы к форме Хессенберга осуществляется процедурой hess(А). Например:

» disp(hess(A)) Более развернутое обращение [P,H] = hess(A) дает возможность получить, кроме матрицы Н в верхней форме Хессенберга, также унитарную матрицу преобразований Р, которая удовлетворяет условиям:

Пример:

» [P,H] = hess(A) Процедура schur (А) предназначена для приведения матрицы к форме Шура. Упрощенное обращение к ней приводит к получению матрицы в форме Шура.

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

Обращение [U,T] = schur(A) позволяет, кроме матрицы Т Шура, получить также унитарную матрицу U, удовлетворяющую условиям:

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

Приведем пример:

» disp(schur(A)) » [U,T] = hess(A) Функция [U,T] = rsf2csf(U,T) преобразует действительную квазитреугольную форму Шура в комплексную треугольную:

» [U,T] = rsf2csf(U,T) Процедура [AA, BB, Q, Z, V] = qz(A,B) приводит пару матриц А и В к обобщенной форме Шура. При этом АА и ВВ являются комплексными верхними треугольными матрицами, Q, Z - матрицами приведения, а V - вектором обобщенных собственных векторов такими, что Обобщенные собственные значения могут быть найдены, исходя из такого условия:

A * V * diag(BB) = B * V * diag(AA).

Необходимость в одновременном приведении пары матриц к форме Шура возникает в многих задачах линейной алгебры - решении матричных уравнений Сильвестра и Риккати, смешанных систем дифференциальных и линейных алгебраических уравнений.

Пример.

Пусть задана система обычных дифференциальных уравнений в неявной форме Коши с одним входом и одним выходом y такого вида:

причем матрицы Q, R и векторы b, c и d равны соответственно 31. 0. d = -0. Необходимо вычислить значения полюсов и нулей соответствующей передаточной функции.

Эта задача сводится к отысканию собственных значений, которые удовлетворяют матричным уравнениям:

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

Ниже приведена совокупность операторов, которая приводит к расчету полюсов:

» [AA, BB] = qz(R,-Q) » diag(AA). /diag(BB) Расчет нулей осуществляется таким образом:

» [AA,BB] = qz(A,B) » diag(AA). /diag(BB) Вычисление собственных значений матричного полинома осуществляет процедура polyeig. Обращение [ R, d ] = polyeig(A0, A1,..., Ap ) позволяет решить полную проблему собственных значений для матричного полинома степени р вида (A0 + *A1 +... + p *Ap)*r = 0.

Входными переменными этой процедуры являются р+1 квадратные матрицы A0, A1,... Ap порядка n. Исходными переменными - матрица собственных векторов R размером (n (n p)) и вектор d собственных значений размером (n p).

Функция polyvalm предназначена для вычисления матричного полинома вида по заданному значению матрицы Х и вектора p = [pn, pn-1,..., p0] коэффициентов полинома. Для этого достаточно обратиться к этой процедуре по схеме:

Пример:

» disp(polyvalm(p,X)) Примечание. Следует различать процедуры polyval и polyvalm. Первая вычисляет значение полинома для любого из элементов матрицы аргумента, а вторая при вычислении полинома возводит в соответствующую степень всю матрицу аргумента.

Процедура subspace(А,В) вычисляет угол между двумя подпространствами, которые "натянуты на столбцы" матриц А и В. Если аргументами являются не матрицы, а векторы A и B, вычисляется угол между этими векторами.

1.4.4. Аппроксимация и интерполяция данных Полиномиальная аппроксимация данных измерений, которые сформированы как некоторый вектор Y, при некоторых значениях аргумента, которые образуют вектор Х такой же длины, что и вектор Y, осуществляется процедурой polyfit(X, Y, n). Здесь n - порядок аппроксимирующего полинома. Результатом действия этой процедуры является вектор длиной (n +1) из коэффициентов аппроксимирующего полинома.

Пусть массив значений аргумента имеет вид:

x = [-0.45 -0.35 -0.25 -0.15 -0.05 0.05 0.15 0.25 0.35 0.45];

а массив соответствующих значений измеренной величины - вид:

y = [-1.1 0.2 0.1 0.8 0.5 0.2 0.4 0.1 0.1 -0.6];

Тогда, применяя указанную функцию при разных значениях порядка аппроксимирующего полинома, получим:

» x = -0.45:0.1:0.45;

» y = [-1.1 0.2 0.1 0.8 0.5 0.2 0.4 0.1 0.1 -0.6];

» polyfit(x,y,1) » polyfit(x,y,2) » polyfit(x,y,3) » polyfit(x,y,4) Это означает, что заданную зависимость можно аппроксимировать или прямой y ( x) = 0,13939x + 0,07, или квадратной параболой y ( x) = 6,1364 x 2 + 0,13939 x + 0,57625, или кубической параболой y ( x) = 7,4592 x 3 6,1364 x 2 0,95338 x + 0,57625, или параболой четвертой степени y ( x) = 35,548 x 4 + 7,4592 x 3 + 1,1509 x 2 0,95338 x + 0,40469.

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

x = -0.45:0.1:0.45;

y = [-1.1 0.2 0.1 0.8 0.5 0.2 0.4 0.1 0.1 -0.6];

stem(x,y); hold on P1=polyfit(x,y,1) P2=polyfit(x,y,2) P3=polyfit(x,y,3) P4=polyfit(x,y,4) x1 = -0.5 : 0.01 : 0.5;

y1=polyval(P1,x1);

y2=polyval(P2,x1);

y3=polyval(P3,x1);

y4=polyval(P4,x1);

plot(x1,y1, '--',x1,y2,':',x1,y3, '.',x1,y4), grid, set(gca, 'FontSize', 12), title('Полиномиальная аппроксимация ');

xlabel('Аргумент'); ylabel('Функция');

legend('исходные','данные','Аппроксимация','линейная','квадратическая','кубическая','четв ертой степени',0) Функция spline(X,Y,Xi) осуществляет интерполяцию кубическими сплайнами. При обращении Yi = spline(X,Y,Xi) она интерполирует значение вектора Y, заданного при значениях аргумента, представленных в векторе Х, и выдает значение интерполирующей функции в виде вектора Yi при значениях аргумента, заданных вектором Xi. В случае, если вектор Х не указан, по умолчанию принимается, что он имеет длину вектора Y и любой его элемент равен номеру этого элемента.

В качестве примера рассмотрим интерполяцию вектора x = -0.45:0.1:0.45;

y = [-1.1 0.2 0.1 0.8 0.5 0.2 0.4 0.1 0.1 -0.6];

x1 = -0.46:0.01:0.46;

y2 = spline(x,y,x1);

stem(x,y); hold on plot (x1,y2, '.'), grid set(gca, 'FontSize',12), title('Интерполяция процедурой SPLINE ');

xlabel('Аргумент'); ylabel('Функция') Результат приведен на рис. 1.19.

Одномерную табличную интерполяцию осуществляет процедура interp1. Обращение к ней в общем случае имеет вид:

Yi = interp1(X,Y,Xi,’’), и позволяет дополнительно указать метод интерполяции в четвертом входном аргументе:

'nearest' ступенчатая интерполяция;

‘spline' кубическими сплайнами.

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

y = [-1.1 0.2 0.1 0.8 0.5 0.2 0.4 0.1 0.1 -0.6];

x1 = -0.46:0.01:0.46;

y1 = interp1(x,y,x1);

y4 = interp1(x,y,x1,'nearest');

y2 = interp1(x,y,x1,'cubic');

y3 = interp1(x,y,x1,'spline');

plot (x1,y1, '--',x1,y2,'.',x1,y3,x1,y4,':'), grid set(gca, 'FontSize',12), legend('линейная','кубическая','сплайновая','ступенчатая',0) title('Интерполяция процедурой INTERP1 ');

xlabel('Аргумент'); ylabel('Функция') Результаты приведены на рис.1.20.

1.4.5. Векторная фильтрация и спектральный анализ В системе MatLAB есть несколько функций для проведения цифрового анализа данных наблюдений (измерений).

Так, функция y = filter(b,a,x) обеспечивает формирование вектора y по заданным векторам b, a, x в соответствии с соотношением:

y(k) = b(1)*x(k) + b(2)*x(k-1) +... + b(nb+1)*x(k-nb) --a(2)*y(k-1) - a(3)*y(k-3) -... - a(na+1)*y(k-na), (1.1) где вектор b имеет такой состав b = [ b(1), b(2),..., b(nb+1)], a вектор а a = [1, a(2), a(3),..., a(na+1)].

Соотношение (1) можно рассматривать как конечно-разностное уравнение фильтра с дискретной передаточной функцией вида рациональной дроби, коэффициенты числителя которого образовывают вектор b, а знаменателя - вектор а, на вход которого подается сигнал x(t), а на выходе формируется сигнал y(t).

Тогда вектор y будет представлять собой значение исходного сигнала этого фильтра в дискретные моменты времени, соответствующие заданным значениям входного сигнала x(t) (вектор х).

Ниже приведен пример применения функции filter.

» y = filter(b,a,x) Функции fft (Fast Fourier Transformation) и ifft (Inverse Fast Fourier Transformation) осуществляют преобразования заданного вектора, соответствующие дискретному прямому и обратному преобразованиям Фурье.

Обращение к этим функциям вида:

приводит к формированию вектора y в первом случае и х - во втором по формулам:

где j - обозначение мнимой единицы; n - число элементов заданного вектора х (оно представляет также размер выходного вектора y).

Приведем пример. Сформируем входной сигнал в виде вектора, элементы которого равняются значениям функции, являющейся суммой двух синусоид с частотами 5 и 12 Гц (рис. 1.21).

t =0:0.001:2;

x = sin(2*pi*5*t) + cos(2*pi*12*t);

figure, plot(t, x); grid, set(gcf,'color','white'), set(gca,'FontName','Arial Cyr','FontSize',16), title('Входной процесс '); xlabel('Время (с)'); ylabel('Х(t)') Найдем Фурье-изображение этого сигнала и выведем графическое представление модуля его Фурьеизображения:

y = fft(x); a =abs(y);

plot(a); grid, set(gca,'FontName','Arial Cyr','FontSize',16), title('Модуль Фурье - изображения ');

xlabel('Номер элемента вектора'); ylabel('abs(F(X(t))') Результат отображен на рис. 1.22.

Теперь осуществим обратное преобразование с помощью функции ifft и результат также выведем в форме графика:

z = ifft(y); plot(t, z);

grid, set(gca,'FontName','Arial Cyr','FontSize',16), title('Обратное Фурье-преобразование ');

xlabel('Время (с)'); ylabel('Z(t)') На рис. 1.23 изображен результат. Нетрудно убедиться, что воспроизведенный процесс совпадает с исходным.

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

а) номер m отвечает моменту времени tm, в который измерен входной сигнал x(m); при этом t б) номер k - это индекс значения частоты fk, которому отвечает найденный элемент y(k) дискретного преобразования Фурье;

в) чтобы перейти от индексов к временной или частотной области, необходимо знать значение h дискрета (шага) времени, через который измерен входной сигнал x(t) и промежуток T времени, на протяжении которого он измеряется; тогда шаг (дискрет) по частоте в изображении Фурье определится соотношением:

а диапазон изменения частоты - формулой так, в анализируемом примере (h =0. 001, T = 2, n = 21):

г) из (2) следует, что индексу k = 1 отвечает нулевое значение частоты (f0 = 0); иначе говоря, первый элемент вектора y(1) является значением Фурье-изображения при нулевой частоте, т. е. - просто суммой всех заданных значений вектора x; отсюда получаем, что вектор y(k) содержит значение Фурье-изображения, начиная из частоты f0 = 0 (которой отвечает k = 1), до максимальной частоты fmax = F (которой отвечает k = n); таким образом, Фурье-изображение определяется функцией fft только для положительных частот в диапазоне от 0 до F; это неудобно для построения графиков Фурье-изображения от частоты; более удобным и привычным представляется переход к вектору Фурье-изображения, определенному в диапазоне частот от (-F/2) до F/2; частота F N = F / 2 получила название частоты Найквиста;

д) как известно, функция e j z является периодической по z с периодом 2 ; поэтому информация об Фурьеизображении при отрицательных частотах расположена во второй половине вектора y(k).

Сформируем для анализируемого примера массив частот, исходя из вышесказанного:

f = 0 : 0.5 : 1000;

и выведем график с аргументом-частотой (рис. 1.24):

plot(f,a); grid, set(gca,'FontName','Arial Cyr','FontSize',16), title('Модуль Фурье - изображения ');

xlabel('Частота (Гц)'); ylabel('abs(F(X(t))') Как следует из рассмотрения рис. 1.24, по нему непросто распознать те частоты (5 и 12 Гц), с которыми изменяется входной сигнал. Это - следствие того обстоятельства, которое было отмечено в примечании г). Чтобы определить частотный спектр входного сигнала, нужно сначала преобразовать полученный вектор y Фурьеизображения с помощью процедуры fftshift.

Функция fftshift (обращение к ней осуществляется таким образом: z = fftshift(y)) предназначена для формирования нового вектора z из заданного вектора y путем перестановки второй половины вектора y в первую половину вектора z. При этом вторая половина вектора z состоит из элементов первой половины вектора y. Более точно эту операцию можно задать соотношениями:

z(1) = y(n/2+1);..., z(k) = y(n/2+k);..., z(n/2) = y(n); z(n/2+1) = y(1);...

..., z(n/2+k) = y(k);... z(n) = y(n/2).

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

f1 = -500 : 0.5 : 500; % Перестройка вектора частот v = fftshift(y); % Перестройка вектора Фурье-изображения plot(f1(970:1030),a(970:1030)); grid, % Вывод графика set(gca,'FontName','Arial Cyr','FontSize',16), title('Модуль Фурье - изображения');

xlabel('Частота (Гц)'); ylabel('abs(F(X(t))') Из графика рис. 1.25 уже становится очевидным, что в спектре входного сигнала есть две гармоники - с частотами 5 и 12 Гц.

Остается лишь то неудобство, что из графика спектра невозможно установить амплитуды этих гармоник. Во избежание этого, нужно весь вектор y Фурье-изображения разделить на число его элементов (n), чтобы получить вектор комплексного спектра сигнала:

N=length(y); a =abs(v)/N;

plot(f1(970:1030),a(970:1030)); grid set(gca,'FontName','Arial Cyr','FontSize',14,'Color','white'), title('Модуль комплексного спектра');

xlabel('Частота (Гц)');ylabel('abs(F(X(t)) / N') Результат приведен на рис. 1.26. Из его рассмотрения вытекает, что "амплитуды" всех составляющих гармоник равны 0.5. При этом нужно принять во внимание, что "амплитуды" распределены между положительными и отрицательными частотами поровну, поэтому они вдвое меньшее действительной амплитуды соответствующей гармоники.

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

1.5.1. Процедура plot Вывод графиков в системе MatLAB - настолько простая и удобная операция, что ее можно использовать даже в режиме калькулятора.

Основной функцией, обеспечивающей построение графиков на экране дисплея, является функция plot. Общая форма обращения к ней такова:

plot(x1, y1, s1, x2, y2, s2,...).

Здесь x1, y1 - заданные векторы, элементами которых являются массивы значений аргумента (х1) и функции (y1), отвечающие первой кривой графика; x2, y2 - массивы значений аргумента и функции второй кривой и т.д.

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

а) тип линии, которая соединяет отдельные точки графика; б) тип точки графика; в) цвет линии. Если переменные s не указаны, то тип линии по умолчанию - отрезок прямой, тип точки - пиксел, а цвет устанавливается в такой очередности:

- синий, зеленый, красный, голубой, фиолетовый, желтый, черный и белый - в зависимости от того, какая по очереди линия выводится на график. Например, обращение вида plot(x1,y1,x2,y2,...) приведет к построению графика, в котором первая кривая будет линией из отрезков прямых синего цвета, вторая кривая - такого же типа зеленой линией и т.д.

Графики в MatLAB всегда выводятся в отдельное графическое окно, которое называют фигурой.

Приведем пример. Пусть нужно вывести график функции y = 3sin(x + /3) на промежутке от -3 до +3 с шагом /100.

Сначала надо сформировать массив значений аргумента х:

потом вычислить массив соответствующих значений функции:

и, наконец, построить график зависимости y(х).

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

» x = -3*pi:pi/100:3*pi;

» y = 3*sin(x+pi/3);

» plot(x,y) В результате на экране появится окно с графиком (см. рис. 1.27).

Если вектор аргумента при обращении к функции plot не указан явно, то система по умолчанию принимает в качестве аргумента номера элементов вектора функции. Например, если ввести команду » plot(y), то результатом будет появление графика в виде, приведенном на рис. 1.28.

Графики, приведенные на рис. 1.27, 1.28, имеют следующие недостатки:

- на них не нанесена сетка из координатных линий, что затрудняет чтение графиков;

- нет общей информации о кривой графика (заголовка);

- неизвестно, какие величины отложены по осям графика.

Первый недостаток устраняется с помощью функции grid. Если к этой функции обратиться сразу после обращения к функции plot:

» x = -3*pi:pi/100:3*pi;

» y = 3*sin(x+pi/3);

» plot(x,y), grid, то график будет снабжен координатной сеткой (рис. 1.29).

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

Заголовок графика выводится с помощью процедуры title. Если после обращения к процедуре plot вызвать title таким образом:

title(‘’), то над графиком появится текст, записанный между апострофами в скобках. При этом следует помнить, что текст всегда должен помещаться в апострофы.

Аналогично можно вывести объяснения к графику, которые размещаются вдоль горизонтальной оси (функция xlabel) и вдоль вертикальной оси (функция ylabel).

Например, совокупность операторов » x = -3*pi : pi/100 : 3*pi;

» y = 3*sin(x+pi/3);

» plot(x,y), grid » title('Функция y = 3*sin(x+pi/3)');

» xlabel('x'); ylabel('y');

приведет к оформлению поля фигуры в виде, представленном на рис. 1.30. Очевидно, такая форма уже целиком удовлетворяет требованиям, предъявляемым к инженерным графикам.

Не сложнее вывод в среде MatLAB графиков функций, заданных параметрически. Пусть, например, необходимо построить график функции y(х), заданной формулами:

x = 4 e-0,05 t sin t; y = 0,2 e-0,1 t sin 2t.

Выберем диапазон изменения параметра t от 0 до 50 с шагом 0.1. Тогда, набирая совокупность операторов t = 0:0.1:50; x = 4*exp(-0.05*t).*sin(t);

y = 0.2*exp(-0.1*t).*sin(2*t);

plot(x,y), grid, set(gcf,'color','white') title('Параметрическая функция x=4*exp(-0.05t)*sin(t);… y= 0.2*exp(-0.1t)*sin(2t) ') получим график рис. 1.31.

1.5.2. Специальные графики Большим удобством, предоставляемым системой MatLAB, является указанная ранее возможность не указывать аргумент функции при построении ее графика. В этом случае в качестве аргумента система принимает номер элемента вектора, график которого строится. Пользуясь этим, например, можно построить "график вектора":

» plot (x), grid, title('График вектора Х') » ylabel('Значение элементов'), xlabel(' Номер элемента').

Результат представлен на рис. 1.32.

Еще более наглядным является представление вектора в виде столбцовой диаграммы с помощью функции bar (см. рис. 1.33):

» bar(x), title('График вектора Х') » xlabel(' Номер элемента'), ylabel('Значение элементов') Если функция задана своими значениями при дискретных значениях аргумента, и неизвестно, как она может изменяться в промежутках между значениями аргумента, удобнее представлять график ее в виде отдельных вертикальных линий для любого из заданных значений аргумента. Это можно сделать, применяя процедуру stem, обращение к которой целиком аналогично обращению к процедуре plot:

stem(x,'k'), grid, set(gca,'FontSize',14), title('График вектора Х') ylabel('Значение элементов'), xlabel(' Номер элемента') На рис. 1.34 изображен полученный при этом график.

Другой пример - построение графика функции в виде столбцовой диаграммы (рис. 1.35):

x = - 2.9 : 0.2 : 2.9; y=exp(-x.*x) bar(x, y), set(gca,'FontSize',14) title('Столбцовая диаграмма функции y = exp(-x^2)') xlabel (' Аргумент х'), ylabel (' Значение функции у') Еще одна полезная инженеру функция - hist (построение графика гистограммы заданного вектора). Стандартное обращение к ней таково:

hist(y, x), где y - вектор, гистограмму которого нужно построить; х - вектор, элементы его определяют интервалы изменения первого вектора, внутри которых подсчитывается количество элементов вектора y.

Эта функция осуществляет две операции:

- подсчитывает количество элементов вектора y, значения которых попадают внутрь соответствующего диапазона, указанного вектором x;

- строит столбцовую диаграмму подсчитанных чисел элементов вектора y как функцию диапазонов, указанных вектором х.

В качестве примера рассмотрим построение гистограммы случайных величин, которые формируются встроенной функцией randn. Примем общее количество элементов вектора этих случайных величин 10 000. Построим гистограмму для диапазона изменения этих величин от -2,9 до +2,9. Интервалы изменения пусть будут равны 0,1. Тогда график гистограммы можно построить с помощью совокупности таких операторов:

x = -2.9:0.1:2.9; y = randn(10000,1);

hist(y,x), set(gca,'fontsize',14) ylabel('Количество из 10000'), xlabel('Аргумент') title('Гистограмма нормального распределения') Результат представлен на рис. 1.36. Из него, в частности, вытекает, что встроенная функция randn достаточно верно отображает нормальный гауссовый закон распределения случайной величины.

Процедура comet(x,y) («комета») строит график зависимости y(х) постепенно во времени в виде траектории кометы. При этом изображающая точка на графике имеет вид маленькой кометы (с головкой и хвостиком), которая плавно перемещается от одной точки к другой. Например, если ввести совокупность операторов:

t = 0:0.1:50;

x = 4 * exp(-0.05*t).* sin(t);

y = 0.2 * exp(-0. 1*t). * sin(2*t);

comet(x,y), то график, приведенный на рис. 1.31, будет построен как траектория последовательного движения кометы. Это обстоятельство может быть полезным при построении пространственных траекторий для выявления характера изменения траектории с течением времени.

MatLAB имеет несколько функций, которые позволяют строить графики в логарифмическом масштабе. К примеру, функция logspace с обращением x = logspace(d1, d2, n) формирует вектор-строку х, содержащую n равноотстоящих в логарифмическом масштабе друг от друга значений, которые покрывают диапазон от 10d1 до 10d2.

Функция loglog полностью аналогична функции plot, но графики по обеим осям строятся в логарифмическом масштабе.

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

В качестве примера рассмотрим построение графиков амплитудно-частотной и фазо-частотной характеристик звена, описываемого передаточной функцией:

Для этого нужно, во-первых, создать полином числителя Pc = [1 4] и знаменателя передаточной функции Pz = [1 4 100]. Во-вторых, определить корни этих двух полиномов:

» P1 = [1 4]; P2 = [1 4 100];

» roots(P1) » roots(P2) -2. 0000e+000 +9. 7980e+000i -2. 0000e+000 -9. 7980e+000i В-третьих, задать диапазон изменения частоты так, чтобы он охватывал все найденные корни:

om0 = 1e-2; omk = 1e2.

Теперь нужно задаться количеством точек будущего графика (например, n = 41), и сформировать массив точек по частоте в логарифмическом масштабе OM = logspace(-2,2,41), где значения -2 и +2 отвечают десятичным порядкам начального om0 и конечного omk значений частоты.

Пользуясь функцией polyval, можно вычислить сначала вектор "сh" комплексных значений числителя частотной передаточной функции, отвечающих заданной передаточной функции по Лапласу, если в качестве аргумента функции polyval использовать сформированный вектор частот ОМ, элементы которого умножены на мнимую единицу. Аналогично вычисляется комплекснозначный вектор "zn" знаменателя ЧПФ.

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

Наконец, для построения графика АЧХ в логарифмическом масштабе, достаточно применить функцию loglog, а для построения ФЧХ удобнее воспользоваться функцией semilogx.

В целом последовательность действий может быть такой:

P1 = [1 4]; P2 = [1 4 100]; OM = logspace(-2,2,40); p=i*OM;

ch = polyval(P1,p); zn = polyval(P2,p);