Методы решения сеточных уравнений
- Людмила Полторацкая 4 лет назад Просмотров:
1 Методы решения сеточных уравнений 1 Прямые и итерационные методы В результате разностной аппроксимации краевых и начально-краевых задач математической физики получаются СЛАУ матрицы которых обладают следующими свойствами: 1) порядок матрицы очень высок и равен числу узлов сетки; 2) матрицы являются разреженными и имеют большое число нулевых элементов; 3) матрицы являются плохо обусловленными: отношение наибольшего собственного значения матрицы к ее наименьшему собственному значению является величиной O( h 2 ). Методы решения СЛАУ получаемых в результате разностной аппроксимации исходной задачи подразделяют на прямые и итерационные. Прямыми методами решения задач называют методы с помощью которых можно получить точное решение задачи за конечное число арифметических операций. Итерационными называют методы с помощью которых можно получить приближенное решение задачи с любой заданной точностью за конечное число арифметических операций. 2 Прямые методы 2.1 Метод монотонной прогонки и метод окаймления Основным прямым методом используемым для решения СЛАУ возникающих в разностных схемах для начально-краевых задач математической физики является метод прогонки и различные его варианты. Классический вариант метода прогонки называемый также монотонной правой прогонкой представляет собой метод Гаусса без выбора главного элемента для систем с трехдиагональной матрицей. Напомним что для системы 1
2 вида c 0 y 0 + b 0 y 1 = f 0 (i = 0) a i y i 1 c i y i + b i y i+1 = f i (i = N 1) a N y N 1 c N y N = f N (i = N) метод правой монотонной прогонки позволяет получить систему y N = β N y i = α i y i+1 + β i i = N 1 N (2.1) где α 0 = b 0 c 0 ; β 0 = f 0 c 0 ; α i = b i c i a i α i 1 i = N 1; β i = f i + a i β i 1 c i a i α i 1 i = N. Для решения системы (2.1) данным методом необходимо затратить Q = 8N + 1 действий. Метод правой монотонной прогонки корректен если выражения в знаменателе для коэффициентов α i β i не обращаются в ноль ни при каких i. Если допустить что прогоночные коэффициенты находятся точно а в y N допущена ошибка ε N то погрешность решения ε i будет удовлетворять однородному уравнению ε i = α i ε i+1 то есть не будет нарастать при выполнении условия α i 1 для всех i. В этом случае можно говорить об устойчивости метода прогонки. Теорема 2.1 Если коэффициенты системы (2.1) удовлетворяют условиям b 0 0 c 0 > 0 a N 0 c N > 0 a i > 0 b i > 0 c i a i + b i i = N 1 (2.2) c 0 b 0 c N a N (2.3) причем хотя бы в одном из неравенств (2.2) или (2.3) выполняется строгое неравенство то есть для матрицы системы (2.1) имеет место диагональное преобладание то c i a i α i 1 0 и α i 1 для всех i = N 1. Доказательство. Из условий теоремы следует что 0 α 0 = b 0 c 0 1. Пусть α i 1 тогда b i+1 α i+1 = c i+1 a i+1 α i b i+1 c i+1 a i+1 b i+1 b i+1 = 1. Следовательно α i 1 для всех i = N 1. Кроме того имеют место неравенства c i a i α i 1 c i a i α i 1 b i + a i (1 α i 1 ) b i > 0 i N 1 откуда получаем что c i a i α i 1 0 при i N 1. 2
3 Остается показать что c N a N α N 1 0. По условию хотя бы в одном из неравенств (2.2) или (2.3) выполняется строгое неравенство. Если c N > a N то c N a N α N 1 0 так как α N 1 1. Если существует 1 i 0 N 1 такое что c i0 > a i0 + b i0 то c i0 a i0 α i0 1 > b i0 откуда следует что α i0 0 так как α N 1 b 0 то неравенство α i 0. Сформулированные условия являются лишь достаточными условиями корректности метода монотонной правой прогонки. Их можно ослабить разрешив некоторым из коэффициентов a i и b i обращаться в ноль. «Частично» пользуясь методом прогонки можно упростить решение системы с разреженной матрицей фрагмент которой является трехдиагональным. На этой идее основаны различные методы окаймления. В качестве примера рассмотрим систему с «почти» трехдиагональной матрицей: a 1 y N c 1 y 1 + b 1 y 2 = f 1 a i y i 1 c i y i + b i y i+1 = f i i = N 1 a N y N 1 c N y N + b N y 1 = f N. (2.4) Такая алгебраическая система возникает при отыскании периодического решения системы трехточечных уравнений: y i+n = y i при условии что: a i+n = a i b i+n = b i c i+n = c i f i+n = f i. Запишем систему (2.4) в матричном виде A N Y N = F N где Y N = (y 1. y N ) T F N = (f 1. f N ) T а матрица A N имеет вид: c 1 b a 1 a 2 c 2 b a 3 c 3 b A N = (2.5) c N 2 b N a N 1 c N 1 b N 1 b N a N c N 3
4 Присутствие ненулевых элементов в правом верхнем и левом нижнем углах матрицы (2.5) не позволяет решать систему (2.4) обычным методом прогонки. Запишем систему (2.4) в виде: A N 1 Y N 1 + U N 1 y N = F N 1 V N 1 Y N 1 c N y N = f N (2.6) где Y N 1 = (y 1. y N 1 ) T F N 1 = (f 1. f N 1 ) T c 1 b a 2 c 2 b a 3 c 3 b A N 1 = c N 2 b N 2 U N 1 = a a N 1 c N 1 b N 1 V N 1 = [ b N a N ]. Решение первого блока уравнений системы (2.6) будем искать в виде: Y N 1 = P N 1 + y N Q N 1 где P N 1 = (p 1 p 2. p N 1 ) T и Q N 1 = (q 1 q 2. q N 1 ) T решения задач: A N 1 P N 1 = F N 1 A N 1 Q N 1 = U N 1. (2.7) Поскольку A N 1 трехдиагональная матрица то системы (2.7) решаются обычной прогонкой: α i = b i c i α i 1 a i β i = f i + a i β i 1 c i α i 1 a i γ i = a iγ i 1 c i α i 1 a i i = N; α 1 = b 1 β 1 = f 1 γ 1 = a 1 ; c 1 c 1 c 1 p N 1 = β N 1 q N 1 = α N 1 + γ N 1 ; p i = α i p i+1 + β i q i = α i q i+1 + γ i i = N 2 N Выразим теперь y N : (2.8) (2.9) V N 1 Y N 1 c N y N = V N 1 P N 1 + y N V N 1 Q N 1 c N y N = f N y N = f N + b N p 1 + a N p N 1 c N b N q 1 a N q N 1 = f N + b N p 1 + a N β N 1 c N b N q 1 a N (α N 1 + γ N 1 ) cn a N α N 1 c N a N α N 1 = 4
5 = β N + α N p 1 1 α N q 1 γ N. Теперь зная y N и коэффициенты p i q i i = 1. N 1 находим решение системы: y i = p i + y N q i i = N 1. Изложенный алгоритм требует Q = 14N 8 действий. Теорема 2.2 Пусть коэффициенты системы (2.4) удовлетворяют условиям: a i > 0 b i > 0 c i a i + b i i = N (2.10) и существует такое число 1 i 0 (N 1) что c i0 > a i0 + b i0. Тогда c i a i α i 1 0 α i 1 и α i + γ i 1 для всех i = N и выполняется условие 1 α N q 1 γ N 0. Доказательство. Если выполнены условия (2.10) то из теоремы 2.1 следует что c i a i α i 1 0 и α i 1 для всех возможных i. Так как a 1 + b 1 c 1 то α 1 + γ 1 1. Предположим что α i + γ i 1. Тогда α i+1 + γ i+1 = b i+1 + a i+1 γ i c i+1 a i+1 α i a i+1 + b i+1 a i+1 (1 γ i ) c i+1 a i+1 α i a i+1 + b i+1 a i+1 α i c i+1 a i+1 α i = 1. Следовательно условие α i + γ i 1 выполнено для всех i. Так как c i0 > a i0 + b i0 то α i0 + γ i0 Метод немонотонной прогонки Рассмотрим систему с трехдиагональной матрицей c 0 y 0 + b 0 y 1 = f 0 (i = 0) a i y i 1 c i y i + b i y i+1 = f i (i = N 1) a N y N 1 c N y N = f N (i = N) для которой не выполнены условия диагонального преобладания. При формальном использовании метода Гаусса без выбора главного элемента на l-м шаге приведения исходной 5
6 системы к системе с верхнетреугольной матрицей получаем: (c l a l α l 1 )y l + b l y l+1 = (f l + a l β l 1 ) (i = l) a i y i 1 c i y i + b i y i+1 = f i (l + 1 i N 1) a N y N 1 c N y N = f N (i = N). (2.11) Если матрица системы имеет диагональное преобладание то c l a l α l 1 0 то есть первое уравнение системы (2.11) можно представить в виде y l = α l y l+1 + β l α l = b l c l a l α l 1 β l = f l + a l β l 1 c l a l α l 1 причем α l 1. Если же диагональное преобладание не имеет места то нельзя гарантировать что c l a l α l 1 0 и даже если это так условие α l 1 может быть не выполнено. В этом случае метод монотонной прогонки уже не применим. Предположим что исходная система имеет единственное решение. Тогда применяя метод Гаусса с выбором главного элемента в строке на l-м шаге получаем: C l y ml + b l y l+1 = F l (i = l) A l y ml c l+1 y l+1 + b l+1 y l+2 = Φ l (i = l + 1) a l+2 y l+1 c l+2 y l+2 + b l+2 y l+3 = f l+2 (i = l + 2) a i y i 1 c i y i + b i y i+1 = f i (l + 3 i N 1) a N y N 1 c N y N = f N (i = N). (2.12) Если l = 0 то имеют место равенства C 0 = c 0 A 0 = a 1 F 0 = f 0 Φ 0 = f 1 и m 0 = 0. При дальнейшем преобразовании системы (2.12) возможно два варианта. В первом случае C l b l. Тогда первое уравнение системы (2.12) можно переписать в виде y ml α l y l+1 = β l α l = b l C l β l = F l C l где α l 1 или же в виде y θl+1 = α l y ml+1 + β l где θ l+1 = m l m l+1 = l + 1. Исключая из уравнения соответствующего i = l + 1 слагаемое A l y ml получаем C l+1 y ml+1 + b l+1 y l+2 = F l+1 (i = l + 1) A l+1 y ml+1 c l+2 y l+2 + b l+2 y l+3 = Φ l+1 (i = l + 2) (2.13) a i y i 1 c i y i + b i y i+1 = f i (l + 3 i N 1) a N y N 1 c N y N = f N (i = N) где m l+1 = l + 1 C l+1 = c l+1 A l α l F l+1 = Φ l + A l β l A l+1 = a l+2 Φ l+1 = f l+2. Во втором случае C l 7 где опять же α l 1 или что то же самое в виде y θl+1 = α l y ml+1 + β l где θ l+1 = l + 1 m l+1 = m l. Полученное уравнение используем для исключения y l+1 из уравнения соответствующего i = l + 1. В результате снова придем к системе (2.13) но теперь m l+1 = m l C l+1 = c l+1 α l A l F l+1 = Φ l c l+1 β l A l+1 = a l+2 α l Φ l+1 = f l+2 + a l+2 β l. Итак алгоритм немонотонной прогонки следующий: 1) C = c 0 A = a 1 F = f 0 Φ = f 1 и m 0 = 0; 2) для всех i = N 1 выполняются действия (прямой ход) а) если C b i то α i = b i C β i = F C C = c i+1 Aα i F = Φ+Aβ i θ i+1 = m i m i+1 = i+1 для i N 2 находим A = a i+2 Φ = f i+2 ; б) если C 8 Для решения системы (2.14) рассмотрим метод матричной прогонки. По аналогии с системой скалярных трехточечных уравнений будем искать решение в виде Y i = α i Y i+1 + β i i = N 1 N где α i матрицы размера M i M i+1 β i векторы размерности M i. При этом формально получаем: α 0 = C 1 0 B 0 α i = (C i A i α i 1 ) 1 B i i = N 1; β 0 = C 1 0 F 0 β i = (C i A i α i 1 ) 1 (F i + A i β i 1 ) i = N; Y N = β N. Если матрицы C 0 и (C i A i α i 1 ) для i = N не вырождены то алгоритм матричной прогонки корректен. Будем говорить что алгоритм устойчив если α i 1 для i = N. Теорема 2.3 Если C i для i = N невырожденные матрицы а A i и B i ненулевые матрицы для i = N 1 и выполнены условия C 1 0 B 0 1 C 1 N A N 1 C 1 i i A i + C 1 B i 1 для i = N 1 причем хотя бы в одном из неравенств имеет место строгое неравенство то алгоритм метода матричной прогонки устойчив и корректен. Если все матрицы A i B i C i квадратные и имеют размер M M а все векторы Y i и F i имеют размерность M то для решения системы (2.14) методом матричной прогонки потребуется Q = O(M 3 N + M 2 N) действий. 3 Итерационные методы 3.1 Итерационные методы решения разностных уравнений как задачи на установление Пусть требуется решить уравнение Au = f (3.1) где A линейный самосопряженный положительно определенный оператор действующий в вещественном гильбертовом пространстве H со скалярным произведением (y v) и нормой y = (y y) f H произвольная функция. 8
9 Уравнению (3.1) можно поставить в соответствие абстрактную задачу Коши: dv + Av = f t > 0 dt v(0) = v 0 (3.2) где v 0 произвольный элемент пространства H v(t) функция со значениями в H. Покажем что при сформулированных условиях на оператор A решение задачи (3.2) стремится по норме к решению задачи (3.1): lim v(t) y = 0. t Для этого введем функцию z(t) = v(t) y. Она будет удовлетворять задаче Коши с однородным уравнением: dz + Az = 0 t > 0 dt z(0) = v 0 y. Умножая уравнение в задаче (3.3) скалярно на z(t) получаем: ( ) dz dt z + (Az z) = 0. Так как ( ) dz dt z = 1 d 2 dt (z z) = 1 d 2 dt z(t) 2 и по условию существует такое δ > 0 что (Az z) δ z 2 то приходим к неравенству из которого следует что d dt z(t) 2 + 2δ z 2 0 e 2δt d dt z(t) 2 + 2δe 2δt z 2 = d dt Интегрируя последнее неравенство получаем: откуда находим e 2δt z(t) 2 z(0) 2 ( e 2δt z 2) 0. v(t) y = z(t) e δt v 0 y 0 при t. (3.3) Следовательно для того чтобы найти приближенное решение задачи (3.1) можно построить разностную схему для задачи (3.2) и вычислять ее решение до тех пор пока не будет выполнено условие dv dt 10 Начальное условие в задаче (3.2) выбирается произвольно. Если Λ линейный самосопряженный положительно определенный разностный оператор аппроксимирующий оператор A то разностную схему для задачи (3.2) можно записать в виде B k y k+1 y k τ k+1 + Λy k = f k = y 0 = v 0 где τ k+1 шаги в общем случае неравномерной сетки по времени B k обратимый оператор. Разностную схему (3.4) можно интерпретировать как итерационный процесс: y 0 = v 0 B k y k+1 = (B k τ k+1 Λ)y k + τ k+1 f k = для нахождения приближенного решения сеточного уравнения (3.4) Λy = f. (3.5) В качестве практического критерия прекращения итераций целесообразно выбирать условие малости невязки: где ε требуемая точность. Λy f 11 с минимальным числом действий Q k для нахождения k-й итерации. Тогда минимальным будет общее число арифметических операций Q(ε) которое нужно выполнить чтобы получить при помощи метода (3.4) решение уравнения (3.5) с заданной точностью ε > 0 при любом выборе начального приближения: n(ε) Q(ε) = Q k. 3.2 Итерационный метод переменных направлений k=1 В качестве первого примера рассмотрим разностную задачу Дирихле для уравнения Пуассона в прямоугольнике Ḡ =
12 Погрешность z +1 = y +1 y удовлетворяет задаче: z z = Λ 1 z Λ2 z x ω h τ (1) +1 z = 0 x γh z +1 z τ (2) +1 z +1 = 0 x γ h z 0 = y 0 y. = Λ 1 z Λ2 z +1 x ω h Если ввести пространство H h сеточных функций со скалярным произведением (y v) = x ω h y(x)v(x)h 1 h 2 все элементы которого ограничены по норме y = (y y) и обращаются в нуль на γ h то операторы A p y = Λ p y p = 1 2 будут самосопряженными положительно определенными и перестановочными причем где λ min p = 4 h 2 p λ min p E A p λ max p E p = 1 2 (3.8) sin 2 πh p 2l p λ max p = 4 h 2 p cos 2 πh p 2l p p = 1 2. Задачу для погрешности z +1 можно переписать в виде: ( ) ( ) E + τ (1) +1 A 1 z = E τ (1) +1 A 2 z ( ) ( ) E + τ (2) +1 A 2 z +1 = E τ (2) +1 A 1 z где z 0 = y 0 y = или исключая z и пользуясь перестановочностью операторов A 1 и A 2 в виде z +1 = S +1 z где S +1 = S (1) +1 S(2) +1 оператор перехода со слоя на слой а операторы S (1) +1 и S(2) +1 имеют вид: ( ) 1 ( ) ( ) 1 ( ) S (1) +1 = E + τ (1) +1 A 1 E τ (2) +1 A 1 S (2) +1 = E + τ (2) +1 A 2 E τ (1) +1 A 2. Таким образом где T n = z n = T n z 0 (3.9) n S разрешающий оператор причем Tn = T n. =1 12
13 Итерационные параметры τ (1) +1 и τ (2) +1 подбираются таким образом чтобы для получения точности ε затратить минимальное число шагов. Для этого необходимо точно знать границы спектра оператора A. Выбор оптимальных (по Жордану) параметров обеспечивает минимум T n. 3.3 Выбор оптимальных параметров итерационной схемы переменных направлений (по Жордану) Из равенства (3.9) следует оценка z n T n z 0. Величина T n зависит от параметров τ (1) и τ (2) = Задача состоит в отыскании таких параметров τ (1) 1 τ (1) (1) 2. τ n и τ (2) 1 τ (2) (2) 2. τ n где число итераций n = n(ε) необходимое для достижения точности ε задано: min <τ (1) τ (2) >T n = q n. Из оценок (3.8) следует что спектр оператора A p принадлежит отрезку: δ p λ(a p ) p p = 1 2. Заменим операторы A 1 и A 2 операторами A 1 и A 2 спектры которых совпадают: ηe A p E p = 1 2 η > 0. Положим A 1 = (qe ra 1 ) 1 (A 1 pe) A 2 = (qe + ra 2 ) 1 (A 2 + pe) где r p q числа подлежащие выбору и введем параметры: ω (1) = τ (1) r q τ (1) p ω(2) = τ (2) + r q + τ (2) p. При этом получим: S = (1) S S (2) где S (1) = (E + ω (1) A 1) 1 (E ω (2) A 1) S(2) = (E + ω (2) A 2) 1 (E ω (1) A 2). Выразим T n через собственные значения операторов A 1 и A 2 : α k1 = λ (1) k 1 (A 1) β k2 = λ (2) k 2 (A 2) k α = N α α = 1 2. Так как A 1 A 2 = A 2 A 1 то они имеют общую систему собственных функций то же что и операторы A 1 A 2 A и T n. Пусть λ k (T n ) собственные значения оператора T n. Поскольку λ( S (1) ) = 1 ω(2) α 1 + ω (1) α λ( S (2) ) = 1 ω(1) β 1 + ω (2) β 13
14 то причем λ(t n ) = n =1 1 ω (2) α 1 + ω (1) α 1 ω (1) β 1 + ω (2) β (3.10) 0 t p > 0. p = κ t κ + t r = ( )p q = r + 1 p 1 (3.13) Пусть задана точность ε > 0 итерационного процесса и известны границы δ α α операторов A α. По формулам (3.13) находим η p q и r. После этого можно определить число итераций n(ε) обеспечивающих заданную точность ε > 0: y n u ε y 0 u. Справедлива приближенная формула: n(ε) 1 π 2 ln 4 ε ln 4 η. (3.14) Введем обозначения: θ = η2 16 ) (1 + η2 2 σ = 2 1 = n. 2n 14
15 Для определения ω имеет место формула: ω = Остается определить τ (1) и τ (2) исходя из формул (1 + 2θ)(1 + θ σ ) 2θ σ/2 (1 + θ 1 σ + θ 1+σ = n (3.15) ) ω = τ (1) r q τ (1) p В частном случае когда δ 1 = δ 2 = δ и 1 = 2 = имеем: p = r = 0 q = 1 η = δ κ = 1 η 1 + η. Преобразование операторов при этом принимает вид: Условие ω (1) = ω (2) дает τ (1) = τ (2) = τ. ω = τ (2) + r q + τ (2) p. (3.16) A 1 = A 1 A 2 = A 2 ω (1) = τ (1) ω (2) = τ (2). 3.4 Итерационный метод на основе эволюционно-факторизованной схемы Рассмотрим краевую задачу для эллиптического уравнения без смешанных производных в прямоугольном параллелепипеде G (прямоугольнике в двумерном случае) либо области G составленной из параллелепипедов: 3 Lu = f(x) x G L = x p=1 p u G = µ(x). ( k p (x) ) k p (x) > 0 x p (3.17) Поставим в области G задачу на установление v = Lv + f(x) x G t > 0 t v G = µ(x) v t=0 = v 0 (x) (3.18) где v 0 (x) произвольная функция желательно удовлетворяющая граничным условиям задачи. Введем в расчетной области пространственную сетку ω h = ω h +γ h где ω h множество внутренних узлов а γ h множество граничных узлов а также сетку по времени с шагом τ и построим для задачи (3.18) эволюционно-факторизованную схему: (E 0.5τΛ 1 ) (E 0.5τΛ 2 ) (E 0.5τΛ 3 )y > <<>t = Λy + f > << y 2 >y 1 15
16 где Λ p разностная аппроксимация L p p = Λ = 3 Λ p. Переход между целыми слоями по времени осуществляется в три шага: (E 0.5τΛ 1 )y 1 = Λy + f y 1 γ 1 h = (E 0.5τΛ 2 )y 2 γ 1 h = 0 (E 0.5τΛ 2 )y 2 = y 1 p=1 y 2 γ 2 h = (E 0.5τΛ 3 )µ t γ 2 h = 0 (E 0.5τΛ 3 )y t = y 2 y γh = µ где γ 1 h граничные узлы по направлению x1 γ 2 h граничные узлы по направлению x2. Пусть λ (i) p собственные значения соответствующих задач Штурма-Лиувилля: Λ p y + λ p y = 0 x ω h p = y = 0 x γ p h. Тогда множители роста ρ для рассматриваемой эволюционно-факторизованной схемы удовлетворяют уравнению: ( τλ 1 )( τλ 2 )( τλ 3 )(ρ 1) = τ(λ 1 + λ 2 + λ 3 ). Одномерный случай: Двумерный случай: Трехмерный случай: ρ = 1 0.5τλ τλ 1. ρ = 1 0.5τλ τλ τλ τλ 2. ρ = 1 τ(λ 1 + λ 2 + λ 3 ) ( τλ 1 )( τλ 2 )( τλ 3 ). Счет на установление по эволюционно-факторизованной схеме сходится значительно быстрее если вместо постоянного шага по времени τ (итерационного параметра) использовать специальный набор τ s s = S. В настоящий момент наиболее эффективным является так называемый логарифмический набор τ s. Идея его построения состоит в том чтобы на каждом шаге по времени гасить одну какую-либо гармонику в численном решении или группу соседних гармоник (то есть добиваться минимума по модулю соответствующих множителей роста). 16
17 В одномерном случае целесообразно выбирать τ min = 2 λ max τ max = 2 λ min так как при этом полностью гасятся первая и последняя гармоники а остальные подавляются частично. Из тех же соображений в двумерном случае выбирают τ min = 2 max <λ max 1 λ max 2 >τ max = 2 min <λ min 1 λ min 2 >. В трехмерном случае уравнение ρ(τ) = 0 либо не имеет положительных корней либо имеет три вещественных корня два из которых положительны а один отрицателен. Если положительных корней у уравнения ρ(τ) = 0 нет то τ min соответствует минимуму ρ(τ) если положить λ p = λ max p : τ min = τ ( 1 b τ = 3 cos ϕ + 2π 3 ) b = λ 1 λ 2 + λ 1 λ 3 + λ 2 λ 3 c = λ 1 λ 2 λ 3 ϕ = 1 3 arccos ( c ( ) ) 3/2 b 3 если ρ(τ ) 0. Если же ρ(τ ) 18 Чебышевский набор: fч = πs S. Наилучшими свойствами в плане скорости сходимости обладает линейная комбинация равномерного и чебышевского наборов: f(s) = C fр(s) + (1 C) fч(s) (3.19) где рекомендуемое значение C = π π + 2. На практике точные границы спектра оператора Λ как правило неизвестны. Для них можно получить только более или менее точные оценки. Поэтому в расчетах вместо границ спектра рекомендуется брать величины λ min и λ max q где значения λ min и λ max получены q оценочно параметр q > 1 то есть производить расширение границ спектра. Практическая рекомендация: q 3. При этом будет иметь место некоторая потеря скорости расчетов но для набора (3.19) она не очень велика. 3.5 Некоторые итерационные методы не сводящиеся к задаче на установление Существует много актуальных задач к которым не применим счет на установление с помощью эволюционно-факторизованной схемы. Например это задачи в криволинейных областях для которых приходится использовать треугольную сетку задачи со смешанными производными и несамосопряженными операторами задачи для уравнений высоких порядков и т.д. При их разностной аппроксимации возникает необходимость решать СЛАУ вида: Au = f где A сильно разреженная матрица размера M M u f векторы размера M 1. Для задач рассматриваемого класса наиболее быстро сходятся многошаговые итерационные методы сопряженных направлений. Они сводятся к минимизации некоторой квадратичной формы Φ(u) и построению соответствующего ортонормированного базиса. Теоретически они позволяют найти точное решение системы но на практике в силу того что M достаточно велико расчеты прекращают по достижении определенной точности. Критерием прекращения итераций выступает малость невязки Метод сопряженных градиентов 18
19 Метод применим для эрмитовых знакоопределенных матриц A. При этом осуществляется поиск экстремума квадратичной формы Φ(u) = (u Au 2f). В качестве начального приближения выбирается произвольное u 1. Вспомогательные величины: r s невязка p s направление очередного спуска q s вспомогательный вектор. Итерационный алгоритм имеет вид: Au s f s = 1 r s = q s 1 r s 1 s = (q s 1 p s 1 ) p s = p s 1 + r s (r s r s ) p 0 = 0 q s = Ap s u s+1 = u s p s (q s p s ). Условие прекращения итераций: r s 20 4 Контрольные вопросы 1) Сформулируйте и докажите теорему о достаточных условиях корректности и устойчивости метода монотонной правой прогонки. 2) В чем состоит идея метода окаймления? Проиллюстрируйте этот метод на примере системы с «почти трехдиагональной»матрицей для периодической сеточной функции. 3) Сформулируйте алгоритм немонотонной прогонки. 4) Сформулируйте алгоритм матричной прогонки и достаточные условия его применимости. 5) В чем состоит идея итерационных методов решения систем сеточных уравнений? Как эти методы связаны с задачами на установление? Запишите каноническую форму двухслойной итерационной схемы. 6) Изложите итерационный метод переменных направлений для приближенного решения задачи для уравнения Пуассона в прямоугольнике. 7) Составьте эволюционно-факторизованную итерационную схему для уравнения Пуассона в прямоугольнике и прямоугольном параллелепипеде. Как выбираются оптимальные итерационные параметры? 20