Управление потоком команд

 

MATLAB имеет пять видов структур управления потоками:

• оператор if

• оператор switch

• циклы for

• циклы while

• оператор break

 

if

Оператор if вычисляет логическое выражение и выполняет группу операторов, если выражение истинно. Необязательные ключевые слова elseif и else служат для выполнения альтернативных групп операторов. Ключевое слово end, кото­рое согласуется с if, завершает последнюю группу операторов. Таким образом, все группы операторов заключены между четырех ключевых слов, без исполь­зования фигурных или обычных скобок.

Простейший условный оператор имеет вид:

If BOOL, EXPR, end

Команды EXPR будут выполнены, если истинно условие BOOL.

Ветка else дает альтернативу:

If BOOL, EXPR, else EXPR_2, end

Напомним, что запятая (или точка с запятой) после условия BOOL необязательна, если команды EXPR расположить на следующей строке. При добавлении еще одно­го служебного слова - elseif - условный оператор становится мощнее:

If BOOL_1, EXPR_1

Elseif BOOL_2, EXPR_2

Elseif BOOL_3, EXPR_3

…

Else EXPR_0

End

Здесь проверки следуют одна за другой. В том случае, если условие BOOL_N истин­но, будет выполнен оператор EXPR_N. Например, при х=4 в результате выполнения условного оператора:

if x<0, 0, elseif x<3, 3, elseif x<5, 5, else 13, end

получим:

ans =

Важно понять, как операторы отношения и оператор if работают с матрицами. Когда операторы отношения применяются к матрицам одного размера, результатом является матрица того же размера, у которой в качестве элементов стоят 0 или 1, в зависимости от соотношения между соответствующими элементами исходных матриц. В этом случае логическое выражение будет оцениваться как Истина только тогда, когда каждый элемент матрицы равен 1 (Истине). Поэтому при исполь­зовании отношений с оператором if большую помощь могут оказать такие функции, как isequal, isempty, all, any.

Когда А и В матрицы, то при проверке их на неравенство вместо конструкции

if A ~= B, <операторы>, end

следует ввести

if any(any(A ~= B)), <операторы>, end

или, что проще,

if A == B, else <операторы>, end

Заметим, что первая конструкция (с A ~= B) почти наверняка не даст того, что нужно, поскольку оператор будет выполняться только если каждый элемент матрицы A будет отличаться от соответствующего элемента матрицы B. Для сведения матричных отношений к вектору или скаляру можно воспользоваться функциями any и all. В предыдущем примере необходимо использование функции any два раза, поскольку эта функция – векторная.

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

If isequal(A,B),...

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

if A>B

'greater'

elseif A<B

'less'

elseif A==B

'equal'

elseif A~=B

'not equal'

Else

Error (' Непредвиденная ситуация ')

End

 

Switch и case

Оператор-переключатель switch предоставляет возможность разветвления и вы­полнения операторов в зависимости от значения переменной VAR (скаляра или стро­ки):

Switch VAR

Case VAR1, EXPR1

case {VAR2, VAR3,...}, EXPR_2

Otherwise EXPR_0

End

Здесь VAR1, VAR2, VAR3,... - различные значения, которые может принимать пере­менная VAR, a EXPR_1, EXPR2 и т. д. - группы операторов. Переключение возможно по единственному значению (выбор VAR1) и по группе значений (выбор из VAR2, VAR3,...). Кроме того, если осознанного выбора не сделано, то выполняются опера­торы группы EXPR_0 (ключевое слово otherwise). Выполняется только первый соответствующий случай, таким образом, нет необходимости в использовании оператора break, как в языке Си.

Приведем пример со строковой переменной, здесь команда disp используется для вывода текстового сообщения:

METHOD = 'Симп';

Switch METHOD

case {'трап','прям'}

Disp('Meтод 2 порядка')

case 'Симп'

Disp('Meтод 4 порядка')

Otherwise disp('Метод неизвестен')

End

В результате получим: Метод 4 порядка

 

For

Перечислительный цикл чаще всего имеет вид:

for N = N0: DN: N1, EXPR, end

Здесь N - счетчик, пробегающий с шагом DN значения от N0 до N1 (учитываются вещественные части чисел); a EXPR обозначает выполняемые в теле цикла ко­манды.

for n = 3:32

r(n) = rank(magic(n));

End

r

Точка с запятой после выражения в теле цикла предотвращает повторения вы­вода результатов на экран, а r после цикла выводит окончательный результат.

Оператор for допускает использовать любую матрицу вместо N0: DN: N1. Пусть A – матрица с n столбцами. Тогда строка

for V = A,..., end

равносильна следующей

for k = 1:n, V = A(:,k);..., end

за исключением того, что в последнем случае в рабочем пространстве появляется ещё и переменная k

 

While

Цикл с условием:

While BOOL, EXPR, end

выполняется до тех пор, пока не будет нарушено условие BOOL.

Ниже приведена полная программа, иллюстрирующая работу операторов while, if, else и end, которая использует метод деления отрезка пополам для на­хождения нулей полинома.

a = 0; fa = -Inf;

b = 3; fb = Inf;

while b-a > eps*b

x = (a+b)/2;

fx = x^3-2*x-5;

if sign(fx) == sign(fa)

a = x; fa = fx;

Else

b = x; fb = fx;

End

End

x

Результатом будет корень полинома х³ -2х-5

х =

2.09455148154233

Для оператора while верны те же предостережения относительно матричного сравнения, что и для оператора if, которые обсуждались ранее.

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

tic, for k = 1:23456, y = sin(k); end, toc

elapsed_time =

0.014780

и расчета с применением векторизации:

tic, k = 1:23456; y = sin(k); toc

e1apsed_time =

0.005615

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

n =1000; tic, for k=1:n, B(k)=k; end; toc

elapsed_time =

0.001883

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

n =1000; A=zeros(1, n); tic, for k = 1:n, A(k)=k; end; toc

elapsed_time =

0.000593

Break и continue

Оператор break позволяет досрочно выходить из циклов for или while. Во вло­женных циклах break осуществляет выход только из самого внутреннего цикла.

Ниже представлен улучшенный вариант предыдущего примера. Почему исполь­зование оператора break в данном случае выгодно?

a = 0; fa = -Inf;

b = 3; fb = Inf;

while b-a > eps*b

x = (a+b)/2;

fx = x^3 - 2*x - 5;

if fx == 0

Break

elseif sign(fx) == sign(fa)

a = x; fa = fx;

Else

b = x; fb = fx;

End

End

x

 

В МATLAB 6 появился оператор continue, означающий пропуск всех следующих до end операторов и переход к следующей итерации цикла, начиная с первого его оператора.

 

Практикум

Лабораторная работа 2. Примеры, иллюстрирующие эффективность MATLAB'а

 

1. Суммирование. Найдем при заданном n частичную сумму ряда s(n) = 1/k^2, k=1:n. Для этого выполним строку

Команда cumsum(f) подсчитывает все частичные суммы s(k) от f(1:k) для каждого k от 1 до n, так что на графике можно наблюдать процесс формирования нужной нам величины. В конце строки выдается численный и точный результаты:

ans = 1.6350 1.6449.

Полагая n=1000, получим

ans = 1.6439 1.6449,

т.е. ошибку в 1 единицу 4-й значащей цифры.

Сходимость не всегда столь очевидна, как на этом графике. Чтобы в этом убедиться, усложним наш пример: при заданных m>1 и n найдем частичную сумму ряда s(m,n) = sum(1/k^m), k=1:n (при m=1 получается уже расходящийся гармонический ряд). Для проведения вычислений отредактируем строку 1:

2; m=2; n=1000; k=1:n;f=k.^(-m); plot(cumsum(f)), sum(f)

и сначала для проверки получим свой старый результат. Но уже при m=1.5 у нас, глядя на график, нет полной уверенности в достижении сходимости. Это тем более так при m=1.2: для n=1000 ans=4.3358, а для n=1e4 ans=4.7991. Факт сходимости ряда при m=1.01 нельзя установить численно из-за низкой скорости его сходимости.

Чтобы лучше запомнить действие команды cumsum, вычислим ò(x/sin(x))dx, xÎ[0, 3]. Подынтегральная функция f=x/sin(x) не имеет в нуле особенности, и поэтому достаточно выполнить строку

3; n=100; h=3/n; x=h/2:h:3-h/2; f=x./sin(x); plot(h*cumsum(f)), grid, h*sum(f)

т.е. аппроксимировать f в серединах интервалов (эти точки x называют полуцелыми в отличие от концов счетных интервалов – целых точек). Сравнение ответа ans = 8.4495 и графика наводит на мысль о том, что пока сходимость еще не достигнута, но при n=1e3 ans = 8.4552, так что при n=1e2 со сходимостью в действительности все в порядке, а возрастание функции h*cumsum(f) на правом конце происходит из-за роста там функции f – это можно увидеть, выполнив

Plot(f)

Для матрицы A команды sum и cumsum работают вдоль столбцов (значит, по первому индексу), а для вектора – вдоль него независимо от того, строка это или столбец. Чтобы провести суммирование для матрицы A вдоль ее строк, нужно выполнить sum(A,2), т.е. указать для выполнения команды второй индекс. Это правило относится ко многим командам MATLAB'a и к многомерным матрицам тоже – по умолчанию имеется в виду первый индекс, а в противном случае нужно всегда указывать, по какому индексу должна работать команда, и это указание не сохраняется для последующих команд.

 

 

2. Произведения. Аналогично суммированию с помощью команд prod и cumprod вычисляются и обрабатываются произведения. Например, найдем Õ(1-1/k^2), k=2:100 (при k®¥ Õ®1/2), выполнив строку

1; n=100; k2=(2:n).^2; a=1-1./k2; cp=cumprod(a); cp(end), plot(cp/.5), grid

Результат cp(end) = 0.5050 говорит о том, что сходимость здесь не очень быстрая. Это видно и из графика, на котором представлена относительная ошибка результата. Обратите внимание на названия переменных k2=k^2 и cp=cumprod(..): при выборе имен переменных всегда нужно стремиться к тому, чтобы эти имена хоть как-то отражали суть дела (это особенно важно при написании больших программ, где много переменных).

При вычислении произведений можно выйти за числовую шкалу. Найдем, например, для каких k можно найти k!. Ясно, что максимально допустимое km вряд ли больше 200, так что строка

должна дать ответ на наш вопрос. Из-за быстрого возрастания kf и ограниченной разрешимости дисплея (это не более 0.5% от максимального значения на графике) мы видим всего одну точку kf(km), перед которой, как нам ошибочно кажется, идут нули и за которой идут числа inf (infinity), вообще никак не представленные на рисунке. Точно так же графика обходится и с переменной NaN (not a number), и это обстоятельство может быть иногда полезным. Переменная NaN возникает в таких ситуациях:

Inf-inf inf/inf

Переменные inf и NaN (они получаются со знаком) можно использовать в программах. Для определения km выполним строку

Sum(isinf(kf))

в которой isinf(kf) выдаст 1 на тех позициях вектора размеров kf, где элементы kf есть inf, и 0 на остальных позициях. Поскольку ans=30, km=n-30=170, что можно было бы получить и сразу, выполнив строку

4; km=sum(isfinite(kf))

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

log (взятие натурального логарифма),

log10 (взятие десятичного логарифма),

abs (взятие модуля),

sign (взятие знака, выдающее 1, 0 и -1).

 

 

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

1. Напишем строку для нахождения общих элементов двух векторов:

x=1:20; y=15:30; [X,Y]=meshgrid(x,y); v=X(X==Y)

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

1; tic, x=0:1/200:1; [X,Y]=meshgrid(x); M=abs(X+i*Y-.5-i*.5)<1/2; s=sum(M(:)), t1=toc

и даст ответ s=31397 точек, t1=0.120 сек, тогда как строка для циклов for-end

2; tic, s=0; w=1:201; for I=w, for J=w, if norm([x(I),x(J)]-.5)<.5,s=s+1; end, end, end, s,t2=toc

дает то же самое s и t2=1.69 сек, так что t2/t1=14. Это лишний раз говорит о том, что нужно разумно подходить к использованию операторов языка программирования.

Упражнения.

1. Из заданной матрицы A размера m*n построить матрицу B с m строками, у которой диагонали с номерами 0, 1, …, n-1 были бы столбцами A с номерами 1:n, а все остальные элементы равнялись бы нулю (решение с помощью циклов).

 

2. Усовершенствуйте, как можете, обе строки последнего примера практикума и посмотрите, как изменится значение t2/t1.

 

3. Проверьте, что при вычислении определителя квадратной матрицы M порядка n методом Гаусса нужно выполнить @2n³/3 операций сложения и умножения. Считая их одинаковыми по времени выполнения, найдите, за какое время t выполняется 1 млн операций, взяв n=500 и выдав время счета det(M).

 

4. Используя предыдущий результат, оцените число арифметических операций для вычисления функций sin и exp.

 

5. Оцените число операций, которое затрачивается в MATLAB'е для организации одного шага в цикле for-end.