Элементы математического моделирования в программных средах MATLAB 5 и Scilab (Андриевский Фрадков)

получилась вектор-строка, в чем убеждаемся, выполнив оператор size (а). В ответ получим: ans = 1 9 . Следовательно, элемент а есть 9-мерная вектор-строка. По принятой в программировании терминологии переменная а является индексированной переменной (именем массива) с двумя измерениями. По первому измерению индекс может принимать значение, равное единице, а по второму - изменяться от одного до девяти. При описании языка MATLAB более естественным является использование математических терминов, что мы обычно и будем делать. В тех ситуациях, когда это

не удается, будем оперировать с понятием массивов				данных
(а не матриц и векторов).	В рассмотренном	выше	примере
оператор-функция size записан без левой части.			Это при-
водит к тому, что результат его выполнения		присваивается
переменной с резервированным именем ans (от англ.				answer
- ответ).
Выполним теперь оператор
> Ъ = а + 2
Мы получили вектор-строку Ъ, каждый элемент			которой
равен значениям элементов вектора а, увеличенных на 2:
Ъ = 3 4 5 6 8 6 5 6 7
Создадим теперь квадратную матрицу А =			2	0
			5	- 1
Для этого выполним команду		10		- 1
> А = [1 2 0; 2 5 -1;	4 10 -1] ;

Заметим, что результат выполнения операции не выводится на дисплей. Это связано с тем, что команда оканчивается символом "точка с запятой".

Пусть теперь матрица В = А (символ "т" обозначает операцию транспонирования). Выполнив оператор

> В=А'

получим

В =
1	2	4
2	5	10
0	-1	-1

31

Произведение этих матриц найдем с помощью оператора > С=А*В;

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

> с=а*Ъ приведет к сообщению об ошибке

??? Error using ==> *

Inner matrix dimensions must agree

- "Неправильное использование *. Внутренние матрицы должны иметь соответствующие размеры".

Можно выполнить оператор > с=а*Ъ'

тогда получим в ответе с = 196, так как мы нашли скалярное произведение векторов а и Ь, т.е. вычислили с = ab . Допускается и выполнение оператора

> с=а'*Ъ;

Тогда в ответе получится матрица размером 9x9 единичного ранга (чтобы проверить это, выполним операторы size (с) и rank(c)). Если же мы хотим выполнить поэлементное умножение массивов а и Ъ, то нам следует использовать оператор

.*, т.е. написать

>с=а.*Ъ

Вответе получим

с = 3 8 15 24 48 24 15 24 35

Этот же оператор применим и к массивам большей размерности, в чем убеждаемся, выполнив действие С=к.*В.1 Для вычисления обратной матрицы используем оператор

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

32

>D=inv(A)

асобственные векторы и собственные числа матрицы А находятся оператором

>[v,p]=eig(A).

Коэффициенты характеристического многочлена матрицы

Авычисляются оператором

>p=poly(А),

выполнив который, получим:

р = 1.00 -5.00 5.00 -1.00 Это значит, что характеристический многочлен

det(AI3 - А) = А3 - 5А2 + 5А - 1.

2.1.2.Действия с многочленами

Выше показано использование функции poly с матричным аргументом. Если аргумент - вектор, то его элементы рассматриваются как корни многочлена, коэффициенты которого следует определить. Например, в результате выполнения оператора

> d=poly( [1,2,3,4]) получим

d = 1 -10 35 -50 24

Это	значит, что d(А) = А4 - 10А3 + 35А2 - 50А + 24 обращается
в ноль	при А € {1,2,3,4}.

Корни многочлена вычисляются функцией roots. Например, выполнив оператор

> r=roots(p) получим:

г =

3.7321

1.0000

0.2679

2 Б. Р. Андриевский и др.

33

Лля умножения многочленов 1 используем функцию свертки векторов conv:

> s=conv(p,[l 1 1 ] )

Вектор s состоит из коэффициентов произведения многочлена р на многочлен А2 + А + 1 и равен

s = 1.00 -4.00 1.00 -1.00 4.00 -1.00. Вычислим корни многочлена 5(A):

> z=roots(s) Получим:

Z =

3.7321

-0.5000 + 0.8660i -0.5000 - 0.8660i

1 . 0 0 00

. 0.2679

Таким образом, результат содержит мнимые значения: пакет MATLAB позволяет выполнять операции и с мнимыми числами. По умолчанию идентификаторы i, j резервированы для мнимой единицы, но пользователь может их переопределить. Отобразим полученные значения на графике. Лля этого выполним операторы

>p l o t ( z , ' * ' )

>xlabelORe \lambda' ) ,ylabel(' Im \lambda') ,grid

Результат показан на рис. 2.1, а. 2

Оператор grid служит для нанесения сетки на график. Заметим, что в системе MATLAB 5 можно отображать и буквы греческого алфавита. Лля этого они в программе должны быть закодированы в формате языка TgX [57].

Построим теперь график функции s(A) для значений А € 6 [0,4] с шагом ДА = 0.05. Выполним операторы

>lambda=0:0.05:4;

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

2В данном издании рисунки, полученные копированием графического окна пакета MATLAB, оставлены в первоначальном виде.

34

>y=polyval(s,lambda);

>plot(lambda,у),grid

>xlabelCUambda') , t i t l e ( 1 s(\lambda) ')

Рис. 2.1. Корни и значения многочлена s(A).

Результат показан на рис. 2.1, 6 Функция polyval служит для вычисления значений многочлена в заданных точках. Обратим внимание на то, что аргумент lambda является вектор-строкой, элементы которой образуют арифметическую прогрессию с нулевым начальным значением и шагом 0.05. Оператор size(lambda) дает ответ ans = 1 81.

Деление многочленов выполняется функцией deconv. Разделим многочлен s(A) на А — 1:

> [q,r]=deconv(s,[1 - 1])

Получим результат - многочлен q(А), коэффициенты которого находятся в массиве

q = 1.00 -3.00 -2.00 -3.00 1.00 и остаток от деления - многочлен

г = 0 0 0 2.2204е-016 0 -2.4425е-015

Функция polyvalm служит для вычисления многочлена от матричного аргумента. Например, выполнив оператор

> polyvalm(p,A)

получим ЗхЗ-матрицу, элементы которой близки к нулевым (не превосходят по модулю 10~13). Так и должно быть по теореме Кэли-Гамильтона, согласно которой каждая квадратная матрица удовлетворяет своему характеристическому уравнению.

35

2.1.3.Действия над функциями

Пользователь имеет возможность определить свою подпро- грамму-функцию, оформив ее в виде т-файла. 1 Имя этого файла и определяет имя функции. Начальная строка подпро- граммы-функции должна быть оформлена по определенным правилам. Она начинается словом function, затем указываются (если несколько, то в квадратных скобках) выходные переменные, потом после знака равенства, в круглых скобках, входные переменные. Имена переменных, используемых в подпрограмме, локализованы внутри нее. Этим подпрограмма отличается от так называемого скрипт-файла, который также является m-файлом, но не содержит заголовка и представляет собой просто последовательность команд. Эффект выполнения скрипт-файла такой же, как если бы эти команды были набраны в командной строке "вручную". Если значения каких-то переменных должны передаваться подпрограмме, минуя аппарат формальных параметров, то они должны быть объявлены глобальными с помощью оператора global. Этот оператор должен присутствовать во всех программах, использующих указанные переменные.

Рассмотрим в качестве примера подпрограмму-функцию, в которой вычисляется значение многочлена 5(A). Создадим гп-файл с именем pol.m, в котором запишем текст 2

function y=pol(lambda)

s = [1.00 -4.00 1.00 -1.00 4.00 -1.00]; y=polyval(s,lambda)

Покажем, как можно найти экстремальные точки этой функции. (Отвлечемся пока от того, что мы знаем аналитическое выражение для s(А)). Для нахождения экстремума функции можно использовать процедуры fmin, fmins.

Выполним оператор

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

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

36

>xl = fminCpol' ,0,5)

Вответе получим xl = 3.0159. При обращении к функции fmin указывается имя функции, минимум которой следует получить (точнее говоря - имя m-файла, содержащего процедуру вычисления данной функции), а также границы интервала минимизации. Исполняя оператор

>yl=pol(xl)

найдем значение функции в точке экстремума у1 = -5.2016е+001

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

>y=pol(lambda);

>plot(lambda,у,xl,yl, ' * ' )

получим график функции, на котором символом * отмечена точка экстремума (см. рис. 2.1, 6). 1 Если внутри заданного интервала имеется точка локального минимума, то она выводится в качестве решения. В противном случае решение находится на границе интервала. Например, при выполнении оператора х2 = fmin('pol' , 1,2) получается х2=2.

Несколько иначе выполняется процедура fmins. Она позволяет находить экстремум функции нескольких переменных. Входным параметром процедуры fmins является начальное приближение к точке минимума (в общем случае - вектор). Диапазон поиска экстремума не ограничен. Поэтому операторы > xl=fmins( 'pol1 ,1) и > xl=fmins( 'pol' , 12) дают одинаковое значение xl = 3.0158е+000. Если функция имеет несколько локальных минимумов, или не имеет глобального минимума, либо обладает большим по норме градиентом, то выбор начального приближения весьма существен. В нашем примере оператор xl=fmins( 'pol' ,13) приводит к появлению предупреждающего сообщения и значению в ответе xl = -8 .2397е+029. Аналогичная ситуация возникает при вызове xl=fmins( 'pol' ,0).

1 Лля вывода графиков функций можно использовать также подпрограмму ezplot.

37

Найти решение уравнения f(x) = 0 можно с помощью функции fsolve. Например, выполнив оператор

xl = f s o l v e ( ' p o l \ 1 0 )

получим один из корней многочлена s(A) xl = 3.7321, а вызовом x2=fsolve( 'pol' ,0) - другой корень х2 = 0.2679. Программа fsolve может использоваться и для решения систем алгебраических уравнений. Вычислить нули функции от одной переменной можно также с помощью программы fzero. Например, выполнение оператора

> x l = f z e r o ( ' p o l ' , 0 ) приводит к результату

Zero	found in the i n t e r v a l : [-0.32, 0.32].
(Ноль найден на интервале [ - 0 . 32, 0.32]).
xl =	2.6795е-001

Значение определенного интеграла J f ( x ) d x можно вычислить подпрограммой quad: оператор

> yi=quad( , pol', - 1,1) выводит значение уi = -4.2667.

2.1.4Построение трехмерных изображений

Выполним операторы

[х,у]=meshgrid(-2:0.25:2); z=x.~2-y.~2;

mesh(x,y,z)

x l a b e l O x ' ) , y l a b e l ( ' y ' ) title('z=x~2 - y~2')

Оператор meshgrid создает двумерные массивы аргументов (х и у) с заданными интервалом и шагом изменения. Далее вычисляются значения z = х2 — у2 и оператором mesh строится поверхность г = z(x,y). Результат показан на рис. 2.2, а. Выполним теперь операторы

[al,bet]=meshgrid([ - pi:pi/20:pi]); x=sin(al) . *(2+sin(bet)); y=cos(al).*(2+sin(bet)); z=cos(bet);

38

Рис. 2.2. Пространственные изображения.

mesh(x,y,z)

a x i s ( [ - 3 3 - 3 3 - 3 3 ] ) axis square

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

полнение первого из них задает диапазоны координатных осей "вручную". (По умолчанию оно происходит автоматически). Второй оператор устанавливает на изображении равные размеры рамки координатных осей.

2.2.Примеры использования системы MATLAB для численного решения задач исследования систем

2.2.1.Применение тулбокса символьных вычислений

Проиллюстрируем рассмотренный выше пример вычислениями в пакете MATLAB 5. Сначала рассмотрим получение расчетных соотношений в общем виде. Для этого используем средства тулбокса символьных вычислений (SYMBOLIC TOOLBOX) [35, 84].

Площадь треугольника ОМК отвечает величине ут. В программе обозначим ее через у_Т, момент окончания движения

t+ - через T, момент переключения	t' - через	t _ l, наиболь-
шую скорость vm (точка М на рис.	1.7, а) -	через v_m, за

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

значения.	Очевидно,	выполнены геометрические		соотноше-
ния: у. =	tgp, =	^ =	tg^2 = ш =	С учетом

39

введенных обозначений напишем программу вычисления неизвестных t',tm)vm :

S=solve(' v ^ / t ^ l s b / m ' , 'v _ m/(T - t _ l)=a/m' ....

' у . Т Ф ч л / г 1 , Ч _ 1 ' , ' Т ' , ' v . m ' ) t_l=simple(S.t_l)

T=simple(S.T) v_m=simple(S.v.m)

Воператоре solve записана исходная система уравнений

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

t . l		=
С		1/(b+a)*2~(l/2)*(m*(b+a)*y_T*a*b)~(1/2)/Ь]
[	-1/(Ь+а)*2~(1/2)*(m*(b+a)*y_T*a*b)~(1/2)/Ъ]
Т =
[		у_Т*т*(b+а)*2~(1/2)/(т*(Ъ+а)*у_Т*а*Ь)~(1/2)]

[ -y_T*m*(b+a)*2~(l/2)/(m*(b+a)*y_T*a*b)~(1/2)] v_m =

[	1/т/(b+a)*2~(1/2)*(т*(b+a)*y_T*a*b)~(1/2)]
[	-l/m/(b+a)*2~(1/2)*(m*(b+a)*y_T*a*b)~(1/2)]

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

граммирования		Фортран. Лля этого используем обращение
f o r t r a n ( t l ) ,	которое выводит строку
t . 1 ( 1 , 1 ) = 1/(b+a)*sqrt(2.ЕО)*sqrt(m*(b+a)*y_T*a*b)/b
Аналогично,	ccode(tl) получает оператор языка Си:
t.lCO][0]	=	l/(b+a)*sqrt(2.0)*sqrt(m*(b+a)*y_T*a*b)/b;

1 Отрицательные решения можно трактовать как движение в противоположном направлении и обратном времени.

40