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

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

Получение краткой справки (help ) в командном окне MATLAB производится с помощью команд:

>> help - вывод сведений о разделах (topics) MATLAB с возможностью

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

>> help - вывод сведений об именах (названиях) функций, входящих в раздел.

>> help <имя функции> - вывод справочных сведений по функции.

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

>> lookfor - вывод справочных сведений по ключевому слову.

>> help demos - вывод списка демонстрационных примеров.

>> hthelp- открывает интерактивное окноMATLAB help .

>> help symbolic– выводит сведения об инструментарии символьной математики (symbolicmathtoolbox)MATLAB.

>> help signal processing toolbox – выводит сведения о функциях пакета Signal Processing Toolbox.Чаще других для справок используется команда help <имя функции>.

Пример .

>> help abs

ABS Absolute value.

ABS(X) is the absolute value of the elements of X. When

X is complex, ABS(X) is the complex modulus (magnitude) of

the elements of X.

See also SIGN, ANGLE, UNWRAP.

Overloaded methods

help iddata/abs.m

Большинство функций имеют несколько вариантов синтаксиса. Имена функций в сообщениях справочной системы выводятся прописными символами, но при вводе имен функций должны использоваться только строчные символы . Вывод справки по необходимой функции сопровождается списком связанных (сопутствующих) функций. Для получения более подробных сведений по нужной функции с примерами вычислений служит командаdoc <имя функции>.

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

Из меню Help c помощью команды Demos можно получить доступ к демонстрационным примерам MATLAB. Эти примеры очень разнообразны и полезны для целей обучения и создания приложений в MATLAB. Доступ к демонстрационным примерам можно получить также с помощью команды demo из командной строки.

Доступ к справочным сведениям в Internet: >>webhttp://www . mathworks . com - загружаетWEB- сайт фирмыMathWorksInc. - производителяMATLAB.

1. Простые вычисления

MATLAB имеет следующие базовые арифметические операции:

Сложение (a+b , 15+23),

Вычитание (a-b , 17-3),

Умножение (a*b , 0.18*6.12),

Деление (a/b , 92.4/15),

Возведение в степень (a ^ b , 7.4^4).

Примеры

Name Size Bytes Class Attributes

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

Для очистки рабочего пространства Workspace, т.е. удаления из него переменных, можно применить команду clear. Команда clc используется для очистки рабочего окна Command Window без очистки рабочей области.

Пакет MATLAB поддерживает также математические функции общего назначения, такие как извлечение квадратного корня sqrt (x ), вычисление прямых и обратных тригонометрических функций, экспоненциальных функций и др. Перечень всех этих функций с возможностью перехода к любой из них можно получить, введя в командную строку help elfun . Все элементарные функцииMATLABявляются функциями, аргументами которых могут быть массивы, т.е. в пакете реализованавекторизация вычислений.

Пример

>> v1 = [ 2 4 sqrt(10)]

2.0000 4.0000 3.1623

0.4161 -0.6536 -0.9998

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

Пример

>> 7*3+5-12/4

>> 7*(3+5-12/4)

Относительная точность выполнения арифметических операций MATLABсоставляет около 16 десятичных цифр в диапазоне чисел от 10 -308 до 10 308 . По умолчанию вMATLABустанавливается формат выводаshort , позволяющий вывести не более 5 значащих цифр числа. Такой формат вывода не всегда достаточен.

Команды установки форматов вывода

>> format short- короткое представление в фиксированном формате (5 знаков),

>> format short e– устанавливает формат научной (экспоненциальной) нотации с 5 десятичными разрядами,

>> format long– формат длинного представления с фиксированной точкой с 15 десятичными разрядами,

>> format long e– формат научной нотации с 15 десятичными разрядами,

>> format bank- денежный формат вывода с двумя десятичными разрядами справа от десятичной точки,

>> format rat- формат вывода в виде рациональной дроби.

Формат вывода может быть также установлен командой меню Preferences .

Следует иметь ввиду, что при вводе чисел в экспоненциальной форме, например, 15.8e-5промежуточные пробелы не допускаются.

Имена переменных MATLABдолжны начинаться с буквы, максимальная длина имени - 31 символ. Имена не должны совпадать с именами функций и процедур и системных переменных. Имена чувствительны к регистру, например,var отличается отVar .

Для создания переменной используется операция присваивания

<имя переменной> = <выражение>;

При этом оператор “;” подавляет эхо-вывод результатов вычисления (присваивания) на экран.

Все объявленные переменные сохраняются в рабочем пространстве (Workspace ) текущей сессииMATLABи доступны для вычислений в данной сессии, кроме случаев, когда переменные специально удаляются из рабочего пространства командойclear .

Примеры

Символьная переменная >> string="hello"

Действительные скаляры (числа)

>> y=5.2*x+15

Для сохранения переменных в файле текущего каталога (по умолчанию папка work) можно использовать команду save

>> save myfile x y

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

Переменные можно удалять из рабочего пространства (Workspace )MATLABкомандойclear

>> clear x y

Undefined function or variable "x".

Undefined function or variable "y".

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

>> load myfile

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

>> sqrt(-3)

Для обозначения мнимой единицы в MATLABзарезервированы переменныеi иj :

3.0000 + 4.0000i

>>y= 2*(1+4*j)

2.0000 + 8.0000i

Специальные функции вычислений с комплексным аргументом:

>> abs(x)% получение модуля числа

>> angle(x)% аргумент (фаза) числа в радианах

>> conj(x)% комплексно сопряженное число

>> imag(x)% мнимая часть числа

>> real(x)% действительная часть числа

Имена предопределенных системных переменных MATLAB нельзя использовать в качестве имен пользовательских переменных. Основные из этих имен:

>> ans имя переменной по умолчанию для результата вычислений.

>> eps переменная машинной точности вычислений, имеет порядок 10 -16 .

>> exit завершение (окончание) работы MATLAB.

>> i или j мнимая единица, т.е. .

>> pi число π.

>> Inf обозначение бесконечности.

>> NaN не числовой результат.

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

>> clear x,y команда удаляет переменные x и y.

>> what вывод списка файлов с расширениями ‘.m’, ‘.mat’, ‘.mex’ из текущего каталога.

>> who отображает переменные текущего рабочего пространства.

>> whos выводит информацию о текущих переменных.

>> dir выводит список файлов текущего каталога.

>> save сохраняет все текущие переменные в файле MATLAB.mat в текущем каталоге.

>> load загружает переменные из MATLAB.mat в текущий сеанс работы.

>> diary сохраняет текст (команды) и результаты вычислений текущего сеанса работы (дневник сессии) в файле с именем diary.

>> diary filename сохраняет текущий сеанс в файле с именем filename.

>> diary off приостанавливает запись в файл.

>> diary on включает запись сессии в файл.

с граничными условиями y (t 0 , t end , p ) = y , где t end , t 0 начальные и конечные точки интервалов. Параметр t (независимая переменная) необязательно означает время, хотя чаще всего решение ДУ ищется во временной области. Система ДУ в форме Коши записывается аналогично (1.1), но под y в этом случае подразумевается вектор-столбец зависимых переменных. Вектор p задает начальные условия.

Для решения ДУ второго и высшего порядка их нужно свести к системе ДУ первого порядка.

Возможны ДУ, не разрешенные относительно производной:

F (t , y , dy /dt ) = 0.

(1.2)

Уравнения (1.2) аналитически к форме (1.1) обычно привести не удается. Однако численное решение особых трудностей не вызывает достаточно для определения f (y , t ) решить (1.2) численно относительно производной при заданных y и t .

Решатели ОДУ

Для решения систем ОДУ в MATLAB реализованы различные численные методы. Их реализации названы решателями ОДУ.

В этом разделе обобщенное название solver (решатель) означает один из возможных численных методов решения ОДУ: ode45, ode23, ode113, ode15s, ode23s, ode23t, ode23tb, bvp4c или pdepe.

Решатели реализуют следующие методы решения систем ДУ:

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

Ode23 одношаговые явные методы Рунге-Кутта 2-го и 4-го порядков в модификации Богацки и Шампина. При умеренной жесткости системы ОДУ и низких требованиях к точности этот метод может дать выигрыш в скорости решения.

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

Ode15s многошаговый метод переменного порядка (от 1 до 5, по умолчанию 5), использующий формулы численного «дифференцирования назад». Это адаптивный метод, его стоит применять, если решатель ode45 не обеспечивает решения и система ДУ жесткая.

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

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

Ode23tb неявный метод Рунге Кутта в начале решения и метод, использующий формулы «дифференцирования назад» 2-го порядка в последующем. Несмотря на сравнительно низкую точность, этот метод может оказаться более эффективным, чем ode15s.

Bvp4c служит для проблемы граничных значений систем ДУ вида y ′ = f (t , y ), F (y (a ), y (b ), p ) = 0 (полная форма системы уравнений Коши). Решаемые им задачи называют двухточечными краевыми задачами, поскольку решение ищется при задании граничных условий как в начале, так и в конце интервала решения.

Все решатели могут решать системы уравнений явного вида y ′ = F (t , y ), причем для решения жестких систем уравнений рекомендуется использовать только специальные решатели ode15s, ode23s, ode23t, ode23tb.

Использование решателей систем ОДУ

tspan вектор, определяющий интервал интегрирования [t 0 t final ]. Для получения решений в конкретные моменты времени t 0 , t 1 , …, t final (расположенные в порядке уменьшения или увеличения) нужно использовать tspan = [ t 0 t 1 … t final ];

y 0 вектор начальных условий;

Options аргумент, создаваемый функцией odeset (еще одна функция odeget или bvpget (только для bvp4c) позволяет вывести параметры, установленные по умолчанию или с помощью функции odeset/bvpset);

p 1, p 2,… произвольные параметры, передаваемые в функцию F ;

T , Y матрица решений Y , где каждая строка соответствует времени, возвращенном в векторе-столбце T .

Перейдем к описанию синтаксиса функций для решения систем ДУ (под именем solver подразумевается любая из представленных выше функций).

[T ,Y ]=solver(@F ,tspan ,y 0) интегрирует систему ДУ вида y ′ = F (t , y ) на интервале tspan с начальными условиями y 0 . @F дескриптор ОДУ-функции (можно также задавать функцию в виде "F "). Каждая строка в массиве решений Y соответствует значению времени, возвращаемому в векторе-столбце T .

[T ,Y ]=solver(@F ,tspan ,y 0 ,options) дает решение, подобное описанному выше, но с параметрами, определяемыми значениями аргумента options, созданного функцией odeset. Обычно используемые параметры включают допустимое значение относительной погрешности RelTol (по умолчанию 1e3) и вектор допустимых значений абсолютной погрешности AbsTol (все компоненты по умолчанию равны 1e6).

[T ,Y ]=solver(@F ,tspan ,y 0 ,options,p 1 ,p 2 …) дает решение, подобное описанному выше, передавая дополнительные параметры p 1 , p 2 , … в m -файл F всякий раз, когда он вызывается. Используйте options=, если никакие параметры не задаются.

Решение ОДУ первого порядка

ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТЫ

· титульный лист;

· исходные данные варианта;

· решение задачи;

· результаты решения задачи.

Пример

Найти решение дифференциального уравнения на отрезке , для которого у (1,7) = 5,3.

Создаем в Command Window функция пользователя

g=@(x,y);

В синтаксисе функции @(x,y) x независимая переменная, y зависимая переменная, x -cos(y /pi ) правая часть ДУ.

Процесс решения осуществляется обращением в Command Window к решателю (солверу) следующим оператором:

Ode23(g,,);

Построение графика с сеткой осуществляется следующими операторами:

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

Рис. 1.2.1. Визуализация численного решения

ЗАДАНИЕ

1. Найдите решения ДУ первого порядка , удовлетворяющего начальным условиям у(х 0 ) = у 0 на промежутке [a , b ].

2. Построить графики функции.

Варианты заданий .

№ варианта		у(х 0 )=у 0	[a , b ]
		y 0 (1,8)=2,6
		y 0 (0,6)=0,8
		y 0 (2,1)=2,5
		y 0 (0,5)=0,6
		y 0 (1,4)=2,2
		y 0 (1,7)=5,3
		y 0 (1,4)=2,5
		y 0 (1,6)=4,6
		y 0 (1,8)=2,6
		y 0 (1,7)=5,3
		y 0 (0,4)=0,8
		y 0 (1,2)=1,4

Лабораторная работа № 2

Решение систем ОДУ

ЦЕЛЬ РАБОТЫ

Сформировать у студентов представления о применении систем ДУ в различных областях; привить умения решать задачу Коши для систем ДУ.

ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТЫ

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

2. Оформите отчет по лабораторной работе, который должен содержать:

· титульный лист;

· исходные данные варианта;

· решение задачи;

· результаты решения задачи.

Пример

Решить систему

с использованием решателя ode23().

Решение:

1. Создать в редакторе m-файл функции вычисления правых частей ДУ.

Пусть имя в редакторе файла sisdu.m, тогда функция может иметь следующий вид:

function z=sisdu(t,y)

z1=-3*y(2)+cos(t)-exp(t);

z2=4*y(2)-cos(t)+2*exp(t);

>> t0=0;tf=5;y0=[-3/17,4/17];

>> =ode23("sisdu",,y0);

>> plot(t,y)

Рис. 1.3.1. Визуализация численного решения, полученного с помощью функции ode23.

1. Что значит решить задачу Коши для системы ДУ?

2. Какие существуют методы решения систем ДУ?

ЗАДАНИЕ

1. Найдите решение системы ДУ

удовлетворяющее начальным условиям на промежутке ;

2. Построить графики функций.

Для примера приводится функция решения 8-го варианта:

function z=ssisdu(t,y)

% вариант 8

z1=-a*y(1)+a*y(2);

z2=a*y(1)-(a-m)*y(2)+2*m*y(3);

z3=a*y(2)-(a-m)*y(3)+3*m*y(4);

z4=a*y(3)-3*m*y(4);

>> =ode23("ssisdu",,);

>> plot(t,100*y)

Рис. 1.3.2. Визуализация численного решения, полученного с помощью функции ode23.

Варианты заданий .

№ варианта	Задания
a	m
	0,1	1,2
	0,2	1,5
	0,3	1,7
	0,4	1,9
	0,5
	0,6	1,9
	0,7	2,3
	0,8	2,7
	0,9
	0,1	1,5
	0,2	1,1
	0,3

Лабораторная работа № 3

1.4Решение ОДУ n -го порядка

ЦЕЛЬ РАБОТЫ

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

ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТЫ

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

2. Оформите отчет по лабораторной работе, который должен содержать:

· титульный лист;

· исходные данные варианта;

· решение задачи;

· результаты решения задачи.

Пример 1.

Решить ДУ второго порядка при данных начальных условиях .

Решение:

Сначала приведем ДУ к системе:

1. Создать m-файл функции вычисления правых частей ДУ.

Пусть имя файла sisdu_3.m, тогда функция может иметь следующий вид:

function z=sisdu_3(x,y)

z2=6*x*exp(x)+2*y(2)+y(1);

2. Выполнить следующие действия:

>> x0=0;xf=10;y0=;

>> =ode23("sisdu_3",,y0);

>> plot(x,y(:,1))

Рис. 1.4.1. Визуализация численного решения, полученного с помощью функции ode23.

ПРИМЕРНЫЕ ВОПРОСЫ НА ЗАЩИТЕ РАБОТЫ

1. Что значит решить задачу Коши для ДУ высших порядков?

2. Как привести ДУ m -го порядка к системе ДУ?

ЗАДАНИЕ

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

2. Построить графики функций.

Варианты заданий .

№ варианта	Задания
Уравнения	Начальные условия

Лабораторная работа № 4 – 5

Динамические системы (ДС)

ЦЕЛЬ РАБОТЫ

Знакомство студентов с основными понятиями ДС, их классификация, фазовое пространство ДС, кинематическая интерпретация системы ДУ, эволюция ДС. Уравнение движения маятника. Динамика осциллятора Ван дер Поля.

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

Состояние ДС описывают набором переменных, выбираемых из соображений естественности их интерпретации, простоты описания, симметрии и т.п. Множество состояний ДС образует фазовое пространство, каждому состоянию отвечает точка в нём, а эволюция изображается (фазовыми) траекториями. Чтобы определить близость состояний, в фазовом пространстве ДС вводят понятие расстояния. Совокупность состояний в фиксированный момент времени характеризуется фазовым объёмом.

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

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

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

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

Работа с системой носит диалоговый характер и происходит по правилу «задал вопрос - получил ответ». Пользователь набирает на клавиатуре вычисляемое выражение, редактирует его (если нужно) в командной строке и завершает ввод нажатием клавиши ENTER. В качестве примера на рисунке показаны простейшие и вполне очевидные вычисления.

Даже из таких простых примеров можно сделать некоторые поучительные выводы:

* для указания ввода исходных данных используется символ >>;

* данные вводятся с помощью простейшего строчного редактора;

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

* если не указана переменная для значения результата вычислений, то MATLAB назначает такую переменную с именем ans;

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

* встроенные функции (например, sin) записываются строчными буквами, и их аргументы указываются в круглых скобках;

* результат вычислений выводится в строках вывода (без знака >>);

* диалог происходит в стиле «задал вопрос - получил ответ».

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

В большинстве математических систем вычисление sin(V) или exp(V), где V - вектор, сопровождалось бы выдачей ошибки, поскольку функции sin и exp должны иметь аргумент в виде скалярной величины. Однако MATLAB - матричная система, а вектор является разновидностью матрицы с размером 1Чn или nЧ1. Поэтому в нашем случае результат вычислений будет вектором того же размера, что и аргумент V, но элементы возвращаемого вектора будут синусами или экспонентами от элементов вектора V.

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

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

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

Работа с редактором массивов

To get started, select "MATLAB Help" from the Help menu.

>> type sin

sin is a built-in function.

>> help sin

SIN(X) is the sine of the elements of X.

Overloaded methods

>> V=

0.8415 0.9093 0.1411 -0.7568

Error using ==> ^

Matrix must be square.

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

Форма вывода и перенос строки в сессии

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

будет дан в виде:

Исключением является вывод векторов столбцов и матриц - тут будет сохранена более наглядная и присущая MATLAB по умолчанию форма вывода.

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

s = 1 - 1/2 + 1/3 - 1/4 + 1/5 - 1/6 + 1/7 ...

1/8 + 1/9 - 1/10 + 1/11 - 1/12;

Максимальное число символов в одной строке командного режима - 4096, а в m-файле - не ограничено, но со столь длинными строками работать неудобно. В ранних версиях в одной строке было не более 256 символов.

Запуск примеров применения MATLAB из командной строки

MATLAB имеет множество примеров применения, часть из которых можно запускать прямо из командной строки. Например, команда

запускает m-файл bench.m демонстрационного примера тестирования системы.

Вычисления и приближение данных в MATLAB

Решение систем линейных алгебраических уравнений и спектральных задач

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

Обзор возможностей MATLAB для решения систем линейных алгебраических уравнений

Среда MATLAB изначально разрабатывалась для работы с матрицами (MATLAB является сокращением от Matrix Laboratory), поэтому арсенал средств MATLAB для решения систем линейных алгебраических уравнений достаточно богат и включает в себя:

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

Решение систем при помощи знака обратной косой черты

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

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

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

X = A\f x = 1 1 1

Вместо знака обратной косой черты можно было вызвать функцию mldivide

X = mldivide(A, f)

Результат будет тем же самым

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

F - A*x ans = 0 0 0

Очень важно не перепутать местами A и f, т.к. при выполнении операции

никакой ошибки не будет, а выведется

X = 0.2707 0.3439 0.3854

что не имеет ничего общего с решением рассматриваемой системы. Разберем, почему так происходит, записав систему с «матрицей» f и «правой частью» A. Неизвестная при этом будет только одна, т.к. f размера 3 на 1:

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

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

(4-7x 1) 2 +(3-11x 1) 2 +(2-12x 1) 2

Несложно убедиться, что минимум доставляет как раз . Воспользуемся, например, средствами Symbolic Math Toolbox

Syms x R = (7*x-4)^2 + (11*x-3)^2 + (12*x-2)^2 dRdx = diff(R) x = solve(dRdx) x = 85/314 double(x) ans = 0.2707

Итак, полученный при записи x = f\A вектор содержит решение трех переопределенных систем, правая часть каждой из которых является соответствующим столбцом матрицы A.

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

k	1	2	3	4	5
x k	1.0	1.5	2.0	2.5	3.0
y k	2.99	2.81	2.89	3.03	3.21

Эта задача сводится к переопределенной системе относительно a и b , если мы потребуем равенства y(x k )=y k для k =1,2,...,5:

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

X = (1:0.5:3)" A = ; f = ; и решаем ее c = A\f получаем c = 2.9903 2.0100 т.е. и. График данных и приближения подтверждает правильность полученных результатов figure fun = @(x)c(1)/x + c(2)*log(x); fplot(fun, ) hold on plot(x, y, "or")

Если в системе больше уравнений, чем неизвестных, то такая система называется недоопределенной (underdetermined system) и решается в MATLAB так же при помощи знака обратной косой черты. У нее будет бесконечно много решений и находится решение, содержащее как можно больше нулевых значений. Рассмотрим такой пример

A = ; f = ; x = A\f x = 1.3333 1.0000 0 1.6667

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

B/A эквивалентно (A"\B")"
где апостроф означает транспонирование. Вместо знака обычной косой черты можно использовать также функцию mrdivide.

Пока мы не будем углубляться в алгоритмы решения систем, заложенные в операции обратной косой черты. Этому посвящен раздел . В частности, для решения систем с квадратной матрицей общего вида используется LU-разложение. Кроме того, проверяется обусловленность матрицы. Классическим примером плохообусловленной матрицы является матрица Гильберта, элементы которой определяются по формуле . Для создания матрицы Гильберта заданного размера в MATLAB служит функция hilb. Решим, например, систему 13-го порядка, правая часть которой такова, что ее решением должны быть все единицы (k -ый элемент правой части есть сумма элементов k -ой строки матрицы). A = hilb(13); f = sum(A, 2); x = A\f

В результате получаем сообщение о плохой обусловленности

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.339949e-018.

(RCOND есть оценка для величины, обратной числу обусловленности) и неверное решение

X = 1.0000 1.0000 1.0004 0.9929 1.0636 0.6552 2.2023 -1.7860 5.3352 -3.4773 3.9431 -0.1145 1.1851

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

Матричные разложения Холецкого, LU и QR

В MATLAB есть функции для следующих разложений:

• разложения Холецкого - функция chol;
• LU-разложения - функция lu;
• QR-разложения - функция qr.

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

Ax=f (k) , k =1, 2,...,N .

Предположим, что матрица A представлена в виде произведения двух матриц, причем системы с матрицами A=BC причем системы с матрицами B и C решаются значительно быстрее, чем с матрицей A . Тогда решения исходных систем получаются следующим образом. Поскольку Ax=f (k) и A=BC , то BCx=f (k) и решая последовательно By=f (k) и Cx=y находим решение k -ой системы. Для заполненных матриц разложение выполняется за O (n 3) операций, где n - размер матрицы (т.е. за то же время, что и решение системы), а решение системы с каждым множителем разложения за O (n 2) операций. Поэтому решение N систем c с предварительным разложением занимает O (n 3)+O (n 2)N . Решение же каждой системы без разложения матрицы потребовало бы O (n 3)N арифметических действий.

Разложение Холецкого для заданной матрицы A заключается в нахождении такой верхней треугольной матрицы R с положительными диагональными элементами, что A=R T R . Известно, что если матрица A симметрична и положительно определена (т.е. для любого вектораx верно: x T Ax0, или, что тоже самое, все собственные числа матрицы положительны), то разложение Холецкого существует и при том единственно.

Рассмотрим в качестве примера разложение Холецкого матрицы

A = ; R = chol(A) R = 2.0000 0.5000 0.5000 0 1.9365 0.3873 0 0 1.8974

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

A - R"*R ans = 1.0e-015 * 0 0 0 0 0 0 0 0 0.4441

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

A = ; R = chol(A) ??? Error using ==> chol Matrix must be positive definite.

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

A = ; R = chol(A) R = 2.0000 0.5000 1.0000 0 1.9365 0.2582 0 0 1.7127 A - R"*R ans = 0 0 0 0 0 0 -1 0 0

Для произвольной квадратной матрицы можно выполнить LU -разложение, т.е. найти нижнюю треугольную матрицу L и верхнюю треугольную матрицу U такие, что A=LU . Более строгое утверждение формулируется следующим образом.

Если Ak - главный минор квадратной матрицы A размера n , составленный из первых k строк и столбцов (т.е. А(1:k, 1:k)) и det(A k )0 для k =1,2,..., n -1, то существует единственная нижняя треугольная матрица L , диагональ которой состоит из единиц, и единственная верхняя треугольная матрица U такие, что A=LU .

При вычислении LU LU-разложения может понадобится перестановка строк для обеспечения численной устойчивости процесса разложения, поэтому функция lu может вернуть матрицу L , которая является нижней треугольной с точностью до перестановки строк, например:

A = ; = lu(A) L = 0 1.0000 0 1.0000 0 0 0.5000 0.5000 1.0000 U = 2 3 1 0 1 1 0 0 4

Произведение LU где L - треугольная с точностью до перестановки строк, и будет равно матрице A (вообще говоря, с учетом погрешности при операциях с вещественными числами).

A-L*U ans = 0 0 0 0 0 0 0 0 0

Если для матрицы A сделано разложение Холецкого или LU -разложение, то решение системы с матрицей A , как уже было сказано выше, сводится к решению двух систем с треугольными матрицами. Это решение может быть выполнено при помощи операции обратной косой черты, поскольку заложенный в ней алгоритм определяет треугольные матрицы и применяет эффективный способ решения системы с ними, требующий O (n 2) арифметических действий. Рассмотрим пример решения системы

с предварительным LU -разложением матрицы.

A = ; f = ;

Выполняем LU -разложение

Lu(A);

и решаем последовательно две системы с треугольными матрицами, сначала с L, затем с U

Y = L\f; x = U\y x = 1 1 1

Решение двух систем можно записать одним выражением

X = U\(L\f)

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

X = U\L\f

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

Рассмотрим последнее матричное разложение - QR -разложение, выполняемое функцией qr. QR -разложение может быть проделано для прямоугольных матриц, именно если A матрица размера m на n , то существует ортогональная матрица Q размера m на m (т.е. такая, что Q -1 =Q T) и матрица R размера m на n со всеми нулевыми элементами под главной диагональю, что A=QR .

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

Qr(A) Q = -0.7428 0.6000 -0.2971 -0.5571 -0.8000 -0.2228 -0.3714 0.0000 0.9285 R = -5.3852 -5.3852 -5.0138 0 -5.0000 0.4000 0 0 6.6108

После того, как сделано QR -разложение, решение системы Ax=f может быть найдено как решение системы с треугольной матрицей R , т.к. поскольку QR=f , то Rx=Q T f :

>> x = R\(Q"*f) x = 1.0000 1.0000 1.0000

Более подробно про матричные разложения, соответствующие функции MATLAB chol, lu и qr и их аргументы написано в разделе Матричные разложения.

Алгоритм решения системы линейных уравнений при помощи знака обратной косой черты

При решении в MATLAB системы линейных алгебраических уравнений при помощи знака обратной косой черты
x = A\b
работает алгоритм, который в зависимости от типа матрицы решает систему при помощи наиболее подходящего метода, реализованного в одной из процедур пакета LAPACK или UMFPACK . Здесь b может быть и матрицей, число строк которой совпадает с числом строк матрицы A. Тогда возвращается матрица x, каждый i-ый ее столбец содержит решение системы Ax (i) = b (i) , i=1,2,...,k, где b (i) - i-ый столбец матрицы b = [ b (1) | b (2) ...| b (k) ].

Алгоритм, реализованный в операции \ состоит из следующих шагов:

• 1. В тривиальном случае, если A разреженная и диагональная (разреженным матрицам в MATLAB посвящен раздел Разреженные матрицы), решение находится по простой формуле x k = b k /a kk , где k=1,2,...n.

2. Если A - квадратная, разреженная и ленточная матрица, то используется солвер для ленточных матриц. При этом могут быть два варианта:

a. Если A - трехдиагональная матрица и b - один столбец из вещественных чисел, тогда система решается гауссовым исключением (его операции совершаются только с элементами на диагоналях). Если в процессе исключения для сохранения устойчивости потребуется сделать перестановки строк, или если матрица A не является трехдиагональной, то работает следующий пункт.
b. Если A - ленточная с плотностью ненулевых элементов больше параметра bandden , по умолчанию равного 0.5, или не выполняются условия предыдущего пункта, тогда, в зависимости от типа A и b (double или single ) вызываются следующие процедуры библиотеки LAPACK:

A и b вещественные типа double - процедуры DGBTRF, DGBTRS
A и b комплексные типа double - процедуры ZGBTRF, ZGBTRS
A или b вещественные типа single - процедуры SGBTRF, SGBTRS
A или b комплексные типа single - процедуры CGBTRF, CGBTRS

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

nz = 1 (на диаг.№-2) + 6 (на диаг.№-1) + 7 (на диаг.№0) + 6 (на диаг.№1) + 1 (на диаг.№2) + 1 (на диаг.№3) = 22,

а число всех элементов в ленте

band = 5 (на диаг.№-2) + 6 (на диаг.№-1) + 7 (на диаг.№0) + 6 (на диаг.№1) + 5 (на диаг.№2) + 4 (на диаг.№3) = 33

и плотность ненулевых элементов будет 2/3. Поскольку 2/3 > 0.5, то будет применен солвер для ленточных матриц
Параметр bandden , как и другие параметры, управляющие алгоритмами для разреженных матриц в MATLAB, устанавливается при помощи функции spparms.



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

• 3. Если A - верхняя или нижняя треугольная матрица, то применяется метод обратной подстановки, или, соответственно, метод прямой подстановки, в котором из последнего (или первого уравнения) находится компонента решения и далее найденные компоненты подставляются в следующие уравнения для нахождения очередных компонент решения системы уравнений.
• 4. Если A - может быть приведена перестановками к треугольной форме, то выполняется соответствующее приведение, а далее система уравнений решается как и в п.3.

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

a. Если А не является разреженной, тогда делается попытка выполнить разложение Холецкого A = R T R с последующим решением систем с треугольными матрицами R T и R: R T y = b и Rx = y. При этом вызываются следующие процедуры библиотеки LAPACK:

A и b вещественные и типа double - DLANGE, DPOTRF, DPOTRS, DPOCON
A и b комплексные и типа double - ZLANGE, ZPOTRF, ZPOTRS, ZPOCON
A и b вещественные и типа single - SLANGE, SPOTRF, SPOTRS, SDPOCON
A и b комплексные и типа single - CLANGE, CPOTRF, CPOTRS, CPOCON

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

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

p = symmmd(A) - в вектор p записаны перестановки

R = chol(A(p, p));

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

и вторая, с множителем разложения и с занесением компонент решения на соответствующие позиции в векторе x
x(p, :) = R \ y

• 6. Если A - хессенбергова матрица, то она приводится к верхней треугольной (с соответствующей модификацией правой части) и затем система решается подстановками

7. Если A - квадратная матрица, не удовлетворяющая условиям пунктов 1-6, то в зависимости от того, разрежена она или нет, выполняются следующие действия

a. Если A не является разреженной матрицей, то при помощи исключения Гаусса с выбором ведущего элемента (для обеспечения устойчивости разложения) выполняется LU-разложение матрицы A = LU, где

U - верхняя треугольная матрица
L - матрица, перестановками строк сводящаяся к треугольной

и решение системы Ax = b находится последовательным решением систем с треугольными матрицами Ly = b, Ux = y.
Для выполнения LU-разложения вызываются следующие процедуры библиотеки LAPACK:

A и b вещественные и типа double - DLANGE, DGESV, DGECON
A и b комплексные и типа double - ZLANGE, ZGESV, ZGECON
A и b вещественные и типа single - SLANGE, SGESV, SGECON
A и b комплексные и типа single - CLANGE, CGESV, CGECON

b. Если A - разреженная матрица, то производится перестановка столбцов с целью уменьшения заполнения множителей L и U в процессе нахождения разложения. Далее перестановкой строк в процессе LU-разложения обеспечивается вычислительная устойчивость, после которого снова решаются системы с треугольными матрицами. Для выполнения LU-разложения разреженной матрицы используются процедуры библиотеки UMFPACK

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

a. Если A не является разреженной матрицей, то выполняется QR-разложение AP = QR, где P - матрица перестановки столбцов, Q - ортогональная матрица (Q T Q = I), а R - верхняя треугольная. Если A размера m на n, то Q размера m на m, а R размера m на n. Далее решение системы находится так:

x = P*(R \ (Q" * b))

поскольку из Ax = b и AP = QR следует: QRx = bP и Rx = Q T bP.

b. В случае разреженной и прямоугольной матрицы A нет смысла вычислять матрицу Q, поскольку она скорее всего окажется заполненной. Поэтому в ходе алгоритма QR-разложения вычисляется c = Q T b (т.е. модифицированная правая часть) и решается система Rx = c с треугольной матрицей для получения решения исходной системы Ax = b.

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

Функция linsolve для решения систем линейных уравнений

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

Функция linsolve, вызванная в самом простом варианте с двумя входными аргументами и одним выходным аргументом
x = linsolve(A, b)
решает систему Ax = b одним из способов, в зависимости от того, квадратная матрица, или нет.

• Если A - квадратная матрица, то предварительно вычисляется ее LU-разложение и затем решаются две системы с треугольными матрицами L и U.
• Если A - прямоугольная матрица, то предварительно вычисляется ее QR-разложение и затем решается система с треугольными матрицами.

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

A = ; b = ; x = linsolve(A,b) Warning: Rank deficient, rank = 1, tol = 4.4686e-015. x = 0 0 2.0000

Для подавления вывода данного сообщения в командное окно следует вызывать функцию linsolve с двумя выходными аргументами
= linsolve(A, b)
тогда в x запишется полученное решение, а в r - или ранг матрицы (если A является прямоугольной матрицей), или обратная величина к оценке для ее числа обусловленности (см. раздел ), например

A = ; b = ; = linsolve(A,b) x = -1.0000e+000 9.9999e-001 3.3307e-006 r = 6.9444e-013

Основные преимущества функции linsolve над операцией \ проявляются при указании типа матрицы системы уравнений. Для этого перед вызовом функции linsolve надо заполнить структуру с информацией о матрице со следующими полями

• SYM - симметричная;
• LT - нижняя треугольная;
• UT - верхняя треугольная;
• UHESS - хессенбергова;
• POSDEF - симметричная, положительно определенная;
• RECT - прямоугольная;
• TRANSA - надо ли решать систему с матрицей, транспонированной к заданной.

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

LT	UT	UHESS	SYM	POSDEF	RECT	TRANSA
true	false	false	false	false	true/false	true/false
false	true	false	false	false	true/false	true/false
false	false	true	false	false	false	true/false
false	false	false	true	true	false	true/false
false	false	false	false	false	true/false	true/false

Если матрица системы положительно определена, то это обязательно учесть при решении, поскольку для положительно определенных матриц решение основано на разложении Холецкого, которое требует меньше операций, чем LU-разложение, используемое при решении систем с квадратными матрицами общего вида. В этом несложно убедиться при помощи следующего примера, в котором создается симметричная положительно определенная матрица (матрица из случайных чисел складывается с транспонированной к ней и на диагональ добавляются большие числа) и система уравнений с этой матрицей сначала решается как система с матрицей общего вида (opts.SYM = false и opts.POSDEF = false), а затем - как с симметричной и положительно определенной матрицей (opts.SYM = true и opts.POSDEF = true).

% задаем все поля структуры ops, кроме SYM и POSDEF opts.TRANSA = false; opts.UHESS = false; opts.RECT = false; opts.UT = false; opts.LT = false; % создаем вектор размеров матриц N = 2.^(8:12); % создаем пустые массивы для записи времен решения TimeGeneral = ; TimePosSym = ; % в цикле создаем матрицы и сравниваем времена решения for n = N % создаем симметричную и положительно определенную матрицу % и вектор правой части A = rand(n); A = A + A" + 2*n*eye(n); b = sum(A, 2); % решаем систему как систему с матрицей общего вида opts.SYM = false; opts.POSDEF = false; Tstart = cputime; x = linsolve(A,b, opts); Tend = cputime; TimeGeneral = ; % решаем систему как систему с симметр и полож. опред матрицей opts.SYM = true; opts.POSDEF = true; Tstart = cputime; x = linsolve(A,b, opts); Tend = cputime; TimePosSym = ; end % выводим графики временных затрат figure loglog(N, TimeGeneral, N, TimePosSym) legend("TimeGeneral", "TimePosSym")

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

Разумеется, если матрица треугольная, то очень важно это указать, поскольку решение системы с треугольной матрицей выполняется за O(n 2) операций, а если к системе с треугольной матрицей применить алгоритм решения, предназначенный для матрицы общего вида, то понадобится порядка O(n 3) операций.

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

A = ; b = ; задать true для поля UT (а все остальные установить в false) opts.UT = true; opts.TRANSA = false; opts.LT = false; opts.UHESS = false; opts.SYM = false; opts.POSDEF = false; opts.RECT = false; то при решении системы функция linsolve будет рассматривать матрицу как верхнюю треугольную и выберет нужные ей элементы из верхнего треугольника A x = linsolve(A,b, opts) При этом получится решение x = 1 1 1 соответствующее матрице A = ;

Влияние обусловленности матрицы на точность решения системы с ней

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

причем предполагается, что каждую из систем мы решаем точно, а основной возникающий вопрос - насколько будет отличаться решение x системы (1) от решения системы (2) с возмущением в правой части δb. Это довольно важный вопрос, поскольку элементы правой части могут быть получены как результат некоторых измерений, соответственно содержащий погрешность δb k в каждой компоненте b k . Хотелось бы, чтобы при измерении одной и той же величины (каждый раз со своими, небольшими погрешностями) соответствующие решения систем с мало отличающимися правыми частями также не очень сильно отличались друг от друга. К сожалению, так бывает не всегда. Причиной этому - характеристика матрицы, называемая ее обусловленностью . Об этом и пойдет дальше речь.

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

Сответственно

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

т.е. спектральная матричная норма равна квадратному корню из максимального собственного числа матрицы A T A. В MATLAB имеется функция norm для вычисления норм матриц и векторов, которая, в частности, умеет вычислять приведенные выше нормы. Почему мы выбрали именно такие нормы векторов и матриц, подробно объясняется в разделе , где приведено немного выкладок и определений. Это связано с оценкой, которой мы будем пользоваться для ошибки в решении системы (вывод этого неравенства также приведен в упомянутом разделе):

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

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

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

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

N = 6; H = hilb(n) H = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909

X = ones(n, 1); b = H*x b = 2.4500 1.5929 1.2179 0.9956 0.8456 0.7365

Видим, что ни в матрице, ни в векторе правой части нет ничего «подозрительного», все числа не сильно отличаются друг от друга. Теперь сформируем возмущенную правую часть b + δb, добавив в вектор b небольшие числа порядка 10 -5 и решим систему с возмущенной правой частью для получения вектора .

Delta_b = 1e-5*(1:n)"; x_tilda = H\(b + delta_b) x_tilda = 0.9978 1.0735 0.4288 2.6632 -1.0160 1.8593

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

Delta_x = x - x_tilda; LEFT = norm(delta_x)/norm(x) LEFT = 1.1475 RIGHT = norm(delta_b)/norm(b) RIGHT = 2.7231e-005

Итак, ошибка в решении порядка единицы, хотя изменения в правой части были порядка 10 -5 . Это полностью укладывается в приведенное выше неравенство для ошибки. Действительно, вычислим число обусловленности cond(H) при помощи функции MATLAB, которая называется cond и по умолчанию вычисляет для спектральной нормы матрицы

C = cond(H) c = 1.4951e+007

LEFT ≈ 1.1475 меньше

C* RIGHT ? 1.4951e+07 * 2.7231e-05 ≈ 407

и неравенство выполнено (даже с некоторым запасом).

При увеличении размера матрицы Гильберта ошибка в решении будет только расти (это несложно проверить, задавая n = 7, 8, …). Причем, при n = 12 выведется сообщение о том, что матрица плохообусловлена и решение может оказаться неверным

Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 2.409320e-017.

В качестве меры обусловленности здесь выбрана величина RCOND равная единице, деленной на оценку числа обусловленности (число обусловленности оценивается по достаточно быстрому алгоритму, который работает намного быстрее, чем cond, о чем более подробно написано в справке по функции rcond). Если значение RCOND мало, то матрица считается плохообусловленной.

Увеличение ошибки в решении при увеличении размера матрицы Гильберта происходит из-за того, что число обусловленности матрицы Гильберта растет очень быстро с ее размером. В этом несложно убедиться при помощи простого цикла и функции semilogy (масштаб по оси ординат - логарифмический):

N = 1:20; C = zeros(1, 20); for n = N H = hilb(n); C(n) = cond(H); end semilogy(N, C) grid on, title("cond(H)"), xlabel("n")

Кстати, поскольку функция cond находит число обусловленности при помощи численного метода (а именно сингулярного разложения для нахождения сингулярных чисел), то число обусловленности после n = 12 вычисляется уже неверно, на самом деле оно должно расти и дальше, в чем можно убедиться при помощи символьных вычислений в MATLAB и операций с заданным числом значащих цифр

N = 1:20; C = zeros(1, N); digits(60) for n = N H = vpa(sym(hilb(n))); % вычисление матрицы Гильберта с точностью до 60-го знака sigma = svd(H); % нахождение сингулярных чисел матрицы Гильберта % вычисление числа обусловленности матрицы Гильберта C(n) = max(double(sigma))/min(double(sigma)); end semilogy(N, C) grid on, title("cond(H)"), xlabel("n")

Рассмотрим теперь три важных момента, касающихся данного примера.

Первый из них - решение системы Hx = b (с вектором правой части, соответствующим решению из единиц) при помощи оператора обратной косой черты не даст точные единицы, погрешность будет уже в десятом знаке (хотя в MATLAB все вычисления по умолчанию выполняются с двойной точностью)

Format long H\b ans = 0.99999999999916 1.00000000002391 0.99999999983793 1.00000000042209 0.99999999953366 1.00000000018389

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

A = rand(n); cond(A) ans = 57.35245279907571

Затем сформируем правую часть, соответствующую точному решению из всех единиц

X = ones(n, 1); b = A*x;

и решим систему, получив в результате хорошую точность

>> A\b ans = 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000

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

>> det(hilb(6)) ans = 5.3673e-018

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

A = ;

Ее определитель очень мал, порядка 10 -121

Det(A) ans = 9.9970e-121

(слово «очень», конечно условно, но он меньше определителя матрицы Гильберта размера 6 на 6, при решении системы с которой у нас возникли проблемы). Однако, малость определителя никак не повлияет на ошибку в решении при возмущении правой части системы, что несложно показать, сформировав систему с известным решением, внеся возмущения в правую часть и решив систему:

X = ; b = A*x; delta_b = 1e-5*; x_tilda = A\(b+delta_b); delta_x = x - x_tilda; LEFT = norm(delta_x)/norm(x) RIGHT = norm(delta_b)/norm(b) LEFT = 2.1272e-005 RIGHT = 2.1179e-005

Итак, относительная ошибка в решении практически равна относительной ошибке в правой части, поскольку число обусловленности приведенной выше матрицы A немного больше 1, а именно:

C = cond(A) c = 1.0303

И, наконец, рассмотрим третий вопрос, касающийся достижения знака равенства в неравенстве для ошибки в решении

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

в котором относительное возмущение правой части равно 10 -5

RIGHT = norm(delta_b)/norm(b) RIGHT = 1.0000e-005

Относительная ошибка в решении системы оказывается 10 5

X = A\b; x_tilda = A\(b+delta_b); delta_x = x - x_tilda; LEFT = norm(delta_x)/norm(x) LEFT = 1.0000e+005

И происходит так потому, что
1) в этом примере число обусловленности матрицы A равно 10 10 ;

C = cond(A) c = 1.0000e+010

2) неравенство для ошибки в решении переходит в равенство.

Если установить формат long, то мы увидим, что LEFT, RIGHT и число обусловленности не в точности 10 5 , 10 -5 и 10 10 , соответственно, но это вызвано ошибками округлений. При решении в точной арифметике, равенство получилось бы точным, что можно показать аналитически (см. раздел ).

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

|| Ax || V = || b || V ⇒ || A || M || x || V ≥ || b || V

где || || M - некоторая матричная норма. Для того, чтобы это можно было сделать, матричная норма || || M и векторная норма || || V должны удовлетворять следующему неравенству

|| A || M || x || V ≥ || Ax || V

для любых матриц A и векторов x. Если это неравенство выполняется, то матричная норма || || M называется согласованной с векторной нормой || || V . Известно, например, что спектральная норма

(квадратный корень из максимального собственного числа матрицы A T) согласована с евклидовой векторной нормой

Именно поэтому в разделе мы проводили все эксперименты с привлечением данных норм.

Получив неравенство || A || M || x || V ≥ || b || V далее заметим, что из Ax = b и следует, что . Поскольку матрица A невырожденная, то отсюда следует, что δx = A -1 δb и || δx || V = || A -1 δb || V . Опять используем согласованность норм и получаем неравенство || δx || V ≤ || A -1 || M || δb || V . Далее в двух полученных неравенствах

|| A || M || x || V ≥ || b || V и || δx || V ≤ || A -1 || M || δb || V

меньшую величину одного из неравенств делим на большую величину другого неравенства, а большую, соответственно, на меньшую

и простым преобразованием окончательно получаем требуемое неравенство

где cond(A) = || A || M * || A -1 || M .

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

согласована с первой векторной нормой

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

согласована с бесконечной векторной нормой

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

• q = norm(A, 1) - максимальная столбцевая норма матрицы A
• q = norm(A, inf) - максимальная строчная норма матрицы A
• q = norm(a, 1) - первая векторная норма a
• q = norm(a, inf) - бесконечная векторная норма a

Число обусловленности матрицы cond(A) = || A || M * || A -1 || M относительно различных матричных норм может быть вычислено при помощи функции cond. Если cond вызывается с одним входным аргументом (матрицей), то вычисляется число обусловленности относительно спектральной матричной нормы. При вызове функции cond с дополнительным аргументом возвращаются числа обусловленности относительно указанной матричной нормы:

• с = cond(A, 1) - число обусловленности относительно максимальной столбцевой нормы матрицы
• с = cond(A, inf) - число обусловленности относительно максимальной строчной нормы матрицы

В качестве примера, демонстрирующего важность согласованности норы матрицы с нормой вектора в неравенстве для ошибки, приведем пример с матрицей A, вектором правой части b и вектором ошибок в правой части delta_b.

A = ; b = [ -5.7373057243726833e-001 -1.5400413072907607e-001 -5.3347697688693385e-001 -6.0209490373259589e-001]; delta_b = [-0.71685839091451e-5 0.54786619630709e-5 0.37746931527138e-5 0.20850322383081e-5];

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

RIGHT = norm(delta_b,1)/norm(b,1); c = cond(A); x = A\b; x_tilda = A\(b+delta_b); delta_x = x - x_tilda; LEFT = norm(delta_x,1)/norm(x,1); format short e disp()

Получаем для левой и правой частей неравенства следующие значения
9.9484e+004 9.9323e+004

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

Разберем теперь, почему число обусловленности матрицы не может быть меньше единицы. Число обусловленности матрицы A определяется как cond(A) = || A || M * || A -1 || M , где || A || M - какая-то матричная норма A. Матричная норма (т.е. правило, по которому каждой матрице сопоставляется число) не может быть произвольной, она должна удовлетворять следующим четырем аксиомам:

A1. || A || ≥ 0 для любой матрицы A и || A || = 0 тогда и только тогда, когда A - нулевая матрица.

А2. || αA || = | α | * || A || для любой матрицы A и числа α.

А3. || A + B || ≤ || A || + || B || для любых матриц A и B

А4. || AB || ≤ || A || * || B || для любых матриц A и B.

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

|| I || = || I*I || ≤ || I || 2 ⇒ || I || ≥ 1.

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

1 ≤ || I || = || AA -1 || ≤ || A || * || A -1 || = cond(A).

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

и построением соответствующих примеров в MATLAB. (Приводимое ниже рассуждение содержится, например, в книге Дж. Форсайт, К Молер. Численное решение систем линейных алгебраических уравнений. М: «Мир», 1969.)

Существенную роль в рассуждениях играет теорема о сингулярном разложении матрицы, согласно которой для любой вещественной матрицы размера на существуют две ортогональные матрицы U и V размера n на n (U T U=UU T и V T V = VV T) такие, что произведение D = U T AV является диагональной матрицей, причем можно выбрать U и V так, что

где μ 1 ≥ μ 2 ≥…≥μ r ≥ μ r+1 =…=μ n =0,

а r - ранг матрицы A. Числа μ k называют спектральными числами матрицы A. Для невырожденной матрицы A верно:

μ 1 ≥ μ 2 ≥ … ≥μ n ≥ 0.

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

|| Ux || = || x ||.

Поскольку умножение на ортогональную матрицу не изменяет спектральной нормы, то

следовательно для числа обусловленности матрицы верно равенство

У нас есть две системы: Ax = b (с точным решением x) и (с точным решением ). Очевидно, что ошибка δx является решением системы, правая часть которой является возмущением δb, т.е. системы Aδx = δb. Пусть D = U T AV есть сингулярное разложение матрицы A, тогда UDV T = A из-за того, что U и V - ортогональные матрицы. Далее

Ax = b ⇒ UDV T x = b.

Введем обозначения

x" = V T x, b" = U T b,

тогда следующие системы эквивалентны

Ax = b ⇔ Dx" = b"

Совершенно аналогично рассмотрим систему относительно ошибки

Aδx = δb ⇒ UDV T δx = δb

Вводим обозначения

δx" = V T δx, δb" = U T δb,

при которых следующие системы оказываются эквивалентны

Aδx = δb ⇔ Dδx" = δb"

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

Выберем теперь специальным образом правые части этих систем, именно, пусть

где β > 0, тогда соответствующее ей решение системы Dx" = b" будет

где μ 1 - максимальное сингулярное число матрицы A. Схожим образом поступим и с системой Dδx" = δ b" , именно, пусть

где γ > 0, тогда соответствующее ей решение системы Dδx" = δb" будет

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

β/μ 1 = || x" || = || V T x || = || x || и γ/μ n = || δx" || = || V T δx || = ||δx ||.

Совершенно аналогично получаем равенства:

β = || b" || = || U T b || = || b || и γ = || δb" || = || U T δb || = || δb ||.

а поскольку

то окончательно получаем

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

Сначала задаем n и получаем две ортогональные матрицы U и V, выполняя QR-разложение матриц n на n из случайных чисел:

N = 4; = qr(rand(n)); = qr(rand(n));

D = diag();

После этого конструируем матрицу A, используя диагональную матрицу D и ортогональные матрицы U и V

A = U*D*V";

Число обусловленности матрицы A будет совпадать с числом обусловленности матрицы D, которое в нашем примере равно 10 10

Beta = 1; gamma = 1e-5;

и строим вектора b" и δb"

B1 = "; db1 = ";

по которым находим вектора b и δb

X = A\b; x_tilda = A\(b+db);

вычисляем левую и правую части неравенства

dx = x - x_tilda; RIGHT = norm(db)/norm(b); LEFT = norm(dx)/norm(x);

и выводим их

Format short e disp()

Получаем равенство

1. Command Window (Окно команд).

Математические выражения пишутся в командной строке после знака приглашения « >> ». Например,

Для выполнения действия нажмем клавишу «Enter».

По умолчанию программа записывает результат в специальную переменную ans.

Для сохранения результата под другим именем используют имя переменной и знак присваивания « = », например

>> z = 1.25 /3.11

Редактировать в Command Window можно только текущую командную строку. Для того чтобы отредактировать ранее введенную команду, необходимо установить курсор в строку ввода и воспользоваться клавишами « » или« ».

Если команда заканчивается «;», то результат её действия не отображается в командной строке.

Командное окно можно закрыть кнопкой « », а кнопка « » служит для отделения окна от интерфейса системы. Вернуться к стандартной форме окна можно с помощью меню:

Главное Меню → Desktop → Desktop Layout → Default .

Очистить командное окно можно с помощью меню:

Главное Меню → Edit → Clear Command Window .

Главное меню системы MatLab.

Главное меню MatLab содержит следующие пункты:

File (Файл) – работа с файлами;

Edit (Правка) – редактирование;

View (Вид) – управление окнами;

Web – связь с фирмой – разработчиком через Интернет;

Window (Окно) – связь с окнами системы;

Help (Справка) – связь со справочной системойMatLab.

Панель инструментов системы MATLAB.

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

New file (Создать) – выводит окна редакторов файлов;

Open file (Открыть) – открывает окна загрузки файлов;

Cut (Вырезать) – вырезает выделенный фрагмент и помещает в буфер обмена;

Copy (Копировать) – копирует выделенный фрагмент в буфера обмена;

Paste (Вставить) – переносит выделенный фрагмент из буфера обмена в строку ввода;

Undo (Отменить) – отменяет результата предыдущей операции;

Redo (Повторить) – повторяет результат предыдущей операции;

Simulink – создает модель Simulink (расширения MatLab );

Help Window (Помощь) – открывает окна справки.

4. Формат вывода результата вычислений.

При вводе вещественных чисел для отделения дробной части используется точка !

>> s = 0.3467819

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

>> format long

0.34678190000000

В списке Numerical Format имеются форматы чисел

short – краткая запись числа;

long – длинная запись числа;

short e – краткая запись числа в формате с плавающей точкой;

long e – длинная запись числа в формате с плавающей точкой;

short g – вторая форма краткой записи числа;

long g – вторая форма длинной записи числа;

Формат отображения числовых данных можно установить в меню File (файл) пункт Preferences (предпочтения). Перейти на вкладку Command Window (окно команд). В опции Text display (отображение текста), в списке Numeric format (числовой формат) установить short g , в опции Numeric display (отображение чисел) установить compact . Эти форматы вывода соответствуют выводу чисел в универсальной форме из пяти значащих цифр и с подавлением пробела между строками.

Основы вычислений в MatLab.

Для выполнения простейших арифметических операций в MatLab применяются операторы:

· сложение и вычитание +, – ;

· умножение и деление *, / ;

· возведение в степень ^ .

Некоторые специальные переменные :

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

eps – относительная погрешность при вычислениях с плавающей точкой;

pi – число ;

i или j – мнимая единица;

Inf – бесконечность ;

NaN – неопределенное значение.

Некоторые встроенные элементарные функции MatLab :

exp(x) – экспонента числа x;

log(x) – натуральный логарифм числа x;

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

abs(x) – модуль числа x;

sin(x), cos(x), tan(x), cot(x) – синус, косинус, тангенс, котангенс числа x;

asin(x), acos(x), atan(x), acot(x) – арксинус, арккосинус, арктангенс, арккотангенс числа x;

sec(x), csc(x) – секанс, косеканс числа x;

round(x) – округление числа x до ближайшего целого;

mod(x,y) – остаток от целочисленного деления x на y с учетом знака;

sign(x) – возвращение знака числа x.

Вычислим значение выражение

>> exp(–2.5)*log(11.3)^0.3 – sqrt((sin(2.45*pi)+cos(3.78*pi))/tan(3.3))

Если оператор не удается разместить в одной строке, то возможно продолжение его ввода в следующей строке, если в конце первой строки указать знак продолжения «…», например,

>> exp(–2.5)*log(11.3)^0.3 – ...

sqrt((sin(2.45*pi)+cos(3.78*pi))/tan(3.3))

Функции для работы с комплексными числами :

Ввод комплексного числа

>> z = 3 + 4i

3.0000 + 4.0000i

Функции abs(z), angle(z) возвращают модуль и аргумент комплексного числа , где , ;

complex(a,b) конструирует комплексное число по его действительной и мнимой части:

conj(z)возвращает комплексно-сопряженное число;

imag(z), real(z) возвращает мнимую и действительную часть комплексного числа z.

6. Векторы .

Ввод, сложение, вычитание, умножение на число.

Вектор в MatLab формируется с помощью оператора квадратные скобки . При этом элементы вектора-столбца разделяют точкой с запятой «;», а элементы вектора-строки разделяют пробелом « » или запятой « , ».

Введем вектор-столбец .

>> x =

Введем вектор-строку .

>> y =

Для транспонирования вектора применяют апостроф «’ »:

>> z = y’

Для нахождения суммы и разности векторов используются знаки « + » и «– »:

>> с = x + z

Умножение вектора на число осуществляется как справа, так и слева при помощи знака « * ».

Векторы могут быть аргументами встроенных функций, например,

>> d = sin(c)

Для обращения к элементам векторов используются скобки (), например,

>> x_2 = x(2)

Последний элемент вектора можно выбрать, набрав команду

>> X_end = x(end)

Из нескольких векторов можно составить один, например

>> r =

1.3 5.4 6.9 7.1 3.5 8.2

Символ двоеточие « : » используется для выделения нескольких элементов из вектора, например

>> w = r(3:5)

Символ двоеточие « : » также позволяет заменять элементы вектора, например,

>> r(3:5)= 0

1.3 5.4 0 0 0 8.2

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

>> h =

1 1.2 1.4 1.6 1.8 2

Шаг может быть отрицательным (в этом случае начальное число должно быть больше конечного).

Шаг, равный единице, можно не указывать

>> k =

Основные функции для работы с векторами .

• length(x) – определение длины вектора x;
• prod(x)– перемножение всех элементов вектора x;
• sum(x) – суммирование всех элементов вектора x;
• max(x)– нахождение максимального элемента вектора x;
• min(x)– нахождение минимального элемента вектора x.

Если вызвать функцию min или max с двумя выходными аргументами = min(x),

то первой переменной присваивается значение минимального (максимального элемента), а второй переменной присваивается номер этого элемента.

7 Матрицы .

Различные способы ввода матрицы .

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

>> A =

2. Матрицу можно вводить построчно, выполняя последовательность команд:

>> A = }