Ортега Дж., Пул У. Введение в численные методы решения дифференциальных уравнений. Под ред. А.А.Абрамова (1986) (1095855), страница 59
Текст из файла (страница 59)
Если Лт сравнительно невелико, скажем Ф= 100, то мы получаем систему с Ф' =10 неизвестными и матрицей коэффициентов размера 10000 Х 10000. В то же время, каково бы ни было Ж, в каждой строке этой матрицы имеется не более пяти отличных от нуля элементов. Следовательно, при больших Я матрица становится очень слабо заполненной. Такие матрицы обычно наэыватотся большими разреженными матрицами, и они возникают не только при численном решении дифференциальных уравнений в частных производных, но и в целом ряде других задач. Именно эта разреженность позволяет современным ЭВМ решать такие большие системы уравнений с относительной легкостью.
Напомним, что, как мы видели в гл. 3, для решения линейной системы размера и Х нпо методу гауссова исключения требуется порядка и' арифметических операций. Следовательно, если бы линейная система размера 104 Х10 была "плотной", т.е. содержала мало нулевых элементов, то для ее решенця методом гауссова исключения потребовалось бы порядка 10'г операций. Даже при использовании самых быстрых современных ЭВМ, обладающих скоростью около 10 операций в секунду, для решения такой системы понадобится несколько часов машинного времени. Тем не менее в этой главе мы покажем, что учет специальной структуры и разреженности систем вида (8.2,11) позволяет решать такие системы сравнительно быстро и легко, несмотря на их большой размер.
Мы закончим этот раздел применением дискретизации, использованной нами-для уравнения Пуассона, к уравнению теплопроводности (8.1.3) с двумя пространственными переменными и, =с(и„,. +и„г). (8.2.12) Дпя простоты изложения мы по-прежнему предполагаем, что областью изменения х и у является единичный квадрат, показанный на рис. 8.3, и что на сторонах этого квадрата заданы условия Дирихле (8.2.3) .
Мы также предполагаем, что задано начальное условие и(О,х, у) =~(х, у). (8.2.13) По аналогии с методом (7.25) для случая одной пространственной переменной мы можем легко построить явный разностный метод дпя уравнения (8,2.12)-: с!1г и".',+' =и'.". + — (и".'. +и".'. +и .+и.т . — 4ит ), П 11 йг Ь/+1 1,/ — 1;+1,! 1 1,! (8.2.15) и, следовательно, мы здесь сталкиваемся с той же самой проблемой: при малости Ь приходится двигаться с очень маленькими шагами по времени. Можно попытаться обойти это ограничение на шаг по времени точно так же, как мы делали это в разеделе 7.3, т.е. используя неявные методы. Так, напрнм .п, неявный метод (7.3.2) здесь принимает вид сЬ| ит.+' =ит + — (ит.+' + и".'.+ + ит ' +ит+1.
— 4ит+' )(8.2.16) 1/ 11 й! 1,/~-1 1,/-! 1+1,!' 1-1,У' су и является безусловно устойчивым. Однако реализация этого метода требуе! решения на каждом шаге по времени системы линейных уравнений и ! + 7 1 ьг т+1 т+1 т+1 т +1 т+1 сЬг( т = — и;., 1, 7 = 1,..., Л~. (8.2.17) сЬ1 Эта система имеет ту же форму, что и система (8.2.7) дпя уравнения Пуассона, с единственным отличием, состбящим в изменении коэффициента т+1 при и, Действительно, левые части уравнений (8.2.17) в точности такие 249 и! = О, 1,..., !',у = 1,..., У. (8.2.14) ЗДЕСЬ ит — ПрнблнжЕННОЕ РЕШЕНИЕ В ТОЧКЕ (1, !' ) ПРОСтраНСтВЕННОй СЕТКИ 1/ т+1 в момент времени п1Ьг, а и1, — приближенное решение на следующем временном слое.
Члены в круглых скобках, стоящие в правой части (8.2.14), в точности соответствуют дискретизации (8.2.7) при Г11 = О. Разностный метод (8.2.14) обладает теми же свойствами, что и его одномерный аналог (7.2.5): он легко реализуем и имеет первый порядок точности по времени и второй порядок точности по пространственным переменным. Этот метод приводит к аналогичному условию устойчивости Д! < Ьзя4с), же, как у конечно-разностньгх уравнений, аппроксимирующих уравнение Гельмгольца и„„+ и„— ои = О, где о — заданная неотрицательная функция от х и у.
Входящий в это дифференциальное уравнение член — ои приводит к появлению в левой части разностных уравнений (8.2.7) члена — //2о//и//, где о// = о (хь у;) . Если взять о = 1! (сЬ/), то мы получим левую часть (82Л 7) (с противоположным знаком) . В случае одной пространственной переменной использование неявного метода типа (7.3.2) не приводит к большим вычислительным трудностям, поскольку решение трехдиагональных систем уравнений находится достаточно быстро.
Однако каждьй шаг по времени в (8.2.17) требует решения задачи типа двумерного уравнения Пуассона. Это представляет более сложную вычислительную проблему, методы решения которой будут рассмотрены в двух последующих разделах. Метод Кранка — Николсона (7.3.12) тоже легко переносится на уравнение (8.2Л2) (см. упражнение 8.2.3), но и он страдает тем же самым недостатком: на каждом шаге по времени приходится решать уравнение типа Пуассона. Мы сейчас, однако, рассмотрим другой класс методов, где основным этапом вычислений снова является решение трехдиагональных систем уравнений. Мы имеем в виду неявные методы переменных направлений. Вероятно, самым простым из этих методов является метод Писмэна— Рэчфорда: и>+1/2 >и 1 с~/ >и+1/2 >и+1/2 и =и- + — — (и +и // >У 2 //2 1+1'/ > ! / (8.2Л8а) т т+1/2 1 ~~~~ т+1 в+1 (8.2Л86) Его можно рассматривать как двухшаговый метод, где первый шаг— т+1/2 (8.2.18а) — состоит в вы/ислении промежуточных величин и,у (г',/ = т+1/2 = 1, ..., У) .
Значения и/ интерпретируются как приближенные значения решения в промежуточный момент времени (пг+ 1/2) >2>/; поскольку здесь шаг по времени равен Ь//2, в правой части (8.2Л8а) появляется множитель 1/2. Вычисления в (8.2.18а) связаны с решением для /' = 1, ..., Л/ трехдиагональных систем уравнений т+1/2 >и+1/2 т+1/2 >и т в (2 + а) и;/ — и;+1 / — и/ 1 / = (/г — 2) и,/ + и, /+, + и/ / 1= 1,...,Л', (8.2.19) где а = 2/1'/(сь/), т.е.
при каждом фиксированном / мы из трехдиагональт+1/2 ной системы (8.2.19) находим значения и,/ (/= 1,..., Ю). Матрица коэффициентов каждой из этих систем есть просто аУ + А, где А — трех- диагональная матрица (7.2.14). Таким образом, матрица о/ + А является строго диагонально доминирующей, и, следовательно, эти системы эффективно решаются методом гауссова исключения без перестановок. 250 т+1/2 После того как промежуточные величины и~; найдены, окончательт+ 1 ные значения и;; получаются из 18.2.18б) решением при 1 = 1,..., Ф трехдиагональных систем т+1 т+1 т+) (8.2.20) Матрицы коэффициентов этих систем также равны а/+ А.
Таким образом, один шаг по времени требует решения 2% трехдиагональных систем с Ф неизвестными. Можно. показать, что этот метод безусловно устойчив. Происхождение термина переменные направления объясняется тем, что мы в некотором смысле сначала с помощью (8.2 18а) аппроксимируем решение в направлении осн х, а затем, с помощью 1'8.2.18б), в направлении оси у. Имеется множество вариантов метода Писмэна — Рэчфорда, в которых применяется та же основная идея перемены направлений.
Этот класс методов является одним нз наиболее широко используемых при решении уравнений параболического типа. Эта же идея применяется и для построения итерационных методов решения эллиптических уравнений, которые будут рассматриваться в разделе 8.4. Дополнительные замечания и ссылки 8.2 Мы в этом разделе ограничились обсуждением двумерного уравнения Пуассона в квадратной области и соответствующего уравнения теплопроводносги, Однако задачи, которые возникают на практике, обычно значительно отклоняются от этих идеальных условий: область может не быть квадратом; уравнение может содержать переменные коэффициенты или вообще быть нелинейным; на одной части границы могут задавать.
ся условия Дирихле, а на другой — условия Неймана; задача может описываться не одним уравнением, а связанной системой уравнений в частных производных и может содержать три и более независимых переменных. Изложенные в этом разделе общие принципы построения разностных схем при этом не изменятся, но каждый из перечисленных факторов может вызвать дополнительные трудности. Одной из классических книг, посвященных конечно-разностным методам решения уравнений эллиптического типа, является книга [89[ (см. также [70[, где рассматриваются задачи тидродинамики и [107[).
Обсуждение и анализ методов переменных направлений, а также родственньгя им методов, таких, как метод дробных шагов, имеется в целом ряде книг (см. например, [7, 64[. В последние несколько лет все более важную роль в решении уравнений эллиптического и параболического типов играют метод конечныв элементов и друп1е проекционные методы. Хотя разработка математических основ метода конечных элементов восходит к 1940-м нэпам, превращение епэ в жизнеспособную вычислительную процедуру, особенно для задач структурнолэ анализа, было осуществлено в 1950-х и 1960-х годах, главным образом инженерами. С тех пор были получены более глубокие и широкие математические результаты и была продемонстрирована применимость метода конечнь1х элементов к эллиптическим и параболическим уравнениям общего вида. Одним из главных достоинств метода является то„что он довольно успешно справляется с криволинейными границами.
Для ознакомления с методом конечных элементов мы отсылаем читателя к книгам [47, 76 ]. УНРАЖНЕНИ,й 8.2 8.2.1. Пусть функция и нужное число раз непрерывно дифференцируема. Разложите функцию и в ряд Тейлора в окрестности точки (х, у) и покажите, что разностные аппроксимации (8.2,4) имеют второй порядок точнос1и, т.е. и„(х, у)— 251 — (1/Ь ') 1и <х + л, у) — 2и (х, у) + и (х — л, у)1 = 0(Л' ), и то же самое для аппроксимации иуу.
82.2. Выпишите в явном виде уравнения (82.7), (8.2.10) и (8.2.11) для случая Л/= 3. 8.2.2. Выпишите формулы метода Крепка — Николсона для уравнения (8.2.12). 8.3. Прямые методы для больших разреженных систем В разделе 8.2 мы пришли к системе линейных уравнений (8.2.10), которую в матричной форме можно, как обычно, записать в виде Ах = Ь. При этом матрица А состоит в основном из нулей — в каждой строке и в каждом столбце имеется не более пяти ненулевых элементов. Кроме того, эти отличные от нуля элементы расположены регулярно — все они находятся на пяти диагоналях. Если организовать вычисления так, чтобы не хранить нулевые элементы и выполнять операции только с ненулевыми элементами, то можно добиться значительной экономии памяти и времени счета.