Содержание к диссертации
Введение
ГЛАВА 1. Обзор и анализ методов обработки матрице задачах экспериментальной физики 16
1.1. Матричные методы редукции в задачах оптической диагностики 16
1.1.1. Основные задачи и схемы измерений в оптической диагностике дисперснофазных струй 16
1.1.2. Обобщенная математическая модель измерения параметров дисперсных веществ и функциональная схема прибора компьютерной диагностики 23
1.1.3. Постановка задачи редукции распределенных параметров и методы решения обратных задач в оптической диагностике 26
1.2. Методы обработки изображений и сигналов 34
1.2.1. Важные особенности постановки задач спектрального анализа сигналов 35
1.2.2. Обычные методы спектрального анализа 36
1.2.3. Методы, основанные на моделях исследуемых процессов 38
1.2.4. Моделирование сигналов на основе сингулярного разложения 44
1.3. Анализ устойчивости определения параметров линейной модели эмпирической зависимости по методу «наименьших квадратов» 47
1.4. Анализ методов определения собственных значений действительных матриц 54
1.4.1. Методы приведения матриц общего вида к форме Хессенберга. Анализ их устойчивости и быстродействия 54
1.4.2. Алгоритмы QR-метода. Анализ их быстродействия и точности 59
1.4.3. Методы определения собственных значений матриц компактной формы с помощью корней характеристического уравнения 64
1.4.4. Апостериорные оценки погрешности нахождения собственных значений действительных матриц общего вида 67
ГЛАВА 2. Устойчивые быстродействующие ортогональные методы приведения действительных матриц к компактной форме и диагонализации симметричных матриц 72
2.1. Модифицированный метод Гивенса 72
2.1.1. Математическое обоснование модификации метода Гивенса 72
2.1.2. Анализ вычислительных погрешностей модифицированного алгоритма Гивенса 77
2.1.3. Анализ результатов численных экспериментов 79
2.2. Алгоритмы диагонализации трехдиагональных симметричных матриц 85
2.2.1. Математическое обоснование алгоритмов диагонализации трехдиагональных симметричных матриц 85
2.2.1.1. Алгоритм MS3DTG90 89
2.2.1.2. Алгоритм MS3DTGJac ...92
2.2.2. Анализ вычислительных погрешностей алгоритмов MS3DTG90, MS3DTGJac 95
2.2.3. Сравнение алгоритмов диагонализации трехдиагональных
симметрических матриц по результатам численных экспериментов 99
2.2.3.1. Анализ точности алгоритмов MS3DTG90, MS3DTGJac и QR-метода 100
2.2.3.2. Анализ быстродействия алгоритмов MS3DTG90, MS3DTGJac и QR-метода 105
2.3. Алгоритмы диагонализации заполненных симметричных матриц 108
2.3.1. Алгоритм MSTG90 110
2.3.2. Алгоритм MSTGJac 1 13
2.3.3. Анализ результатов численных экспериментов 119
ГЛАВА 3. Методы "коррекции линейной интерполяции" и "коррекции кратного корня" для нахождения собственных значений матриц специального вида 121
3.1. Математическое обоснование методов 121
3.2. Анализ быстродействия и точности метода "коррекции линейной интерполяции" (МКЛИ) и метода "коррекции кратного корня" (МККК) 129
3.2.1. Анализ скорости сходимости методов 129
3.2.2. Анализ точности методов 134
3.2.3. Численные примеры 135
ГЛАВА 4. Метод диакоптики спектра собственных чисел действительных матриц компактного вида 139
4.1. Алгоритм "затухающего маятника" метода "диакоптики" спектра собственных чисел трехдиагональных симметрических матриц 139
4.1.1. Математическое обоснование алгоритма 139
4.1.2. Анализ вычислительных погрешностей и быстродействия алгоритма "затухающего маятника" метода диакоптики спектра собственных чисел трехдиагональных симметрических матриц 149
4.2. Алгоритм "затухающего маятника" метода диакоптики спектра собственных чисел матриц в форме Хессенберга 149
ГЛАВА 5. Физические приложения ортогональных методов диагонализации и диакоптики действительных матриц 158
5.1. Измерение геометрических параметров изделий с помощью оптического анализа поверхностных неоднородностей 158
5.2. Интегральный метод измерения скорости частиц двухфазного потока 159
5.3. Интегральный метод определения температурного распределения частиц дисперснофазных струй 164
5.3.1. Метод редуцирования температурного распределения частиц по их интегральному тепловому спектру 165
5.3.1.1. Постановка задачи и аналитическая модель измерения 165
5.3.1.2. Физическая модель измерения и методика калибровки 167
5.3.2. Автоматизированный спектрофотометр "Диагностик-Т" на базе интегральной МДП-фотодиодной линейки 173
5.4. Диагностика дисперсности в процессе впрыска топлива 176
5.5. Метод "диакоптики" спектра энергий и электронной структуры в кластерных расчетах неупорядоченных систем 181
Выводы 187
Литература 189
Приложение 197
- Обобщенная математическая модель измерения параметров дисперсных веществ и функциональная схема прибора компьютерной диагностики
- Математическое обоснование алгоритмов диагонализации трехдиагональных симметричных матриц
- Алгоритм "затухающего маятника" метода диакоптики спектра собственных чисел матриц в форме Хессенберга
- Автоматизированный спектрофотометр "Диагностик-Т" на базе интегральной МДП-фотодиодной линейки
Введение к работе
Актуальность исследований. В настоящее время существо проблем в технике физического эксперимента в большей степени связано с совершенствованием техники регистрации и обработки сигналов, автоматизацией физического эксперимента с повышением его технологичности, обеспечивающей увеличение информативности с понижением трудоемкости затрат на него. Сложные дорогостоящие физические эксперименты (например, в физике атомного ядра и элементарных частиц) диктуют необходимость быстродействующей обработки больших массивов экспериментальных данных для достижения цели проведения эксперимента в реальном масштабе времени. Поэтому актуальным современным направлением является компьютерный физический эксперимент с использованием «интеллектуальных инструментов», другими словами, метакомпьютера - виртуального сверхкомпьютера в виде глобальной компьютерной сети распределенных вычислений для проведения анализа, обработки и визуализации данных, поступающих с датчиков, и требующих специальных средств взаимодействия программных компонентов и вычислительных систем, физически размещенных в научных центрах разных стран.
Обеспечение возможности работы измерительного оборудования и высокопроизводительных систем обработки информации в реальном масштабе времени достигается за счет применения параллельных вычислительных многопроцессорных систем на базе параллельных матричных процессоров в виде СБИС, встроенных в датчики измерительных устройств и архитектура которых использует структуру методов матричной алгебры. Например, методы решения множества линейных алгебраических задач, матричные методы цифровой фильтрации сигналов определили архитектуру «систолических» матричных процессоров конвейерных вычислений (встречно-поточных процессоров), процессоров «волновой обработки», высокопараллельных процессоров с перестраиваемой конфигурацией.
При решении многих физических задач интегро-дифференциальные уравнения, описывающие математическую модель измерений физического эксперимента (либо физической системы), сводятся к матричным уравнениям, которые необходимо редуцировать с определенной точностью. Данная проблема связана с применением и разработкой методов регуляризации, использующих
ортогональные методы приведения больших матриц к компактному виду (почти треугольной, треугольной, трехдиагональной, блочно-диагональной, диагональной), диагонализацией с определением их «энергетического» спектра, либо с их обращением. Например, в оптических обратных задачах, в статистической обработке экспериментальных данных, а также, и при разработке специализированных матричных процессоров на базе СБИС, устойчивость обращения матриц достигается ортогональными преобразованиями, реализующими разложение по сингулярным (собственным) значениям (разложение по энергетическому спектру). Хорошо известна, например, роль трехдиагональных матриц в теории ортогональных многочленов, к которым сводятся современные алгоритмы сплайновых аппроксимаций эмпирических зависимостей, в разностных методах решения задач математической физики, в расчетах энергетического спектра и электронной структуры «цепочных» молекул. Для них необходимы численно устойчивые алгоритмы определения собственных чисел с гарантированной точностью.
Состояние вопроса. Недостаточно развитое состояние в области создания ортогональных методов, которые характеризовались бы широкой универсальностью применения, нечувствительностью к свойству «плохой обусловленности» матриц, побуждает предпринимать попытки эффективно решить эту проблему. По данной проблеме широко известны работы зарубежных и отечественных научных школ под руководством Дж.Х. Уилкинсона, Б. Парлетта, В.В. Воеводина, С.К. Годунова, В.Н. Фадеевой, Д.К. Фадеева, А.Н. Тихонова и многих других авторов. Практически отсутствуют какие-либо эффективные методы декомпозиции задачи на собственные числа.
В качестве примеров актуальных задач, в которых и до настоящего времени существует потребность в совершенствовании методов диагонализации, декомпозиции и приведения к компактному виду больших матриц, можно привести следующие:
1) методы статистического анализа экспериментальных данных и
планирования эксперимента, линейного прогнозирования и моделирования
процессов измерений в технике физического эксперимента;
2) современные методы обработки сигналов: калмановская и винеровская
фильтрация, адаптивная трансверсальная фильтрация, спектральный анализ с
использованием сингулярного разложения, цифровая обработка изображений и
распознавание образов с помощью двумерных унитарных преобразований,
кодирование сигналов в системах передачи информации, проблемы устойчивости систем управления и связи, обработка сигналов в фазированных антенных решетках и другие;
3) квантовомеханические расчеты атомных и молекулярных систем,
кластерные методы в задачах физики твердого тела;
4) теория колебаний и теория распространения волн в различных средах и
многие другие.
Например, в квантовомеханических расчетах молекулярных систем и кластеров твердого тела практически до 70 - 80% «машинного» времени затрачивается на диагонализацию матриц и для «больших» систем с большим порядком матриц (до десятка тысяч и выше) практически не реально за один непрерывный счет на ЭВМ получить результаты (проблема «метакомпьютинга»). Поэтому в настоящее время актуальным является решение проблемы декомпозиции задачи диагонализации матриц «по частям». То есть задача сводится к разработке такого метода, осуществляющего инвариантное преобразование большой матрицы к блочно-диагональному виду, в результате которого появилась бы возможность каждый матричный блок более низкого порядка диагонализировать «по частям» в рамках «однопроцессорной системы» по очереди, либо параллельно в рамках «многопроцессорной системы».
Цель исследований заключается в расширении возможностей и
совершенствовании методов экспериментальной физики за счет создания
устойчивых, с гарантированной точностью и высоким быстродействием методов
приведения к компактному виду, диагонализации и декомпозиции матриц
большой размерности, представляющих массивы многомерных
экспериментальных данных в матричных измерительных уравнениях.
Задачи исследования:
1. Определение широкого класса задач экспериментальной физики, использующих модель измерительного уравнения в матричном виде, с анализом погрешностей вычисления собственных значений матриц большой размерности;
2. Разработка быстродействующих методов диагонализации симметричных
трехдиагональных и заполненных матриц, обладающих абсолютной сходимостью
и гарантированной точностью независимо от свойств «обусловленности» матриц;
3. Создание быстродействующего метода декомпозиции «компактной»
формы матриц большой размерности.
4. Разработка методов для решения обратных спектральных задач экспериментального исследования температурного распределения частиц в гетерогенных потоках.
Научная новизна результатов исследований:
Получены «апостериорные» оценки «ошибки смещения», «стандартного отклонения» и максимальной погрешности вычисленных различными методами собственных чисел действительных матриц, использующих инварианты матриц: след (S) и евклидову норму (Е).
Разработан «модифицированный» метод Гивенса, который за счет использования свойств рекуррентности пересчета определенной части элементов матрицы сокращает количество операций умножений в 4/3 раза, а, значит, имеет более высокое быстродействие, чем традиционный алгоритм Гивенса.
Разработаны алгоритмы диагонализации симметричных трехдиагональных и заполненных матриц, сходимость которых в отличие от наиболее эффективных, например, алгоритмов QR-метода является абсолютной и точность не зависит от свойств «обусловленности» исходных матриц.
Разработаны быстродействующие метод «коррекции линейной интерполяции» и метод «коррекции кратного корня» для определения корней «характеристического» уравнения трехдиагональной симметричной матрицы, позволяющие вычислять собственные числа матриц с высокой точностью и в случае «патологически близких» корней (с учетом ограниченности «машинной» точности, практически кратных корней).
Разработаны алгоритмы «затухающего маятника» метода «диакоптики» спектра собственных чисел симметричных трехдиагональных матриц и матриц в форме Хессенберга, осуществляющие декомпозицию матриц «по частям», характеризующиеся замедлением сходимости процесса диакоптики лишь для «плохо обусловленных» матриц.
Новизна технических решений. Разработан принципиально новый способ определения температуры частиц конденсированной фазы движущихся гетерогенных объектов. Новизна технического решения подтверждена патентом № 2107899 на изобретение, зарегистрированным в Госреестре изобретений (РОСПАТЕНТ), Москва, 27.03.1998.
Методы исследования.
В диссертации использованы методы экспериментальной и теоретической физики, методы оптической диагностики, теория взаимодействия светового
излучения с веществом, вычислительные методы линейной алгебры, численное моделирование и методы решения обратных задач, теория вероятности и случайных процессов, методы статистической обработки данных.
На всех этапах исследований проводилось сопоставление теоретических выводов и оценок с результатами компьютерного и физического эксперимента.
Практическая ценность работы:
Полученные апостериорные оценки погрешностей вычисленных собственных значений действительных матриц, которые являются оценками «снизу», наряду с «традиционными» априорными оценками «сверху», которые часто могут оказаться «завышенными», позволяют более достоверно оценивать погрешности вычисленных собственных значений.
Разработанный «модифицированный» метод Гивенса, практически не уступающий в отношении быстродействия методу Хаусхолдера как наиболее эффективному методу приведения действительной матрицы общего вида к компактной форме, характеризуется гарантированной численной устойчивостью в отличие от метода Хаусхолдера и поэтому может быть рекомендован к широкому применению.
Разработанные алгоритмы диагонализации симметричных
трехдиагональных и заполненных матриц могут быть рекомендованы в различных прикладных задачах на собственные значения для матриц большой размерности (я > 1000).
Разработанные быстродействующие метод «коррекции линейной интерполяции» и метод «коррекции кратного корня» для определения корней «характеристического» уравнения трехдиагональной симметричной матрицы позволяют вычислять собственные числа матриц с высокой точностью и в случае «патологически близких» и кратных корней;
Разработанные алгоритмы «затухающего маятника» метода «диакоптики» спектра собственных чисел трехдиагональных матриц и матриц в форме Хессенберга, осуществляющие декомпозицию матриц «по частям», позволяют для «сверхбольших» матриц решать задачу на собственные числа.
Реализация результатов.
Для всех предлагаемых к широкому использованию методов и алгоритмов, изложенных в диссертационной работе, разработано соответствующее программное обеспечение на алгоритмических языках ФОРТРАН и ПАСКАЛЬ.
Проведено их детальное тестирование в ходе численных экспериментов, в результате которого получено подтверждение всех теоретических оценок и выводов.
На новый способ определения температуры частиц конденсированной фазы движущихся гетерогенных объектов получен патент РФ № 2107899, использующий предлагаемые в работе матричные методы.
Разработанные методы, алгоритмы и соответствующие им программы использовались при решении конкретных прикладных задач, результаты которых отражены в публикациях.
На защиту выносятся следующие основные научные результаты:
Методика «апостериорных» оценок погрешности вычисляемых различными методами собственных значений действительных матриц;
Модифицированный метод Гивенса для приведения матрицы общего вида к компактной форме и быстродействующие методы диагонализации симметричных трехдиагональных и заполненных матриц, основанные на устойчивых элементарных вращениях, обладающие абсолютной сходимостью и гарантированной точностью независимо от свойств «обусловленности» матриц;
Быстродействующие алгоритмы «затухающего маятника» метода диакоптики спектра собственных чисел трехдиагональных симметричных матриц и матриц в форме Хессенберга, реализующего процедуру декомпозиции матрицы путем последовательного приведения б л очно-диагональной формы к диагональной и форме Шура соответственно.
4. Метод редукции спектральной задачи экспериментального исследования
температурного распределения частиц в гетерогенных потоках.
Публикации. Основные результаты диссертационной работы опубликованы в 24 печатных работах, получен один патент на изобретение.
Апробация работы. Основные положения и результаты диссертационной работы докладывались и обсуждались на следующих всероссийских и международных конференциях и совещаниях: Всесоюзная научная конференция «Современное состояние теории атомов и молекул», г. Вильнюс, 1979г„ Всесоюзная научная конференция «X Сибирское совещание по спектроскопии», г. Томск, 1981 г., Всесоюзная научная конференция «Планарные дефекты в упорядоченных сплавах и интерметалл ид ах», г. Барнаул, 1987г., Всесоюзная научная конференция «Координатно-чувствительные фотоприемники и оптико-электронные устройства на их основе», г. Барнаул, 1989г., Всесоюзная научная
конференция «Оптические сканирующие устройства и измерительные приборы на их основе», г. Барнаул, 1990 г., 10-ое Всесоюзное координационное совещание по квантовой химии, г. Казань, 1991г., Международная конференция по алгебре памяти А.И. Ширшова, г. Новосибирск, 1991г., Первая международная конференция «Нанотехнология, наноэлектроника и криоэлектроника», г. Барнаул, 1992 г., Третья международная конференция памяти М.И. Каргаполова, г. Красноярск, 1993г., Международная конференция «Всесибирские чтения по математике и механике», г. Томск, 1997г., Международная научно-техническая конференция «Совершенствование быстроходных ДВС», Барнаул, 1999 г., 8-ая международная конференция «Математические методы в электромагнитной теории», г. Харьков, Украина, 2000г.
Структура и объем работы. Диссертационная работа состоит из введения, 5 глав, заключения, списка литературы и приложения. Работа изложена на 197 страницах машинописного текста, содержит 34 рисунка, 11 таблиц и список литературы из 102 наименований.
Краткое содержание работы
Во введении обоснованы актуальность, научная и практическая значимость проблемы, сформулированы цель и задачи исследований, их научная и практическая новизна, изложены основные выносимые на защиту положения, приведена краткая характеристика работы.
В первой главе диссертации проведен обзор широкого класса задач экспериментальной физики, в которых возникает проблема обработки многомерных данных, представляющих собой матрицы большой размерности. К ним сводятся задачи статистической обработки экспериментальных данных с использованием, например, метода «главных компонент» или метода «наименьших квадратов» (МНК), задачи планирования физического эксперимента, обработка сигналов, квантовомеханические расчеты кластерных неупорядоченных систем, оптические и электрические измерения и другие. Их общим признаком является наличие матричного уравнения, или сводящегося к нему интегрального или интегро-дифференциального уравнения, связывающего экспериментальные данные с физическими параметрами объекта исследования.
Показано, что решение таких измерительных уравнений в экспериментальной физике сводится к некорректным обратным задачам, в которых актуальной проблемой является разработка устойчивых
быстродействующих методов обращения, приведения к компактному виду и диакоптики матриц большой размерности.
В научной литературе, главным образом, приведены оценки погрешностей вычисления собственных значений в виде «априорных» оценок «сверху», которые не редко оказываются «завышенными». Для более объективной оценки реальной погрешности в конце первой главы предложены «апостериорные» оценки «снизу»: оценка ошибки «смещения» , минимальная оценка
«стандартного отклонения»
Выяснено, что в настоящее время применяемые для этой цели методы Гаусса-Жордана, Хаусхолдера, QR - метод и другие не гарантируют необходимой точности и абсолютной сходимости. Анализ причин, приводящих к этим недостаткам, позволил сделать вывод, который обосновывает выбор направления разработки методов решения таких задач, основанных на определении энергетического спектра матриц большой размерности.
Во второй главе представлен к рассмотрению модифицированный метод Гивенса, который в результате «рационализации» стандартного алгоритма Гивенса и с учетом свойств рекуррентности пересчета элементов сокращает число умножений в 4/3 раза.
Показано, что ошибки накопления в стандартном и модифицированном алгоритмах одного порядка. Численные эксперименты подтверждают такой вывод, а также вывод о том, что метод Хаусхолдера не обладает гарантированной устойчивостью и точностью преобразований. Численные эксперименты показали практически одинаковое быстродействие модификации метода Гивенса и метода Хаусхолдера.
Далее рассмотрены предлагаемые в диссертации быстродействующие алгоритмы диагонализации симметричных трехдиагональных и заполненных матриц, абсолютная сходимость и точность которых гарантирована независимо от «обусловленности» матриц.
Для трехдиагональных матриц в предлагаемых двух алгоритмах на каждом этапе, который представляет собой итерационный процесс, в верхней части матрицы происходит отделение двух собственных чисел для первого алгоритма (его обозначение MS3DTG90), а для второго алгоритма (обозначение
lVrS3DTGJac) - трех собственных чисел. То есть, в верхней части матрицы на каждом этапе происходит преобразование трехдиагональной формы в диагональную. Итерационные процессы характеризуются абсолютной сходимостью.
Для алгоритмов MS3DTG90 и MS3DTGJac получены оценки погрешности на собственные числа, пропорциональные величине пъ'2. Численные эксперименты показали, что реальные погрешности несколько ниже. Сравнение по быстродействию с алгоритмами QR - метода в численных экспериментах показало, что для широкого класса матриц с условием на ее порядок п > 1000 наступает преимущество предлагаемых алгоритмов MS3DTG90 и MS3DTGJac над QR — методом.
Для заполненных матриц каждый из двух предлагаемых алгоритмов является обобщением, соответственно, алгоритмов MS3DTG90 и MS3DTGJac, используемых для диагонализации трех диагональных матриц. Первый алгоритм (его обозначение MSTG90) также в результате одного этапа в верхней части матрицы отделяет два собственных числа, а второй алгоритм (обозначение MSTGJac) - три собственных числа. Алгоритмы MSTG90 и MSTGJac характеризуются абсолютной сходимостью, не зависящей о і cbomclb «обусловленности» заполненных матриц, устойчивостью и гарантированной точностью вычислений собственных чисел, а также высоким быстродействием. Поэтому эти алгоритмы могут быть рекомендованы к широкому использованию в задачах на собственные числа.
В третьей главе предлагаются два метода, которые для матрицы компактной формы, например, трехдиагональной симметричной формы, позволяют вычислять с высокой точностью корни (собственные числа) характеристического полинома этой матрицы даже и в случае «патологически близких (кратных)» корней (то есть, для плохо обусловленных матриц). Метод «коррекции кратного корня» (МККК) позволяет разделять «спектр собственных чисел» на простые корни и в случае «патологически близких» корней. Если «патологическая близость» корней находится в пределах «машинного нуля», тогда МККК позволяет локализовать «почти кратные» (и кратные) корни с заданной точностью. Второй метод - метод «коррекции линейной интерполяции», ускоряя сходимость, уточняет локализованные простые корни.
Алгоритм, использующий комбинацию методов МККК и МКЛИ, позволяет определять собственные числа трехдиагональной матрицы с высокой точностью (до 16 значащих цифр) и со скоростью, достаточно близкой к скорости QR-метода. Проведенные численные эксперименты подтвердили эффективность такого алгоритма.
В четвертой главе рассмотрены два алгоритма диакоптики спектра собственных чисел матриц компактной формы (трехдиагональной и формы Хессенберга - почти треугольной). Блочные структуры матрицы А[к~ц перед к — ым шагом, который будем называть к - ой «большой итерацией» (БИ), и ее матричных «клеток» А(2\~1) и А^~[) для случая трехдиагональной симметричной
формы имеют вид
л(к-\)
1_Л21
А(Ь-1)
о ь(к-])' о о
j(A-l) Л 2
=«-п)г =
ьа-и
О О
Через Ь(к Х) обозначен элемент «связи» а^Л трехдиагональных
симметричных клеток 4, ' и А\{ > гДе m-порядок клетки А\х~ \ На диагонали
клетки А{2\~]) слева от 6'*_1) находится нулевая строка, ниже - нулевой столбец.
Оставшаяся нулевая клетка имеет размерности на единицу меньше, чем А^"'\ По
отношению к структуре клетки А^~[) в структуре клетки А\кг~]) нулевые строка и
111 " "22
клетка Л(2\~]) имеет структуру, аналогичную структуре А2]
столбец поменялись местами, а также элемент b{k~l) с нулевой клеткой. В случае формы Хессенберга клетки А{\к~]) и А22'^ имеют также форму Хессенберга,
(*Ч) трехдиагональной
симметричной формы Atk~]). Поэтому элемент b(k~1} в клетке ^Г'1 в этом случае
выполняет ту же роль элемента «связи». Для клетки Д(2_1> свойство
АІ2 =(^2i~>}У не выполняется и она, как правило, является заполненной. В обоих алгоритмах соблюдается правило чередования «прямой» и «обратной» БИ, которые состоят из последовательности вращений, и обеспечивают инвариантность каждой из компактных форм.
Рекуррентные схемы в обоих алгоритмах обеспечивают сходимость элемента связи диагональных клеток к нулю при к~>оо по принципу «затухающих колебаний» маятника за счет «сжимающего» оператора, определяющегося синусами вращения.
В обоих алгоритмах для плохо обусловленных матриц каждой из компактных форм значения «сжимающего» оператора по модулю могут быть
близкими к единице, что ведет к замедлению сходимости. Объем вычислений, соответственно, и время вычислений для трехдиагональных матриц оказывается порядка &, п log2 п, для почти треугольных — порядка к2 п2, что примерно в п
раз меньше, чем для наиболее эффективных методов. Поэтому применение алгоритмов диакоптики спектра собственных чисел для больших и сверхбольших матриц актуально.
В пятой главе рассмотрены конкретные приложения вынесенных на защиту методов, использованных в некоторых задачах экспериментальной и прикладной физики, и показана перспективность применения этих методов в квантовомеханических расчетах кластерных неупорядоченных систем.
В заключении диссертации сформулированы основные выводы и результаты, а также проблемы, требующие дальнейшего решения.
Личный вклад автора заключается в формулировке и постановке задач, выборе методов их решений, разработке конкретных методов определения энергетического спектра матриц большой размерности и анализе их устойчивости, точности и быстродействия и непосредственном участии в экспериментальных исследованиях.
Автор выражает глубокую благодарность всем соавторам. Особую благодарность хочется выразить научному руководителю доктору физико-математических наук, профессору, заслуженному деятелю науки РФ Евстигнееву В.В., а также доктору технических наук, профессору Гуляеву П.Ю за постоянное внимание и помощь при подготовке диссертационной работы. Отдельно следует отметить важную роль кандидата физико-математических наук, доцента Аникеева B.C. |в становлении и формировании научного мировоззрения автора.
Обобщенная математическая модель измерения параметров дисперсных веществ и функциональная схема прибора компьютерной диагностики
Спектральная задача нахождения функции распределения частиц по размерам C(z) в этом элементарном объеме представляется как задача по обращению интегрального уравнения Фредгольма первого рода [45]: где F(x, г)-ядро уравнения, известное из теории рассеяния света на отдельной частице [45,46], (х) - экспериментально определяемая функция параметра х. Ядро уравнения F(x,z) может быть индикатрисой рассеяния монодисперсной системы частиц радиуса z, тогда Ч (х) - индикатриса полидисперсной системы; либо F(x, z) — показатель рассеяния, тогда 4у(х) —показатель рассеяния при длине волны А = х, и т. д. Поиск метода расчета функции распределения C(z) по известным функциям F(x, z) и Ч (х) составляет задачу теории обращения.
Обращение интегрального уравнения первого рода (1.1.1) через замену его алгебраической системой, выраженной в матричной форме: где -столбец значений экспериментально определяемой функции H ix); F- матрица, составляющая ядро уравнения (1.1.1); С-столбец значений функции распределения по размерам, приводит к принципиальным трудностям, свойственным задачам такого вида. Погрешности в измерении значений Ц и неточности в расчетах ядра F, приводят к серьезным ошибкам в определении значений функции С, причем правильное математическое решение уравнения (1.1.2) может содержать отрицательные корни, то есть физически не реализуемые. Такие задачи являются некорректными, причем, чем больше размерность уравнения (1.1.2), тем более чувствительно к ошибкам ее решение. Однако, в некоторых случаях задачу удается решить в матричном виде (1.1.2).
В одном из таких случаев экспериментально находят оптический коэффициент пропускания прямо прошедшего излучения г(Л,) или оптическую толщину среды через слой рассеивающих частиц толщины L на различных длинах волн А . Этот метод разработан К.С. Шифриным с сотрудниками [20,45] и известен как метод спектральной прозрачности. Метод состоит в экспериментальном измерении показателя рассеяния ст(Я.) частиц в некотором объеме среды на различных длинах волн. Показатель ослабления слоя рассеивающих частиц определяется по закону Бугера: е(А ) = 1п(1 / г(Л .)) / Z,. Определение функции распределения частиц C(z) производится согласно: 1([3) определяется в результате опыта, у(Р) вычисляется по (1.1.9) с известным 1(р). Угол Р определяется размерами фокального пятна, pmax=5 -6 . 1.1.2. Обобщенная математическая модель измерения параметров дисперсных веществ и функциональная схема , прибора компьютерной диагностики диагностики параметров ДФС интегральными методами контроля рассматривается на примере импульсного ел абозапыл енного потока частиц (рис. 1.3), поочередно и в случайный момент времени пересекающих измерительный объем. Очевидно, что такая модель так же справедлива для контроля неподвижной ДФС и ТИС с известным законом сканирования измерительного объема, заполненного случайным образом частицами, размер которых не меньше величины оптической разрешающей способности прибора. В нашем случае сканирование частиц обеспечивается движением потока относительно измерительной оптической системы (ОС). Каждой частице соответствует свое значение контролируемого физического параметра z(, связанного с интенсивностью оптического излучения Q \х) законом физической оптики A(znx), например дифракционного рассеяния, спектром теплового излучения, лучевой проекции и т.п. При этом л:-регистрируемый параметр оптического излучения (угол дифракционного рассеяния, длина волны, координата проекции и т.п.). Фотоприемник ТИС и тракт ОС вносят искажения в оптический сигнал, которые описываются соответствующими аппаратными функциями А:{ (х) и к2 (х) (спектральная чувствительность, дифракционный предел разрешающей способности, астигматизм, аберрация и т.п.). Преобразование параметра z, в выходной электрический сигнал ТИС gi(x) = A(zrx)- к{(х)- к2(х), может быть описано с помощью обобщенной аппаратной функции прибора K(z,x). Входной оптический сигнал, регистрируемый за полное время сканирования частиц, в интегральном виде выглядит так: где f(z) -dil dz- искомая функция распределения контролируемого параметра ДФС. В случае неподвижных ДФС, (х) может иметь физический смысл индикатрисы рассеяния полидисперсной среды, теплового спектра разнородно нагретых тел и т.п. В рассматриваемом случае, регистрации движущегося потока частиц, С І ( ) может проявляться в виде треков или "смаза" изображения. Решение задачи «интегрального метода контроля (ИМК)» состоит в определении функции f(z) по выходному сигналу ТИС, когда оптический параметр х взаимно однозначно определяется координатой точки сканируемого изображения, например, когда каждому элементу телевизионного изображения соответствует своя длина волны, угловой или линейный параллакс световых лучей рассеянных или излучаемых частицами ДФС. В простейшем случае проекционной ОС параметр х может быть просто координатой изображения частицы. Прямую задачу измерения можно записать в интегральном виде: g(x) = кх (х)к2 (хХ(х) = Я , {Ф2 {x)A(z, x)]f(z)dz = \K(z, x)f(z)dz, которая в дискретном виде соответствует операторному уравнению g - А f -для идеального прибора, и g = К f - для реального прибора. Рассмотрим задачу определения величины / прибором А по результатам косвенных измерений g. Конкретный вид оператора А определяется физическим методом измерения и аппаратной функцией прибора к(г), но предполагается выполнение основного интегрального измерительного уравнения диагностики: На практике истинное значение g никогда не известно, так как оно всегда содержит некоторую экспериментальную ошибку п. Поэтому измеренные данные gd можно представить в виде где gd -известная матрица размера Mxl, а п матрица экспериментальных ошибок размера Mxl. Решение задачи диагностики сводится к процедуре подбора искомой функции/и минимизации следующей положительной величины: где индекс + обозначает комплексное сопряжение и транспонирование, j\ и Уг положительные константы, /0 -"пробная" функция, а В -матрица размера МхМ, описывающее некоторое сглаживание /. Первый член в (1.1.12) является мерой точности /. Если экспериментальных ошибок нет, то gd = g, у] = у2 = О, U - 0 минимальное значение U и решением является f = А ]g. Второй член в (1.1.12) указывает на отклонение / от пробной функции /0, а третий член характеризует отклонение / от идеального сглаживания, соответствующего (Bf) - 0. Обычно (Bf) описывает первую или вторую производную [1-3]. Рассмотрим случай у2=0. Дифференцируя U по /, находим где / - квадратная матрица размера N N. Выбор ух должен быть сделан так, чтобы обеспечить разумный компромисс между первым и вторым членами в (1.1.12). Если у{ мало, то преобладает первый член, и решение является сильно осциллирующим. Если j\ велико, то преобладает второй член, и решение получается сильно сглаженным. Дальнейшие подробности можно найти в работе [19].
Математическое обоснование алгоритмов диагонализации трехдиагональных симметричных матриц
Следует отметить, что поскольку усредняющее ядро (1.1.35) должно быть подобно дельта-функции, для подавления A(rQ, г) в точках rx, отличных от rQ, требуется некоторое число членов N. Это число N может несколько превышать число осцилляции функции К(г). Влияние числа измерений и их ошибок на разрешение рассмотрено в [3]. Компьютерный прибор можно характеризовать тремя паспортными параметрами как обычный физический прибор: Зависимости Я = Н(є, 5),G = G{s,S),q = q(s, S) задают оперативную характеристику комплекса, которая служит паспортом, определяющим его "приборные" возможности. Оптимизация параметров прибора производится стандартными методами регуляризации в соответствии с алгоритмом (рис. 1.5).
Кроме рассмотренных в подразделах 1.1.2, 1.1.3.1, 1.1.3.2 методов редукции матричных уравнений в задачах обработки экспериментальных данных необходимо отметить методы пространственной фильтрации дискретных изображений и сигналов в задачах обработки различных снимков: рентгеновских, электронно-микроскопических, аэрокосмических, томографических, термограмм, радиолокационных и звуколокационных карт, телевизионных изображений и многих других. К этим методам относятся: метод регрессии, винеровское оценивание, методы двумерной рекурсивной и нерекурсивной фильтрации [50]. В задачах реставрации и улучшения изображений снимков важную роль играют двумерные унитарные преобразования Фурье, Хаара, Адамара-Уолша, Карунена-Лоэва и другие [50]. В задаче регистрации движущегося потока частиц, рассмотренной в подразделе 1.1.2, при использовании оптических систем "смаз" изображения можно характеризовать пространственно-инвариантной функцией рассеяния, для которой модель искажений имеет вид где G—искаженное изображение; Q- истинное изображение.
Операторы смазывания А и В, соответственно, столбцов и строк изображения можно считать «циркулянтами» [50]. Тогда имеет место представление где F - матрица Фурье-преобразования, диагонализирующая матрицы А и В с результатами А]А и Л1 2. Реставрированное изображение получается путем устранения разделимого пространственно-инвариантного смаза за счет инверсной Фурье-фильтрации согласно (1.2.1) и (1.2.2) в виде Инверсная Фурье-фильтрация по своей сути связана с методом редукции «псевдообращения на базе разложения по сингулярным (собственным) значениям» (ПО+РСЗ), который играет важную роль в задачах спектрального анализа сигналов. Кратко рассмотрим основные методы и постановки задач спектрального анализа сигналов [48,49].
Спектральный анализ составляет основу большинства методов обработки сигналов, используемых обычно для выделения представляющих интерес сигналов, а также для получения информации из соответствующих данных. Задача классического спектрального анализа состояла в оценке энергетического спектра по конечному числу зашумленных отсчетов дискретного стохастического процесса или по нескольким значениям оценки его корреляционной функции. Задача современного спектрального анализа в некоторых практических случаях, например, в радио- и гидролокации, при обработке сигналов в фазированных антенных решетках, где спектр является линейчатым, состоит в оценке положения спектральных линий. Как результат этого понятия смещения и дисперсии, которые прежде относились к оценке формы спектра, теперь относятся к оценке частоты. Другой качественной характеристикой классического спектрального анализа является разрешающая способность по частоте, которая определяет степень детальности исследования спектра [48,49]. В современной постановке разрешающая способность есть способность различать и идентифицировать спектральные линии, которые расположены близко по частоте [49]. Во многих современных случаях применения спектральный анализ приходится производить, имея в распоряжении короткие записи данных (в радиолокации, например, в каждом отраженном импульсе может содержаться всего несколько отсчетов), но, тем не менее необходимы оценки с высоким разрешением, малыми смещением и дисперсией.
Разрешающая способность по частоте в обычных методах, основанных на преобразовании Фурье, примерно обратно пропорциональна объему выборки. Поэтому для улучшения разрешающей способности необходимо ввести некоторые дополнительные ограничения (или априорную информацию). С этой целью в современных методах данные моделируются как процесс на выходе линейной системы, на вход которой подан белый шум. Когда такая модель соответствует реальному сигналу, эти методы приводят к улучшению спектральных оценок. Цена, которую приходиться платить за это улучшение, как правило, состоит в увеличении сложности вычислений по сравнению с таковой при обычных методах, основанных на использовании быстрого преобразования Фурье (БПФ) в виде спецпроцессора на базе СБИС.
Рассмотрим обычные методы анализа энергетического спектра, а затем сосредоточим внимание на современных методах, которые подходят для разрешения узко полосных сигналов, в том числе и на методе ПО+РСЗ.
Алгоритм "затухающего маятника" метода диакоптики спектра собственных чисел матриц в форме Хессенберга
Если шаг разбиения h достаточно мал, то есть nh«\, то элементы / будут иметь значения на уровне погрешностей вычислений (либо последние элементы окажутся равными нулю), то есть постановка обратной задачи согласно (1.3.11) с увеличением п (объема экспериментальных наблюдений) имеет тенденцию к «вырождению». Это означает, что решение ни с помощью метода Гаусса, как это показано выше в (1.3.20)-(1.3.21), ни с помощью метода Хаусхолдера, согласно (1.3.15)-(1.3.19), может быть не найдено, так как обратная подстановка производится по одной и той же процедуре (1.3.18) и (1.3.19). Как видно из (1.3.18), (1.3.19) и (1.3.21), компоненты вектора решения, начиная с последнего, последовательно домножаются на обратные элементы диагонали треугольной формы, которые могут оказаться очень малыми. Также необходимо отметить, если значение шага h, напротив, окажется порядка 10 и более, то с увеличением п в вычислениях может произойти «переполнение порядка». Ситуация мало изменится, если даже шаг h будет не постоянен. Из анализа, проведенного выше, можно сделать вывод о том, что «корректной» постановкой обратной задачи МНК является «классическая» постановка задачи в виде (1.3.2)-(1.3.10), связанная с решением «нормальной» системы уравнений (1.3.7)-(1.3.10), матрица А которой является симметричной. В случае неортогонального базиса, являющегося полной системой функций, симметричная матрица А скалярных произведений является невырожденной. Поэтому требуется лишь применять устойчивые методы решения. Как уже известно из многих источниках [52,54,66,69], методы Гаусса и Хаусхолдера не являются абсолютно устойчивыми и для симметричной невырожденной матрицы (1.3.7). Как видно из (1.3.18) (1.3.19), компоненты вектора решения, начиная с последнего, в процессе обратной подстановки последовательно домножаются на обратные элементы диагонали треугольной формы, и это произведение согласно (1.3.21) накапливается «факториально». Поэтому решение системы (1.3.7)-(1.3.10) предлагается в виде С-А ]-В. Обратную матрицу А предлагается определять в виде A =QT-A j-Q, используя для симметричной матрицы А представление A = Q -AA-Q, где -диагональная матрица спектра точных собственных чисел матрицы А и Q - матрица точного ортогонального преобразования, для которой выполняется условие Q = QT. В реальных вычислениях спектр собственных чисел A-A=Q-A-Q находится с определенной точностью. Запишем А =D-Q -A A Q,TO есть C=D-B будет зависеть от обратных значений собственных чисел линейно в силу ортонормированности матриц QT и Q. Известный своей эффективностью метод Якоби по существу носит итерационный характер [52-54]. Методы Гивенса и Хаусхолдера [52], которые излагаются ниже, приводят с помощью унитарных преобразований матрицу общего вида к компактной форме - к почти треугольной форме (форме Хессенберга), а симметрическую матрицу - к трехдиагональной форме. Такое приведение может быть выполнено не итерационными методами и, следовательно, оно требует намного меньше вычислений. Собственные числа матриц формы Хессенберга (трехдиагональной формы) определяются простыми и весьма эффективными методами, например QR-методом и другими [52]. Приведение по методу Гивенса состоит из {п - 2) основных шагов. Перед m-ым шагом первые (т-1) столбцов исходной матрицы А уже приведены к верхней форме Хессенберга, так что, например, при п-в и т = 3 имеем Можно сказать, что главный ведущий минор порядка (т -1) уже приведен к форме Хессенберга, и эта часть не меняется при дальнейших преобразованиях. Через АтЛ обозначим матрицу, полученную из AQ после приведения столбцов с номерами 1,2,....,т-1. Основной т-ый шаг состоит из (п-т-1) вспомогательных шагов, в которых последовательно аннулируются элементы в позициях {т + 2,т + 3,...,п) столбца с индексом т. Нуль в позиции {m + \+i), где г = 1,2,...,«-т-1, получается при помощи вращения в плоскости (т + \,т + 1 + і). В (1.4.1) подчеркнуты элементы, которые нужно аннулировать. Основной т-ый шаг можно описать следующим образом: Для всех / от 1 до (и - т - 1) выполнить шаги от 1) до 4): 1) вычислить 6, =\bl a2m Jn\ К=атЯ„, 2) вычислить с( = 6М /6., s: amrUiii) (br Записать bt на место атЛт и нуль на место am+l+i w; 3) для всех значений j: = т +1, т + 2,..., п вычислить ат+] j-c,+ am+Ul s, и а„,, после чего записать их на место а „., , и а ,,,,,; 4) для всех значений j = 1,2,..., п вычислить [- ,+ +, -5, и й,,++; с, «y,m+i s/ после чего записать их на место ajrtl+l и aym+1+j; Полное количество умножений приблизительно равно а число операций сложений (вычитаний) Для симметричных матриц с учетом использования лишь симметрии наддиагональных элементов число умножений для полного приведения к трехдиагональной форме требуется около Апъ /3 и число операций сложений (вычитаний) приблизительно равно 2пг 13.
Автоматизированный спектрофотометр "Диагностик-Т" на базе интегральной МДП-фотодиодной линейки
Гарантированные оценки для преобразований отражений с использованием накопления скалярных произведений не могут быть лучше оценок типа (1.4.8). Однако вероятностный анализ и гарантированные оценки, если при реализации преобразований отражения не используют накопления скалярных произведений, приводят к оценкам порядка п1 [52,57]. Если вычисления ведутся с одинарной точностью, то влияние ошибок округления в типичных последовательностях преобразований вращений будет в ып или даже п раз по порядку меньше, чем в преобразованиях отражений, решающих ту же задачу [57]. В [52] в отношении преобразований отражений указывается на важность правильного выбора знака sm, соответствующего знаку а( 1 , так как в случае неправильного выбора знака относительная ошибка преобразований становится большой. Такая ситуация в реальных вычислениях на ЭВМ встречается не редко, когда для чисел, близких к машинному нулю, знак а\ \ определяется неправильно. В разделе 2.1.3 приведены конкретные примеры неустойчивого поведения метода Хаусхолдера. В [69] метод Хаусхолдера также подвергается критике и приводится способ, который чаще всего дает малые величины ошибок округлений. В [69] обсуждается вариант метода Хаусхолдера, использующий стратегию перестановок строк (столбцов) на каждом этапе приведения, которая, вероятнее всего, не может гарантировать абсолютную устойчивость. Это предположение можно проанализировать на следующем матричном примере. Рассмотрим матрицу А-[а ], где i,j = 1,2,...,п. Порядок матрицы п может принимать любое целое значение (п 3). Допустим, что элементы аи} = 1 для всех / и j.B результате применения точных преобразований по методу Гивенса и методу Хаусхолдера уже после 1-го этапа приведения преобразуемая матрица принимает трехдиатональную форму, На главной диагонали первый элемент равен 1, второй элемент равен значению (п -1) и остальные элементы главной диагонали оказываются равными нулю. Первые элементы обоих кодиагоналей (побочных диагоналей) оказываются равными значению -Jn-l и остальные элементы кодиагоналей равны нулю. Таким образом, в преобразованной матрице верхняя клетка 2-го порядка не содержит нулевых элементов, а остальные элементы матрицы равны нулю. Собственные значения этой клетки равны значению п и нулю. Следовательно, исходная матрица рассмотренного примера имеет собственное число со значением п и нулевое собственное число кратности («-1). В реальных вычислениях на ЭВМ после одного этапа приведения исходной матрицы, элементы которой равны 1 (или матриц, отличающихся от нее на некоторое возмущение малого порядка), вместо нулевых элементов будем иметь малые элементы на уровне ошибок округления. В этом случае все оставшиеся не приведенными строки (столбцы) будут равносильны между собой в отношении стратегии их перестановок. Можно сделать вывод о том, что для матриц, подобных свойствам матриц рассмотренного примера, стратегия перестановок строк не гарантирует численную устойчивость. Устойчивость же метода Гивенса обеспечивается устойчивым вычислением "косинуса и синуса" каждого элементарного вращения с гарантированной точностью. Численные эксперименты подтверждают этот вывод. Резюмируя выше сказанное, можно сделать следующее заключение. Существует необходимость в разработке метода приведения больших матриц к компактной форме, основанного на устойчивых элементарных вращениях, быстродействие которого было бы не хуже, чем у метода Хаусхолдера. 1.4.2. Алгоритмы QR-метода. Анализ их быстродействия и точности Широко известные алгоритмы QR-метода [52,54,55-58,63-65] позволяют приводить произвольную действительную матрицу к треугольному (диагональному) виду с помощью унитарных преобразований. Для матриц в форме Хессенберга (в трехдиагональной форме) эти алгоритмы во многих отношениях являются наиболее эффективными средствами определения собственных значений [52,54,67], Идея, заложенная в QR-методе, состоит в использовании разложения матрицы в произведение унитарной матрицы Q и верхней треугольной матрицы R. Процесс определяется соотношениями A„=QmRm, A =Q:n-AmQm=QlQmRmQm=RmQm, (1.4.9) где Q m - эрмитово сопряженная матрица к матрице Qm на т - ом шаге унитарного преобразования подобия. Если Ат - симметричная матрица, то Qm = Qm -результат транспонирования матрицы Qm, которая является ортогональной матрицей, то есть Q m = Q = Q m. В результате итерационного процесса унитарных преобразований по формуле (1.4.9), начиная с некоторого т, матрицы Ат будут близки к треугольной (диагональной) матрице, на диагонали которой содержатся собственные числа. Доказательство сходимости такого процесса содержится в [52,63,64]. Достоинством QR-метода является то, что метод сохраняет форму Хессенберга (трехдиагональную форму) и "ленточную" форму в процессе преобразований [52,54,67]. Для ускорения QR метода используется так называемый сдвиг sm с восстановлением. Кублановская В.Н. доказала [63,64], что в итерационном процессе последний внедиагональный элемент убывает к нулю со скоростью {Яп I Яп )т, где Щ \Я ЯП_, \Яп\. Процесс, использующий сдвиги с восстановлением, имеет вид Am-sm-I = Qm-Rm, Rm-Qm+sm-I = Am+l. (1.4.10) Тогда скорость убывания к нулю последнего внедиагонального элемента будет оцениваться величиной [(А„ -sm)!(X sm)] ". Величина сдвига sm выбирается в соответствии с определенной стратегией выбора, обеспечивающей близость сдвига к Я„. Процесс, не использующий сдвиги sm, может сходиться медленно, если некоторые или все собственные числа плохо разделены. Процесс (1.4.10) будет сходиться тем быстрее, чем сильнее различаются ближайшие друг к другу собственные числа. В случае очень близких друг к другу собственных чисел Яп_{ и Яп выбор sm, который в каждой стратегии определяется через элементы матрицы Ат, может оказаться следующим:
