сходимость, не должны существенно зависеть от дискретизации, что было также подтверждено эмпирически [51].

Будко и Самохин [148] обобщили результат Рахолы на произвольные неоднородные и анизотропные рассеиватели, описав область на комплексной плоскости, содержащую весь спектр интегрального оператора рассеяния. Эта область полностью определяется значениями, которые принимает m внутри частицы, и не зависит от x. Они показали, что для вещественных m или m, имеющих очень малую мнимую часть, эта область проходит близко от начала координат, поэтому спектр может содержать очень малые собственные значения для частиц больше длины волны, что приводит к большому числу обусловленности матрицы A и медленной сходимости. В подтверждение этих теоретических выводов, в недавней работе Айранчи (Ayranci) и др. [149] было показано, что Niter действительно уменьшается с увеличением Im(m),

хотя авторы ограничились размерами x < 8. Анализируя спектр интегрального оператора для частиц много меньших длины волны, Будко и др. [150] предложили эффективный итерационный метод для этого случая.

В заключение, существует несколько современных итерационных методов [КМН(КС), Би-СГ(КС) и Би-СГСтаб], эффективность которых была доказана применительно к МДД. Но невозможно установить однозначное превосходство одного из этих методов над другими – необходимо сравнивать их применительно к конкретной задаче. Более того, предобусловливание матрицы взаимодействия в МДД практически не изучено, за исключением простейших вариантов, в то время как оно явно необходимо для больших x и m, когда все существующие методы сходятся очень медленно или расходятся. Нам кажется, что следующий основной прорыв в численных аспектах МДД будет достигнут именно за счёт эффективного предобуславливателя.

Большое количество диполей требует больших вычислительных ресурсов, что достигается обычно использованием параллельных компьютеров [35]. В данной диссертации не обсуждается параллельная эффективность, но для итерационных методов она обычно близка к единице [34]. Однако это не всегда так для предобуславливателей, что следует принимать во внимание при использовании «тяжёлых» предобуславливателей, требующих существенное время, в параллельной реализации МДД.

2.1.4.2. Разложение по порядкам рассеяния

Приближение РДГ [14] состоит в предположении, что E(r) равно Einc(r). F(n) тогда получается напрямую из формулы (26). Обобщением РДГ является итерационное

решение интегрального уравнения (2), которое может быть переписано в виде

49

E(r) = Einc (r) + ΛE(r) ,

(80)

где Λ это линейный интегральный оператор, описывающий рассеиватель. Итерационная схема получается путём подстановки текущей (l-ой) итерации электрического поля E(l)(r) в правую часть уравнения (80) и вычисления следующей итерации в левой части:

E(l+1) (r) = Einc (r) + ΛE(l ) (r) .

(81)

Начальное значение берётся таким же как и в РДГ, E(0)(r) = Einc(r), что приводит к следующей общей формуле для решения:

			∞
E(r) = ∑Λl Einc (r) ,								(82)
			l=0
что является прямым следствием известного ряда Неймана:
							∞
(I − Λ)−1 = ∑Λl ,								(83)
							l=0
где I это единичный оператор. Необходимым и достаточным условием сходимости
ряда Неймана является
			Λ				<1.	(84)

Физический смысл данного метода состоит в последовательном вычислении взаимодействия между различными частями рассеивателя. Нулевой порядок (или РДГ) совсем не учитывает взаимодействие, первый порядок учитывает один раз влияние рассеянного поля каждого диполя на остальных, и т.д. Условие (84) утверждает, что взаимодействие частей рассеивателя между собой должно быть малó, но не настолько малó, как требуется для применимости РДГ ( Λ <<1). В общих задачах рассеяния,

особенно в квантовой физике, формула (82) называется разложением Борна.

Хотя теоретически это разложение простое, оно не может быть непосредственно применено [151], так как каждая последующая итерация требует аналитического вычисления многомерных интегралов всё возрастающей сложности, что становится невыполнимо даже для рассеивателей простейшей формы. Последним достижением в этом направлении является, вероятно, работа Аквиста (Acquista) [151], который вычислил разложение Борна для однородного шара вплоть до второго порядка. Поэтому практическое применение разложения Борна требует дискретизации интегрального оператора, на чём собственно и основан МДД.

Разложение по порядкам рассеяния (РПР) в МДД было независимо разработано Чиапеттой (Chiappetta) [152] и Сингхэмом и Бореном [54,153], применяя ряд Неймана к уравнению (19). В этом случае Λ является матрицей Λij = Gij αj , каждый элемент

50

Gij = G ji = Gi′−j .

которой это тензор, выраженный матрицей 3×3. Явная проверка условия (84) для конкретного рассеивателя численно очень затруднительна, однако де Хоп (de Hoop) [154] вывел достаточное условие для скалярных волн:

2π(kR0 )2 max		χ(r)		<1 ,	(85)
r

где R0 это наименьший радиус описанной вокруг рассеивателя сферы. Хотя она, строго говоря, не применима к светорассеянию, можно использовать формулу (85) в качестве оценки.

РПР сходится в ограниченном диапазоне x и m [153], и даже там, где оно сходится, более современные итерационные методы сходятся быстрее (см. параграф 2.1.4.1). Но РПР обладает ясным физическим смыслом и может использоваться для изучения многократного рассеяния.

2.1.4.3. Блочно-топлицева структура

Квадратная матрица A называется топлицевой (Toeplitz), если Aij = ai − j, т.е. элементы матрицы не меняются вдоль любой линии параллельной главной диагонали [140]. В блочно-топлицевой (БТ) матрице порядка K в качестве элементов ai выступают не числа, а другие квадратные матрицы:

a0		a1	K aK −1
	a−1	a0	O	M
A =					.	(86)
	M	O	O	a1
		K a−1		a0
a−K +1

Двухуровневая БТ матрица имеет БТ матрицы в качестве компонентов ai. Продолжая рекурсивно, определяется многоуровневая БТ (МБТ) матрица.

Рассмотрим прямоугольную решётку nx×ny×nz, пронумерованную как

i = ny nz (nx −1)ix + nz (ny −1)iy + nziz ,

(87)

где iμ {1,...,nμ} обозначают декартовы координаты элемента решётки. Определим векторный индекс i = (ix,iy,iz), тогда легко проверить, что матрица взаимодействия из уравнения (22), определённая формулой (15), удовлетворяет следующему условию:

(88)

Само по себе это условие позволяет сильно уменьшить используемую итерационными методами память за счёт косвенной адресации. Дополнительные упрощения связаны с тем, что формула (88) определяет симметричную трёхуровневую БТ матрицу (порядки уровней: nx, ny и nz), чьими базовыми элементами являются 3×3 матрицы (тензора) Gij .

51

С одной стороны, прямоугольная решётка не является большим ограничением, так как любой рассеиватель может быть помещён внутрь подходящей решётки. С другой стороны, при этом добавляются «пустые» диполи, чтобы достроить решётку до полного параллелепипеда. Более того, положение и размер диполей не может выбираться произвольно, например, для более точного описания формы рассеивателя, что особенно актуально для очень пористых частиц или кластеров многих частиц, в которых размер мономера порядка одного диполя. Для большинства задач эти ограничения малы по сравнению с большим ускорением вычислений, основанном на БТ структуре. Произведение матрицы на вектор в этом случае преобразуется к виду свёртки, которая вычисляется с помощью быстрого преобразования Фурье (БПФ) за время порядка O(nlnn) (см. параграф 2.1.4.4). Следует отметить, однако, что существуют альтернативные подходы, не требующие регулярную решётку (см. параграф 2.1.4.5).

БТ структура также может быть использована для ускорения прямых методов. Флато и др. [155] применили алгоритм для обращения симметричной БТ матрицы, имеющий сложность O(n3/nx) и требующий O(n2/nx) памяти, поскольку хранятся только два блочных столбца обратной матрицы. При этом ось x выбирается вдоль самого большого размера частицы. Недавно Флато [156] рассмотрел частный случай одномерного МДД, когда все диполи эквидистантно расположены на прямой, при этом системы уравнений для каждой декартовой компоненты разделяются. В этом случае матрица взаимодействия для каждой компоненты является симметричной топлицевой и может быть обращена с использованием быстрого алгоритма, который сперва решает систему уравнений для двух правых частей (например, каким-нибудь итерационным методом), после чего умножение обратной матрицы на любой вектор (т.е. решение системы для любой правой части) требует лишь O(nlnn) операций. Но там же Флато отметил существенное ограничение для всех методов быстрого обращения матрицы: они применимы только, если поляризуемости всех диполей одинаковы, так как иначе первый член в правой части уравнения (22) разрушает БТ структуру на диагонали матрицы. Тем самым, все эти методы применимы лишь к однородным прямоугольным рассеивателям. К счастью, это не затрудняет умножение матрицы на вектор, так как можно вычислить диагональные члены независимо и добавить к конечному результату.

2.1.4.4. Быстрое преобразование Фурье

Гудмен и др. [157] показали, что умножение матрицы взаимодействия для прямоугольной решётки (см. параграф 2.1.4.3) на вектор может быть преобразовано к дискретной свёртке

52