моделирования светорассеяния частицей с большим отношением поперечных размеров, при этом оставаясь совместимым с методиками на основе БПФ (параграф 2.1.4.4).
Хлебцов [129] предложил упрощённый вариант МДД, основанный на предположении, что все вектора поляризации параллельны падающему электрическому полю. Таким образом количество неизвестных уменьшается в три раза, но страдает точность, в частности, пропадают все деполяризационные эффекты.
Маркел [130] аналитически решил уравнения МДД для случая бесконечной одномерной периодической цепочки диполей, что похоже на подход, используемый при выводе ДСР [46].
Шомэ и др. [131] обобщили МДД на периодические структуры и далее на дефекты в периодической структуре на поверхности [132]. Идея использования сложной функции Грина, соответствующей некой внешней геометрии, в стандартном МДД была описана в обзоре Мартина [133].
Янг (Yang) и др. [134] вычисляли с помощью МДД поверхностное электрическое поле и определяли интенсивности комбинационного рассеяния света малыми металлическим частицами произвольной формы.
Лемар и Бассрей (Bassrei) [135] показали, что можно восстановить форму объекта по измеренной индикатрисе рассеяния. Эта процедура основана на обращении зависимости между поляризуемостями диполей и рассеянным полем, которая берётся из формулировки МДД. Похожая идея используется в современных работах по оптической томографии [136-138].
Зубко и др. [139] пробовали различные изменения в тензоре Грина при МДД моделировании светорассеяния назад частицами неправильной формы, и показали, что дальнодействующая часть тензора Грина ответственна как за пик яркости при рассеянии назад, так и за отрицательную ветвь поляризации.
2.1.4.Численные соображения
Вэтом подразделе обсуждаются численные аспекты МДД. При этом надо иметь ввиду, что итоговое время вычислений зависит не только от численного метода, но и от его конкретной реализации в компьютерной программе. В частности, в разделе 2.5 сравниваются четыре программы на основе МДД, которые используют практически одинаковые численные методы, но их скорость вычисления различается в разы. Оптимизация компьютерных программ рассматривается в разделах 2.4 и 2.5.
45
2.1.4.1. Прямые и итерационные методы
Есть два общих класса методов решения систем линейных уравнений Ax = y, где x это неизвестный вектор, а A и y – известные матрица и вектор соответственно: прямые и итерационные [140]. Прямые методы приводят к результату за определённое количество шагов, в то время как количество итераций Niter, требуемое итерационными методами, в общем случае априори не известно. Типичным примером прямого метода является LU разложение, которое позволяет быстро решать систему для многих векторов y, как только проведено разложение. Обычно итерационные методы быстрее, более численно устойчивы и требуют меньше памяти, но нельзя утверждать, что итерационные методы заведомо превосходят прямые, так как эффективность первых сильно зависит от конкретной проблемы [141].
Для произвольной n×n матрицы (в МДД n = 3N) LU разложение требует время O(n3) и памяти O(n2), в то время как время одной итерации O(n2) [141]. В общем случае итерационные методы сходятся за O(n) итераций, а иногда могут и не сойтись совсем. Но обычно удовлетворительная точность достигается уже после много меньшего чем n количества итераций, тогда итерационные методы становятся существенно быстрее прямых, особенно при больших n. Многие итерационные методы используют матрицу A только для вычисления произведения матрицы на вектор (иногда также с транспонированной матрицей), что позволяет использовать специальные методики вычисления этих произведений. Такие методики уменьшают требуемый объём памяти, поскольку нет необходимости хранить всю матрицу, особенно для матриц специального вида (см. параграф 2.1.4.3). Специальная структура матрицы позволяет также ускорить вычисление произведения матрицы на вектор с O(n2) до O(nlnn) (см. параграфы 2.1.4.4 и 2.1.4.5). Подобное ускорение, однако, имеет место и для прямых методов (см. параграф 2.1.4.3).
Применительно к МДД обычно использовались итерационные методы (за исключением работ, описанных в параграфе 2.1.4.6) – сначала они использовались для ускорения вычислений [44], а потом позволили использовать больше диполей [35,142], поскольку необходимость хранения всей матрицы A в памяти ограничивает прямые методы. В МДД широко используются методы на основе пространства Крылова, такие как [141] метод сопряжённых градиентов (СГ), СГ в применении к нормализованному уравнению с минимизацией нормы ошибки или невязки (СГНО и СГНH соответственно, CGNE, CGNR), Би-СГ, Би-СГ стабилизированный (Би-СГСтаб, BiCGSTAB), СГ в квадрате (СГК, CG squared), обобщённый метод минимальной невязки (ОММН, GMRES), метод квазиминимальной невязки (КМН, QMR), КМН без
46
транспонирования (КМНБТ, transpose free QMR), и обобщённые методы производного типа основанные на Би-СГ (ОПБи-СГ, GPBi-CG) [143].
Важной составляющей итерационного метода является предобусловливание, которое эффективно уменьшает число обусловленности матрицы A и тем самым ускоряет сходимость. Это, однако, требует дополнительных вычислений и при инициализации итерационного метода, и во время каждой итерации. Предобусловливание исходной системы может быть описано [141] как
M1AM−21 (M2x) = M1y , (77)
где M1 и M2 называются соответственно левым и правым предобуславливателем. Они должны быть либо быстро обращаемыми, либо встроенными в итерационный процесс. Простейший пример первого типа это предобуславливатель Якоби (точечный), что есть просто диагональная часть матрицы A, а примером второго типа является полиномиальный предобуславливатель Неймана порядка l:
l |
|
M = ∑(I − A) j . |
(78) |
j=1
Как КМН, так и Би-СГ могут использовать комплексную симметричность (КС) матрицы взаимодействия в МДД для сокращения в два раза количества произведений матрицы на вектор [144] (и следовательно времени вычислений). Лумме и Рахола [77]
первыми применили КМН(КС) в МДД и сравнили его с СГНH. При m от 1.6 + 0.1i до
3 + 4i и x от 1.3 до 13.5, что соответствует N от 136 до 20336, КМН(КС) был в 2-4 раза быстрее СГНH.
Рахола [51] далее исследовал КМН(КС) и сравнил его с СГНH, Би-СГ(КС), БиСГСтаб, СГК и ОММН (как полный, так и с ограниченной размерностью рабочего подпространства). Сходимость методов изучалась для “типичной небольшой задачи” (к сожалению, авторы не указали конкретные параметры), и лучшие результаты показали КМН(КС) и Би-СГ(КС). Хотя полный ОММН требовал меньше итераций, ограничение размерности рабочего подпространства хотя бы до 40 делало его медленнее КМН(КС).
Флато [145] описал использование различных итерационных методов в МДД и проверил многие из них в комбинации с несколькими предобуславливателями,
моделируя светорассеяния однородным шаром с x = 0.1 и m от 1.33 до 5 + 0.0001i, а
также с x = 1 и m от 1.33 до 1.33 + i и 3 + 0.0001i. Использовались левый (Л) и правый (П) предобуславливатель Якоби, а также полиномиальный предобуславливатель Неймана первого порядка. К сожалению, не было указано количество диполей N, что затрудняет сравнение с другими работами. Для малых частиц СГ(Л) показал наилучшие
результаты для всех показателей преломления, СГ и СГ(П) были к нему близки, а 47
СГНH(Л) и Би-СГСтаб(Л) были медленнее примерно в 4 раза. Для x = 1 наилучшие результаты показал Би-СГСтаб(Л), немного хуже были Би-СГСтаб(П) и СГК(Л),(П), а КМНБТ (и с предобуславливателями Якоби, и без) был медленнее в 3-4 раза. В целом предобуславливатель Неймана первого порядка показал себя неудовлетворительно, а самым подходящим методом для МДД был признан Би-СГСтаб(Л).
Недавно Фэн (Fan) и др. [146] сравнили ОММН, КМН(КС), Би-СГСтаб, ОПБи-СГ и Би-СГ(КС) применительно к рассеивателям размером порядка длины волны (x до 10)
с m до 4.5 + 0.2i и заключили, что ОММН с размерностью рабочего подпространства 30 наиболее быстрый, но он требует в четыре раз больше памяти чем другие методы. При этом учитывалось только время вычисления произведений матрицы на вектор, в то время как остальные части итерационного метода также могут занимать значительное время, особенно для ОММН(30). Среди менее требовательных к памяти методов КМН(КС) и Би-СГ(КС) показали лучшую сходимость чем Би-СГСтаб и ОПБи-СГ,
особенно при |m| > 2. Более того, авторы указали недостатки в работе Флато [145], которые делают его выводы спорными.
В разделе 2.4 сравниваются КМН(КС), Би-СГ(КС) и Би-СГСтаб при моделировании светорассеяния шарами в интервале x вплоть до 160 и 40 для m = 1.05 и 2 соответственно.
Рахола [147] показал, что спектр интегрального оператора рассеяния для любого однородного рассеивателя представляет собой отрезок в комплексной плоскости от 1 до m2, за исключением небольшого количества точек, соответствующих резонансным показателям преломления для конкретной формы. Спектр матрицы A примерно такой же, поскольку эта матрица получается дискретизацией интегрального оператора (см. также [51]). Предполагая, что спектр A полностью лежит на указанном отрезке, была
выведена следующая оценка для оптимального фактора уменьшения* γr: |
|
|||||||||
γr = |
|
|
|
1 |
|
|
|
|
. |
|
|
|
|
2 |
|
|
(79) |
||||
|
|
|||||||||
|
|
|
|
2 |
|
|
|
|||
|
|
1+ |
|
+ |
|
|
|
|
|
|
m2 −1 |
m2 −1 |
|
||||||||
Формула (79) приближённая, применимая только для малых частиц, которые обладают лишь несколькими резонансами. Тем не менее даже в общем случае спектр A близок к спектру интегрального оператора, который определяется формой, размером и показателем преломления рассеивателя. Поэтому спектр, а следовательно и
* Норма невязки уменьшается на этот множитель каждую итерацию.
48