книги / Математическое моделирование дискретные подходы и численные методы

Чтобы увидеть разделение доступных веществ по классам, применяется команда “ChemicalData["Classes"]”. Далее имя класса можно использовать для получения списка веществ этого класса, например, выполнив команду “ChemicalData["Nanomaterials"]”, получим список химических соединений, молекулы которых используются в нанотехнологиях как наночастицы. Используя их названия, можно узнать их химическую формулу, ее структурное представление, физические свойства. Например, для получения химической формулы надо выполнить команду

“ChemicalData["YttriumIronOxide", "HillFormulaDisplay"]”,

в результате получим Fe5O12Y3. Структурная формула этой молекулы имеет вид, получаемый с помощью функции “ChemicalData["YttriumIronOxide"]” (рис. 1.15, а). Аналогично можно получить формулу молекулы кофеина:

“ChemicalData["Caffeine", "HillFormulaDisplay"]”,

как в виде C8H10N4O2 или структурной формулы (рис. 1.15, б) с помощью функции “ChemicalData["Caffeine"]”, так и в виде трехмерного изображения молекулы(рис. 1.15, в) спомощью функциисдополнительнымпараметром

“ChemicalData["Caffeine", "MoleculePlot"]”.

Переменные и функции. Для создания символьной переменной выполняется команда, например, “a=5”, которая автоматически создает в ядре пакета переменную a с текущим значением 5, и во всех выражениях, содержащих a, которые будут выполнены позже, вместо a будет подставлено число “Integer[5]”. То есть появляется глобальная переменная с заданным значением. Отметим, что после выполнения команды присваивания переменная a стала изображаться черным цветом, а не синим, как до выполнения команды. Для того чтобы очистить переменную a, например, для ее использования в символьных преобразованиях, выполняется команда “Clear[a]” или “a=.”. После этого цвет переменной a на экране снова становится синим.

Пример: a = 63.84 ---> 63.84

b = a2 + a – 1 ---> 4138.39

Последовательность символов “--->” будем использовать, чтобы показать результат выполнения той или иной команды в пакете. Итак, b сохранилось в ядре не как выражение “a2 + a – 1”, а как действительное число 4138.39. Если изменить a, то b не изменится:

81

a = 19 ---> 19

b ---> 4138.39

Выясните, что пакет «знает» о величине b: “? b”.

В «Mathematica» может использоваться несколько форм выполнения некоторой операции. В частности, для присваивания может использоваться функциональная форма “Set[]” – “Set[a,19]”.

Если придать какой-либо переменной новое значение, а затем выполнить записанные ранее выражения с этой переменной, то будут получены другие ответы. Часто бывает необходимо изменить значение переменной с использованием ее предыдущего значения, например “a = a + 7”. Для упрощения аналогичных команд в пакете реализованы операции: “+=”, “–=”, “*=”, “/=”. Если требуется изменять значение переменной на единицу (например, при реализации циклов), то можно использовать удобные короткие формы: “++” и “––”, которые используются как “i++”, “i––”, “++i”, “––i”.

Выполните приведенную последовательность команд и объясните, как получаются такие результаты.

«Mathematica» позволяет вычислять выражении без присваивания переменной некоторого значения. Для этого используются так называемые подстановки. При этом в конец выражения, которое требуется вычислить, дописывается последовательность символов вида “/.x→ 5”, в которой вместо числа 5 может задаваться любое значение переменной x, в том числе и другая переменная. Набор символов “x→ 5”, идущий после “/.”, называется правилом замены переменной или правилом подстановки. При преобразованиях выражений последовательно можно применять любое число правил замены. Для ввода стрелки “→ ” используется последовательность символов “–” и “>”.

Если в вычислениях возникает необходимость многократного использования некоторой последовательности операций – это хороший повод создать собственную функцию. Синтаксис задания функции требует указания ее аргумента, который определяется знаком подчеркивания после имени переменной, и задания тела функции, например: “myFunc[x_]: =x2+x+1”. Имя собственной функции пользователя может начинаться с любой буквы,

82

не обязательно заглавной. Подчеркивание указывает, что переменная с именем x является локальной, и при вызове функции на ее месте может оказаться любое выражение или число. При вводе приведенной команды видно, что пакет «Mathematica» автоматически изображает переменные курсивом зеленого цвета. Операция “:=” называется отложенным присваиванием и означает, что правая часть будет вычисляться только при вызове функции, то есть тогда, когда функции f[x] будет передано конкретное выражение для x. А обычное присваивание “=” по аналогии можно назвать немедленным присваиванием. Итак, в пакете мы встретили уже три различных «знака ра-

венства»: “==”, “=”, “: =”.

Отметим, что если в конце текста команды ставится символ “;”, то выходная ячейка не выводится. Это удобно делать при выполнении присваивания, поскольку нет необходимости еще раз видеть на экране то выражение, которое присваивается переменной. Также символ “;” позволяет разделить несколько команд и записать их в одну строчку. При задании функции с помощью отложенного присваивания на экран также ничего не выводится.

Рассмотрим обычные ошибки, допускаемые при задании функций.

Пример: bad[x]: = x2 + x – 1 bad[x] ---> – 1 + x + x2

bad[y] ---> bad[y]

То есть пакет «не знает», что такое “bad[y]”.

Внекоторых случаях при задании функции может использоваться

изнак немедленного присваивания “=”, но с ним надо обходиться осторожно. Приведем примеры.

Пример: x = 5; bad2[x_] = x2 + x – 1;

bad[y] ---> 29

Здесь при инициализации функции произошло вычисление ее правой части x2 + x – 1 при заданном значении x, то есть была введена функция с постоянным значением 29. В таком случае использование немедленного присваивания “=” приводит к ошибке.

Рассмотрим другой пример. Для нахождения производной функции по некоторой переменной используется встроенная функция “D[]”, например “D[x2,x]”. Создадим функцию, которая будет возвращать значение производной выражения x2 + x – 1 при различных значениях переменной x.

83

Пример: bad3[x_]: = D[x2 + x – 1, x]

bad3[5] ---> 5 is not a valid variable.

Произошло не то, что ожидалось. Из-за отложенного присваивания в правую часть тела функции D[x2 + x – 1, x] вместо x

не существует! В этом примере при задании функции необходимо использовать немедленное присваивание, чтобы сначала была найдена производная, а потом функция возвращала ее значения при различных значениях переменной x:

x =.; g[x_] = D[x2 + x – 1, x]; g[5] ---> 11

Помимо функций, задаваемых выше с помощью явного указания имен функций и переменных и поэтому называемых поименованными функциями, пакет «Mathematica» позволяет работать с так называемыми анонимными функциями. Они удобны при однократном применении некоторого выражения. Вместо имени функции при использовании анонимной функции указывается конструкция, оканчивающаяся знаком “&”. Поскольку эта конструкция заменяет имя функции, то после нее можно ставить квадратные скобки с указанием аргумента (или нескольких аргументов). Конструкция, о которой идет речь представляет собой тело функции, в котором вместо имени переменной используется символ “#” для функции одного аргумента или набор вида “#1”, “#2” и так далее для функции нескольких аргументов. Например, для вычисления значения выражения x2 + x – 1 может быть создана анонимная функция “ (#2 + # – 1) &”. Ее применение выглядит как “ (#2 + # – 1) &[x]”. Если в теле функции используются логические операции, то результатом будет оценка истинности получаемого высказывания, например, “ (#12 > #22) &” будет сравнивать квадраты двух передаваемых функций аргументов. Для вычисления производной с помощью анонимной функции необходимо заставить пакет предварительно выполнить символьное преобразование производной. Сделать это можно с помощью функции “Evaluate[]”. Например, в ситуации, аналогичной рассмотренной в приведенном выше примере, запись “Evaluate[D[#2 + # – 1, #]]&” дает корректный результат.

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

84

ваются элементами второго уровня, а список – многоуровневым. Заголовок этого объекта называется “List”, и задавать список можно функцией “List[]”, перечисляя в качестве ее аргументов все элементы списка

“List[a,b,c,3,–1,x]”.

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

“atomChain={{1,0,0},{2,0,0},{3,0,0},{4,0,0}}”.

Функция “Length[]” возвращает длину списка, что в нашем случае соответствует количеству атомов:

“ nAtoms=Length[atomChain]”.

Для обращения к элементу списка используется имя списка и номер элемента, который указывается в двойных квадратных скобках, например “atomChain[[3]]”. Если в двойных квадратных скобках указано отрицательное число, например “atomChain[[–1]]”, то номер элемента будет отсчитываться с конца списка (в приведенном примере будет взят последний элемент). Для выбора случайного атома необходимо получить случайное число от 1 до nAtoms. Случайные числа в пакете генерирует функция “Random[]”. Ее первый аргумент указывает тип возвращаемого числа, второй – интервал, из которого выбирается число. В нашем случае нужны целые случайные числа:

“Random[Integer,{1, nAtoms}]”. “atomChain[[Random[Integer,{1, nAtoms}] ]]”.

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

Ранее список задавался перечислением его элементов. В «Mathematica» существуют функции для генерации списков. Одним из наиболее универсальных способов создания списков является использование функ-

85

ции “Table[expr,iterator]”. Первым аргументом стоит выражение expr, которое в результате повторных исполнений будет задавать элементы списка, второй аргумент – итератор, который определяет характер повторения expr. Итератор всегда задается в фигурных скобках и может иметь следующий вид: “{n}” – выражение expr будет повторено n раз (если оно содержит функцию, генерирующую случайные числа, то результаты исполнения будут различными), “{i,n}” – итератор содержит переменную итерирования i, которая принимает значения от 1 до n, “{i,n1,n2}” – переменная i принимает значения от n1 до n2, “{i,n1,n2,st}” – переменная i принимает значения от n1 до n2 с шагом st. Для задания цепочки атомов заданного их количества n выполним

“Table[{i,0,0},{i,n}]”,

предварительно задав значение n, например, “n=20”. Итераторов может быть несколько – для одного итератора получаемый список соответствует некоторому условному вектору, для двух итераторов – плоской матрице, для трех итераторов будет получаться трехмерный массив и так далее.

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

Функция “Join[]” соединяет списки, не выполняя упрощений. Функция “Union[]” объединяет любое количество списков в теоретико-множествен- ном смысле, удаляя повторяющиеся элементы и сортируя оставшиеся элементы списка в алфавитном порядке (применение к одному списку приводит к удалению повторов и сортировке его элементов). Функция, которая задает теоретико-множественное пересечение списков, имеет название и синтаксис применения вида “Complement[list,list1,lit2,…]”, где из списка list удаляются все элементы, встречающиеся в списках list1, lit2 и так далее.

Функции “Append[list,el]” и “Prepend[list,el]” добавляют элемент el к списку list. Первая добавляет новый элемент в конец списка, вторая – в начало. Функция “Delete[list,k]” удаляет элемент с номером k из списка list. Если число k является отрицательным, то отсчет номера удаляемого элемента ведется с конца списка. Функция “Flatten[list]” удаляет все внутренние фигурные скобки, то есть делает многоуровневый список одноуровневым. Дополнительный аргумент задает номер уровня, до которого необходимо проводить упрощение списка: “Flatten[list,L]”. Напри-

мер, генерируя с помощью функции “Table[a{i,j,k},{i,n},{j,n},{k,n}]” на-

86

бор атомов с межатомным расстоянием a, соответствующий простой кубической решетке, получим список сложной структуры:

“n=2; atomList = Table[a{i,j,k},{i,n},{j,n},{k,n}]” ---> “{{{{a,a,a},{a,a,2a}},{{a,2a,a},{a,2a,2a}}},

{{{2a,a,a},{2a,a,2a}},{{2a,2a,a},{2a,2a,2a}}}}”.

Применениефункции“Flatten[atomList]” приводитктакомурезультату:

“{a,a,a,a,a,2a,a,2a,a,a,2a,2a,2a,a,a,2a,a,2a,2a,2a,a,2a,2a,2a}”,

когда структура списка полностью утрачивается. Применение функции с дополнительным аргументом “Flatten[atomList,1]” приводит к списку

“{{{a,a,a},{a,a,2a}},{{a,2a,a},{a,2a,2a}},

{{2a,a,a},{2a,a,2a}},{{2a,2a,a},{2a,2a,2a}}}”,

в котором наборы координат атомов сгруппированы по парам. Применение “Flatten[atomList,2]” приводит к двухуровневому списку

“{{a,a,a},{a,a,2a},{a,2a,a},{a,2a,2a},

{2a,a,a},{2a,a,2a},{2a,2a,a},{2a,2a,2a}}”,

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

“n=4; atomList = Flatten[Table[a{i,j,k},{i,n},{j,n},{k,n}],2]”.

Для работы с векторами и матрицами в пакете «Mathematica» реализованы операции скалярного умножения – десятичная точка “.”: “a.b”, где списки a и b могут быть векторами одной длины или матрицами с соответствующими размерностями. Функциональная форма этой операции – функция “Dot[a,b]”. Векторное произведение двух векторов из двумерного или трехмерного пространства задается с помощью команды “a×b”, где символ “×” вводится с помощью последовательности “#Cross#” (“#” появляется после нажатия клавиши Esc). Соответствующая функция– “Cross[a,b]”. Для квадратных матриц может быть найден определитель: “Det[A]”, обратная матрица: “Inverse[A]”, транспонированная матрица любой размерности получается при выполнении команды: “Transpose[A]”. В ряде преобразований надо задавать единичную матрицу, для этой цели в пакете предусмотрена функция “IdentityMatrix[n]”, где n – ее размерность.

87

Матрица компонент тензора поворота задается с помощью функции “RotationMatrix[θ]” в двумерном случае, где θ – угол поворота вокруг оси, ортогональной плоскости, либо с помощью “RotationMatrix[θ,e]”, где переменная e задает ось поворота (здесь это не обязательно единичный вектор, при вызове указанной функции этот вектор автоматически нормируется). Для получения системы собственных чисел тензора второго ранга используется функция “Eigenvalues[A]”, для получения системы собственных векторов, соответствующих этим числам предусмотрена функция “Eigenvectors[A]”. полная система собственных чисел и векторов находится с помощью функции “Eigensystem[]”.

Для задания тензорного умножения применяется частный вариант функции “Outer[f,list1,list2,…]”, которая возвращает список значений функции f нескольких аргументов на всех возможных комбинациях аргументов, когда первый берется из первого списка list1, второй аргумент – из второго списка list2 и так далее. Например, выполним команду

“Outer[f,{a,b,c},{x1,x2}]” --->

“{{f[a, x1], f[a, x2]}, {f[b, x1], f[b, x2]}, {f[c, x1], f[c, x2]}}”.

Здесь для получения конкретного результата вычислений функция f должна быть задана заранее. В частности, для получения тензорного произведения двух и более векторов в качестве f надо указать имя функции “Times”, реализующей умножение своих аргументов. Функция “Outer[f,list1,list2,L]” имеет дополнительный аргумент L, задающий для многоуровневых списков уровень, элементы которого воспринимаются как единое целое и подставляются целиком в аргументы функции f.

Какой операции тензорного анализа соответствует команда “Outer[D,v,{x,y,z}]”, где v представляет собой векторное поле?

Для нахождения суммы, произведения элементов списка или результата применения какой-либо другой функции ко всему списку сразу используется функция “Apply[f,list]” или в компактной форме “f@@list”. Эта функция заменяет заголовок “List” списка на имя f. То есть для нахождения суммы элементов списка {a,b,2a,–1,c,3,b} выполняется команда

“Plus@@{a,b,2a,–1,c,3,b}”.

При этом происходит замена заголовка списка и вычисляется значение функции вида “Plus[a,b,2a,–1,c,3,b]”, что приводит к результату “2+3a+2b+c”. Для получения произведения элементов надо применить имя функции “Times”.

88

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

“cMass=Plus@@ atomList/Length[atomList]”.

Для смещения массива атомов в пространстве так, чтобы его центр масс совпал с началом координат, необходимо от радиус-вектора каждого атома отнять вектор, задающий положение центра масс. Команда “atomList – cMass”, очевидно не даст желаемого результата. В данном случае к каждому элементу списка atomList необходимо применить анонимную функцию смещения “ (# – cMass) &”. Применение некоторой функции f к каждому элементу списка list реализуется с помощью команды “Map[f,list]”. В нашем примере:

“Map[ (# – cMass) &, atomList]”.

Короткий вариант применения функции “Map[]” задается симво-

лами “/@”:

“ (# – cMass) & /@ atomList”.

Функция “Map[f,list]” также может использоваться для визуализации сгенерированного кристаллического образца. Действительно, команда

“Graphics3D[{Specularity[White,80],ColorData["Atoms","Lu"], Sphere[#, 0.35]} & /@ (atomList /. a → 1), Lighting → "Neutral"]”

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

Часто возникают задачи выбора элементов из списка по некоторому критерию. Например, выясним, как выбрать атомы, принадлежащие одной из граней кубического образца с простой кубической решеткой (их координаты содержатся в переменной atomList). Для выбора элементов списка list по заданному критерию crit применяется функция “Select[list, crit]”. Критерий представляет собой имя функции одной переменной. Выбираются те элементы, для которых вычисленная функция критерия возвращает значение “True”. Рассмотрим грань, ортогональную оси Ox1 для набора атомов atomList, центр масс которого еще не смещен в начало координат. На ребре образца находится n атомов, следовательно, абсцисса атомов нужной грани равна a n, тогда команда

89

“rightFace = Select[atomList, #[[1]] == n a&]”

даст требуемый список атомов одной грани (рис. 1.16, б). Получить рисунок с выделенной боковой гранью с использованием функции “Map”.

Рис. 1.16. Образец с простой кубической решеткой:

а – структура кристалла; б – цветом выделена одна боковая грань

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

“Complement[atomList, rightFace]”.

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

rightFace, Complement[atomList, rightFace],1],1]]”.

Элементы программирования. Рассмотрим основные подходы к программированию, реализованные в пакете «Mathematica». Каждый из этих подходов определяет лишь структуру кода, но не дает дополнительных вычислительных возможностей. В различных случаях выбор подхода может

90