автор - Виктор Кон (Курчатовский институт), URL: http://kohnvict.narod.ru
Содержание
01. Вычисление и показ зависимостей
02. Суммирование матрицы и график
03. Комплексные вычисления
04. Вычисление точки прохождения функции через ноль
1. Вычисление и показ зависимостей
Самой простой задачей, на мой взгляд, является вычисление каких-то функций, то есть зависимостей каких-то величин от аргумента. Аргумент изменяется с постоянным шагом от начального значения и нас интересует конечное число точек в этой цепочке. Очевидно, для таких расчетов надо организовать цикл, в котором описать процедуру вычисления функций при постоянном изменении аргумента. Далее надо каким-то образом запоминать результаты расчета, чтобы их потом можно было показать в каком-то виде. Ну и надо иметь средства, чтобы реально их показать. В качестве функций возьмем функции Бесселя целых индексов 0, 1 и 2.
На языке ACL данная задача имеет несколько решений. Рассмотрим первое решение. Прежде всего, функции Бесселя J0(x) и J1(x) в ACL считаются элементарными и вычисляются непосредственно. А функцию J2(x) можно выразить через первые две с помощью рекуррентной формулы J2(x)=(2/x)J1(x)-J0(x). Нам надо задать начальное значение x0, шаг x1 и число точек n. Затем необходимо организовать цикл. Результаты в цикле будем конвертировать сразу в текст и запоминать в текстовом массиве. После окончания цикла записанный текст надо показать в окне редактора в режиме показа. Сложность возникает только в первой точке, если начинать от нуля. Ее надо выделить отдельно и написать ответы в этой точке в явном виде. Вот полный текст программы.
Этот текст можно скопировать в редактор интерпретатора vkACL.jar и выполнить, нажав клавишу F12. В результате появится окно, в котором будет показана вот такая таблица. А затем она же будет записана в файл bessel.txt.
Рассмотрим более подробно, что же написано в коде программы. Сначала задаются параметры расчета. Среди них надо отметить команду s(3)=1; . Элемент s(3) служебного массива задает позицию курсора текстового массива. Команда печати в массив #pri начинает печатать именно с этой позиции. При печати она автоматически перемещает курсор на новое место за напечатанным текстом. Таким образом, следующие команды продолжают заполнять текстовый массив последовательно. Чтобы узнать сколько же символов было напечатано, достаточно из последнего значения параметра s(3) вычесть его первое значение. Именно это и было сделано при показе текста в текстовом редакторе, а затем при его записи в файл.
Команда печати использует форматы и специальные символы. Формат типа \G диктует интерпретатору каким именно образом надо конвертировать число в текст, а символ \n указывает на перенос текста в следующую строку. Каждая из используемых команд описана в полном техническом описании. Здесь я пытаюсь объяснить саму стратегию языка. В данном решении задачи текст сразу показывается в текстовом редакторе, а потом записывается в файл. Но его можно было и не записывать. Текст можно просто скопировать в файл через буфер обмена. Хотя, конечно, записать автоматически удобнее для пользователя.
Но возможно и другое решение задачи, когда результаты в цикле запоминаются в реальном массиве, в виде чисел. Затем эти числа записываются в файл по соответствующему формату. А уже потом, из файла показываются в текстовом редакторе. Удобство такого подхода в том, что числа можно записывать не только в текстовом формате, но и в формате данных. В таком формате они записываются в файл намного компактнее, но зато стандартные тестовые редакторы их не показывают. Чтобы их посмотреть надо их снова прочитать в формате данных, а затем конвертировать в текст. Текст программы, выполняющей второе решение выглядит так
А результат ее выполнения выглядит так
Легко заметить, что во втором случае появилась буква Е, разделяющая само число от показателя степени десяти. Так принято писать в других языках программирования и этот формат сделан как раз для совместимости с ними. При использовании команды #pri этот формат записывается как \B. Форматов записи чисел сразу в файл намного меньше, чем форматов записи в текстовый массив. Но иногда они бывают полезными. Если посмотреть на текст программы, то легко заметить, что теперь роль курсора играет переменная j. Она задает индекс реального массива, а результат записывается в последовательные элементы этого массива.
При расчете зависимостей каких-то функций от аргумента бывает полезно не только просматривать таблицы чисел, но и строить графики этих зависимостей, которые показывают их более наглядно. Так как графики имеют очень много параметров, то команды ACL, которые выполняют их рисование также имеют много параметров. График какого-то определенного типа можно записать не одной командой, а серией команд, то есть подпрограммой. Эту подпрограмму можно оформить в виде суперкоманды, у которой уже меньше параметров, а остальные параметры жестко заданы. Использование такой суперкоманды проще, и часто ее вполне хватает. Примером такой суперкоманды является ##smau. Ее название является сокращением от английских слов show memory with argument universally. Ниже показан код программы, которая вычисляет те же зависимости, но показывает их на графике с помощью указанной суперкоманды.
График я показывать не буду, каждый может скопировать текст в окно редактора программы, выполнить его и посмотреть на графики. Описание самой суперкоманды дано в отдельной главе данной книги. Я здесь повторю, что она делает. Так как запись чисел в массив шла таким образом, что все функции в первой точке, затем все -- во второй, и так далее, то матрицу надо транспонировать. Это делает команда #ma [op=mtr;]. После ее выполнения все значения первой функции идут первыми, затем все значения второй и так далее. Для суперкоманды надо задать начальный элемент массива, где записана первая функция в переменную A, затем число точек в переменную B, затем число функций в переменную C. В переменные D и E задаются начальное и конечное значения аргумента. В переменной F задается номер цвета для рисования первой функции, остальные рисуются номерами со сдвигом на 1, наконец переменная G задает режим показа графика и масштабирование. Так значение 1100 означает показывать (первая единица) и масштабировать 100% (младшие разряды, там 100). Режим 100% показывает график так как задано в суперкоманде. Суперкоманда рисует только функции, поэтому начальный индекс массива дан со сдвигом, исключающим строку аргументов.
Наконец, надо задать два текста. Первый -- это имя файла для записи картинки без расширения, расширение всегда png, второй - заголовок графика. Эти тексты разделяются символом вертикальной черты. Вот так, очень просто, можно не только сделать вычисление каких-то зависимостей, но и показать их на графике.
2. Суммирование матрицы и график
В этом примере показана конкретная программа, которая решает задачу преобразования данных. Одна из моих научных программ, представленных онлайн в интернете, показывает в окне текстового редактора матрицу целых чисел, которые получены умножением реальных чисел на 1000. Эти числа можно скопировать в файл anen.txt. Задача состоит в том, чтобы прочитать эти числа, просуммировать их по второму индексу, промасштабировать на реальные значения и затем записать их в файл, а также показать на графике. Ниже я сначала покажу всю программу целиком, а потом объясню как она работает.
Итак, начинаем разбирать программу. Из входных данных расчета следует, что матрица имеет размер 124*76, поэтому сначала определяем эти значения в переменных и вычисляем полное число чисел. Команда во второй строке прочитывает все числа из файла в массив r(), начиная с первого элемента. Команда третьей строки суммирует матрицу по второму индексу. Результат записывается сразу за матрицей. В четвертой строке выполняется нормировка. В пятой строке числа записываются в файл. Наконец, в шестой строке записан код для графика по этим числам. Более подробно об этом можно узнать в описании суперкоманды ##smau.
3. Комплексные вычисления
В языке ACL нет такого типа данных, как комплексные числа. Каждое комплексное число просто считается вектором из двух элементов, сначала действительная часть, потом мнимая часть. Вектор из n комплексных чисел на самом деле представляет собой матрицу (2,n) и к ней применимы все операции, какие работают с матрицами. В частности можно сделать транспонирование #ma [op=mtr; b=1; nx=2; ny=n;] и в результате в том же месте массива сначала будут записаны все реальные части вектора, а потом все мнимые части вектора. Так же точно двумерная матрица комплексных чисел -- это просто трехмерная матрица, но при транспонировании двумерную матрицу можно рассматривать как одномерную.
Дело в том, что двумерных матриц ведь в ACL тоже нет. Есть только один большой массив-вектор r(). Но команды языка ACL сами знают как надо обрабатывать элементы части этого массива, если нужно выполнять заданную операцию. Среди таких команд есть удобные команды, которые часто необходимы. Вот две их них.
#ma [op=cne; b=; n=; nx=; trx=;] и #ma [op=cen; b=; n=; trx=;]
Первая из двух преобразует комплексный вектор из нормальной в экспоненциальную форму (conversion normal to exponent). Параметр [b] всегда указывает начало вектора, [n] указывает число комплексных чисел, [trx] -- смещение индекса от [b] для записи результата. Оно может быть как положительным, так и отрицательным и при этом по модулю больше, чем 2*n, чтобы новый и старый массивы не перекрывались. Экспоненциальная форма комплексного числа -- это модуль и фаза. Но фаза определяется только с точностью до 2*пи*k. Программа устраняет скачки и делает фазу гладкой и непрерывной функцией. Но для ее работы надо указать точку (номер комплексного числа), с которой выполнять выпрямление фазы. Эта точка задается параметром [nx]. Она должна быть больше [b], но меньше [b]+[n]. Если это не так, то ее значение автоматически исправляется на середину интервала.
Вторая операция делает обратное преобразование из экспоненциальной в нормальную форму. Это бывает необходимо, если комплексный вектор сразу задан своей фазой при единичном модуле или модулем и фазой. Такое задание часто встречается в оптике вообще, и в рентгеновской оптике, в частности. В данном случае в параметре [nx] нет необходимости, а все остальные параметры имеют тот же самой смысл.
В операциях расчета свертки комплексных волновых функций приходится умножать (а в обратной операции делить) комплексные векторы поэлементно. Для таких операций тоже есть готовые команды
#ma [op=mcc; b=; nx=;] и #ma [op=dcc; b=; nx=;]
У них совсем мало параметров, потому что они требуют жесткого расположения массивов. Так в первой операции [b] указывает начало комплексного вектора первого множителя, а [nx] число элементов этого вектора. Второй множитель должен быть расположен сразу за первым. Такое расположение программист должен организовать в самом начале написания программы. Комплексные векторы перемножаются поэлементно и результат умножения записывается на место первого вектора. При этом первый вектор изменяет свои значения, а второй -- нет. Я всюду пишу вектор, но это может быть и многомерная матрица, просто в данном случае ее внутренняя структура не имеет значения, и она может рассматриваться как вектор.
Вторая команда обратная первой, но деление производится умножением на комплексно-сопряженное значение и делением на квадрат модуля. При этом, во избежание неприятности с делением на ноль (иначе будет выдано сообщение об ошибке и счет остановится), к квадрату модуля прибавляется маленькое число, которое надо задать в переменной E.
Все указанные операции выполняют внутренние циклы поэлементных вычислений, но в ACL как бы выполняется операция над сложным объектом в виде вектора (или матрицы). И все записывается очень компактно и просто. Есть, однако, еще одна операция, в которой существенно используется векторная природа объекта. Это преобразование Фурье комплексной волновой функции. Хотя сама операция называется FFT, и она действительно использует алгоритм быстрого преобразования Фурье, но реально выполняется преобразование из системы точек в X пространстве в систему точек в Q пространстве и наоборот. При этом X пространство определено как реальное пространство конкретной размерности, а Q пространство является сопряженным и полностью определяется X пространством. Вызов команды очень простой
#ma [op=ftt; b=; nx=;]
Но для ее успешного выполнения необходимо понимать как она работает. Точки в X пространстве должны быть заданы по формуле x[n]=D*(n+(1-N)/2). При этом N -- это полное число точек сетки, которое равно N = 2^[nx]. Имеется в виду модуль параметра [nx]. Реально оно может быть как положительным, так и отрицательным. Знак этого параметра задает знак аргумента экспоненты. А знак зависит от системы, которая используется. Так я лично при преобразовании из X в Q ставлю знак минус, а при обратном преобразовании ставлю знак плюс. Но можно делать и наоборот, только следить за тем, чтобы при прямом и обратном преобразовании были разные знаки. Число D действительно задается в переменной D и оно имеет размерность длины (мкм, см и так далее)
В обратном пространстве Q будет столько же точек, но их координаты будут равны q[n]=d*(n+(1-N)/2), где d = 2п/X, X = D*N, а полная область в q пространстве равна Q = 2п/D. Сам комплексный вектор в нормальной форме должен быть записан начиная с индекса [b]. Ответ получается в том же месте. То есть исходный вектор преобразуется в ответ.
Операции поэлементного сложения и вычитания комплексных векторов нет, так как их можно выполнить рассматривая комплексный вектор как матрицу, то есть реальный объект, а для таких объектов указанные операции есть. Точнее есть операции поэлементного сложения и умножения вектора на число
#ma [op=vva; b=; le=; trx=;] и #ma [op=vmс; b=; le=;]
В первой операции [b] -- начало, [le] -- длина вектора, [trx] -- смещение начала для записи результата. Желательно, чтобы первый вектор (и там же результат) находились отдельно от второго вектора, то есть у них не было общих элементов. Хотя это не является ошибкой и иногда может приводить к интересным эффектам. Чтобы сделать вычитание, надо предварительно второй вектор умножить на -1. Умножение вектора на константу делает вторая операция, причем сама константа должна быть задана в переменной C.
Команды для выполнения двумерного преобразования Фурье от комплексной матрицы в чистом виде тоже нет. На самом деле ее можно заказать как частный случай другой более сложной операции. Но для двумерного преобразования Фурье в принципе можно написать программу на ACL на основе одномерного преобразования. И такая программа не будет считать заметно дольше, потому что главное время расчета все равно отнимает процедура одномерного преобразования, которая делается быстро. Такая программа пока не написана.
Интересно, что такие расчеты выполняются в другой команде #ma [op=mcg;] которая выполняет свертку или развертку реальной матрицы с гауссианом. Такая операция часто используется для уменьшения или увеличения четкости изображения. В физических задачах она используется для учета конечного разрешения детектора и любых факторов, приводящих к сглаживанию изображения, то есть двумерного распределения интенсивности излучения.
Иногда необходимо выполнить достаточно сложные вычисления с комплексными числами непосредственно на ACL. Так как такого типа данных нет, то нельзя просто так писать знаки умножения, деления и так далее. С другой стороны расписывать всю комплексную арифметику через реальные числа очень громоздко. Я решил проблему следующим образом. Вычисления можно организовать по принципу компьютерного кода для машин самого первого поколения. Вместо имен комплексных переменных будем использовать индекс реального массива r(). А операцию будем записывать числом. Команда, которая выполняет такие операции, аналогична команде реальных вычислений и ее аргументом является набор инструкций таких вычислений. На самом деле это тоже операция com команды ma. То есть ее синтаксис такой
#ma [op=com;] K A B C K A B C . . .
где числа разделены на группы по 4, в каждой группе первое число -- это код операции, затем индексы массива r() как первый операнд, второй операнд и результат операции. Рассмотрим пример, нам надо выполнить такое выражение с комплексными числами a=a+b*exp(c)/(d+e). Вводим еще два числа f и g и все числа размещаем в реальном массиве так что a=[r(1),r(2)], b=[r(3),r(4)], . . ., g=[r(13),r(14)]. Тогда выражение можно вычислить следующей системой команд #ma [op=com;] 5 5 0 11 3 11 3 11 1 7 9 13 4 11 13 11 1 1 11 1 . Хотя такая запись выглядит достаточно абстрактной, но она очень компактная и позволяет быстро и просто выполнить комплексную арифметику. Реализовано 12 операций, которые имеют следующие коды: 1 -- сложение, 2 -- вычитание, 3 -- умножение, 4 -- деление, 5 -- exp(), 6 -- cos(), 7 -- sin(), 8 -- tan(), 9 -- pow(,), 10 -- log(), 11 -- J0(A), 12 -- J1(A). Для функций второй аргумент не используется. Исключением является функция pow(), для которой используется реальная часть второго аргумента, то есть комплексное число возводится только в реальную степень. J0 и J1 -- это функции Бесселя нулевого и первого порядков. Есть и альтернативная форма такой записи. Более подробную информацию можно получить при чтении описания команды #ma.
4. Вычисление точки прохождения функции через ноль
Рассмотрим пример решения так называемой обратной задачи. Она состоит в том, что нужно определить значение аргумента, при котором функция равна нулю. Мы умеем вычислять функцию, она может быть любой, но нам надо узнать обратную зависимость, точнее просто значение аргумента, при котором функция имеет заданное значение, равное нулю. Это возможно только методом последовательных приближений. Решение такой задачи разумно делать в виде процедуры с входными данными. В частности надо иметь процедуру вычисления функции, пусть она имеет название (func).
Она вычисляет значение функции в переменную f по значению аргумента в переменной x. И еще нужно указать границы интервала, на котором мы хотим искать решение. Пусть они заданы в переменных xa (первое значение) и xb (последнее значение). Наконец, надо указать точность, с которой мы хотим знать ответ. Пусть она задана в переменной e. Тогда код процедуры может быть таким
Осталось объяснить как это работает. Мы используем 3 разные точки на оси аргумента x1, x2, x3 и три значения функции f1, f2, f3. Сначала вычисляем значение функции в первых двух точках, которые равны границам интервала. В первой точке значение запоминаем также в переменной f0. Затем открываем цикл, который работает при значении переменной &, равном 1. Аппроксимируем функцию по двум точкам линейной фукнцией f = Ax +B и находим такое значение x3, где эта функция равна нулю. Вычисляем истинное значение f3 функции в этой точке. Затем первую точку заменяем на вторую, а вторую на третью.
И проверяем условие, что f3 по модулю меньше, чем e*f0. Если это условие выполняется, то & будет меньше, чем 1, и цикл прервется. Если нет, то & будет снова равен 1 и цикл повторится снова, но аппроксимация функции будет происходить по другим точкам, более близким к решению. Поэтому рано или поздно условие все же нарушится, если такая точка существует. В данной процедуре не предусмотрен обратный случай. Это легко исправить, если ввести ограничение на число попыток и прерывать цикл по истечении этого числа. Как это сделать -- домашнее задание для читателя.