Числовой инструментарий, часть вторая |
В первой части мы разобрали основные возможности пакетов GNU/Octave 2.1.34, Scilab 2.6 и Tela 1.32. Сегодня мы поговорим о матрицах, рассмотрим некоторые предопределенные функции и научимся создавать свои собственные, введем операторы управления ходом выполнения программы. В завершении главы мы кратко обсудим имеющихся в пакетах возможности ввода/вывода. МатрицыВекторы хороши, когда данные зависят от одного параметра. Разные значения параметра отображаются в разные индексных значениях. Но если данные зависят от двух параметров, использовать векторы неудобно и требуется структура более общего характера, которая позволяет использование двух независимых индексов. Такая структура называется матрицей. Матрица напоминает только что купленый ящик пива: ячейки образуют прямоугольник и все бутылки -- пардон, элементы -- на месте. Матрицы могут быть составлены из скалярных величин, как в этом примере (GNU/Octave): octave:1> # температура осадки солнечные octave:1> # гр. F дюймы часы octave:1> weather_data = [ 73.4, 0.0, 10.8; ... > 70.7, 0.0, 8.5; ... > 65.2, 1.3, 0.7; ... > 68.2, 0.2, 4.1] weather_data = 73.40000 0.00000 10.80000 70.70000 0.00000 8.50000 65.20000 1.30000 0.70000 68.20000 0.20000 4.10000 Здесь можно увидеть три новых идеи. Во-первых, мы используем комментарии
для того, чтобы дать обозначения колонкам нашей матрицы. Комментарий
начинается с символа "решетки" " Подобно векторам, матрицы могут состоять не только из скаляров, но и из векторов или других матриц.Например, если у нас есть несколько переменных, хранящих ежедневные данные по погоде: weather_mon = [73.4, 0.0, 10.8] weather_tue = [70.7, 0.0, 8.5] weather_wed = [65.2, 1.3, 0.7] weather_thu = [68.2, 0.2, 4.1] мы можем определить weather_data = [weather_mon; weather_tue; weather_wed; weather_thu] или, если мы получаем данные от каждого измерительного инструмента: temperature = [73.4; 70.7; 65.2; 68.2] rain = [0.0; 0.0; 1.3; 0.2] sunshine = [10.8; 8.5; 0.7; 4.1]
weather_data = [temperature, rain, sunshine] Общее правило: запятые разделяют колонки, символы "точка с запятой" отделяют строки. Доступ к скалярным значениям, содержащимся в матрице octave:10> weather_data(3, 2) ans = 1.3000 Можно изменять элементы, присваивая им новые значения: octave:11> weather_data(3, 2) = 1.1 weather_data = 73.40000 0.00000 10.80000 70.70000 0.00000 8.50000 65.20000 1.10000 0.70000 68.20000 0.20000 4.10000 Теперь, когда мы смогли определить матрицу данные_о_погоде_в_тропическом_лесу = weather_data + 2.1 данные_по_зимней_погоде_в_сибири = weather_data / 3.8 не слишком осмыслены, хотя компьютер и не будет жаловаться. В первом
случае он покорно прибавляет Положим, мы хотим сделать с octave:16> temp = [weather_data(1, 1); ... > weather_data(2, 1); ... > weather_data(3, 1); ... > weather_data(4, 1)] temp = 73.400 70.700 65.200 68.200 Очевидно, что индексы столбцов temp = weather_data([1, 2, 3, 4], 1) Вообще говоря, любой вектор может быть использован, как вектор индексов.
Единственно, следите за тем, чтобы не выходить за пределы индексирования.
Порядок индексов имеет значение (напрмер В нашем примере вектор индексов мог бы создаваться специальным встроенным
оператором диапазона [range] " low:high Например: octave:1> -5:2 ans = -5 -4 -3 -2 -1 0 1 2 Пример с погодой упрощается до вида: temp = weather_data(1:4, 1) Поскольу извлечение целого столбца или строки -- очень частая операция,
имеется возможность "срезать угол" еще круче. Если мы опустим и low,
и high в операторе octave:17> temp = weather_data(:, 1) temp = 73.400 70.700 65.200 68.200 Воспользуемся вновь обретенным знанием и извлечем число солнечных часов во вторник, среду и четверг: octave:19> sunnyhours = weather_data(2:4, 3) sunnyhours = 8.50000 0.70000 4.10000 а теперь -- сотояние погоды во вторник: octave:20> tue_all = weather_data(2, :) tue_all = 70.70000 0.00000 8.50000 Преобразования количества выпавших осадков из дюймов в миллиметры
теперь тривиально: умножим вторую колонку
octave:21> rain_in_mm = 25.4 * weather_data(:, 2) rain_in_mm = 0.00000 0.00000 27.94000 5.08000 Мы уже видели, что векторы можно использовать со скалярными величинами: 1.25 + [0.5, 0.75, 1.0] или [-4.49, -4.32, 1.76] * 2 Точно также, скаляры можно использовать с матрицами. octave:1> 1.25 + [ 0.5, 0.75, 1.0; ... > -0.75, 0.5, 1.25; ... > -1.0, -1.25, 0.5] ans = 1.75000 2.00000 2.25000 0.50000 1.75000 2.50000 0.25000 0.00000 1.75000 octave:2> [-4.49, -4.32, 1.76; ... > 9.17, 6.35, 3.27] * 2 ans = -8.9800 -8.6400 3.5200 18.3400 12.7000 6.5400 В обоих случаях результатом будет применение скаляра к каждому элементу вектора или матрицы. А как векторы и матрицы? Очевидно, что выражения подобные [7, 4, 9] + [3, 2, 7, 6, 6] [2, 4; 1, 6] - [1, 1, 9, 4] бессмыслены. В первом случае не согласуется размер векторов (3 и
5 элементов соответственно), во втором случае -- различна форма (матрица
2х2 и четыре столбца). Для того, чтобы операция имела смысл, векторы
и матрицы, учавствующие в сложении или вычитании, должны иметь одинаковую
форму, т.е. одинаковое число строк и столбцов. Технический термин
"форма" [shape] в данном контексте означает размерность. Узнать размерность
любого объекта можно с помощью встроеной функции octave:22> size(weather_data) ans = 4 3 octave:23> size(sunnyhours) ans = 3 1 Ответ -- вектор, первый элемент которого соответствует числу строк,
а второй -- числу столбцов в аргументе Умножение и деление матриц можно задать двумя способами, в наших "числодробилках" реализованы оба.
Подробности
Различия
Встроенные функции матричных вычисленийОх -- это сколько же всего надо упомянуть! Наши числовые инструменты предлагают десятки предопределенных функций. Здесь я могу только подогреть интерес читателя.
Замечание по поводу производительности: все три программы -- интерпретаторы. Это означает, что каждое выражение сначала разбирается [parsed], затем интерпретатор выполняет необходимые вычисления и наконец вызывает внутри выражения функции -- в сравнении с откомпилированной программой это относительно медленный процесс. Но упомянутые ранее функции используются в откомпилированном виде! Они выполняются почти с максимальной скоростью. Интерпретатор в данном случае просто передает матрицу целеком в откомпилированную функцию, написанную на Fortran'е, C или C++. Функция делает свое дело, после чего интерпретатор забирает результат. Мы, таким образом, вывели фундаментальное правило успешного использования "числодробилок": отдавайте предпочтение откомпилированным функциям, а не интерпретируемому коду. Что значительно ускоряет процесс счета. Функции, определяемые пользователемНе зависимо от того, как много предопределенных функций предоставляет программа, их никогда не не будет достаточно. Всегда требуются специализированные функции, решающие конкретные проблемы конкретного пользователя, или же требуется просто сгруппировать часто повторяющиеся предопределенные операторы. Другими словами, всегда есть потребность в функциях, определяемых пользователем Пользовательские функции лучше определять в файлах для того, чтобы
ими можно было воспользоваться в последующих сессиях. В GNU/Octave
файлы с функциями имеют расширение .m, а загружаются либо
автоматически, либо командой
GNU/Octave и Scilab function [res1, res2, ..., resM] = foo(arg1, arg2, ..., argN) # function body endfunction Tela function [res1, res2, ..., resM] = foo(arg1, arg2, ..., argN) { // function body }; где arg1 по argN соответсвуют аргументам функции (также известные, как параметры), а res1 по resN соотсветствуют возвращаемым значениям. Да, ваши глаза не лгут, допускаются множественные возвращаемые значения, что может удивить большинство читателей, знакомых с наиболее распространенными языками программирования. Это, впрочем, является необходимостью, поскольку функции не могут изменять значений входных аргументов. Хватит теории! Давайте напишем функцию, которая на входе принимает матрицу и возвращает матрицу такой же размерности, но с элементами, отмасштабированными в интервале (0, 1). ### Octave function y = normalize(x) ## Возвращает матрицу X, масштабированную в интервале (0, 1). minval = min(min(x)) maxval = max(max(x)) y = (x - minval) / (maxval - minval) endfunction А теперь в Scilab определим функцию, возвращающую спектральный радиус
[spectral radius] матрицы. Мы воспользуемся встроенной функцией // Scilab function r = spectral_radius(m) // Возварщает спектралный радиус R матрицы M. r = max(abs(spec(m))) endfunction Наконец, в Tela напишем функцию, которая вычесляет Евклидову (Фробениуса) норму [Frobenius norm] для данной матрицы. // Tela function x = frobenius(m) // Возвращает Фробенисову норму X of matrix M. { x = sqrt(sum(abs(m)^2)) }; "Автомагическая" загрузка файлов функций в GNU/Octave работает следующим
образом: когда Octave сталкивается с функцией без определения, просматривается
список директорий, задаваемый встороенной переменной Операторы управления ходом выполнения программыВесь код, который мы написали до сих пор, выполняется строго последовательно, мы не использовали таких операторов, как операторы условного перехода или операторы цикла. Прежде, чем мы начнем манипулировать с передачей управления, нам нужно рассмотреть логические выражения, поскольку на них опирается провека условий в операторах условного перехода и в циклах. Логические выражения образуются из (1.) чисел, (2.) операций сравнения и (3.) логических выражений, объединенных в одно целое с помощью логических операторов.
Теперь мы готовы для первого условного оператора, а именно оператора
условного перехода # Octave // Scilab // Tela if ( условие ) if условие then if ( условие ) { # если да // если да // если да else else } else { # иначе // иначе // иначе endif end }; условие является логическим выражением в описанном выше смысле. оператор цикла # Octave // Scilab // Tela while ( условие ) while условие while ( условие ) { # тело цикла // тело цикла // тело цикла endwhile end }; условие, конечно же, является логическим выражением. В Octave и Scilab оператор # Octave // Scilab // Tela for var = expr for var = expr for (init; cond; step) { # работа // работа // работа endfor end }; Вот несколько примеров. Octave function n = catch22(x0) ## Знаминитая функция catch-22 (двадцать второй капкан:) ## или уловка-22 (по названию знаменитого романа Дж. Хеллера, ## подробностей переводчик не знает:) : ## невозможно определить заранее, завершатся ли вычисления ## для каждого конкретного значения аргумента. ## Возвращает число итераций. n = 0 x = x0 while (x != 1) if (x - floor(x/2)*2 == 0) x = x / 2 else x = 3*x + 1 endif n = n + 1 endwhile endfunction Scilab function m = vandermonde(v) // Возвращает M: матрицу Вандермонде, // основанную на векторе V. [rows, cols] = size(v) m = [] // empty matrix if rows < cols then for i = 0 : (cols-1) m = [m; v^i] end else for i = 0 : (rows-1) m = [m, v^i] end end endfunction Tela function vp = sieve(n) // Решето Эратосфена; возвращает вектор VP всех // простых чисел VP, которые строго меньше 2*N. // 1 не рассматривается, как простое число. { vp = #(); // пустой вектор if (n <= 2) { return }; vp = #(2); flags = ones(1, n + 1); for (i = 0; i <= n - 2; i = i + 1) { if (flags[i + 1]) { p = i + i + 3; vp = #(vp, p); for (j = p + i; j <= n; j = j + p) { flags[j + 1] = 0 } } } }; Ввод/ВыводМы долго-долго работаем с одним из наших цифровых инструментов и в определенный момент нам уже хочется закончить наш долгий рабочий день. Но нам и не хочется потерять все наработанные за день результаты. Функции мы уже сохранили в файлах с исходным кодом. Пришло время узнать, как дать возможность существования вне программы нашим данным. Простые операции ввода/выводаВо всех трех программах используется одинаковая модель ввода/вывода (I/O), в значительной степени заимствованная из языка программирования C. В этой модели есть возможность тонкого управления нюансами чтения и записи. Однако нередко нет необходимости прямого задания формата записываемого файла. Если все, что нужно -- сохранить переменные для того, чтобы ими можно было бы воспользоваться в следующий раз, то вполне подойдут и упрощенные команды I/O.
Ввод/вывод, ориентированный на матрицыПоскольку матрицы используются очень часто, то для их загрузки и сохранения матриц целиком имеются специальные команды. Считывание матрицы одной командой очень практично и удобно для загрузки результатов эксперемента или данных, созданных другими программами. Предположим, что у нас есть ASCII файл datafile.ascii, содержащий строки # run 271 # 2000-4-27 # # P/bar T/K R/Ohm # ====== ====== ====== 19.6 0.118352 0.893906e4 15.9846 0.1 0.253311e5 39.66 0.378377 0.678877e4 13.6 0.752707 0.00622945e4 12.4877 0.126462 0.61755e5 и находящийся в рабочей директории. Пять первых строк файла -- не числовые. Наши программы пропускают их при чтении, но они могли бы помочь пользователю разобраться с данными. Я специально взял не слишком аккуратно отформатированные данные, как это обычно и бывает в жизни. Функции загрузки матриц разбивают читаемый поток по пробельным символам, а не по конкретным колонкам, так что файл datafile.ascii им понравится. octave:1> data = load("datafile.ascii") data = 1.9600e+01 1.1835e-01 8.9391e+03 1.5985e+01 1.0000e-01 2.5331e+04 3.9660e+01 3.7838e-01 6.7888e+03 1.3600e+01 7.5271e-01 6.2294e+01 1.2488e+01 1.2646e-01 6.1755e+04 или в Scilab'е -->data = fscanfMat("datafile.ascii") data = ! 19.6 0.118352 8939.06 ! ! 15.9846 0.1 25331.1 ! ! 39.66 0.378377 6788.77 ! ! 13.6 0.752707 62.2945 ! ! 12.4877 0.126462 61755. ! а так в Tela >data1 = import1("datafile.ascii") >data1 #( 19.6, 0.118352, 8939.06; 15.9846, 0.1, 25331.1; 39.66, 0.378377, 6788.77; 13.6, 0.752707, 62.2945; 12.4877, 0.126462, 61755) Во всех трех примерах data будут содержать матрицу 5 на 3, заполненную значениями из datafile.ascii. Парные команды для сохранения отдельной матрице в файле формата ASCII: save("data.ascii", "data") # GNU/Octave fprintfMat("data.ascii", data, "%12.6g") // Scilab export_ASCII("data.ascii", data) // Tela Обратите внимание, что Ни одна из описанных выше команд, естественно, не сохраняет оригинальный заголовок: т.е. строки файла datafile.ascii, начинающиеся с символа "решетки". Для того, чтобы записать и их, необходимо прибегнуть к "низкоуровневым" функциям ввода/вывода в стиле языка C, которые присутствуют в каждой из рассматриваемых программ. Ввод/вывод в стиле языка CДля точного управления вводом и выводом предлагается модель C I/O. Во всех трех приложениях реализована функция printf(format, ...) Более того, GNU/Octave и Tela следуют схеме именования языка C: handle = fopen(filename) fprintf(handle, format, ...) fclose(handle) в то время, как Scilab добавляет к именам функций префикс " handle = mopen(filename) mprintf(handle, format, ...) mclose(handle) Не зависимо от того, называются ли функция В следующем месяце: Графика, графики функций и отображение данных. Кристоф Шпиль [Christoph Spiel]Крис управляет консалтинговой компанией по открытому программному обеспечению в верхней Баварии, Германия. По своей специальности - физике -- у него есть докторская степень от Munich University of Technology -- его основные интересы лежат в области теории чисел, гетерогенных программных сред и программного инжиниринга. С ним можно связаться по адресу [email protected]. Copyright (C) 2001, Christoph Spiel.
|
Вернуться на главную страницу |