Самарский А.А. Гулин А.В. - Численные методы (1078412), страница 76
Текст из файла (страница 76)
Тогда получим 91-х+2р,— С(у~ ~+у ~~)+уыз= — (г,-,+сыч), откуда, учитывая, что и,, + у„, = Сн,— Е„ придем к уравнению й-з — (С вЂ” 2Е) у~+ уме = — (Ез, + СЕ;+ Е;,), (3) связывающему значения искомого вектора в узлах одинаковой четности. В частности, если 1 — четные, то проведено исключение нечетных узлов. Далее этот процесс исключения можно прадолжнть /с=О 0 1 ' 5 С 5 о 7 о .910111а151а151о Л=1 и . С Л Е ГП ГЯ 1 1Л Рнс, 16. Порядок исключения неизвестных в методе редукции аналогичным образом. При этом необходимо предположить, что число узлов М является степенью двойки, Л'=2'". Прежде чем переходить к случаю произвольного лт, рассмотрим для наглядности случай лт= 4, т. е.
У =16. Обозначим через я номер атака исключения неизвестных. ))ри я=и система уравнений совпадает с исходной и содержат значения неизвестных во всех внутренних узлах. На рис. 16 это соответствует верхней горизонтальной черте, где кружочками отмечены номера неизвестных ун входящих в систему.
На следующем этапе 41В у;, — С (')у; + у;, = — Йо, ( = 2, 4, 8,..., 2 — 2, (4) где С(о = (С(н~) и 2В С(н) = С У(з) =У, „+СУ(+Г(„. (5) (6) По индукции легко доказать, что на й-и этапе исключения, Й=1, 2,..., т, получаем систему У, ...— С У(+У„...= — У( (и-() й-н) (7) (=2' ', 3-2н ',..., 2'" — 2" ', Ун = ((1 Ук = рн где матрицы С'"-" и векторы Р( ') находятся из рекуррентных со- отношений С"'=С, у(н) 7((н-н), ~ (нэо)7((н-н), н(н-1) л(н) — л. 1=2н 2.2( 3.2н 2н — 2н (8) (9) Таким образом, весь процесс решения состоит из прямого хода и обратного хода. Прямой ход заключается в нахождении матриц С") и векторов Р(" по формулам (8), (9).
Обратный ход состоит в нахождении векторов у, нз системы (7), начиная с й=т. Метод редукции в том виде, как он здесь изложен, не применяется в реальных вычислениях по двум причинам. Во-первых, он не- экономичен из-за того, что на каждом этапе приходится обращать матрицу Спч общей структуры. Во-вторых, вычисление правых частей по формулам (9) неустойчиво.
В следующих пунктах будет показано, как можно устранить указанные недостатки метода редукции. 420 (1=1) происходит исключение неизвестных с нечетными номерами, в результате чего получаем систему вида (3), содергкащую значения неизвестных только в четных узлах. Этап й=1 изображен на рис. 16 второй сверху горизонтальной чертой. Стрелки указывают, какие неизвестные были исключены. На втором этапе ()(=2) остается каждый четвертый узел и на заключительном этапе (А=З) остается только одно уравнение, связывающее у„у„у„. Поскольку у, и у„заданы (см.
(2)), из последнего уравнения можно найти у,. Тем самым начинает осуществляться обратный ход в методе исключения. Зная у„можно найти у, и у„, далее — все неизвестные с четными номерами и, наконец, все остальные неизвестные. Вернемся к общему случаю, когда Л)=2'". Согласно (3), в результате первого этапа исключения (1=1) получаем систему урав- нений 2. Обращение матриц. Соотношение (8) представляет собой нелинейное разностное уравнение первого порядка для матриц. Его решение можно найти в янном виде. Рассмотрим сначала числовой аналог уравнения (8), а именно разностное уравнение У»=У»,— 2, й=1,2,..., у,=х, (10) где х — заданное число.
Покажем, что решение у„=у„(х) уравнения (1О) выра>кается через многочлен Чебышева, а именно у» (х) = 2Т» ( — ~, А = 1, 2,..., (11) где соз(аагссозх), если )х! (1, Т„(х) = — 1(х+1/х' — 1)" +(х — 1/х' — 1)"1, если )х) «1. 2 Из выражения для Т„(х) следует, что Т,„(х) =2(Т„(х))» — 1. Отсюда при а =2'-' получаем Т,,( — ') =2 (»Т... ( — '))' — 1, т.
е. функция (1! ) удовлетворяет уравнению (10). Корнями многочлена (11) являются числа хс»=2 соэ „„—, 1=1,2,..., 2. (2! — 1) н Поэтому функция у„(х), представляющая собой многочлен относи- тельно х степени 2' со старшим коэффициентом 1, разлагается в произведение »» (21 — 1) к» у»= Ц ~х — 2 соз 1=» Приведенные выше выводы имеют место и для матричного уравнения (8). По индукции легко доказывается, что решением Снч уравнения (8) является многочлен (относительно матрицы С) степени 2" со старшим коэффициентом 1, Для этого многочлена справедливо представление, аналогичное (!1), т. е.
С'»' = 2Т,» ~ — ), и справедливо разложение на линейные множители »» >=» (12) 42! Наличие такого разложения позволяет упростить провед)ру обращения матриц. В случае разностных аппроксимаций уравнений эллиптического типа матрица С является трехдиагоиальной, в то время как Сич — матрицы общей структуры.
Благодаря разло кению (12) обращение матрицы С'"' сводится к последовательному обращению трехднагональных матриц Сзн=С вЂ” 2соэ з Е, 1=1,2,...,2. (13) зз Действительно, пусть требуется решить уравнение Сзмо=зг, (14) где з См'= з1Сь . !=з Рассмотрим для наглядности сначала случай, когда 1=2. Тогда (14) принимает вид С,,С„С,,С,,о= (15) Обозначим о,=гч о,=С,,СззС,зо, о,=С,зС, о, о,=С„о, оз=о. Тогда получим, что решение системы (15) сводится к последовательному решению четырех систем уравнений Сз зоз=аз, Сззоз=аз, Сззоз=оь Сз,оз=оз, где о,=з(х о,=о.
Точно так же решение системы (14) в общем случае сводится к последовательному решенню систем уравнений С,,о,=о, „1=1, 2, ..., 2", (16) где оз=йх ом =о. Если матрицы С„, трехдиагональные, то кал дую нз систем (16) можно решить методом прогонки. Указанный вьцпе способ обращения матрицы Сзм предполагает определенный порядок выполнения промежуточных этапов: сначала надо обратить матрицу С,„затем — матрицу С„л и т. д. Однако, поскольку в матрице С'"' все сомножители перестановочны, можно использовать и другие способы введения промежуточных значений о„1=1, 2,..., 2' — 1. Так, систему (15) можно записать в виде С,,С,,С,,С,, = р и заменить системой уравнений Сззоз=оз Сззоз=из, Сззоз = оз, Сззо з = оз где оз=зр, о,=о.
Теоретически прн любом порядке выполнения промежуточных этапов мы должны получить одно и то же решение о. Однако, если число промежуточных этапов велико, то при ре:пении систем уравнений (16) на ЭВМ будет происходить накопление погрешностей округления. Рост погрешностей округления зависит от порядка выполнения промежуточных этапов. Здесь имеет место 422 примерно та же ситуация, что н в итерационном методе с чебышевским набором параметров (см. п. 2 $ б гл. 2 ч. П). Поэтому прн реальном решении системы (14) рекомендуется обращать внимание на порядок выполнения промежуточных этапов. Более подробно этот вопрос рассмотрен в книге (35).
3. Вычисление правых частей. В методе редукции правые части г'," " уравнения (7) должны удовлетворять соотношению (9). Будем искать г( в виде (») (17) где векторы р( ), ()~~ подлежат определению. Подставляя (17) в (9), получим , И) (Ь (М г (»-1) (» — И, (»-1) Р( 1 ((( = ' Р( »- 1 ()( »-1 1 +С" п(С" " '» и+ '" и) +С"" (~о + '.~п . Отсюда, учитывая, что С'"=(С'"-о)' — 2Е, придем к уравнению (с'»-о)2(!() — Р(»-")+ '," =2Р',"+ '.—",, +,".—,и + 1-2 1 (22 +С'"п(Р('-,",, + Р',",—,,", + Е( ')).
Выберем теперь вектор ()( таким, чтобы выполнялось соотно(») шение а, =2Р +9,,» +9,„1»ги (») (») (»-1) (»-1) (18) Тогда предыдущее уравнение после умножения на (С'"-о)-' примет внд +;„.,"- + (19) = Р(»-1+ Р(„»-1+(1( где а (А — О (А> (А-!) Тогда придем к уравнению (»-)) (»-н 8-и С 11 =д( +р,. „,+д„,, „ (20) где 1( =- у( — р; . Из этого уравнения, обрацая матрицу С( (»-1) (»-1) (»-1) находим 1(» " и затем вычисляем р(= р(» о+ 1';)' ".
422 Таким образом, в,(числение векторов Р,'») можно заменить нахо>кдением векторов р(,()( из система уравнений (17) — (19). Сначала, (») и) обращая матрицу С" ', находим вектор а(' ", затем в»числяем р';"= = р,'» " + а(» и и затем по формуле (18) начисляем ()('. Правую часть Р)(~ " уравнения (7) надо заменить согласно (17) в (ражением 2 Ф-и (»-1) (»-1) (»-1) '( ='- Р( + ((( 4.
Формулировка и обсуждение алгоритма. Сформулируем более детально алгоритм метода редукции. Вычисления ведутся в цикле по индексу й. Сначала осуществляется прямой ход, т. е. для Й=1, 2, ..., т — 1 решаются уравнечия С()1-1) и-1) (4-1) '1, м-1), и-1) = Р;,1- т Р(„(- +()) (21г и находятся векторы И) (1 — 1) (4-1) Р = р;)-5; (М М-1) И 1) (), =2р) +д,,л +4, 11-1, где )=2", 2.2', 3 2', ..., 2" — 2". Вычисления ведутся, начиная с 1=1, причем при я=! задаются начальные значения )л) б (1) Р Решение систем (21) сводится, как это было показано в п.
2, к ре(пению более простых систем уравнений Сл 1,(п(,(= и( 1 ь 1= 1, 2,, 2 (22г ) = 2', 2 2',..., 2 — 2', где М-1) М-П (Л-1) пал = Р(,1-1 + Р(, 4-1 + Ч1 (Л 1! П(1 1,.=5; 11 Системы (22) решаются при фиксированном 1 для каждого 1=2', 2 2', ..., 2'" — 2'. Поскольку матр ща С... не зависит от ), вычисления целесообразно организовать так, чтобы матрица С... обращалась один раз. Подсчитаем число действий, необходимых для решения всех систем вида (22). Предположим, что для решения системы (22) с фиксированными я, 1, ( гоебуется () арифметических действий.
Тогда при фиксированных й, 1для решения систем (22) с (=2', 2.2',... — 1 ..., 2'" — 2' требуется (1=(2 — 1)() арифметических дей- 2 етний. При фиксированном й и для 1=1, 2, ..., 2' ', 1=2", 2 2", ... ..., 2'" — 2', требуется 2' '(2"' ' — 1)()= (2" ' — 2' 1)() действ й. Наконец, для всех и= 1, 2,..., т — 1 необходимо выполнить 1й — 1 ()'Я(2 ' — 2' ')=(! +(т — 2) 2'" 1)(1 Л'=1 арифметических действий.
Самое существенное здесь то, что при больших Х число действий является величиной 0(()М!оц,Л(). 1М 1) После того как все векторы р) ', д) найдены, осуществляется обратный ход метода редукции, т. е. начиная с я=т решаются 424 уравнения (23) и вычисляются искомые векторы Й-о и-в у'=Р +1 где )г = т, т — 1,..., 1, 1= 2' ', 3 2'-',..., Л' — 2'-', Л'= 2". Система (23) при фиксированных 1, я решается, как и ранее, путем последовательного решения систем уравнений С. с.шс,=ш, сь 1=1, 2,..., 2" ', (24) 1=2' ', 3 2' ' 5 2" ',..., 2"" — 2'-', ) %-1) Й-о где ."агь =у + у, е- + у„,,ь-ь ш,ь,,= 1,, Решение системы (23) также требует 0(дЛГ!оп, Л') действий. Рассмотрим применение метода редукции к решению разностного уравнения Пуассона (4), (5) из $5. Предположим, что сетка содержит Л~,=2"' точек по направлению х, и Л', точек — по х,.