Методы решения СЛАУ в задачах вычислительной динамики Текст научной статьи по специальности «Математика»

Литература

1. Кручинин В,В, Разработка компьютерных учебных программ. Tcí'-;k, 1998.

2. Поликахин А.В. Гипертекст: сущность, состояние, проблемы, перспективы / А,В. Поликахин, А.Ю, Совин. М„ 1993,

3. Агеев В.Н. Примеры гипертекстовых и гипермедиа систем (обзор) II Компьютерные технологии в высшем образовании. М., 1994,

4. Palmer W. Agnew, Anne S. Keiierman, Jeanine Meyer. Multimedia in the Classroom Aliyn & Bacon; 1 edition, 1998,

5. Сваровский И.Н. Анализ технических средств для дистанционного обучения / И,Н. Сваровский, О.Б. Журавлева II Тез. докл. научно-методической конференции «Дистанционное образование. Состояние, проблемы, перспективы». Томск, 1997,

6. Park 0. Hypermedia: Functional Features and Research Issues II Educational Technology. 1991. August,

7. Бондарь B.A. Проблемы приема экзаменов в дистанционном образовании и новые подходы к их решению / В.А, Бондарь, 0,Ю. Исакова, В,В. Кручинин II Тез, док. Международной методической конф, «Новые информационные технологии в университетском образовании!!, Кемерово, 2002.

8. Исакова О.Ю, Основные направления совершенствования контроля знаний в Томском межвузовском центре дистанционного образования I О.Ю. Исакова, В.В, Кручинин II Единая образовательная информационная среда: проблемы и пути развития: Материалы II Всероссийской научно-практической конф. Томск, 2003,

9. Исакова О.Ю. Проблемы контроля знаний в дистанционной технологии и пути их решения / О.Ю. Исакова, В.В, Кручинин II Международный конгресс конференции «Информационные технологии в образовании»: Сбооник трудов участников конференции. Часть V. М., 2003.

10. Исакова О.Ю, Опыт организации контроля знаний в Томском межвузовском центре дистанционного образования / О.Ю, Исакова, В.В. Кручинин, А.Ф. Уваров II Тез, докл. Международной научно-методической конф. «Инновационные технологии организации обучения в техническом вузе: на пути к новому качеству образования», Пенза, 2004.

11. Исакова О.Ю, Автоматизация синтеза вопросов в компьютерных учебных программах I О.Ю. Исакова, В,В, Кручинин //Доклады III научно-практической конф. «Современные средства и системы автоматизации - гарантия высокой эффективности производства», Томск, 2002

УДК 519.612:537.8

С. П. Куксенко, Т.Р, Газизов

МЕТОДЫ РЕШЕНИЯ СЛАУ В ЗАДАЧАХ ВЫЧИСЛИТЕЛЬНОЙ ЭЛЕКТРОДИНАМИКИ

Томский государственный университет систем управления и радиоэлектроники

Введение

Необходимость в решении системы линейных алгебраических уравнений (СЛАУ) возникает при использовании широкого класса моделей и подходов, используемых при автоматизирован ном проектировании радиоэлектронной аппаратуры. В частности, решение задач излучения или рассеяния электромагнитной волны сложными объектами, являющихся одними из основных задач электромагнитной теории, может быть получено с помощью интегральных уравнений, сводящихся методом моментов к СЛАУ с плотными матрицами [1, с. 5-19].

При компьютерном моделировании основные вычислительные затраты состоят из суммы затрат на формирование матрицы и затрат непосредственно на решение СЛАУ. Следовательно. выбор наиболее эффективного способа решения СЛАУ позволит снизить общие временные затраты. Следует отметать, что

в большинстве публикаций по решению СЛАУ большого порядка рассматриваются разреженные матрицы. Поэтому разработка и исследование новых математических методов и подходов для решения СЛАУ с плотными матрицами весьма актуальны. В известных авторам публикациях, [2, р. 158-178; 3-4] не приведено сравнение имеющихся методов решения СЛАУ в приложении к вышеописанным задачам.

Цель данной работы ~ сравнение методов решения СЛАУ.

Рассматривается вопрос о решении системы вида

а ! | х! + а 1 з т 2 + ■ • ’ ~г а \ д< х ы Ъ: ’

а-цХI аггХг + • •• + атнХи ~~ Ьг>

, (])

. а N і х і+ а N і х і+ ■■■+ а а* х >■ “ Ь ы * или иначе, векторно-матричного уравнения Ах = Ь, (1а)

где Ь - вектор свободных членов, х - вектор неизвестных (вектор-решение) размера М,

А - МхАт матрица коэффициентов данной системы, Эффективность способа решения системы (1), с заданной точностью, во многом зависит от структуры матрицы А: размера, обусловленности, симметричности, специфики расположения ненулевых элементов в матрице и др.

Все методы решения СЛАУ можно разбить на два класса: прямые и итерационные. Прямыми называются методы, которые приводят к решению за конечное число арифметических операций. Если операции реализуются точно, то и решение будет точным. Итерационными являются методы, в которых точное решение может быть получено лишь в результате бесконечного повторения единообразных действий.

Точные методы Метод Гаусса

Наиболее известным и популярным способом решения линейных систем вида (1) является метод последовательного исключения неизвестных -метод Гаусса (СЕ). Систему (1) приводят к треугольному виду последовательно, исключая сначала х\ из второго, третьего, ..., Л^-го уравнений, затем х~1 из третьего, четвертого, .... ДТ'-го уравнений преобразованной системы, и т.д. Па первом этапе заменяют второе, третье, /У-е уравнение на уравнения, получающиеся сложением этих уравнений с первым, умноженным соот-

вететвенно на . Еи _£и яу> . Результатом

аи 1./у

этого этапа преобразований будет эквивалентная (1) система. На втором этапе проделываются такие же операции, как и на первом, с подсистемой полученной системы, без учета первого уравнения. Продолжая этот процесс, на (Л-1)-м этапе так называемого прямого хода метода Гаусса данную систему (1) приводят к треугольному виду. Очевидно, что треугольная структура системы позволяет последовательно одно за другим вычислять значения неизвестных, начиная с последнего. Этот процесс последовательного вычисления значений неизвестных называют обратным ходом метода Гаусса.

Метод Ш-разложения матриц (Ш)

Пусть А = («Доданная Агх/У-матрица, а

соответственно нижняя (левая) и верхняя (правая) треугольные матрицы. Справедливо следующее утверждение. Если все главные миноры квадратной матрицы А отличны от нуля, то существуют такие нижняя Ь и верхняя и треугольные матрицы, что А = Ш, Если элементы диагонали одной из матриц 1Ь или II фиксированы (ненулевые), то такое разложение единственно. Обычно разложение матриц осуществляют при фиксировании диагонали нижней треугольной матрицы (/,; = 1 при /' -- /'). Далее находят I,у при I > у (/,, = 0 при / <у) и и,, при / < /' (м,/= 0 при I >/').

Если матрица А исходной системы (1) разложена в произведение треугольных Ь и 11, то вместо (1а) можно записать эквивалентное уравнение Шх = Ь.

Введя вектор вспомогательных переменных у, последнее выражение можно переписать

(Ьу = Ь,

в виде системы <

(их = у.

Таким образом, решение исходной системы с квадратной матрицей коэффициентов свелось к последовательному решению двух систем с тре-у гол ьны м и м ат р и цам и козф фициен то в.

Итерационные методы

Точные методы понятны и просты для программной реализации, однако, вычислительные затраты этих методов, которые пропорциональны Аг3, серьезно ограничивают круг решаемых задач. Поэтом)' в последнее время широко применяются итерационные методы. Их вычислительные затраты пропорциональны N¿,■N1, если для сходимости потребуется Ы„ итераций.

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

Рассмотрим систему (1) и сформируем для нее следующую задачу. Пусть заданы некоторые два подпространства Кг— и Е~ .

Требуется найти такой вектор х^К, который обеспечивал бы решение исходной системы, «оптимальное относительно подпространства Ь», то есть чтобы выполнялось условие \/ 1е I,: (Аж, /) = (Ь, /), называемое условием Петрова-Галеркина. Сгруппировав обе части равенства по свойствам скалярного произведения и заметив, что Ь-Ах = гх, это условие можно переписать в виде V /е Ь: (гх, /) = 0, то есть гх=Ь-Ах 1 Такая задача называется задачей проецирования х на подпространство К ортогонально к подпространству I.

При построении и реализации проекционных методов важную роль играют так называемые подпространства Крылова, часто выбираемые в качестве К. Подпространством Крылова размерности т, порожденным вектором V и матрицей А, называется линейное пространство К„(V, А) = зрап{у, Ау, А2у,---,А">~1у}.

В качестве вектора V обычно выбирается невязка начального приближения го; тогда выбор подпространства X и способ построения базисов подпространств полностью определяет вычислительную схему метода [3, 4].

Предобусловливание

Вычисление итерационными методами зависит от обусловленности матрицы А, оцениваемой числом обусловленности С0ПС1А = ||А!Н|А^|Н|/ц„от|/|Хм/;,|.

С ростом сопс!А обусловленность ухудшается, и для ряда проблем сходимость может оказаться очень медленной и поэтому итерационный процесс может застопорится или даже оборваться,

О «гибельном» влиянии числа обусловленное™ на решение СЛАУ точными методами можно узнать из [6, р. 847-853]. Однако применение так называемого предобусловливания улучшает сходимость к требуемому решению.

Пусть М - некоторая невырожденная матрица размера N. Домножив (1) на матрицу М"1, получим систему

МГ'Аж = ТУГ1!}, (2)

которая в силу невырожденности М имеет то же точное решение х„ Введя обозначения А=М_|А и Ь == М 'ь, запишем (2) в виде А х = Ь . (3)

Хотя (3) алгебраически эквивалентна (1), спектральные характеристики матрицы А отличаются от характеристик матрицы А, что, вообще говоря, ведет к изменению скорости сходимости решения (3) итерационными методами по сравнению с (1) в конечной арифметике.

Процесс перехода от (1) к (3) с целью улучшения характеристик, матрицы для ускорения сходимости к решению называется предобу-словливанием, а матрица М-1- матрицей предо-бусловливателя. Из (2) сразу же вытекает важное требование: матрица М должна быть близка к матрице А (чтобы сопс! А -> 1). Действительно, выбор М = А сразу же приводит (1) к виду 1х = А' ‘Ь (I - единичная матрица), однако не имеет практического смысла, так как требует нахождения А’1, что по существу .сводится к решению (1). Вторым естественным требованием является требование легкой вычислимости матрицы М.

Невязкой системы (I), соответствующей вектору х, называется вектор г = Ь - Ах = А(х,- х), где х, - точное решение. Невязка г системы (3) связана с невязкой системы (1) очевидным соотношением

М г = г, (4)

которое справедливо и для матрично-векторных произведений г = Ац и г = А ц:

г = А ч = М.Л'Ц! => М г = Ац = г. (5)

Это позволяет вместо явного перехода от (Г) к (3) вводить в схемы методов корректирующие шаги для учета влияния предобусловливающей матрицы (в некоторых случаях, когда матрица М имеет простую структуру, переход от А и Ь

к А и Ь может выполняться и явно).

Из (5) следует еще одно условие: структура матрицы предобусловливателя должна допускать легкое и быстрое решение «обратных к предобусловливателю» систем вида М г = г.

Таким образом, матрица М должна быть: по возможности близка к матрице А, легко вычислима, легко обратима.

Для заданной матрицы А потребуем представить ее в виде А = Ш + К, где Т и II — как и прежде нижнее- и верхнетреугольные матрицы, соответственно, а матрица К является матрицей ошибки. Тогда приближенное представление А«Ш называется неполным ИТ разложением матрицы А или коротко ее 1ЬП-разложением. Таким образом, при использовании в качестве матрицы М неполного ЬП-разложения матрицы А получим выполнение вышеописанных требований.

Для систем с разреженными матрицами хорошим предобусловливанием служит НАДО) [3]. Оно заключается в применении ЬП-разложения к матрице А, но если а(/= 0, то полагается /,,= О или и-ч= 0. Поскольку в задачах электромагнитной совместимости матрица А является плотной, то применение 1Ы)(0) не целесообразно.

Но по аналогий с ILU(O) матрицу М предлагается формировать так: перед преобразованием каждой строки вычислить ее норму, затем' умножением значения нормы на задаваемый допуск обнуления т (drop tolerance) получить значение порога (threshold value). Далее проверить, если какой-либо элемент преобразуемой строки меньше порога, то этот элемент положить нулем. После проверки всей строки применить к ней преобразования из ILU(O), После преобразования всех строк получают разреженную матрицу М такую, что матрица М"!А обусловлена лучше, чем матрица А.

Далее приведен описанный выше алгоритм. ДляN w = а№

Вычислить wvv = ||wj|9-x Для s - 7V

Если w.v < ww

то w.v = О

Увеличить s Для к = 1,..., i - 1 и wk ■/- О wk= wk/ukk Для j = к + N

W = W - Wk- lib.

Увеличить j Увеличить к hj= wh для у = 1, i-l

щ= wj, для j = /, ...,N

Увеличить /

Далее приведены, пожалуй, три самых популярных на данный день итерационных метода (крыловского типа) различающихся способом построения базиса в пространстве Крылова. Необходимые пояснения можно найти в [3-5; 7; 8, р. 129-142]. Параметр ТЫ, использованный в алгоритмах, определяет точность итерационного процесса, a iV/'ax - максимальное задаваемое число итераций.

Метод обобщенных минимальных невязок (GMRES(m)) с вращениями Гивенса

Построить матрицу нредобусловливателя М

Выбрать m < N, произвольное начальное прибли-

max

жение Хо и Nj, г~Ъ..Ахо

Найти го из системы Мго= г

P4WI, •

Создать ((/■« + 1) х т) - матрицу II

Создать вектор (от + 1) - вектор q Для N„ = 1,2,... до сходимости или до Н = 0; q = pe¡ v, = rc/p Для /= 1,..., от г = Av;

Найти w из системы Mw = г Для к= ],..., i

hk¡~ (w, v/() w = w — hkj vk

увеличить к

h¡ f i,/= liw||,,

V,, = W//?, v

Для k= i~ 1

(/?!,■, ..., h, + ij)' — G/,(h\j, .... hj-,. | ,У увеличить к

Построить G, так, что [G,-/?..,-],. i _ 0 4 = G,q увеличить i

Решить треугольную СЛАУ H, :: my == q¡.

m

xo=xo+ Zj/V,

/=1

2 = Ь - ÂX()

Найти r0 из системы Mrn = z

P = !Ni2

Продолжать пока p/jlbjl, > Toi

Основной проблемой в GMRES является выбор размерности т пространства К; хранение базисных векторов этого пространства определяет' основные затраты памяти при программной реализации метода, составляющие mN чисел. Наиболее распространенной модификацией GMRBS является перезапускаемый GMRES или GMRES(m). Выбирается некоторая размерность т < JV, определяемая, по существу, лишь доступным объемом памяти, И после обновления решения X = Х0+ У/пУ/д оно принимается за новое начальное приближение, после чего процесс повторяется.

Метод бисопряженных градиентов (BiCG) Построить матрицу нредобусловливателя М

;ПЙХ

Выбрать начальное приближение х<> и N„

Го = Ь-Ахо

Выбрать вектор r¡, удовлетворяющий

условию (r0, r¡ ) Ф 0 ( Гц = Го)

Найти 20 из системы Mz0- г0 Найти %{) из системы М'~ г0

р0= Щ

Р(, = 4

тах

Для] = 1, 2, ... до сходимости или до N¡1 а1~ (г/. Г/ У(Ар/, р* )

х/,, = х/+о/р/

Г/,1 = Гу- О/Ару

г‘и = г* - ауАгр*

Если ||^)||2/||Ь|12 < Го1

то КОНЕЦ (хун - полученное решение)

Иначе

Найти 2/+| из системы Мг/+| = Найти г/+1 из системы

М х г

т *\/+{ *_/+!

[3; = (1у+,/ Жу (, )(г;/ Жу)

Р/И = г/-г^ ¡3/Р/

Р*+1 = *-*+)+Эур’

увеличить/

Стабилизированный метод бисопряженных градиентов (В!СС31аЬ)

Построить матрицу предобусловливателя М

тах

Выбрать начальное приближение х0 и А7,-,

Го = Ь - Ах0

Выбрать вектор г , удовлетворяющий

условию (Г(Ь г ) * 0 (г = Го)

П • _ 1 /1 „г ,ПаХ

Для I =- 1,2,.. .до сходимости или до Ни

Им = (г*, Го)

Если I = I

Р/= Г,_,

Иначе

0ы = («/. 1/и/..2)(а,.,,/ю,..1)

Р/= Г/_. | + Р/-)(р;_1 — С.О,_|У,_))

* ,(с

Найти р из системы Мр = р,

V, = Ар

ос,- = о,ч/(г , у,)

8 = Г/..,, СХ/ V /

Если ||8Ц2/||Ь||2 < То!

то КОНЕЦ (ж,-= + а,- р -

полученное решение)

Иначе

% * Найти 8 из системы М* = 8

I = Ав

СО, = (Ч, 8>/(1, 1)

& * х|~х,_1 + а, р +о>,8 г, = я - ю,-1

увеличить г

Вычислительный эксперимент

Система (1) решалась вышеописанными методами (СЕ, Ш, В|СО, В1СС381аЬ, СМПЕ8(т)) на примере определения токов в проводной антенне (ТУ = 721). При вычислениях полагалось, что Го/ = 1(Г8, Л7,"1“ - 500. Результаты приведены в табл. I, где х - параметр, используемый для построения матрицы предобусловливателя, Т - время решения СЛАУ, Т/и/ -время неполного Ьи - разложения, Ы-, - число итераций, а столбец «Решение» показывает: совпадает ли полученное решение соответствующим методом с решением, полученным методом Г аусса.

Таблица 1

Результаты работы разных методов решения СЛА У

Метод т Г, с Тку, с «я Решение

вЕ _ 30,26 п _ - 4-

Ш - 30,81 - _ 4-

- 135,12 _ 400 4~

5-10-г 14,89 0,14 37 4-

10-2 10,88 0,28 27 +

В|СС 5-10-3 9,06 0,39 22 4-

Ю-з 5,77 0,44 13 +

5.10-4 5,32 0,66 11 4-

Ю-4 7,80 4,06 8 4-

5-10-5 16,10 12,25 7 4-

131,33 - 500 -

5-10-2 11,59 0,14 40 4-

10-2 7,25 0,28 24 4-

В!СС51аЬ 5-10-3 7,91 0,39 26 4-

Ю-з 3,74 0,44 11 4"

5-10-'* 3,35 0,66 9 4-

Ю-4 6,31 4,06 6 +

5 -10-а 14,77 12,25 5 4-

- 912,15 - 500 -

5-10-2 284,40 0,14 141 +

10-2 57,78 0,28 29 4-

01ШЕ5(1О) 5-10-5 36,85 0,39 18 +

10-5 4,78 0,44 2

5-10-1 4,72 0,66 2 4-

Ю-4 6,43 4,06 1 +

5-10-5 14,94 12,25 1 4-

1926,40 - 500 -

5-Ю-2 49,38 0,14 12 +

10"2 12,74 0,28 3 +

GMRES(20) 510-3 8,73 0,39 2 +

10-3 4,78 0,44 1 4*

5-10“5 5,00 0,68 1

10-4 8,84 4,08 1 н-

5-10-s 17,79 12,25 1 •+-

5-10-2 18,73 0,14 3

GMRES(30) 102 6,70 0,28 1 +

5-Ю-3 6,76 0,39 1 -ь

Из табл. 1 видно, что из всех рассмотренных методов выгоднее использовать итерационные методы с предобуеловливаиием при оптималь-

ном значении параметра т. Так, требуемое решение методом В1СС81аЬ при т = 5- КГ4 было получено в 9 раз быстрее (наибольший выигрыш), чем при использовании метода Гаусса. Использование метода ЕИСС^аЬ позволяет получить результат в 1,6 и 1,4 раза быстрее (при оптимальных значениях т), чем при использовании методов ЕНСС и СМК.Е5(т), соответственно.

Как видно из таблицы, для каждого итерационного метода существует оптимальное значение параметра т. Исследования по выяснению зависимости оптимального значения т от матрицы (структуры объекта), а также пояснения наличия этого факта можно найти в [9, с. 26-28]. В работе [10, с. 54-57] показано, что можно получить ускорение решения СЛАУ итерационными методами за счет снижения точности вычисления (числа итераций) до достаточной для получения заданных характеристик (диаграммы направленности антенны).

Литература

1, Харрингтон Р.Ф. Применение матричных методов к задачам теории поля IIТИЭЭР. 1967. №2.

2, Lee J„ Zhang J., Lu С. Incomplete LU preconditioning for scale dense complex linear systems from electromagnetic wave scattering problems. Computational Physics, 2003. Vol. 185.

3, Saad Y. Iterative methods for sparse linear systems. PWS Publishing Company. 1996.

4, Баландин М.Ю., Шурина Э.П. Методы решения СЛАУ большой размерности. Новосибирск, 2000.

5, Ильин В.П, Методы неполной факторизации для решения линейных систем. М., 1995,

6, Тарап К. Sarkar at al. Survey of Numerical Methods for Solution of Large Systems of Linear Equations for Electromagnetic Field Problems, IEEE Trans, on Antennas and Propagat, Vol. AP-29. Nov, 1981, № 6.

7, http://www.cerfacs.fr/algor/softs/

8, Guennouni A. EL., Jbilou K„ Sadok H, A block version of BICGSTAB for linear systems with multiple right-hand sides. Electronic Transactions on Numerical Analysis. 2003, Vol. 16.

9, Газизов T.P., Куксенко С.П. Оптимизация допуска обнуления при решении СЛАУ итерационными методами с предобуеловливаиием в задачах вычислительной электродинамики II Электромагнитные волны и электронные системы, 2004, № 8,

10, Куксенко С.П., Газизов Т.Р. Ускорение решения СЛАУ в задачах вычислительной электродинамики: Сборник научных трудов VII Всероссийской научно-практической конф, «Проблемы информационной безопасности государства, общества и личности (16-18 февраля 2005)». Томск, 2005.

УДК 512.643:512.644

И.С. Костарев, С.П. Куксенко, Т.Р. Газизов

ПОВЫШЕНИЕ ЭФФЕКТИВНОСТИ РЕШЕНИЯ СИСТЕМЫ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ ИТЕРАЦИОННЫМИ МЕТОДАМИ

Томский государственный университет систем управления и радиоэлектроники

Введение

Необходимость в решении системы линейных алгебраических уравнений (СЛАУ) возни-

кает при использовании широкого класса моделей и подходов, используемых при автоматизированном проектировании радиоэлектронной аппаратуры. В частности, решение задач излучения