Ортега Дж., Пул У. Введение в численные методы решения дифференциальных уравнений. Под ред. А.А.Абрамова (1986) (1095855), страница 33
Текст из файла (страница 33)
Мы поступим точно так же, как в разделе 3.2. Разобьем отрезок ~а, Ь] узловыми точками 0 =хе <х, «... х„<х„„= 1 с шагом й и в каждой внутренней точке х,,...,х„аппроксимируем вторую производную центральными разностями. Подставляя эти аппроксимации в (4.4.1), приходим к следующей системе. уравнений (соответствующей (3.2.12) ); ц,~, --2и;+и;, =й~у(х;,и;), ю'=1,...,и, (4.4,3) где и, = а и и„+, = !3 известны из граничных условий (4.4.2). Если функция е нелинейна по и, то эта система и уравнений е и неизвестными также будет нелинейной.
Решение и;,..., и,*, системы (4.4.3), если оно существует, является приближением в узловых точках х,,..., х„к соответствукнцему решению и задачи (4.4.1) — (4.4.2) . Систему (4.4.3) можно следуюшим образом переписать в матричновекторной форме. Пусть А — трехдиагональная матрица (4.4.10), все элементы главной диагонали которой равны 2, а все элементы двух прилежащих диагоналей равны — 1, и пусть Н вЂ” вектор-функция от ~, ч = (и,,... ..., и„), опрецеленная как вектор-функция Н нз (4.4.5) определяются равенствами зи, +Ьг + 10из Зиг + 4Ь' + 10иг 2 — 1 — 1 2 Н(ч) — л . (4.4.10) А= 3и„+пгниг +10и~ 2 — 1 — 1 2 Рассмотрим теперь некоторые численные методы решения системы (4.4.5).
Итерации Пикара, которые мы обсуждали в разделе 4.3, в данном случае имеют вид Ач»+' = — Н(ч»). (4.4.1 1) Время проведения одной итерации почти полностью зависит от сложности вычисления Н, потому это решение трехдиагональной линейной системы находится очень быстро. Более того, в этом случае ЬУразложение мат- рицы А достаточно получить только один раз, не повторяя его на каждой итерации. Будут ли итерации (4.4.11) сходиться, целиком зависит от свойств Н. Давайте теперь рассмотрим применение к (4.4.5) метода Ньютона.
Как мы видели в разделе 4.3, в этом случае матрица Якоби Е (ч) = А + Н '(ч). (4.4.12) Таким образом, й-я итерация метода Ньютона состоит в следующем, 1. Решить 1А + Н'(ч»)] у" = — [Ач» + Н(ч")1, 2. Положить ч»+1 = ч» +у». Так как 1-я координата Н~ вектор. функции Н Н;(ч) = Ьг~(х;, и;) (4,4.13) зависит только от и~, то дН;(ди~ =О, / чь ~'.
(4.4.14) Следовательно, матрица Якоби Н' (ч) диагональна, так что линейные.системы, которые приходится решать на итерациях Ньютона (4.4.13), оказываются трехдиагональными и требуют сравнительно небольшого числа операций, даже если для нх решения используется алгоритм с перестановкой строк. В данном случае матрица Якоби А +Н'(ч) имеет вид д 2+юг — »(х,, и,) д 2+ я~ ' — у(хг иг) ди г д 2+ Ь~ — х(х„, и„) (4.4.15) 134 Таблица 43 Метод Ньютона для разностных уравненнв (4.4.9) х а=0,1 «=0,01 Й = 0,001 — 0,0058 — 0.0118 — 0,0176 — 0,0230 — 0,0276 — 0,0304 — 0,0305 -0,0266 — 0,0171 — 0,0058 — 0,0118 — 0,0176 — 0,0230 — 0,0276 — 0,0304 -0,0305 — 0,0266 — 0,0171 — 0,0058 — 0,0116 — 0,0174 — 0,0223 — 0,0274 — 0,0302 — 0.0303 — 0,0265.
-0,0170 0,1 0.2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 135 Обратите внимание, что для реализации метода Ньютона необходимо вы- числять частные производные от 8, входящие в (4.4.15). Если функция У достаточно сложная, то для получения явных выражений для этих частных производных может оказаться полезной техника символьных вычислений, рассмотренная в гл. 1. Эти выражения были бы затем включены в програм- му вычисления матрицы Якоби (4.4.15) .
Давайте вернемся к краевой задаче (4.4.б) и соответствующим разност- ным уравнениям (4.4.9). Здесь 8 определяется формулой (4.4.7), так что де(х, и) — =3+30и~, ди и 1'-й диагональный элемент матриць1 Якоби равен просто 2+ 71~ (3+ 30и~). Так как матрица А в (4.4.10) диагонально доминирующая, добавление по- ложительных членов йа (3+ 30и,'.) только увеличивает степень преоблада- ния диагональных элементов.
И в общем случае, если д8(х,и)/ди > О, 0 ~.~' < 1, и А есть матрица из (4.4.10), то А + Н (т) — диагонально доминирующая матрица, (4.4.1б) А + Н (т ) положительно определенная матрица. (4.4.17) Как мы видели в разделе 3.3, любого из этих свойств достаточно для того, чтобы трехдиагональные системы (4.4.13) метода Ньютона можно было бы численно устойчиво решать алгоритмом гауссова исключения без какой- либо перестановки строк. Можно также доказать (но это выходит за рамки настоящей книги), что любое из условий (4.4.1б) или (4.4.17) гарантирует, что сис тема (4.4.5) имеет единственное решение.
В табл. 4.3 приведены результаты решения методом Ньютона разност- ных уравнений (4.4.9). Значения даны в узлах сетки 0,1; 0,2; ...; 0,9 при Л = 0,1; 0,01 и 0,001 (л = 9, 99 и 999) . Во всех случаях в качестве на- чального приближения для метода Ньютона брался вектор чо = О, а итера- ции прекращались, когда все координаты вектора поправкиу" из (4.4.13) становились меньше по абсолютной величине, чем 10 ~. Предположим теперь, что уравнение (4.4.1) имеет более общую форму: (а(х)и') =8(х, и), 0<х~1, (4.4.18) а граничные условия оставим в виде (4.4.2) . Если для левой части этого уравнения воспользоваться разностной аппроксимацией (3.2.29), то разностные уравнения примут вид а; 1~э в| 1 — (а„1/2 +а; 1/2)ьс ьа;+1/з ь';,, =й у(х„и;), (4.4.19) с=1,...,л„ где а;, 1/з = а (х; + й/2) и по-прежнему о 0 = а, с „, = Р.
Если мы теперь обозначим через А трехдиагональную матрицу, соответствую:цую (4.4.19) (матрица А выписана явно в (3.2.30)), то система уравнений снова запишется в виде (4.4.5), итерации Пикара в виде (4.4.11) и итерации Ньютона в виде (4.4,13). Матрица Якоби будет опять симметричной и трехдиагональной. Кроме того, если будут выполнены условия ду(х, с) >О, а(х)>О, (4.4.20)„ ди 0<х(1, (Мы несколько изменили обозначения раздела 3.1: в качестве независимой переменной вместо г используется х.) Если узлы сетки по-прежнему определены формулами (4.4.8), то а«,,2 =(/+ 1/2)'й', и, следовательно, разностные уравнения (4.4.19) принимают вид (/ — 1/2)'и; 1 — (2Р + 1/2) и; + (1 + 1/2)' пс+, = с ос(й)' /(и; + с/), (4 4 24) где / = 1,..., и.
Из граничного условия и (1) = Р имеем и„„= ф, Однако другое граничное условие, и (0) = О, теперь не определяет значения и„. Вместо этого аппроксимируем и (0) оцносторонней разностью: и'(0) = (и, — в0)/а. Отсюда, так как и (0) = О, мы получаем дополнительное условие в0 = в,, подстановка которого в первое (с = 1) уравнение (4.4.24) приводит к соотношению 9 — — (и1 — и1) = с и, й~ /(и, + с/). 4 (4.4.25) Это уравнение вместе с уравнениями (4.4.24) при / =2,...,и образует систему и уравнений для определения и „..., в„. Так каку(х, ь) =сих'/(и+с/),то ду/до = сс/х'/(и + с/)1 (4.4.26) 1за то останутся в силе свойства (4.4.16) и (4.4.17) .
Вид разностных уравнений (4.4.19), в принципе, сохраняется и цля уравнения (х'и ) =сох-/(и+с/), О~;х..-1, (4.4.21): которое возникло нз задачи о диффузии, рассмотренной в разделе 3.1. В этом случае вместо (4.4.2) задаются граничные условия и (0)=0, а(1) = Д и матрица Якоби для этой системы имеет вид 4сДйг 9+ + ~)г 4с с~(2)~)г 34+ + у)2 1 4 4сд(пзбг)г (8п +2)+ —— + ~)г — (2п — 1) Лополнигельные замечания и ссылки 4.4 Локазательство того, что при условиях (4.4.16) или (4.4,17) система (4.4.5) имеет единственное решение, можно найти, например, в книге ) 54), раздел 4.4.
Вопросы численного решения уравнения диффузии (4.4.21) проанализированы Келлером )30), и мы здесь следуем его результатам. Как отмечалось в дополнительных замечаниях к разделу ЗЛ, это уравнение имеет особенность в начале координат, так как в этой точке коэффициент при е' обращается в бесконечность. УПРАЖНЕНИЯ 4.
4 4.4.1. Выпишите раэностные уравнения (44.3) и соответствующую матрицу Якоби для; а) х (х, е) = с + с', б) е (х, с) =хе'. 4.4.2. Выпишите явные формулы итераций Ньютона (4.4.13) для разностных уравнений (4.4.3) с функциями х из упражнения 4.4.1. 4.4.3. Выполните снова упражнение 4.4.1 для разностных уравнений (4.4.19) при: а) а(х) =е; б) а(х) = 1+х'.
4.4.4. Покажите, что в предположениях (4.4.20) справедливы утверждения(4.4.16) и (4.4.17) . 4.4.5. Рассмотрите двухточечную краевую задачу с"= е + 2 — е, 0 <х < 1, с (О) = и хг = О, с (1) = 1; а) выпишите в матрично-векторной форме стандартные конечно-разностные уравнения для этой задачи прин = 0,01; б) проанализируйте подробно, как бы вы стали решать систему уравнений пункта а на ЭВМ; анализ должен включать четкое описание метода, исследование того, какие проблемы могут возникнуть при применении этого метода, сколько потребуется машинного времени и т.д.; в) подробно опишите использование метода стрельбы для решения этой краевой задачи и раэберите все трудности. которые могут возникнуть при реализации этого метода на ЭВМ. Мы видим, что матрица Якоби является симметричной.
Более того, так как по предположению константы с н д положительны, нз (4.4,26) следует, что ду/д ц > О, так что, как и прежде, матрица Якоби оказывается положительно определенной. Следовательно, мы .и в этом случае можем применять метод Ньютона, не делая при решении трехдиагональных систем никаких перестановок. Глава 5 ЕСТЬ ЛИ ЧТО-НИБУДЬ ЕЩЕ, КРОМЕ КОНЕЧНЫХ РАЗНОСТЕЙ? 5.!. Введение в проекционные методы (р (х) и') + ~у (х) и = ~(х), О < х < 1, (5.1.1) и (О) = О, и (1) = О, (5.1.2) где для простоты взят отрезок [О, 1] и нулевые граничные условия.
Будем искать приближенное решение задачи (5.1.1) — (5.1.2) в виде и(х)= 2' с;р;(х), 1=1 (5.1.3) где базисные функции ~р1 удовлетворяют граничным условиям задачи, т.е. 4у(О)=~р1(1)=О, 7= 1,...,и. (5.1.4) В этом случае приближенное решение и, задаваемое формулой (5.1.3), будет удовлетворять граничным условиям. Классическим примером набора базисных функций, удовлетворяющих условиям (5.1.4), являются функции р; (х) = а1п у ах, 1' = 1,..., и. (5.1.5) Другим примером может служить набор полнномов р; (х) = х' (1 — х), у = 1,..., и. (5.1.6) !ЗВ В трех предыдуших главах мы довольно подробно рассмотрели применение конечно-разностных методов для построения приближенных решений дифференциальных уравнений. В этой главе мы остановимся на другом подходе, реализация которого имеет несколько вариантов, называемых методом конечных элементов, методом Галеркина, методом Рэлея — Ритца н т.д.