Федоренко Р.П. Введение в вычислительную физику (1185915), страница 52
Текст из файла (страница 52)
Итак, опишем многообразие Я(!) в форме Ж) х() =И ). (12) где б есть матрица и- 1В, р — й-вектор. При этом мы должны позаботиться о том, чтобы описание (12) было эквивалентно описа- Повторим (в матричной форме) вывод уравнений для Ь(!) и Ь(С). Дифференцируя по ! соотношение 2.х = Ь, получаем е,х + 2,х = Ь. За- мена х на Ах + а дает (з. + 2,А)х + Еа = Ь. Примем для з. и Ь урав- нения В гв1 жесткие лииейные кРАеВые зАдАчи нию (9), а строки С(г) образовывали хорошо обусловленный базис. Лучше всего, чтобы они были взаимно-ортогональными.
(Требования к 13 сформулируем несколько позже.) Будем искать с(1) в форме с(г) = м(г) Й(г), где м(г)— матрица х -Р х, пока ие определенная. при любой матрице М(1) строки С суть некоторые линейные комбинации строк Е. Если М вЂ” невырожденная матрица (т.е. из Мг = 0 следует г = О), то линейные подпространства, натянутые иа строки з. и С, совпадают, или, другими словами, з.
и С отображают л-мерное подпространство в одно и то же х-мерное, определенное лишь разными базисами. Требование равномерной по г хорошей обусловленности базиса из строк С можно оформить в виде требования постоянства матрицы С(г) С (г). Это — матрица к- х, элементами которой являются всевозможные скалярные произведения строк С. Предполагая, что строки 1,, входящие в левое краевое условие, предварительно ортонормированы, будем считать, что М(0) = Е (так будем обозначать единичную матрицу /с- х).
Постараемся определить М(г) таким образом, чтобы СС' при всех г оставалась единичной, т,е, чтобы строки С(г) образовывали ортонормированный базис. Для этого нужно обеспечить равенство — ",(СС') = СС*+ СС' = СС'+ (СС')* = 0. (1З) Из (10) имеем С = (МЦ, = МС + Мт. = М~.
— МаА. Используя выражение для С, а также соотношение Ь = М 'С, преобразуем СС'. СС = (МА'. — МА,А)С' = (ММ 'С вЂ” ММ 'СА)С' = = ММ 'СС' — САС'. Определим М таким образом, чтобы это выражение было равно нулю. (Второй член (13) при этом тоже автоматически обратится в нуль, и будет обеспечено С(г)С'(1) = сопзг.) Итак, после несложных формальных преобразований для М получаем задачу Коши: М = САС'(СС') ~М, М(0) = Е. (14) Это уравнение решать не придется. Мы используем его при выводе уравнения для С: С = (МС),= МЕ, + М)- = МЕ, — МЕА = САС'(СС') М1. — МСА. лгивлижвнныв методы вычислительной ьизики сч. п И наконец, учитывая, что МЕ=б, бб' = Е, получаем С= САС'Π— бА, 0(0) = ь(0). (15) Теперь приступим к выводу уравнения для р(с).
Из (9) имеем ь(с) х(с) = Ь(с), где Ь(с) — решение задачи Коши (10). Умножение на М дает мс. х(с) = м ь(с), т.е. 0(с) х(с) = м(с) ь(с). Однако брать р(с) = м(с) ь(с) нельзя, так как мы не собираемся вычислять М, Е, Ь. Выведем уравнение для р(С), используя то соотношение, которое мы хотим получить: 0(С) х(с) = р(С), где С(С) — уже известная матрица и-~ Ь, х(с) — решение краевой задачи. Дифференцируя по с, имеем р= Сх+ Сх = (С+ СА)х+ Оа. Учитывая (15) и связь р = Ох, преобразуем первый член: (О + ОА) х = СА6 Ох = ОАО'р. Окончательно для р(с) получаем задачу Коши: (16) р = ОАО" р + ба, р(0) = Ь (0), Отметим важное обстоятельство.
Так как р(с) =б(с) х(с), а строки 0(с) образуют ортонормированную (хотя и не полную в лмерном пространстве) систему, то р(С) — величина того же поридка, что и х(с). Если задача вычислительно корректна и искомое решение х(с) ограничено в смысле (3), то р(с) — такая же величина. Но пока мы получили только часть уравнений для х(С), порожденную переносом (»прогонкой») левого краевого условия на весь интервал 10, Т1. Точно так же можно перенести правые краевые условия, интегрируя задачи Коши справа-налево.
Сами же уравнения для б+ и р~ по форме в точности совпадают с уравнениями для С и р (выше мы выводили их в общей форме (15), (16)). Итак, в каждой точке мы получаем систему л уравнений С-(с) х(с) = 0-(с), или 0(с) х(с) = р(с), 0+(с) х(с) = р" (с), где теперь 0(с) — матрица л- и, р(с) — л-вектор. Заметим, что Ц р(с) Ц = О(йх(с) Ц). Строки матрицы С, вообще говоря, ортонормированной системы не образуют, так как ее части 0 и С+ получены независимо. Можно исправить и этот недостаток, если сначала решить задачу для С (с) и при ортогонализации векторов 1ь входящих в прз- й 13! жесткие линеЙные кРАеВые 3АдАчи вые краевые условия, потребовать еще и ортогональности к строкам матрицы С (Т).
Если краевая задача не вырождена, это требование выполняется. В принципе при ортогонзлизации может возникнуть большая потеря точности, если пространства, натянутые на «правые» векторы 1,. (! = lс + 1, /с + 2, ..., И), и строки С (Т) образуют очень малый угол, хотя и остаются еще формально независимыми, т.е. дают в сумме все л-мерное пространство. Такая ситуация возникает тогда, когда исходная краевая задача почти вырождена, т.е. среди точек спектра однородной задачи х — Ах = Лх с однородными краевыми условиями (2) имеется точка, близкая к нулю. Если нуль принадлежит спектру, краевая задача вырождена, теряются стандартные свойства существования и единственности решения. Такую задачу называют «задачей на спектре» (мы здесь ее не рассматриваем).
Однако чем меньше по модулю ближайшее к нулю собственное значение спектральной задачи, тем хуже обусловлена исходная задача, тем больше постоянная С в оценке (4). Это почти очевидно. Если у(1) — собственная функция (!!у)! = 1), соответствующая малому собственному числу Л, то функция х(!) = х(г) + у(г) удовлетворяет уравнению х = Ах + а, а = а + Лу, !!х — х!! = О(! ), с малым возмущением 0(! Л)) правой части (х, очевидно, удовлетворяет невозмущенным краевым условиям), Таким образом, постоянная С в оценке (4) не может быть меньше 1/! Л!. Что касается обоснования алгоритма, то оно пока носит качественный характер: при его реализации мы не имеем дела с решениями типа ес', как это было бы, если бы мы попытались решать задачу обычным методом «фундаментальных решений» (см.
6 8). Однако полной ясности все-такн нет. Не очень ясен механизм преодоления влияния большого параметра. Он как был в исходной задаче, так н остался в уравнениях (!5), (! 6), в которых присутствует матрица А. В з 9 для простейшего случая было проведено достаточно полное исследование и, в частности, было обнаружено, что метод прогонки требует интегрирования задач Коши с большим параметром, но устойчивых.
Здесь, видимо, механизм носит несколько иной характер. Ниже в более, простой и прозрачной ситуации мы попробуем прояснить этот вопрос. Он состоит в следующем: как происходит накопление вычислительных погрешностей при численном решении задач Коши (15), (16) по стандартным схемам типа Рунге — Кутты, например? Ведь если ориентироваться на общие оценки, то для (15) имеем !!(САС С)о!! !)САС'!! ж !)А!! ПРИбЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ !ч. и (это тривиальная оценка, не учитывающая возможных тонких «компенсаций» больших величин при вычислении правой части (15)).
Как было сказано выше, выбор шага численного интегрирования т из условия !!А)!т.ес! считается здесь вполне приемлемым, Однако почему при оценке погрешности численного интегрирования не возникает стандартной и в общем случае неулучшаемой величины тге! «! г (р — порядок аппроксимации), не очень понятно. Тригонометрическая прогонка. Рассмотрим часто встречающуюся в приложениях задачу Штурма — Лиувилля — р(г) ф + д(г) у(г)+у(г) =О, оагкт, (17) с краевыми условиями, которые запишем в удобной для дальнейшего форме: у(0) соз а, + у(0) в!и а,=Ь,, !=О, у(Т) сов аз+ у(Т) в!и а = Ьз, Г = Т, (18) где р(г), а(г), У(г), а,, аз, ь„ьп т — заданные функции и числа. запишем (17), (18) в стандартной форме — в виде системы двух уравнений первого порядка; суть дела от этого не меняется.
Обозначая х, = у, х = ру (разумеется, предполагается р(г) ж р > О), получаем систему 0О + (19) Краевые условия имеют стандартную форму с векторами 1, = (сов ап з!и а,). Алгоритм прогонки состоит в том, что соотношение, имекяпее форму краевого условия, продолжается (в силу уравнения) на весь интервал !О, Т1 Рассмотрим «прогоночное соотношение слева»ч х,(г) соз р(е) + хз(г) в!и р(!) = р(г). (20) х, соз р — х, в!п р ф+ х з!и р + х соз р ф = р. Преобразуем это выражение, используя (20) и исключая хы хьс (1/р)ха соз р — х, з!и у ф + дх, в!п у + хз сов у ф = р — У в!п ~р.
ДлЯ Р(!) и Р(г) имеем Данные Коши: Р(0) = ап Р(0) = Ьи Выве- дем уравнения для них, дифференцируя (20) по ! и заменяя произ- водные от х из (19): жесткие линейные кгАевые з»д»чи в !в) Приводя подобные члены, имеем х1 яп ф (в — ф) + хз сов ф (1/р+ р) = р — / з1п ф. ( ) Умножая выражение (21) на я!п р и вычитая из него (20), умно- женное на (1/р+ ф) сов ~р, исключим члены с хз.' х,[(я —, р) зш р — (1/р+ р) соз р[— = Д в!и р — /3!пз р — Щ1/р+ ф) соз р. Оно выполняется при любом х„если потребовать одновременного обращения в нуль коэффициента при х1 и свободного члена.
В итоге мы получаем уравнения для ф и р: ф= ввйпз р — (1/р) совз р, р(0) = ан (22) Д= ~(1/р+ ф) созз ~р/3!пз <р+/з!и ~р, ЯО) бн (23) Деление на 3!и р не приводит к неприятностям, так как из (22) следует (1/р+ р) = (д + 1/р) щп ~р. Уравнение для р лучше использовать в виде (24) р = / ей и р + р з!и р соз р ( д + 1/р). Заметим, что нас интересуют задачи с большим параметром. Большой величиной считаем о, причем знак д(!) не определяем однозначно. При в ) 0 однородное уравнение (17) имеет решения экспоненциального типа (одно быстро убывающее, другое быстро растущее вправо).
При в< 0 уравнение (17) имеет решения колебательного характера, причем, если [г/[Т:в. ! (1000, например), на [О, Т[ укладываются десятки колебаний. В некоторых задачах (см. в 15) в(!) может иметь разные знаки на разных частях интервала [О, Т[. Что можно сказать в этих условиях об уравнении (22) для р и о возможности достаточно аккуратного численного интегрирования его на большом интервале времени? А интервал действительно большой, так как величина Т [!! з!пз р+ (1/р) совз р[ ж 2То 3!и ф соз р в принципе может быть достаточно большой, если только функция р(!) не «застревает» надолго в окрестности таких значений, где ейп р соз р = О, " — 1833 пеиелижениые методы вычислительной оизики 1ч. и 258 Пытаясь разобраться в характере уравнения (22), заметим, что ) ф~ ж а. (Мы пользуемся не очень строгими оценками; речь идет о грубом анализе, в котором величиной (1/р) созе р О(1) можно пренебречь.) Следовательно, все траектории (22) проходят в конусе с раствором, определяемым величиной ~д), т.е.
в целом траектории (22) не имеют экспоненциального роста, траекторий типа е~о1' .среди решений (23) нет, хотя на малых частях интервала [О, Т1 такие решения могут и появиться. В з 7 специально подчеркивалось, что процесс накопления погрешностей при численном интегрировании задачи Коши существенным образом зависит от устойчивости вычисляемой траектории.