Решение систем линейных алгебраических уравнений

Глава 1 Решение систем линейных алгебраических уравнений В данной главе мы будем рассматривать следующую задачу — требуется решить систему линейных алгебраических уравнений (в дальнейшем СЛАУ), задаваемых в виде матричного уравнения: Ax = F, где A = (aij ) — матрица размерности n×n, причем det A 6= 0; x = (x1 , x2 , . . . , xn )T — вектор-столбец неизвестных, F = (f1 , f2 , . . . , fn )T — заранее заданный вектор-столбец правой части. Далее мы будем предполагать, что у всех рассматриваемых задач решение существует и единственно (в данном случае это следует из того, что определитель матрицы A не равен нулю). Рассмотрим несколько методов решения СЛАУ. Методы решения делятся на прямые, где мы получаем x за конечное число действий, и итерационные — получаем x как предел некоторой сходящейся последовательности {xk } : x = lim xk . k→∞ В случае прямых методов получается точное (до погрешности аппаратуры) решение. В итерационных методах вводится точность решения ε > 0. Если |xk −x| < ε (x — точное решение) — заканчиваем процесс вычислений, иначе — вычисляем очередное xk+1 и т. д. К прямым методам относятся формулы Крамера, метод Гаусса. К итерационным методам — метод Зейделя, метод верхней релаксации и т.д. Основным показателем при оценке эффективности метода является его сложность. В прямых методах это чаще всего количество арифметических операций, необходимых для вычисления x, а в итерационных — количество итераций, необходимых для достижения заданной точности ε (его записывают как k0 (ε)). К примеру, метод Крамера (модифицированный) имеет сложность O(n4 ), а метод Гаусса — O(n3 ). Выбор метода, в основном, зависит от размерности матрицы. При решении систем большой размерности чаще используют итерационные методы. 1.1 Прямые методы решения СЛАУ. Метод квадратного корня В этой главе мы рассмотрим несколько прямых методов для решения системы Ax = F. (1.1) В самом общем случае используется метод Гаусса1 , в котором мы будем выделять прямой ход — приведение матрицы коэффициентов к верхнетреугольному виду (мы думаем, что все читатели знают, 1 Впервые описан К. Гауссом в 1849г., правда, без обратного хода. Этот же метод часто называют методом ГауссаОстроградского. 4 1.1. ПРЯМЫЕ МЕТОДЫ РЕШЕНИЯ СЛАУ. МЕТОД КВАДРАТНОГО КОРНЯ 5 как это делается), и обратный ход — приведение верхнетреугольной матрицы к диагональной и нахождение самого решения. Однако существует целый класс методов, использующих специфику конкретной матрицы — например, метод прогонки. Метод прогонки имеет сложность O(n), но он применим только к трехдиагональным матрицам. Вычисление элементов обратной матрицы Чтобы решить систему (1.1), можно вычислить A−1 , тогда решение будет представимо в виде x = A−1 F. Для определения элементов матрицы A−1 рассмотрим уравнение AA−1 = E, где матрица A — задана, а требуется найти A−1 . Обозначим за aij элементы матрицы A, а за zij — элементы матрицы A−1 . Умножая строки матрицы A на столбцы A−1 (по определению произведения матриц), получим n2 алгебраических уравнений: n X ail zlj = δij , i, j = 1, n. l=1 Эти уравнения можно объединить в группы с фиксированным индексом j, то есть у нас получится по n уравнений для определения каждого столбца матрицы A−1 . Это будут уравнения вида Az j = ej , где z j — j-й столбец матрицы A−1 , а ej — j-й столбец матрицы E. Получаем n систем по n уравнений в каждой. Далее используем метод Гаусса решения систем уравнений, при этом сложность данного метода будет равна O(n4 ). Кроме того, можно применить прямой ход метода Гаусса сразу ко всем системам (т. к. матрица 3 A — общая для всех систем), и сложность метода станет ≈ n3 . Наложим ограничения на матрицу A из уравнения (1.1). Определение. Матрица A называется положительно определенной, если выполнено одно из условий (они эквивалентны друг другу): 1. все главные миноры больше нуля; 2. все собственные числа матрицы A (далее обозначаются как λ(A)) положительны; 3. ∀ x 6= 0 hAx, xi > 0. Примечание. Из пункта 1 (а также из 3) определения следует необходимое условие положительной определенности — все диагональные элементы положительны. Достаточным условием положительной определенности является условие диагонального преn X обладания: если элементы aii > |aij |, то матрица A = (aij ) положительно определена. j=1 j6=i Метод квадратного корня Будем рассматривать симметричные (т. е. A = AT ) положительно определенные матрицы. Для таких матриц справедливо представление (разложение Холецкого): A = LLT , (1.2) 6 Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ где L — нижнетреугольная матрица:  l11   L=   .. . ... lij ..  ..  .  .  0  . lnn Рассмотрим задачу (1.1). Подставив разложение (1.2) для A в уравнение (1.1), получим: LLT x = F. Обозначим LT x = y, тогда получим систему уравнений относительно новой переменной y : Ly = F. Учитывая то, что матрица L нижнетреугольная, для решения этой системы требуется только обратный ход метода Гаусса (решение уравнений построчно сверху вниз). Аналогично, так как y = LT x, то нужно решать только СЛАУ с верхней треугольной матрицей. Поэтому основной вклад в сложность 3 даст вычисление элементов матрицы L — и общая сложность метода будет приближенно равна n6 . Теперь покажем, что матрицу L можно найти. Обозначим L = (lij ), LT = (lij ) — из определения транспонированной матрицы вытекает, что lij = lji . Найдем элементы этих матриц, для чего распишем разложение (1.2):    l11 0 . . . 0 l11   ..   .. ..    . lij . .   = (aij ).  0 LLT =      . . . . . .   . lij . 0   . lnn 0 . . . 0 lnn По правилу умножения двух матриц получаем: aij = n X min(i, j) lil llj = {учитывая треугольный вид матриц} = l=1 X lil llj . l=1 Найдем элементы lij матрицы L. Сначала вычислим l11 : a11 = l11 l11 =⇒ 2 a11 = l11 =⇒ l11 = √ a11 (a11 > 0). Получив l11 , найдем элементы первого столбца матрицы L: ai1 = li1 l11 = li1 l11 =⇒ li1 = ai1 , l11 i = 2, n. Рассмотрим второй столбец матрицы A : ai2 = li1 l12 + li2 l22 =⇒ li1 l21 + li2 l22 = ai2 , i = 2, n. (1.3) 2 2 2 2 — а l21 уже найдено. = a22 − l21 = a22 =⇒ l22 + l22 Пусть i = 2, тогда l21 Извлекая корень, получим: s s s q 2 2 a a a11 a22 − a221 2 = l22 = a22 − l21 a22 − 221 =⇒ l22 = a22 − 21 = . l11 a11 a11 Корректность вытекает из того, что a11 > 0, тельной определенности матрицы A. a11 a22 − a212 = a11 a21 a12 > 0 по условию положиa22 1.1. ПРЯМЫЕ МЕТОДЫ РЕШЕНИЯ СЛАУ. МЕТОД КВАДРАТНОГО КОРНЯ 7 Подставив в (1.3) только что вычисленное значение l22 , получим представление для остальных элементов второго столбца: ai2 − li1 l21 li2 = , i = 3, n. l22 Рассмотрим произвольный ( k-й ) столбец матрицы A : li1 l1k + li2 l2k + . . . + li(k−1) l(k−1)k + lik lkk = aik ⇐⇒ ⇐⇒ li1 lk1 + li2 lk2 + . . . + li(k−1) lk(k−1) + lik lkk = aik , i = k, n. (1.4) Пусть i = k, тогда из (1.4) получаем: 2 2 2 2 + lkk = akk . lk1 + lk2 + . . . + lk(k−1) 2 Переносим все слагаемые, кроме lkk , в правую часть и извлекаем корень: lkk = q 2 − . . . − l2 akk − lk1 k(k−1) . (1.5) Можно показать, что под корнем получается отношение положительного минора k-го порядка, к положительной сумме: s |...| > 0 . lkk = (. . .) > 0 Подставив (1.5) в (1.4), получим: lik = aik − li1 lk1 − . . . − li(k−1) lk(k−1) , lkk i = k + 1, n. Мы полностью определили элементы матрицы L, доказав тем самым работоспособность нашего метода. Заметим, что мы использовали симметричность матрицы A, когда брали в уравнениях для нахождения lij только те ее элементы, которые лежат на главной диагонали или ниже ее. Модифицированный метод квадратного корня Теперь несколько модифицируем метод квадратного корня, а точнее, распространим его на более широкий класс задач. Итак, мы снова рассматриваем уравнение Ax = f при условии, что det A 6= 0. (1.6) Предположим, что матрица A симметрична (A = AT ) — положительной определенности требовать не будем. Покажем, что для таких матриц справедливо представление A = LDLT , (1.7) где L — нижнетреугольная матрица с единицами на главной диагонали, а D — диагональная матрица:    L=   1 .. . lij  ... 0 ..  . ;  .. . 0 1    D=   d11 .. . .. . ... .. ... .  ..  .  .  0  dnn 8 Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ Для наглядности (1.7) можно представить так (LT = (lij )):  1      .. . lij  ... 0 d11 ..    .  0  . .. . 0  .. 1 .. . ... .. ... .  1 ..   ..  . .   0  . 0   .. dnn 0 ...  lij    = A.   1 .. . Обозначим произведение матриц DLT за новую матрицу B = (bij ) и заметим, что lij = lji . Теперь найдем элементы матриц L и D. Согласно определению произведения матриц, элемент aik при i > k определяется так: aik = li1 b1k + li2 b2k + . . . + li(k−1) b(k−1)k + lik bkk — мы отбросили нулевые последние слагаемые, обусловленные видом матрицы L. А элемент aik при i = k определяется так: akk = lk1 b1k + lk2 b2k + . . . + lk(k−1) b(k−1)k + lkk bkk . Теперь воспользуемся тем, что bjk = djj ljk = djj lkj , и перепишем два последних уравнения: aik = li1 (d11 lk1 ) + li2 (d22 lk2 ) + . . . + li(k−1) (d(k−1)(k−1) lk(k−1) ) + lik bkk ; (1.8) akk = lk1 (d11 lk1 ) + lk2 (d22 lk2 ) + . . . + lk(k−1) (d(k−1)(k−1) lk(k−1) ) + lkk bkk . (1.9) При k = 1 (1.9) выглядит так: 2 d11 = a11 . l11 Зная, что l11 = 1, получим формулу для d11 : d11 = a11 . Теперь рассмотрим случай k = 1, i = 2, n. Из (1.8) получим: li1 (d11 l11 ) = ai1 =⇒ {l11 = 1} =⇒ li1 = ai1 . d11 Мы получили уравнения, определяющие первый столбец матрицы L. Аналогично, при k = 2 (1.9) и (1.8) перепишутся так: l21 (d11 l21 ) + l22 (d22 l22 ) = a22 li1 (d11 l21 ) + li2 (d22 l22 ) = ai2 =⇒ =⇒ {l22 = 1} li2 = =⇒ 2 d11 ; d22 = a22 − l21 ai2 − li1 d11 l21 . d22 И так далее. Вот как будут выглядеть общие формулы для элементов матриц D и L : 2 2 dkk = akk − (lk1 d11 + . . . + lk(k−1) d(k−1)(k−1) ); lik = (aik − li1 d11 lk1 − . . . − li(k−1) d(k−1)(k−1) lk(k−1) ) . dkk Заметим, что формулы корректны — элементы dkk не могут быть равны нулю, так как в этом случае определитель матрицы A обращался бы в ноль. 1.2. ЛИНЕЙНЫЕ ОДНОШАГОВЫЕ ИТЕРАЦИОННЫЕ МЕТОДЫ 9 Итак, разложение (1.7) получено. Подставим его в исходное уравнение (1.6): LDLT x = f =⇒ {обозначим DLT x = y} =⇒ Ly = f. Мы получили уравнение относительно вектора y. В левой части стоит нижнетреугольная матрица L, что позволяет нам сильно сократить вычисления, применяя только обратный ход метода Гаусса. После нахождения y вспомним, как мы его определяли: DLT x = y =⇒ {обозначим LT x = z} =⇒ Dz = y — это система уравнений, определяемая диагональной матрицей D. Решая ее, получим z. Теперь у нас есть простая система для нахождения x : LT x = z. Наличие в левой части верхнетреугольной матрицы L, опять же, сильно упрощает вычисления. Сложив количество операций, необходимых для получения разложения (1.7) и для поиска x, по3 лучим сложность данного метода, равную ≈ n6 . Итерационные методы решения СЛАУ Это следующая группа методов поиска решений систем линейных уравнений, особенно эффективная при решении больших систем (с числом неизвестных порядка тысячи и более). В общем случае сначала задаются некоторым вектором x0 , называемым начальным приближением. От него строится последовательность x1 , x2 , . . . , xk и так далее, где число k называют номером итерации. В общем случае (k + 1)-е приближение зависит от всех предыдущих: xk+1 = Fk+1 (x0 , x1 , . . . , xk ). От последовательности, естественно, ожидается сходимость к вектору x, который будет являться решением исходной системы. Определение. Итерационный метод называется m−шаговым, если каждое последующее итерационное приближение строится лишь по m предыдущим: xk+1 = Fk+1 (xk−m+1 , . . . , xk−1 , xk ) — на практике наиболее часто используется m = 1 и m = 2. Определение. Если Fk+1 — линейная функция, то соответствующий итерационный метод называется линейным. 1.2 Линейные одношаговые итерационные методы Согласно определению, общее выражение для xk+1 в линейных одношаговых итерационных методах таково: xk+1 = S k+1 xk + ψk+1 , (1.10) где S k+1 — некоторая матрица, а ψk+1 — некоторый вектор. Логично потребовать, чтобы вектор x = A−1 f (то есть искомое решение) при подстановке вместо xk+1 и xk обращал (1.10) в тождество: A−1 f = S k+1 A−1 f + ψk+1 , 10 Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ верное для всех k. Преобразовав, получим: (A−1 − S k+1 A−1 )f = ψk+1 =⇒ ψk+1 = Qk+1 f, где Qk+1 = A−1 − S k+1 A−1 . Потребуем, чтобы эта матрица была обратима, и формулу для Qk+1 домножим на A : S k+1 = E − Qk+1 A. Это выражение подставим в (1.10): xk+1 = xk − Qk+1 Axk + Qk+1 f. Домножим это выражение на (Qk+1 )−1 и перенесем слагаемые: Qk+1
−1 (xk+1 − xk ) + Axk = f. Теперь домножим и разделим первое слагаемое на некоторую константу τk+1 : Qk+1
−1 τk+1 xk+1 − xk + Axk = f. τk+1
−1 Обозначив Bk+1 = Qk+1 τk+1 , получим так называемую каноническую форму записи одношагового итерационного метода: xk+1 − xk + Axk = f. (1.11) Bk+1 τk+1 Примечание. Мы могли бы и не вводить константу τk+1 , однако впоследствии будет видно, что с матрицей Bk+1 и константой τk+1 удобнее работать по отдельности. Примечание. Далее мы иногда будем отождествлять итерационный метод и его каноническую форму — так что не пугайтесь. Примечание. Матрица B обязана быть обратимой — иначе вывод канонической формы записи итерационного метода будет некорректен. Определение. Если Bk+1 = E, то соответствующий метод называется явным (находим xk+1 , не решая системы уравнений), в противном случае — неявным. Определение. Если Bk+1 = B и τk+1 = τ, то метод называется стационарным, в противном случае — нестационарным. Определение. Вектор z k = xk − x (отклонение от точного решения) называется погрешностью на k-й итерации. k→∞ Определение. Итерационный метод называется сходящимся, если kz k k −→ 0 для некоторой выбранной нормы || · ||. Понятно, что предел может не быть достигнут за конечное число шагов. В таком случае мы задаемся некоторой точностью ε — достаточно малым числом ( ∼ 10−3 , 10−5 ) — и будем останавливать итерационный процесс при выполнении неравенства ||xk − x|| 6 ||x0 − x|| 1 ε , или ||xk − x|| 6 ε||x0 − x|| для некоторого k. В этом случае мы будем говорить, что решение получено с точностью ε. 1.2. ЛИНЕЙНЫЕ ОДНОШАГОВЫЕ ИТЕРАЦИОННЫЕ МЕТОДЫ 11 Кроме того, естественно потребовать, что если это неравенство выполняется для k0 = k0 (ε), то оно должно выполняться и для любого k > k0 (ε). Число k0 (ε) называется минимальным числом итераций, необходимых для достижении заданной точности ε. Понятно, что чем меньше k0 (ε), тем лучше метод. Примеры одношаговых итерационных методов Для иллюстрации нижеприведенных методов сначала более детально распишем уравнение Ax = f :  a11 x1 + a12 x2 + . . . + a1n xn = f1 ;       a21 x1 + a22 x2 + . . . + a2n xn = f2 ;  . . . . . . . . . . . . . . . . . . . . .      an1 x1 + an2 x2 + . . . + ann xn = fn . Расставляя итерационные индексы над компонентами вектора X, можно получить различные итерационные процессы. Метод Якоби Расставим индекс следующего приближения (k + 1) над диагональными элементами, а индекс k — над остальными:  + a12 xk2 + . . . + a1n xkn = f1 ;  a11 xk+1 1      a21 xk + a22 xk+1 + . . . + a2n xk = f2 ; 1       n 2 . . . . . . . . . an1 xk1 + an2 xk2 + . . . + ann xk+1 = fn . n : Из этой системы легко получить общее уравнение для xk+1 i fi − = xk+1 i i−1 X aij xkj − j=1 n X aij xkj j=i+1 , i = 1, n. (1.12) aii Чтобы от этой расчетной формулы перейти к канонической форме записи метода, представим матрицу A в виде суммы трех матриц — нижнетреугольной, диагональной, и верхнетреугольной:    A1 =     .. . aij ..   ;   . A = A1 + D + A2 ,   a11   ..   . 0  ; D=  ..   . ann    A2 =     .. . aij .. .   .   (1.13) Из (1.12) легко вывести, что A1 xk + Dxk+1 + A2 xk = f. Для приведения к канонической форме (1.11) прибавим и вычтем Dxk : D(xk+1 − xk ) + Axk = f. Отсюда получаем формальное определение метода Якоби через матрицу Bk+1 и константу τk+1 :  Bk+1 = D; τk+1 = 1. Легко видеть, что он неявный и стационарный. 12 Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ Метод Зейделя2 Пусть теперь итерационные индексы (k + 1) стоят не только на диагонали, но и под ней:  = f1 ; + . . . + a1n xkn a11 xk+1 + a12 xk2  1      a21 xk+1 + a22 xk+1 + . . . + a2n xk = f2 ; n 1 2       . . . . . . . . . . an1 xk+1 1 + an2 xk+1 2 + . . . + ann xk+1 n = fn . Относительно xk+1 эта система решается так: сначала из первого уравнения находится xk+1 , потом 1 из второго — xk+1 и т.д. 2 Для записи метода в канонической форме представим матрицу A в виде суммы матриц A1 , D, A2 (как и в методе Якоби): A = A1 + D + A2 , где A1 — нижнетреугольная, D — диагональная, а A2 — верхнетреугольная матрица. Тогда написанная выше итерационная схема будет выглядеть так: A1 xk+1 + Dxk+1 + A2 xk = f. Преобразуя это уравнение, приведем его к каноническому виду: (A1 + D)xk+1 + A2 xk = f ⇐⇒ ⇐⇒ (A1 + D)(xk+1 − xk ) + Axk = f. Отсюда получаем формальное определение метода Зейделя:  Bk+1 = A1 + D; τk+1 = 1. Как несложно заметить, он является неявным и стационарным. Метод релаксации Формально определим метод следующим образом:  Bk+1 = D + ωA1 ; τk+1 = ω, где ω — некий числовой параметр. Нетрудно заметить, что метод Зейделя является частным случаем этого метода при ω = 1. Кроме того, метод релаксации является стационарным и неявным. Метод простой итерации И здесь никакой подробной системы не будет:  Bk+1 τk+1 = E; = τ. Подставив эти равенства в (1.11), получим каноническую форму метода простой итерации: xk+1 − xk + Axk = f. τ На этот раз метод является явным и стационарным. 2 Предложен Л. Зейделем в 1874 г. 1.3. СХОДИМОСТЬ ОДНОШАГОВЫХ СТАЦИОНАРНЫХ МЕТОДОВ 13 Метод Ричардсона В данном случае матрица Bk+1 — единичная, а константа τk+1 определяется по некоторым расчетным формулам. Выражение для τk+1 приводить не будем, запишем только каноническую форму записи метода: xk+1 − xk + Axk = f. τk+1 1.3 Сходимость одношаговых стационарных методов Теорема 1.1 (Достаточное условие сходимости одношагового стационарного итерационного метода). Пусть матрица A — симметрическая положительно определенная, B — положительно определена, (xk+1 − xk ) + Axk = f сходится для любого начального и τ > 0, тогда итерационный процесс B τ τ приближения x0 , если B − A > 0. 2 Доказательство этой теоремы было дано в курсе «Введение в численные методы», кроме того, его можно найти в [1]. Утверждение 1.1. Пусть матрица A — симметрическая положительно определенная, и выполнено n X условие диагонального преобладания: aii > |aij |. Тогда метод Якоби является сходящимся для j=1 i6=j любого начального приближения. Доказательство. Так как матрица A — положительно определена, ∀ x 6= 0 Распишем скалярное произведение поподробнее: hAx, xi = n X n X aij xi xj 6 i=1 j=1 6 1 2 n X n X n X n X i=1 j=1 |aij |x2i + i=1 j=1 |aij ||xi xj | 6 hAx, xi > 0.   x2i + x2j 6 |xi xj | 6 2 n n n X n n n   X X X 1 XX x2i |aii | + |aij | < |aij |x2i = |aji |x2j = {i ↔ j} = 2 i=1 j=1 i=1 j=1 i=1 j=1 i6=j < {так как матрица A — матрица с диагональным преобладанием} < n n X X |aii |xi xi = 2 hDx, xi . < 2 x2i |aii | = 2 i=1 i=1 Таким образом, получаем: 0 < hAx, xi < 2 hDx, xi . Откуда следует, что 2 hDx, xi − hAx, xi > 0 =⇒ 2D − A > 0. То есть матрица (2D − A) положительно определена. В этом случае выполняются все условия теоремы 1.1 при B = D и τ = 1 — это и дает сходимость метода Якоби. Утверждение 1.2 (Сходимость метода релаксации). Пусть B = D + ωA1 , τ = ω, A — симметрическая положительно определенная матрица и ω ∈ [0; 2], тогда метод релаксации является сходящимся для любого начального приближения. 14 Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ Доказательство. Проверим выполнение достаточного условия сходимости стационарного одношагового итерационного метода (теорема 1.1): B− τ A > 0, 2 подставив формулы для B и τ из условия теоремы. Таким образом, для доказательства теоремы достаточно показать, что 2D + 2ωA1 − ωA > 0 (1.14) Представим матрицу A в виде (1.13), тогда: hAx, xi = hA1 x, xi + hDx, xi + x, AT2 x = { x, AT 2 x = hA1 x, xi} = = hDx, xi + 2 hA1 x, xi . Теперь перепишем (1.14): h(2D + 2ωA1 − ωA)x, xi > 0. Подставим сюда разложение матрицы A и раскроем скобки: 2 hDx, xi + 2ω hA1 x, xi − 2ω hA1 x, xi − ω hDx, xi > 0. В результате получаем: (2 − ω) hDx, xi > 0, где hDx, xi > 0 в силу того, что hAx, xi > 0 и примечания к определению положительно определенной матрицы, а множитель (2 − ω) больше нуля в силу условия теоремы. Утверждение 1.3 (необходимое и достаточное условие сходимости метода простой итерации). Пусть матрица A — симметрическая положительно определенная, а параметр τ больше нуля. Тогда для 2 , где λ(A) — сходимости метода простой итерации необходимо и достаточно, чтобы τ < λmax (A) 3 некоторое собственное значение матрицы A. Доказательство. Достаточность. Пусть, согласно требованию теоремы, τ< 2 λmax (A) . Тогда это верно для любого собственного значения матрицы A : τ< 2 λ(A) ⇐⇒ 1− τ λ(A) > 0 2 ⇐⇒ λ(E − τ A) > 0. 2 Это значит, что все собственные значения матрицы E − τ2 A положительны, то есть она положительно определена. Так как B = E, то по теореме 1.1 метод простой итерации является сходящимся. Необходимость. Возьмем в качестве начального приближения x0 = x + µ, где µ — собственный вектор A : Aµ = λmax (A)µ. Рассмотрим наш итерационный процесс: xk+1 − xk + Axk = f. τ (1.15) 3 Далее такое обозначение будет использовано для собственных значений других матриц; иногда (например, в неравенствах) оно может обозначать весь спектр. 1.3. СХОДИМОСТЬ ОДНОШАГОВЫХ СТАЦИОНАРНЫХ МЕТОДОВ 15 Заменим xk на погрешность решения z k : z k = xk − x, где x — точное решение. Тогда для погрешности получим такое выражение: z k+1 − z k + Az k = 0 τ =⇒ z k+1 = (E − τ A)z k ∀k. Выразим z k через предыдущие значения: z k = (E − τ A)z k−1 = . . . = (E − τ A)k z 0 . Заметим, что z 0 = µ : (E − τ A)k µ = (E − τ A)k−1 (µ − τ Aµ) = (E − τ A)k−1 (1 − τ λmax (A))µ = . . . = (1 − τ λmax (A))k µ. Таким образом: z k = (1 − τ λmax (A))k µ. Посчитаем норму z k : ||z k || = |1 − τ λmax (A)|k ||µ||. k→∞ Так как итерационный процесс сходится ( ||z k || −→ 0 ), то |1−τ λmax (A)| < 1. Раскрыв знак модуля, получим:  2 1 − τ λmax (A) < 1; . =⇒ 0 < τ < 1 − τ λmax (A) > −1. λmax (A) Утверждение доказано. Сходимость стационарных методов Теперь установим необходимое и достаточное условие сходимости для стационарных одношаговых итерационных методов. Для этого перейдем в итерационном процессе B xk+1 − xk + Axk = f τ B z k+1 − z k + Az k = 0. τ (1.16) к погрешности z k = xk − x : Отсюда получим формулу z k+1 = (E − τ B −1 A)z k , и обозначим S = E − τ B −1 A. Матрица S называется матрицей перехода. Теперь выражение для z k+1 выглядит так: z k+1 = Sz k . (1.17) Теорема 1.2 (Критерий сходимости одношагового стационарного итерационного метода). Итерационный метод (1.16) сходится для любого начального приближения x0 тогда и только тогда, когда для всех собственных значений λ(S) выполнено |λ(S)| < 1. Доказательство. Необходимость. Пусть µ — собственный вектор S, соответствующий собственному значению λ(S). Рассмотрим вектор начального приближения x0 = µ + x, тогда z 0 = x0 − x = µ. Все погрешности связаны соотношением (1.17). Выразим из него z k z k = Sz k−1 = S 2 z k−2 = . . . = S k z 0 , где z 0 = µ, тогда z k = S k µ = S k−1 λ(S)µ = . . . = λk (S)µ. 16 Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ В силу сходимости верно, что k→∞ ||z k || −→ 0 =⇒ k→∞ ||λk (S)µ|| −→ 0 =⇒ k→∞ |λ(S)|k ||µ|| −→ 0. Отсюда следует, что |λ(S)| < 1. Достаточность. Пусть матрица S — матрица простой структуры, причем |λ(S)| < 1. Из этого следует, что существуют собственные вектора, образующие ортонормированный базис {ξi }ni=1 . n X Пусть z 0 = ci ξi — разложение z 0 по базису. Тогда получим следующее выражение для погрешi=1 ности: zk = S k z0 = S k n X ci ξi = S k−1 n X ci Sξi = S k−1 ci λi ξi = . . . = n X ci λki ξi . i=1 i=1 i=1 i=1 n X Рассмотрим норму погрешности: ||z k || = n X i=1 в силу того что n X ci λki ξi 6 n X |ci | · |λi |k ||ξi || 6 {ρ = max |λi (S)| < 1} 6 ρk i=1 i n X k→∞ |ci | · ||ξi || −→ 0, i=1 |ci | · ||ξi || ограничено. i=1 Примечание. 1. При доказательстве достаточности предполагалась, что матрица S — матрица простой структуры. Теорема верна и без этого предположения, которое сделали для упрощения доказательства. 2. На первый взгляд кажется, что область применения этого утверждения широка, но на практике найти весь спектр матрицы перехода непросто (иногда сложнее, чем решить систему прямым методом). Пусть для погрешности итерационного метода справедливо неравенство: ||xk − x|| 6 q k ||x0 − x||, где q ∈ (0, 1). (1.18) Определение. Итерационный метод, погрешность итерационного приближения которого удовлетворяет (1.18), сходится со скоростью геометрической прогрессии со знаменателем q. Потребуем, взяв ε > 0, чтобы q k < ε. Тогда для погрешности итерационного метода будет выполнена оценка: ||xk − x|| < ε||x0 − x||. 1 А это означает, что к k-й итерации погрешность начального приближения уменьшится в раз. ε Таким образом, число итераций, необходимых для достижения требуемой точности ε, будет & ' ln 1ε k > k0 (ε) = . ln 1q Определение. Число ln 1q называется скоростью сходимости итерационного метода. Примечание. Вообще говоря, число q из вышеприведенных неравенств определяется неоднозначно. Для формальности можно считать, что это минимальное из всех q, удовлетворяющих (1.18). 1.4. ОЦЕНКА ПОГРЕШНОСТИ ОДНОШАГОВЫХ СТАЦИОНАРНЫХ МЕТОДОВ 1.4 17 Оценка погрешности одношаговых стационарных методов Далее мы будем активно использовать матричные неравенства, например A > B. Оно означает, что для всех x hAx, xi > hBx, xi . Введем норму вектора x, порожденную симметричной положительно определенной матрицей A : p ||x||A = hAx, xi. Теорема 1.3 (Оценка погрешности стационарных одношаговых итерационных методов). Пусть матрицы A, B симметричны и положительно определены, и существуют такие положительные константы γ1 , γ2 , что γ1 B 6 A 6 γ2 B. (1.19) Тогда итерационный метод, задаваемый уравнением B xk+1 − xk + Axk = f, τ где τ = 2 , γ1 + γ2 (1.20) сходится для любого начального приближения x0 со скоростью геометрической прогрессии: ||xk − x||A 6 q k ||x0 − x||A , где q = γ1 1−ξ , ξ= . 1+ξ γ2 Доказательство. Перейдем от приближений xk к погрешности на k-й итерации: z k = xk − x. Как уже было показано, уравнение (1.20) можно переписать так: B z k+1 − z k + Az k = 0, τ и из него следует, что z k+1 = Sz k , S = E − τ B −1 A. (1.21) Теперь установим справедливость такого неравенства: ||z k+1 ||A 6 q||z k ||A . (1.22) 1 Известно, что для матрицы A = AT > 0 существует такая матрица, обозначаемая A 2 , что 1 1 (A )2 = A, причем A 2 = (A 2 )T > 0. Домножим (1.21) слева на эту матрицу: 1 2 1 1 A 2 z k+1 = A 2 Sz k 1 ⇐⇒ 1 1 1 1 A 2 z k+1 = A 2 SA− 2 A 2 z k . 1 1 Обозначив ω k = A 2 z k и S = A 2 SA− 2 , получим, что ω k+1 = Sω k . 1 1 1 1 1 1 Заметим, что S = A 2 SA− 2 = A 2 (E − τ B −1 A)A− 2 = E − τ A 2 B −1 A 2 . 1 1 Обозначив C = A 2 B −1 A 2 , получим, что S = E − τ C. (1.23) Легко проверить, что матрицы C и S будут симметричны, а C — еще и положительно определена (это понадобится нам чуть позже). Теперь преобразуем ||z k ||A : rD q E rD 1 E q 1 1 1 k k k k k ||z ||A = hAz , z i = A2 A2 z , z = A 2 z k , A 2 z k = hω k , ω k i = ||ω k ||. 18 Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ Таким образом, от доказательства неравенства (1.22) можно перейти к доказательству того, что ||ω k+1 || 6 q||ω k ||. (1.24) Преобразуем ω k+1 : E D 2 ||ω k+1 ||2 = ω k+1 , ω k+1 = Sω k , Sω k = S ω k , ω k . 2 Теперь будем искать такие q, что S 6 q 2 E — как видно из последних преобразований, этого будет достаточно для доказательства (1.24). 2 Итак, предположим, что S 6 q 2 E. Тогда, по свойству симметричных матриц, и в силу того, что C положительно определена, это эквивалентно тому, что −qE 6 S 6 qE ⇐⇒ {подставим (1.23)} ⇐⇒ −qE 6 E − τ C 6 qE ⇐⇒   τ C 6 (1 + q)E; τ E 6 (1 + q)C −1 ; −1 ⇐⇒ ⇐⇒ {домножим на C } ⇐⇒ τ C > (1 − q)E. τ E > (1 − q)C −1 . Из этой системы вытекает такое двойное неравенство: 1 + q −1 1 − q −1 C 6E6 C . τ τ 1 1 1 1 Из того, что C = A 2 B −1 A 2 , следует, что C −1 = A− 2 BA− 2 . Таким образом, предыдущее двойное неравенство эквивалентно такому: 1 1 1 − q −1 1 + q −1 A 2 BA− 2 6 E 6 A 2 BA− 2 . τ τ 1 Умножив его слева и справа на A 2 , получим 1+q 1−q B6A6 B. τ τ (1.25) Очень похоже на одно из условий теоремы. Покажем, что 1−q = γ1 , τ 1+q = γ2 . τ (1.26) Действительно: 1 1 − γγ22 −γ 1−q (γ2 + γ1 ) − (γ2 − γ1 ) +γ1 = = γ1 ; = 2 τ 2 γ2 +γ1 1 1 + γγ22 −γ 1+q (γ2 + γ1 ) + (γ2 − γ1 ) +γ1 = = γ2 . = 2 τ 2 γ2 +γ1 Итак, мы показали, что неравенство ||z k+1 ||A 6 q||z k ||A эквивалентно матричному неравенству γ1 B 6 A 6 γ2 B из условия теоремы. Таким образом, показано, что ||z k+1 ||A 6 q||z k ||A для любого k. Тогда, переходя к более ранним членам последовательности, получим: ||z k ||A 6 q||z k−1 ||A 6 . . . 6 q k ||z 0 ||. Из этого напрямую следует, что ||xk − x||A 6 q k ||x0 − x||A . Теорема доказана. 1.4. ОЦЕНКА ПОГРЕШНОСТИ ОДНОШАГОВЫХ СТАЦИОНАРНЫХ МЕТОДОВ 19 Замечание 1. В случае, когда ξ мало, можно получить такую оценку для скорости сходимости:     1+ξ 1 2ξ ln = ln = ln 1 + ≈ 2ξ. q 1−ξ 1−ξ Замечание 2. Из теоремы 1.3 следует, что выбор чисел γ1 , γ2 напрямую влияет на скорость сходимости. Для выяснения их возможных значений рассмотрим произвольный собственный вектор µ матрицы B −1 A : B −1 Aµ = λ(B −1 A)µ. Это эквивалентно тому, что Aµ = λ(B −1 A)Bµ (мы просто домножили на B ). Как известно, неравенство γ1 B 6 A 6 γ2 B (1.27) означает, что γ1 hBx, xi 6 hAx, xi 6 γ2 hBx, xi для любого x. Положим x = µ и преобразуем это двойное неравенство: γ1 hBµ, µi 6 hAµ, µi 6 γ2 hBµ, µi ⇐⇒ ⇐⇒ γ1 hBµ, µi 6 λ(B −1 A) hBµ, µi 6 γ2 hBµ, µi =⇒ =⇒ γ1 6 λ(B −1 A) 6 γ2 . Так как собственный вектор мы выбирали произвольно, то получаем, что γ1 6 λmin (B −1 A), γ2 > λmax (B −1 A). Таким образом, наиболее точными константами, с которыми выполняется неравенство (1.27), являются константы γ1 = λmin (B −1 A), γ2 = λmax (B −1 A). В этом случае параметр τопт = 2 λmin (B −1 A) + λmax (B −1 A) называется оптимальным итерационным параметром. Кроме того, так как мы берем q = 1− 1+ γ1 γ2 γ1 γ2 , то наилучшим вариантом будет как раз λmin (B −1 A) = λmax (B −1 A). Все, мы закончили с главной теоремой этого параграфа, пора построить пример использования наших методов. Примечание. Выполнение неравенства (1.27) при γ1 = λmin (B −1 A), γ2 = λmax (B −1 A) следует из возможности построения в данном линейном пространстве базиса из собственных векторов. Модельная задача. Сравнение скорости сходимости различных итерационных методов Рассмотрим краевую задачу:  −u00 (x) = f (x), 0 < x < 1; u(0) = u(1) = 0. Найдем ее решение, используя численные методы. Для этого сначала разделим отрезок [0; 1] на N равных промежутков длины h = N1 , обозначив границы отрезков как xi : 0 = x0 < x1 < x2 < . . . < xN = 1; h = xi+1 − xi . 20 Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ Воспользуемся тем, что u00 (xi ) можно приблизить 2-й разностной производной: u00 (xi ) ≈ u(xi+1 ) − 2u(xi ) + u(xi−1 ) h2 def Пусть y(x) — приближение нашей функции. Потребуем, чтобы yi = y(xi ) = u(xi ). Тогда, зная, что u (xi ) = f (xi ), и используя разностное приближение, получим такую систему уравнений относительно значений y в узлах сетки: 00 − y(xi+1 ) − 2y(xi ) + y(xi−1 ) = f (xi ), h2 i = 1, N − 1. Обозначим fi = f (xi ). Из краевых условий можно получить, что y0 = yN = 0. Тогда получим такую систему:   − yi+1 − 2yi + yi−1 = fi , i = 1, N − 1; h2 (1.28)  y0 = yN = 0. Понятно, что решив ее, мы получим приближенные значения u(xi ) = yi . Теперь обозначим Y = (y1 , y2 , . . . , yN −1 )T ; F = (f1 , f2 , . . . , fN −1 )T . Тогда систему (1.28), использовав значения y0 и yN , можно переписать в ния AY = F,   2 1  −  h2  2 −1 h2    1   .. 2 −   .   2 −1 2 . . . h2  h      1    .. .. .. .. .. где A =  = 2 . . . . .   h          0 ... 0 . . . . . . − h12        1 2  − 2 h h2 виде матричного уравне(1.29)        .. . .    .. . −1  −1 2 Легко видеть, что A = AT . Покажем, что A > 0 — для этого обоснуем положительность всех собственных значений. Утверждается, что уравнение для поиска собственных значений Aξ = λξ в каком-то смысле эквивалентно задаче Штурма-Лиувилля:  X 00 + λX = 0, 0 < x < 1; X(0) = X(1) = 0, решением которой будут собственные функции Xm (x) = sin πmx. Отсюда делается вывод, что для собственных значений матрицы A справедлива формула: λm 4 = 2 h  πmh sin 2 2 , m = 1, N − 1. Таким образом, все λm положительны и матрица A положительно определена. Легко проверить, что справедлива следующая цепочка неравенств: 0 < λ1 < λ2 < . . . < λN −1 , 1.5. ПОПЕРЕМЕННО-ТРЕУГОЛЬНЫЙ ИТЕРАЦИОННЫЙ МЕТОД 21 2 2   4 4 πh πh причем λ1 = λmin = 2 sin и λN −1 = λmax = 2 cos . ҷ Будем решать уравнение (1.29) h 2 h 2 методом простой итерации: xk+1 − xk + Axk = f. τ Для него справедлива теорема 1.3 о скорости сходимости, так как A = AT > 0 (по доказанному) и B = B −1 = E — симметричная и положительно определенная матрица. Тогда, согласно замечанию 2 к теореме 1.3 мы можем взять γ1 = λmin (A), γ2 = λmax (A) и положить τ= 2 2 h2 = . = 4 λmin (A) + λmax (A) 2 2 πh 2 ( πh )) (sin ( ) + cos 2 2 h2 При этом для погрешности верна оценка ||xk − x||A 6 q k ||x0 − x||A , где q = λmin 1−ξ , ξ= . 1+ξ λmax Тогда, если в качестве условия выхода из итерационного процесса использовать неравенство ||xk − x||A 6 ε||x0 − x||A , то число итераций, необходимых для достижения точности ε, равно & ' ln 1ε k0 (ε) = . ln 1q В нашем случае ξ достаточно мало, поэтому верна такая оценка для ln λmin (A) 1 = 2 tg2 ln ≈ 2ξ = 2 q λmax (A)  πh 2  ≈ 1 : q 1 π2 2π 2 h2 = {h = } = . 4 N 2N 2 Отсюда следует, что k0 (ε) = N2 1 2N 2 1 ≈ ln . ln π2 ε 5 ε К примеру, для ε = 0, 5 · 10−4 ≈ e−10 имеем: k0 (ε) = 2N 2 . Таким образом, например, если разбить отрезок на 10 частей, получив систему уравнений с числом уравнений N = 10, требуется выполнить 200 итераций, в случае разбиения отрезка на 100 частей (соответственно N = 100 уравнений) уже требуется выполнить 20000 итераций для получения решения с заданной точностью ε. Впоследствии мы увидим, что по сравнению с другими методами метод простой итерации является очень медленно сходящимся. Впоследствии мы увидим, что по сравнению с другими методами метод простой итерации является очень медленно сходящимся (хотя это понятно и так, 20000 итераций — это же убить себя можно...). 1.5 Попеременно-треугольный итерационный метод В этом пункте мы рассмотрим еще один итерационный метод решения СЛАУ. Этот метод работает быстрее метода простой итерации, но и параметры B и τ в нем выбираются не так тривиально. 22 Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ Рассмотрим одношаговый итерационный метод: B xk+1 − xk + Axk = f. τ Возьмем матрицу B как произведение верхнетреугольной и нижнетреугольной матриц особого вида: B = (E + ωR1 )(E + ωR2 ), (1.30) где матрицы  r11   R1 =     .. . rij ..   ,   .  r11   R2 =     .. . rij ..      . rnn rnn такие, что их сумма R1 + R2 = A, а их элементы на их диагонали равны половине соответствующих элементов A. Теорема 1.4 (Достаточное условие сходимости ПТИМ). Пусть A — симметрическая положительно τ (просто согласование ω определенная матрица, а матрица B задается формулой (1.30), где ω > 4 с τ ). Тогда соответствующий попеременно-треугольный итерационный метод сходится. Доказательство. Преобразуем матрицу B к виду, показывающему выполнение достаточного условия сходимости. Раскроем скобки в представлении (1.30) B = (E + ωR1 )(E + ωR2 ) = E + ω(R1 + R2 ) + ω 2 R1 R2 = = {так как R1 + R2 = A} = E + ωA + ω 2 R1 R2 . Заметим, что в силу симметричности матрицы A, матрицы R1 и R2 связаны следующим образом: R1 = R2T . Вычтем и добавим ωA : B = E − ωA + ω 2 R1 R2 + 2ωA = (E − ωR1 )(E − ωR2 ) + 2ωA. Покажем, что выполняется достаточное условие сходимости из теоремы 1.1. Для этого распишем такое скалярное произведение: hBx, xi = h(E − ωR1 )(E − ωR2 )x, xi + 2ω hAx, xi . Перенесем скобку (E − ωR1 ) в правую часть скалярного произведения, а потом воспользуемся тем, что (E − ωR1 )∗ = E − ωR1T = E − ωR2 , получим: hBx, xi = h(E − ωR2 )x, (E − ωR2 )xi + 2ω hAx, xi = ||(E − ωR2 )x||2 + 2ω hAx, xi . В силу неотрицательности нормы hBx, xi > 2ω hAx, xi , что равносильно B > 2ωA, а, так как в τ τ условии теоремы мы взяли ω > , то выполняется достаточное условие сходимости: B > A. Теорема 4 2 доказана. Теорема 1.5. Пусть матрица A — симметрическая положительно определенная, и пусть существуют такие константы δ > 0 и ∆ > 0, что A > δE, ∆ A > R1 R2 4 (R1 и R2 — такие же, как в теореме 1.4). 1.5. ПОПЕРЕМЕННО-ТРЕУГОЛЬНЫЙ ИТЕРАЦИОННЫЙ МЕТОД 23 Тогда ПТИМ сходится при любом начальном приближении со скоростью геометрической прогрессии, причем при значениях параметров  √ √ √  δ ∆δ δ∆ 2   τ= √ , γ2 = ·√ ; , где γ1 =  γ1 + γ2 2 4 ∆+ δ  2    ω=√ . ∆δ √ δ 1− ξ k+1 k √ , ξ= . скорость сходимости наилучшая: ||x − x||A 6 q||x − x||A , где q = ∆ 1+3 ξ δ Доказательство. (1) Для корректности последующих формул сначала убедимся, что 6 1. ∆ ∆ Из условий теоремы, так как A > δE и A > R1 R2 , следует выполнение неравенств: 4 δ||x||2 = δ hx, xi 6 hAx, xi , hR1 R2 x, xi 6 ∆ hAx, xi . 4 (1.31) Левую часть второго неравенства можно переписать в виде: hR1 R2 x, xi = hR2 x, R2 xi = ||R2 x||2 . (1.32) Из (1.31) очевидно следует, что 2 δ||x||2 6 hAx, xi = hAx, xi . hAx, xi Кроме того, в силу представления матрицы A hAx, xi = hR1 x, xi + hR2 x, xi = hR2 x, xi + x, R1T x = 2 hR2 x, xi , поэтому 2 2 hAx, xi 4 hR2 x, xi 4||R2 x||2 · ||x||2 = 6 {| hu, vi | 6 ||u|| · ||v||} 6 . hAx, xi hAx, xi hAx, xi Распишем это с помощью (1.31) и (1.32): 4||R2 x||2 · ||x||2 4∆ hAx, xi ||x||2 6 = ∆||x||2 . hAx, xi 4 hAx, xi В итоге получаем, что δ||x||2 6 ∆||x||2 для любого x 6= 0, или, что то же самое, δ 6 ∆ =⇒ ξ = δ 6 1. ∆ (2) Теперь фиксируем некоторое ω и посмотрим, что происходит. Из доказательства теоремы 1.4 следует, что B > 2ωA , откуда 1 A6 B. 2ω 1 некоторой константой γ2 , и перейдем к доказательству симметричного неравенства Обозначим 2ω для γ1 . Преобразуем условия нашей теоремы:  (   E 6 1 A; A > δE; δ =⇒ ∆ ∆  A > R1 R2 .  R R A. 1 2 6 4 4 24 Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ Использовав эти неравенства в очевидном (хотя и ранее доказанном) равенстве B = E + ωA + ω 2 R1 R2 , получим, что 1 ∆ A + ωA + ω 2 A. δ 4 Отсюда легко получается выражение для A : −1  ω2 ∆ 1 +ω+ B. A> δ 4 −1  1 ω2 ∆ Обозначив γ1 = +ω+ , получим, что выполняется условие теоремы 1.3, то есть δ 4 B6 γ1 B 6 A 6 γ2 B. Кроме того, мы можем выбирать различные γ1 и γ2 , варьируя параметр ω. Теперь воспользуемся теоремой 1.3, и получим: ||xk − x||A 6 q||xk−1 − x||A , 1− γ1 γ2 γ1 γ2 2 . Здесь γ1 и γ2 (а, значит, и q ) неявно зависят от ω. +1 1+ При уменьшении q скорость сходимости возрастает. Выберем ω так, чтобы q было минимально. γ2 γ2 (ω) Очевидно, для этого надо минимизировать отношение > 1, и получим для . Обозначим g = γ1 γ1 (ω) него такую формулу: 1 ∆ + ω + ω2 δ 4 = 1 + ω∆ + 1 . g= 2ω 2 8 2ωδ Чтобы найти экстремум, возьмем производную от g по ω : где q = =1− γ2 γ1 g0 = ∆ 1 . − 8 2δω 2 2 (заметим, что вторая произПриравняв производную к нулю, получим точку минимума ω = √ δ∆ водная будет больше нуля). Теперь найдем значения γ1 и γ2 в этой точке.  √ √ 1 δ δ∆ 1    √ ; ·√ = 2 =  γ1 = 1 2 ω2 ∆ √ 2 ∆+ δ δ + δ∆ δ +ω+ 4 ω= √2 δ∆ √     γ = 1 = δ∆ . 2 2ω 4 Воспользовавшись найденными значениями γ1 и γ2 , подсчитаем q : √ γ2 − γ1 q= = γ2 + γ1 δ∆ 4 √ δ∆ 4 √ − + δ 2 √ δ 2 · · √ √ √ δ∆ ∆+ δ √ √ √ δ∆ ∆+ δ = q √ √ √ √ √ √ √ √ δ 1 − ∆ δ∆ + δ δ∆ − 2 δ δ∆ ∆− δ ∆ √ √ √ √ √ = q . =√ √ =√ δ ∆ δ∆ + δ δ∆ + 2 δ δ∆ ∆+3 δ 1+3 ∆ Итак, мы получили, что наибольшая скорость сходимости достигается при параметрах, указанных в условии теоремы. Это, в общем-то, и требовалось доказать. 1.6. ЧЕБЫШЕВСКИЙ НАБОР ИТЕРАЦИОННЫХ ПАРАМЕТРОВ 25 Оценим скорость сходимости этого метода при q 6= 0. По определению q   r r δ δ 1 + 3 4∆ δ δ 1 ∆   q q ≈4 ln = ln = ln 1 + 4 + . q ∆ ∆ 1− δ 1− δ ∆ ∆ Применив попеременно-треугольный итерационный метод к решению системы уравнений (1.29) с трехдиагональной матрицей A, получим оценку на k0 (ε) : ln 1 k0 (ε) = qε . δ 4 ∆ δ λmin (B −1 A) Как уже показывалось для метода простой итерации, = = tg2 ∆ λmax (B −1 A) откуда получаем оценку на k0 (ε) : k0 (ε) ≈ ln 1ε 4 πh 2 =  πh 2  ≈ π 2 h2 , 4 ln 1ε N = {ε = 5 · 10−5 ≈ e−10 } ≈ 1, 6N. 2π Это уже неплохой результат (гораздо меньше, чем в случае метода простой итерации), например, взяв N = 100, мы должны будем сделать 160 итераций вместо 20000 в МПИ. 1.6 Чебышевский набор итерационных параметров В данном разделе мы будем применять для уже порядком заколебавшего нас уравнения Ax = f такую итерационную схему: xl − xl−1 B + Axl−1 = f. (1.33) τl Теперь фиксируем число итераций (скажем, k, а то m — маловато...), и постараемся выбрать τl так, чтобы погрешность на k-й итерации была минимальной: ||z k || −→ inf . τl Чтобы получить оценку на погрешность, перейдем от l-го приближения к погрешности z l = xl − x, как это уже делали раньше сто раз. Тогда итерационный процесс будет выглядеть просто прекрасно: B z l − z l−1 + Az l−1 = 0. τl Отсюда без проблем получается выражение для z l : z l = (E − τl B −1 A)z l−1 , l = 1, k. Рекурсивно применяя эту формулу для z k , получим очень-очень длинную формулу: z k = (E − τk B −1 A)(E − τk−1 B −1 A) . . . (E − τ1 B −1 A)z 0 . (1.34) Для тех, кто в танке, повторяем: необходимо так подобрать τl , чтобы минимизировать z k . Оказывается, что эта задача не только имеет решение, но нам даже разрешили его опубликовать (тоже мне, ценность великая...). Это решение, а, точнее, набор τl называется чебышевским набором итерационных параметров. Прямо корпускулярно-волновая теория света, мдаа... А теперь приколитесь — следующая теорема пойдет без доказательства (кому неймется, может посмотреть в [1]). 26 Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ Теорема 1.6. Пусть матрицы A и B симметричны и положительно определены, τl считаются по следующей формуле (k — фиксировано):  2   ;  τ0 = λ −1  A) + λmax (B −1 A) min (B     τ0 1−ξ λmin (B −1 A) τl = , где ρ0 = , ξ= ;  1 + ρ0 tl 1+ξ λmax (B −1 A)        tl = cos (2l − 1)π , l = 1, k, 2k а xl вычисляется по формуле (1.33). Тогда погрешность ||xk − x||A будет минимальной среди всех возможных, и для нее справедливо неравенство: ||xk − x||A 6 qk ||x0 − x||A , √ ρk1 1− ξ √ . где qk = , ρ = 1 1 + ρ2k 1+ ξ 1 Замечание 1. Подробную инструкцию о том, что надо выкурить, чтобы догадаться до таких формул, можно найти в [2]. Замечание 2. Найдем, при каких k выполняется условие остановки итерационного процесса: ||xk − x||A 6 ε||x0 − x||A . (1.35) Из теоремы 1.6 следует, что оно заведомо (т.е. всегда (т.е. даже ночью)) выполняется, если qk 6 ε. Это, в свою очередь, эквивалентно тому, что ρk1 6 ε ⇐⇒ ε(ρk1 )2 − ρk1 + ε > 0. 1 + ρ2k 1 Корнями этого квадратного трехчлена будут ρk1 = 1± √  k ρ1 = ε; 1 − 4ε2 1 ± (1 − 2ε2 ) ≈ =⇒  1 2ε 2ε ρk1 = + ε. ε Отсюда следует, что qk 6 ε при  ρk1 6 ε;  ρk1 > 1 + ε. ε В теореме 1.6 мы выбирали ρk1 так, что оно было меньше единицы. Поэтому второй случай нам не подойдет ( 1ε  1), из чего следует, что qk 6 ε ⇐⇒ ρk1 6 ε. Таким образом, (1.35) верно при & k > k0 (ε) = ln 1ε ln ρ11 ' . В данном случае за скорость сходимости принимается ln ρ11 . Оценим этот параметр: √   p p 1+ ξ 1 √ = ln 1 + 2 ξ + O(ξ) ≈ 2 ξ. = ln ln ρ1 1− ξ (1.36) 1.6. ЧЕБЫШЕВСКИЙ НАБОР ИТЕРАЦИОННЫХ ПАРАМЕТРОВ 27 Посмотрим, как работает этот метод (1.33) (для упрощения сделаем матрицу B единичной) применительно к нашей модельной задаче (1.28). Напомним, матрица системы в примере была вида:   2 −1     −1 2 . . . 0      1   .. .. .. A= 2 . . . .   h     .. ..   . . −1   −1 2 Собственные значения для нее выглядят так: λmin = πh 4 sin2 ( ), h2 2 λmax = 4 πh cos2 ( ), h2 2 min поэтому ξ = λλmax = tg2 ( πh 2 ). И снова, в который раз (и не в последний, по-моему) зададимся ε = 0, 5 · 10−4 ≈ e−10 . В этом случае из (1.36) можно поиметь следующую формулу: k0 (ε) = 1 10 10 10 = {h = } = N ≈ 3, 5N. ≈ πh πh N π 2 tg( 2 ) Тогда для N = 10 нам понадобится 35 итераций, а для N = 100 — 350 итераций, то есть уже можно себя не убивать... Откровенно говоря, этот метод немеряно круче (на порядок) метода простой итерации, да и то, что схема явная, не может не радовать. С другой стороны, этот метод медленнее, чем попеременно-треугольный, но там-то надо было матрицы обращать, а здесь этого иногда можно и не делать... Упорядоченный набор чебышевских параметров Вообще говоря, в описанном в предыдущем разделе процессе не предполагалось монотонности убывания погрешности. Обычно ее и нет, что на практике при реализации метода с чебышевским набором параметров часто приводит к возникновению неприятных ситуаций (переполнение регистров и прочий отстой). Для предотвращения таких дел набор параметров можно особым образом упорядочить, чтобы как-то скомпенсировать увеличение погрешности за счет ее уменьшения на соседних итерациях. К примеру, в нашем методе отличия на разных итерациях заключались в изменяющемся параметре ). По-другому это можно записать так: τl , который в свою очередь зависит от tl = cos( (2l−1)π 2k   π  (2l − 1)π tl = cos = cos θlk , θlk = {1, 3, . . . , 2k − 1}. 2k 2k Одним из способов упорядочивания параметров τl является простая перестановка значений θlk . Мы рассмотрим пример такой перестановки, когда k является степенью двойки: k = 2p . Тогда значения будут считаться по рекуррентным формулам:  θ11 = 1;  2m θ = θim , i = 1, m; (1.37)  2i−1 2m 2m θ2i = 4m − θ2i−1 , i = 1, m. — последние две формулы необходимо применять для всех m = 1, 2, 4, . . . , 2p−1 . Проще это будет понять на примере: k = 24 = 16 : 28 Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ θ11 2 θ1,2 4 θ1,...,4 8 θ1,...,8 m=1: m=2: m=4: 1; 1, 4-1=3; 1, 8-1=7, 3, 8-3=5; 1, 15, 7, 9, 3, 13, 5, 11; — потом по θlk строятся tl и τl , при этом получается более устойчивый метод. В следующей теореме, которую мы опять примем без доказательства, будет показано, как можно соединить попеременно-треугольный итерационный метод с упорядоченным набором чебышевских параметров. Теорема 1.7. Пусть уравнение Ax = f решается по следующей итерационной схеме (число шагов — k — фиксировано): xl − xl−1 B + Axl−1 = f, l = 1, k, τl где B = (E + ωR1 )(E + ωR2 ), R1 + R2 = A (матрицы R1 , R2 имеют тот же смысл, что и в теореме 1.4). Пусть также существуют такие положительные константы δ и ∆, что A > δE, Тогда, если выбирать    τ0 =          τl =            tl = ∆ A > R 1 R2 . 4 итерационные параметры τl так: 2 ; γ1 + γ2 √ δ ∆ √ ; √ = γ1 τ0 1−ξ 2( ∆ + δ) , ξ= , , ρ0 = √  1 + ρ0 tl 1+ξ γ2  δ∆   γ2 = , 4 π  cos θk , θlk считаются по формулам вида (1.37). 2k l      γ1 то для достижения точности ||A(xk − x)||B −1 6 ε||A(x0 − x)||B −1 ln 2ε достаточно выполнить k0 (ε) = √ q итераций. δ 2 2 4 ∆ Доказательство. Доказательство можно найти в [1]. Теперь посмотрим, какие улучшения даст этот метод в модельной задаче (1.28). Напомним, что константы δ, ∆ можно взять такими:  4  2 πh   δ = λmin (A) = 2 sin ( ); h 2  4 πh   ∆ = λmax (A) = cos2 ( ). h2 2 Тогда при ε = 0, 5 · 10−4 ≈ e−10 получим, что √ √ ln 2 + 10 ln 2 + 10 1 (ln 2 + 10) N √ k0 (ε) = √ q ≈ √ q = {h = } = ≈ 3 N. N 2 π 2 2 4 tg2 ( πh 2 2 πh 2 ) 2 1.7. ОДНОШАГОВЫЕ ИТЕРАЦИОННЫЕ МЕТОДЫ ВАРИАЦИОННОГО ТИПА 29 Таким образом, для N = 10 нам понадобится 10 итераций, а для N = 100 — 30 итераций. Превосходство в скорости над ранее рассмотренными методами очевидно, однако недостатки тоже налицо: необходимо уметь обращать матрицу B и много знать о спектре матрицы A. Заметим, что это общие требования у быстрых итерационных методов — но, как мы увидим далее, они не являются необходимыми. 1.7 Одношаговые итерационные методы вариационного типа В данном разделе мы вновь будем работать с одношаговыми методами. Запишем итерационную схему метода: xk − xk−1 B + Axk−1 = f. τk Как и раньше, перейдем от xk к погрешности z k = xk − x : B z k − z k−1 + Az k−1 = 0. τk Отсюда выразим (k + 1)-ю погрешность: z k+1 = (E − τk+1 B −1 A)z k . (1.38) Теперь зафиксируем некоторую матрицу D : D = DT > 0. Как уже говорилось, в этом случае 1 существует такая матрица D 2 , что 1 1 D 2 D 2 = D, 1 1 D 2 = (D 2 )T . Далее мы будем подбирать τk+1 так, чтобы минимизировать ||z k+1 ||D (считаем, что z k уже задана) — такой способ построения итерационного процесса называется локальной минимизацией. Согласно определению ||z k+1 ||D , имеем: rD q E ||z k+1 ||D = hDz k+1 , z k+1 i = 1 1 D 2 z k+1 , D 2 z k+1 = ||y k+1 ||, 1 где y k+1 = D 2 z k+1 . То есть, минимизация ||z k+1 ||D эквивалентна минимизации ||y k+1 ||. Преобразуем немного выражение для z k+1 в (1.38): 1 1 z k+1 = (E − τk+1 B −1 A)D− 2 D 2 z k . 1 Домножим обе части равенства на D 2 : 1 D 2 z k+1 1 ⇐⇒ D 2 z k+1 1 1 1 = D 2 (E − τk+1 B −1 A)D− 2 D 2 z k ⇐⇒ = 1 1 1 (E − τk+1 D 2 B −1 AD− 2 )D 2 z k . Согласно принятым обозначениям, имеем: 1 1 y k+1 = (E − τk+1 D 2 B −1 AD− 2 )y k . 1 1 Переобозначив C = D 2 B −1 AD− 2 , получим более короткую запись: y k+1 = (E − τk+1 C)y k . Вспомним, что нам надо уменьшать ||y k+1 ||. Более удобно работать с квадратом этого выражения: ||y k+1 ||2 = y k+1 , y k+1 = (E − τk+1 C)y k , (E − τk+1 C)y k . 30 Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ Из линейности скалярного произведения получаем: 2 ||y k+1 ||2 = ||y k ||2 + τk+1 ||Cy k ||2 − 2τk+1 y k , Cy k . Потребуем дополнительное условие: положительную определенность C. В этом случае квадрат нормы ||y k+1 || примет удобный для нас вид: # Cy k , y k = ||y || = ||y || + Cy , Cy − 2τk+1 hCy k , Cy k i #2 " 2 Cy k , y k Cy k , y k k 2 k k . − = ||y || + Cy , Cy τk+1 − hCy k , Cy k i hCy k , Cy k i " k+1 2 k 2 k k 2 τk+1 Очевидно, что наименьшее значение ||y k+1 || достигается, когда второе слагаемое обращается в нуль, то есть, когда: Cy k , y k τk+1 = hCy k , Cy k i — это выражение больше нуля, так как C > 0. Теперь преобразуем эту формулу, используя ранние 1 1 1 обозначения (y k = D 2 z k , C = D 2 B −1 AD− 2 ): D τk+1 = D 1 1 1 1 1 D 2 B −1 AD− 2 D 2 z k , D 2 z k 1 1 E 1 1 1 D 2 B −1 AD− 2 D 2 z k , D 2 B −1 AD− 2 D 2 z k E= DB −1 Az k , z k . hDB −1 Az k , B −1 Az k i (1.39) Получилось, что для получения оптимального выражения для τk+1 необходимо знать z k . Далее с помощью варьирования матрицы D сведем это к вычислению известных нам величин. Для этого введем некоторые дополнительные понятия, связанные с погрешностью итерационного метода. Рассмотрим произведение Az k = A(xk − x) = Axk − Ax = Axk − f = rk — полученная разность вычислима на каждом шаге итерационного процесса и называется невязкой4 . Запишем схему итерационного процесса B xk+1 − xk + Axk = f, τk+1 и выразим xk+1 : xk+1 = xk − τk+1 B −1 (Axk − f ) = xk − τk+1 B −1 rk = xk − τk+1 ω k . Число ω k = B −1 rk называется поправкой. Перепишем (1.39), используя новые определения: τk+1 = Dω k , z k . hDω k , ω k i (1.40) Проще ли нам от этого стало? Да пока не очень, потому что величина z k по-прежнему неизвестна, и вычислить ее мы не можем. . . Воспользуемся возможностью варьировать D. 4 невязка(от греч. ηεϑαςυς ) — дисбаланс. 1.8. ПРИМЕРЫ ИТЕРАЦИОННЫХ МЕТОДОВ ВАРИАЦИОННОГО ТИПА 1.8 31 Примеры итерационных методов вариационного типа Метод скорейшего спуска Пусть матрица A, задающая систему уравнений Ax = f, симметрична и положительно определена (A = AT > 0), и выберем матрицу D = A, тогда τk+1 = Aω k , z k ω k , Az k ω k , rk = = hAω k , ω k i hAω k , ω k i hAω k , ω k i — теперь, так как нам известно rk на каждом шаге, то τk+1 становится вычислимым. Получим необходимое условие на матрицу B для применимости итерационных методов вариационного типа. Используя то, что 1 1 1 1 1 1 1 1 C = D 2 B −1 AD− 2 = {D = A =⇒ D 2 = A 2 } = A 2 B −1 AA− 2 = A 2 B −1 A 2 , и положительную определенность матрицы C, получим D 1 E D E 1 1 1 0 < hCx, xi = A 2 B −1 A 2 x, x = B −1 A 2 x, A 2 x . 1 Обозначив A 2 x = y, получим B −1 y, y > 0, или, еще раз переобозначив B −1 y = u, hu, Bui > 0, то есть матрица B должна быть положительно определена — это и есть необходимое условие применимости метода вариационного типа. Возьмем B = E : rk , rk . xk+1 = xk − τk+1 rk , где τk+1 = hArk , rk i Этот набор (e τ ) соответствует методу скорейшего спуска. Пояснение названия метода. F (x) = hAx, xi − 2 hf, xi = Введем функцию F (x) : n X i=1 (Ax)i xi − 2 n X fi xi = i=1 n X n X aij xj xi − 2 i=1 j=1 n X fi xi , где x ∈ Rn . i=1 Пусть A = AT , найдем градиент функции F (x), чтобы определить направление максимального убывания: T  ∂F ∂F , ..., , grad F (x) = ∂x1 ∂xn где n X n n n X X X ∂F ∂ = aij (xj , xi ) − 2fl = ail xi + alj xj − 2fl . ∂xl ∂xl i=1 j=1 i=1 j=1 Используя симметричность матрицы A, получим: n X ∂F =2 ali xi − 2fl = 2((Ax)l − fl ) = 2rl . ∂xl i=1 Таким образом, мы показали, что градиент функции F (x) равен удвоенному вектору невязки: grad F (x) = 2r. То есть, учитывая формулу xk+1 = xk − τk+1 rk , и представление для невязки, получается, что каждый раз мы строим xk+1 от xk в направлении антиградиента функции F (x). 32 Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ При этом остается свободным параметр τk+1 . Будем его выбирать так, чтобы минимизировать значение функции F (xk+1 ) на k + 1 приближении.
min F (xk+1 ) = min Axk+1 , xk+1 − 2 f, xk+1 = τk+1 τk+1
= min A(xk − τk+1 rk ), xk − τk+1 rk − 2 f, xk − τk+1 rk . τk+1 Раскроем скобки, и, в силу свойств скалярного произведения, получим   2 min F (xk+1 ) = min Axk , xk − 2τk+1 rk , Axk + τk+1 Ark , rk − 2 f, xk + 2τk+1 f, rk . τk+1 τk+1 Чтобы найти минимум по τk+1 , возьмем производную по этой переменной от функции F (x) и приравняем ее к 0. ∂F (xk+1 ) = −2 rk , Axk + 2τk+1 Ark , rk + 2 f, rk = 0. ∂τk+1 Откуда получаем выражение для τk+1 : τk+1 = rk , rk rk , Axk − rk , f = . hArk , rk i hArk , rk i Как нетрудно заметить оно совпадает с записанным ранее τk+1 в определении метода скорейшего спуска. Таким образом, этот метод гарантирует построение итерационной последовательности, которая строится в направлении максимального уменьшения функции F (x). Из курса вариационного исчисления известно, что минимизация функции F (x) приводит к наиболее быстрой сходимости xk к точному решению системы. Рассмотрим скорость сходимости этого метода. Будем подбирать τk+1 так, чтобы минимизировать z k+1 в методе скорейшего спуска. Понятно, что при выборе τk+1 отличным от оптимального, гарантирующего минимум, будут получены величины не меньше. Таким образом, получаем оценку сверху: ||z k+1 ||A 6 ||(E − τ0 B −1 A)z k ||A . В параграфе об оценке сходимости одношаговых стационарных методов была получена верхняя оценка ||z k+1 ||A 6 q||z k ||A . Очевидно, метод скорейшего спуска имеет оценку не хуже, чем ОСИМ, и, соответственно, итерационные методы вариационного типа имеют погрешность не хуже, чем соответствующие (по матрице B ) ОСИМ с оптимальным выбором параметра τ0 . Отсюда можно сделать заключение, что вариационные методы сходятся не медленнее стационарных методов. Кроме того, для метода вариационного типа следует обратить внимание на то, что данный результат (оценку погрешности) можно получить, практически не исследуя матрицу B −1 A. Метод минимальных невязок В качестве матрицы D возьмем D = AT A. Из этого следует, что D = DT > 0. Тогда, учитывая представление (1.40), параметр τk+1 будет считаться так: τk+1 = Dω k , z k . hDω k , ω k i Перепишем τk+1 так, чтобы его можно было вычислить. Распишем D : τk+1 = AT Aω k , z k Aω k , rk Aω k , rk = = hAT Aω k , ω k i hAω k , Aω k i ||Aω k ||2 1.8. ПРИМЕРЫ ИТЕРАЦИОННЫХ МЕТОДОВ ВАРИАЦИОННОГО ТИПА 33 — теперь параметр τk+1 вычислим, так как невязка и поправка нам известны на каждом шаге. Для сходимости итерационных методов вариационного типа мы требовали положительной определенности матрицы C. Посмотрим к чему это условие приведет в методе минимальных невязок. Итак, найдем ограничения на исходную систему. D 1 E 1 0 < hCx, xi = D 2 B −1 AD− 2 x, x = D 1 E 1 1 = {обозначим D− 2 x = y} = D 2 B −1 Ay, D 2 y = DB −1 Ay, y = AT AB −1 Ay, y . Сделаем еще две замены: Ay = u и u = Bv, тогда получим hCx, xi = AB −1 u, u = hAv, Bvi = B T Av, v . То есть, надо следить за тем, чтобы матрица B T A была положительно определена. Если в качестве B возьмем единичную матрицу, то метод минимальных невязок будет применим к системам уравнений с положительно определенными матрицами. Объяснение названия метода. Параметр τk+1 мы подбираем таким образом, чтобы минимизировать норму погрешности: q q q min ||z k+1 ||D = min hDz k+1 , z k+1 i = min hAT Az k+1 , z k+1 i = min hrk+1 , rk+1 i = min ||rk+1 || τk+1 τk+1 τk+1 τk+1 τk+1 — таким образом, минимизируя погрешность, мы минимизируем невязки. Метод минимальных поправок В качестве матрицы D возьмем D = AT B −1 A, и наложим ограничения на матрицу: D = DT > 0, т.е. матрица B должна быть симметрической положительно определенной. Возьмем B = E =⇒ rk = ω k . При таком выборе матриц B и D получим следующее выражение для параметра τk+1 : AT B −1 Aω k , z k B −1 Aω k , rk Aω k , ω k Dω k , z k = = = . hDω k , ω k i hAT B −1 Aω k , ω k i hB −1 Aω k , Aω k i hAω k , Aω k i Итак, формула для итерационного параметра становится реализуемой. Исходя из ограничений наложенных на матрицы C выше (C > 0 ) опишем класс систем, которые можно решать с помощью метода минимальных поправок. Раскроем условие положительной определенности: D 1 E 1 0 < hCx, xi = D 2 B −1 AD− 2 x, x . τk+1 = 1 Как и в методе минимальных невязок, обозначим D− 2 x = y, Ay = u, B −1 u = v, тогда получим: E D 1 1 0 < D 2 B −1 Ay, D 2 y = DB −1 Ay, y = AT B −1 AB −1 Ay, y = = B −1 AB −1 Ay, Ay = B −1 AB −1 u, u = B −1 Av, Bv = hAv, vi — то есть, метод минимальных поправок применим для систем с положительно определенной матрицей. Объяснение названия метода. Как и в методе минимальных невязок, будем минимизировать норму погрешности ||z k+1 ||2D = Dz k+1 , z k+1 = AT B −1 Az k+1 , z k+1 = = B −1 Az k+1 , Az k+1 = B −1 rk+1 , rk+1 = ω k+1 , Bω k+1 = ||ω k+1 ||2B . Таким образом, минимизируя норму погрешности, мы минимизируем норму поправки. 34 Глава 1. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ Метод минимальных погрешностей В этом примере матрица D берется равной некоторой матрице B0 : B0 = B0T > 0, которая должна быть к тому же легко обратимой. Матрица B определяется так: B = (AT )−1 B0 . Как и ранее, покажем корректность формулы для τk+1 : τk+1 = Dω k , z k B0 B −1 rk , z k = = hDω k , ω k i hB0 B −1 rk , ω k i T = {B−1 = B−1 0 A } = B0 B0−1 AT rk , z k B0 B0−1 AT rk , ω k = rk , rk rk , Az k = . hrk , Aω k i hrk , Aω k i — как уже показывалось, rk (невязку) и wk (поправку) мы умеем вычислять на каждом шаге. Другим ограничением была положительная определенность матрицы C. Проверим это: E D 1 E D 1 1 1 1 hCx, xi = D 2 B −1 AD− 2 x, x = {обозначим y = D− 2 x} = D 2 B −1 Ay, D 2 y = = DB −1 Ay, y = {D = B0 , B = (AT )−1 B0 } = B0 B0−1 AT Ay, y = hAy, Ayi . Очевидно, hAx, Axi всегда больше нуля, если матрица A — невырожденная. Таким образом, метод минимальных погрешностей верен для любых невырожденных матриц A. Название метода произошло от требования минимизации погрешности ||z k+1 ||D = ||z k+1 ||B0 на (k + 1)-м шаге. Замечание. Все эти методы возникли из стратегии локальной минимизации погрешности от шага к шагу. При таких требованиях достигается высокая скорость сходимости, причем нам не требуется дополнительная информация, к примеру, о спектре матрицы A. Однако нет предела совершенству... 1.9 Двухшаговые итерационные методы вариационного типа Данные методы используются при решении систем большой размерности (тысячи, десятки тысяч уравнений). Решаем мы по-прежнему систему Ax = f. Действуя так же, как и при работе с одношаговыми итерационными методами, можно показать, что каноническая форма двухшаговых методов такова: B xk+1 − xk + (1 − αk+1 )(xk − xk−1 ) + Axk = f, αk+1 τk+1 k = 1, 2, . . . (1.41) Для того, чтобы задать итерационный процесс, мы определяем два начальных приближения: x0 и x . Остальные вычисляются по этой формуле, где B, αk+1 , τk+1 — параметры, определяющие метод. Также неявным параметром будет матрица D. Она определяет норму ||z k+1 ||D , которую мы, как и раньше, будем минимизировать на каждой итерации. Для этого необходимо брать αk+1 и τk+1 такими:  Dω k , z k   , k = 1, 2, . . . ; τ =  k+1   hDω k , ω k i   α1 = 1; !−1   k k  Dω , z  τ k+1   1− · , k = 1, 2, . . .  αk+1 = τk αk hDω k−1 , z k−1 i 1 — подробный вывод этих формул можно найти в [2]. Как и раньше, при надлежащем выборе мы можем избавиться от z k в скалярном произведении и корректно применять эти формулы при подсчете приближений. Поясним это на примерах. 1.9. ДВУХШАГОВЫЕ ИТЕРАЦИОННЫЕ МЕТОДЫ ВАРИАЦИОННОГО ТИПА 35 Метод сопряженных градиентов В данном методе берется D = A, поэтому матрица A должна быть симметричной и положительно определенной. Как уже показывалось ранее, в этом случае для τk+1 можно получить такую формулу: τk+1 = ω k , rk . hAω k , ω k i В свою очередь, для αk+1 будет справедливо такое выражение: !−1 ωk , rk τk+1 αk+1 = 1 − · τk αk hω k−1 , rk−1 i — в состав этих формул входит уже не погрешность, а невязка и поправка, которые мы вычислять умеем. Метод сопряженных невязок В данном случае берется D = AT A, и несложно проверить, что D = DT > 0. Получающиеся формулы для τk+1 и αk+1 мы опускаем (вывод предоставляется читателю). Кроме того, в данном методе и во всех последующих вводится дополнительное требование: DB −1 A = (DB −1 A)T > 0. (1.42) В данном случае оно эквивалентно тому, что B T A > 0. Это небольшое ограничение на матрицу B. Метод сопряженных поправок Здесь мы берем D = AT B −1 A. Легко проверить, что D = DT > 0, если B = B T > 0. Ограничение (1.42) приводит к требованию положительной определенности матрицы A. Формулы для τk+1 и αk+1 также становятся корректными. Метод сопряженных погрешностей Здесь все тоже просто: D = B0 , где B0 = B0T > 0, причем матрица B0 должна быть легко обратимой. Матрица B задается строго: B = (AT )−1 B0 . Ограничение (1.42) приводит к требованию невырожденности матрицы A. Общие замечания (1) Можно показать, что любой из выше перечисленных методов сходится не медленнее, чем соответствующий ему одношаговый итерационный метод с чебышевским набором параметров. (2) Кроме того, если количество шагов итерации превысит размерность системы, то последовательность итерационных приближений выйдет на точное решение — то есть вышеперечисленные методы фактически являются прямыми. Другое дело, что нужная точность будет достигнута намного раньше. (3) Учитывая, что α1 = 1, можно не угадывать x1 (первое приближение), а просто подставить в (1.41) k = 0 : x1 − x0 B + Ax0 = f α1 τ1 — это простая одношаговая схема для поиска x1 . Глава 2 Задачи на собственные значения Сначала о терминологии. Пусть A = (aij ) — матрица размера n × n, а x = (x1 , x2 , . . . , xn )T — вектор неизвестных. Тогда поиск таких констант λ и векторов x, что Ax = λx, называется задачей на собственные значения. Нетрудно заметить, что эта задача эквивалентна поиску таких x и λ, для которых (A − λE)x = 0. При этом λ называются собственными значениями, а соответствующие им вектора x — собственными векторами. Из курса линейной алгебры известно, что если det(A − λE) = 0, то решение существует. Для его нахождения надо найти такие λ, чтобы определитель |A − λE| обращался в ноль. Этот определитель является полиномом от λ с коэффициентами из A — поэтому, если n мало, то корни этого полинома легко найти. Это также просто сделать, если матрица A является диагональной или верхнетреугольной: в этом случае определитель будет равен произведению диагональных элементов. Далее мы будем рассматривать различные методы нахождения собственных значений. Определение. Метод решает полную проблему собственных значений, если он находит их все. В противном случае (например, надо найти только границы спектра матрицы A) говорят, что решается частичная проблема. 2.1 Поиск собственных значений методом вращений Данный метод, предложенный К. Якоби, позволяет найти все собственные значения (решает полную проблему поиска). Для этого требуется симметричность матрицы A : A = AT . Известно, что в этом случае для A справедливо такое представление (оно единственно): A = QT Λ Q, (2.1) где Q — некая ортогональная матрица (QT = Q−1 ), а Λ — диагональная матрица, причем на диагонали у нее стоят именно собственные значения матрицы A : Λ = diag (λ1 , . . . , λn ) . Если домножить равенство (2.1) на (QT )−1 = Q слева, а на Q−1 = QT — справа, то получим такую формулу: Λ = QAQT (2.2) — то есть для нахождения собственных значений матрицы A нам необходимо найти соответствующую матрицу Q и провести два матричных умножения. 36 2.1. ПОИСК СОБСТВЕННЫХ ЗНАЧЕНИЙ МЕТОДОМ ВРАЩЕНИЙ 37 Матрицу Q будем находить с помощью ортогональных преобразований матрицы A, постепенно уменьшая абсолютные значения ее внедиагональных элементов: A1 = Vij1 A (Vij1 )T ; A2 = Vij2 A1 (Vij2 )T ; ··· Ak = Vijk Ak−1 (Vijk )T , и так далее. Здесь Vijk — некие ортогональные матрицы, а индексы в них говорят о номере преобразования — (k ) и об индексе уменьшаемого элемента из Ak−1 — ( ij ). Произведение ортогональных матриц дает ортогональную матрицу, поэтому, если мы на некотором шаге придем к диагональной матрице, то это будет означать получение преобразования (2.2). Матрицы Vijk будут задаваться так (ϕ — пока свободный параметр):   1 .  ..    1   cos ϕ − sin ϕ     1   k   . Vij =  . . 0 0  0 1   sin ϕ cos ϕ     1  ..  . 1 — здесь на диагонали стоят единицы везде, кроме vii и vjj , где стоят косинусы, а в позициях vji и vij стоят синусы. Все остальные элементы — нулевые. Легко проверить, что все эти матрицы ортогональны. Индексы i и j задаются на каждом шаге заново. Они обозначают индекс максимального по модулю внедиагонального элемента, то есть |aij | = max |aklm |, l,m l6=m где aklm — элементы матрицы Ak (k-е приближение к Λ ). Если i и j можно выбрать несколькими способами, то берется произвольная пара. Если же все внедиагональные элементы — нулевые, то процесс прекращается. Назовем количественной характеристикой диагональности матрицы Ak такое число: t(Ak ) = n X n X (akij )2 i=1 j=1 j6=i k→∞ Понятно, что если t(Ak ) −→ 0, то Ak сходится к диагональному виду. Таким образом, через t(Ak ) можно наглядно оценивать скорость сходимости данного метода. Теперь фиксируем максимальный элемент aij и установим соотношение между матрицами Ak+1 и Ak . Ak+1 мы задавали так: Ak+1 = Vijk+1 Ak (Vijk+1 )T . Обозначим Bk = Ak (Vijk+1 )T = (bkms ) и (Vijk+1 )T = (v k+1 lm ) — обозначение для элементов матрицы k+1 T (Vij ) . Тогда, согласно определению матричного произведения,  k a , s 6= i, j;   n  ms X akmi cos ϕ − akmj sin ϕ, s = i; bkms = akmp v k+1 (2.3) ps =   p=1  k k ami sin ϕ + amj cos ϕ, s = j. 38 Глава 2. ЗАДАЧИ НА СОБСТВЕННЫЕ ЗНАЧЕНИЯ один ненулевой при p 6= i, j, и — эта сложная формула вытекает из того, что среди элементов v k+1 ps два ненулевых при p = i, j. Теперь, согласно принятым обозначениям, Ak+1 = Vijk+1 Bk . Выведем отсюда формулу для элементов матрицы Ak+1 :  k m 6= i, j; b ,   n  ms X k+1 k bkis cos ϕ − bkjs sin ϕ, m = i; (2.4) ak+1 vmp bps = ms =   p=1  k bis sin ϕ + bkjs cos ϕ, m = j. k+1 — элементы матрицы Vijk+1 ). — идея построения системы та же (vmp Напомним, что aij — максимальный по модулю элемент из A. Положив в (2.4) m = i и s = j, получим такую формулу для него: k+1 aij = bkij cos ϕ − bkjj sin ϕ. Теперь распишем bkij и bkjj через (2.3), взяв m = i, j и s = j : ak+1 ij = (akii sin ϕ + akij cos ϕ) cos ϕ − (akji sin ϕ + akjj cos ϕ) sin ϕ = = (akii − akjj ) sin ϕ cos ϕ + akij (cos2 ϕ − sin2 ϕ) = {A = AT ⇒ Ak = AT k} = (akii − akjj ) sin 2ϕ + akij cos 2ϕ. 2 Напомним, мы пытаемся уменьшить внедиагональные элементы матрицы A при ортогональном преобразовании. Потребуем равенство нулю для ak+1 ij . Таким образом, написанное выше преобразование превращается в уравнение относительно ϕ. Решая его, получим: ϕ= 2akij 1 arctg k . 2 aii − akjj Вычислим количественную характеристику диагональности получившейся матрицы Ak+1 : t(Ak+1 ) = n X n X 2 (ak+1 ms ) . m=1 s=1 s6=m Согласно формулам, элементы в Ak при умножении на (Vijk+1 )T изменяются только в i-м и j-м столбцах. Аналогично, элементы в Ak+1 изменяются относительно B k только в i-й и j-й строках. То есть Ak+1 отлична от Ak только в i-х и j -х строках и столбцах. Выделим в сумме неменяющиеся элементы. Распишем сумму, раскроем скобки, и проведем все преобразования: t(Ak+1 ) = = n X n X m=1 s=1 m6=i,j s6=i,j,m n n X X m=1 s=1 m6=i,j s6=i,j,m (akms )2 + (akms )2 + n n X X   k+1 2   k 2 2 (ais ) + (ak+1 = (bmi ) + (bkmj )2 + js ) m=1 m6=i,j n X s=1 s6=i,j  k 2 (ami ) cos2 ϕ + (akmj )2 sin2 ϕ − 2akmi akmj sin ϕ cos ϕ+ m=1 m6=i,j n X   k 2 +(akmi )2 sin2 ϕ + (akmj )2 cos2 ϕ + 2akmi akmj sin ϕ cos ϕ + (bis ) cos2 ϕ + (bkjs )2 sin2 ϕ− s=1 s6=i,j  −2bkis bkjs sin ϕ cos ϕ + (bkis )2 sin2 ϕ + (bkjs )2 cos2 ϕ + 2bkis bkjs sin ϕ cos ϕ = 2.2. СТЕПЕННОЙ МЕТОД ПОИСКА СОБСТВЕННЫХ ЗНАЧЕНИЙ = n X n X n n X X  k 2   k 2  (ami ) + (akmj )2 + (ais ) + (akjs )2 + 2(akij )2 − 2(akij )2 = (akms )2 + s=1 m=1 m6=i,j s6=i,j,m 39 m=1 m6=i,j s=1 s6=i,j = t(Ak ) − 2(akij )2 . Из всей этой последовательности формул следует, что t(Ak+1 ) < t(Ak ) — как видно, характеристика диагональности монотонно убывает с ростом индекса, причем уменьшается каждый раз на 2(akij )2 — удвоенный квадрат максимального элемента. Из этого следует, что описанный итерационный процесс приводит к диагональной матрице. Запишем оценку на скорость сходимости процесса. Пусть akij — максимальный внедиагональный элемент. Простым подсчетом элементов матрицы можно получить такое неравенство: t(Ak ) 6 n(n − 1)(akij )2 . t(Ak ) для n > 2. Подставим это неравенство в ранее полученное n(n − 1) соотношение между t(Ak ) и t(Ak+1 ) : Отсюда следует, что (akij )2 > t(Ak+1 ) = t(Ak ) − 2(akij )2 6 t(Ak ) − 2 t(Ak ) = qt(Ak ), n(n − 1) где q = 1 − 2 < 1. n(n − 1) Применив эту операцию k раз, получим: t(Ak ) 6 q k t(A). Из последнего неравенства видно, что процесс нахождения матрицы Q сходится к диагональной матрице Λ со скоростью геометрической прогрессии со знаменателем q. Можно немного оптимизировать способ выбора «плохого» элемента: сначала выбираем «плохую» n X строку (например, выбираем строку s, где (asi )2 имеет максимальное значение), а потом из этой i=1 строчки выбираем максимальный по модулю элемент. При этом, немного неоптимальный выбор ask компенсируется скоростью нахождения «плохого» элемента (2n вместо n2 сравнений). Теперь будем решать частичные проблемы. Как пример частичной проблемы, можно привести задачу нахождения границ спектра в итерационном методе. 2.2 Степенной метод поиска собственных значений Будем рассматривать задачу нахождения максимального по модулю собственного значения симметрической матрицы A. Алгоритм поиска Возьмем любой вектор x0 , отличный от нуля, и построим последовательность векторов xk такую, что xk+1 = Axk , k = 0, 1, . . . (2.5) По этой последовательности построим числовую последовательность {Λk1 } по следующему правилу: Λk1 = xk+1 , xk . hxk , xk i Пусть собственные числа занумерованы так, что первым стоит максимальное по модулю — искомое: |λ1 | > |λ2 | > |λ3 | > . . . > |λn |. 40 Глава 2. ЗАДАЧИ НА СОБСТВЕННЫЕ ЗНАЧЕНИЯ Примечание. Обратите внимание, что максимальное по модулю собственное значение единственно — мы будем этим пользоваться. Покажем, что в этом случае последовательность {Λk1 } будет сходиться к λ1 . Так как матрица A симметрическая, то для нее существует ортонормированный базис {ξi } из собственных векторов. Разложим начальное приближение x0 по этому базису: x0 = n X αi ξi . (2.6) i=1 Из (2.5) получаем выражение для xk через x0 : xk = Axk−1 = . . . = Ak x0 . Используя разложение x0 по базису, запишем следующую оценку на xk : xk = Ak n X αi ξi = Ak−1 i=1 n X αi Aξi = Ak−1 i=1 n X αi λi ξi = . . . = i=1 В силу того, что |λ2 | > |λi |, где i = 3, n, получим, что n X αi λki ξi = α1 λk1 ξ1 + n X i=2 i=1 n X αi λki ξi αi λki ξi = O(|λ2 |k ). Подсчитаем два i=2 скалярных произведения для нахождения Λk1 . xk , xk = α1 λk1 ξ1 + O(|λ2 |k ), α1 λk1 ξ1 + O(|λ2 |k ) k k = α12 λ2k 1 + O(|λ1 | · |λ2 | ); xk+1 , xk = α1 λk+1 ξ1 + O(|λ2 |k+1 ), α1 λk1 ξ1 + O(|λ2 |k ) 1 + O(|λ1 |k+1 · |λ2 |k ). = α12 λ2k+1 1 Зная скалярные произведения, вычислим члены последовательности Λk1 :    k λ2 2 2k+1 α1 λ1 1 + O λ1 α2 λ2k+1 + O(|λ1 |k+1 · |λ2 |k )    = λ1 1 + O = Λk1 = 1 21 2k k α1 λ1 + O(|λ1 |k · |λ2 |k ) λ2 2k 2 α1 λ1 1 + O λ1 λ2 λ1 k !! −→ λ1 при k → ∞. То есть видно, что последовательность Λk1 сходится к λ1 — искомому максимальному собственному значению. Определим, к чему же сходится последовательность xk . Для этого рассмотрим ! k xk α1 λk1 ξ1 + O(|λ2 |k ) α1 λk1 ξ1 + O(|λ2 |k ) λ2    = ±ξ1 + O p = = k ||xk || λ1 hxk , xk i λ2 k |α1 ||λ1 | 1 + O λ1 Откуда видно, что xk сходится к собственному вектору ξ1 по направлению (чем больше k, тем k ближе ||xxk || по направлению к ξ1 ). Примечание. В разложении (2.6) в общем случае коэффициент α1 может быть равен нулю, и итерационный процесс может не сходиться. Для устранения неопределенности обычно прогоняют алгоритм для нескольких случайно выбранных начальных приближений. Примеры 2.3. МЕТОД ОБРАТНОЙ ИТЕРАЦИИ 41 Поиск максимального собственного значения. В общем случае максимальное собственное значение может не быть максимальным по модулю. В этом случае переходят от матрицы A к матрице D = A + cE, где c — некоторая положительная константа. Легко заметить, что λ(D) = λ(A) + c и λmax (D) = λmax (A) + c. Очевидно, что можно подобрать такое c, что все собственные значения D будут положительны. Теперь мы пользуемся вышеприведенным алгоритмом, находим λmax (D), а потом и λmax (A). Для поиска минимального собственного значения можно брать константу c < 0 и проводить аналогичные рассуждения. Поиск собственного значения, близкого к заданному числу Λ. λ, что |λ − Λ| = min |λi − Λ|. Это задача о поиске такого i Аналогично предыдущему примеру, переходим от матрицы A к матрице D = E −c(A−ΛE)2 , c > 0. В этом случае λ(D) = 1 − c(λ(A) − Λ)2 , и максимальное собственное значение D будет соответствовать собственному значению A, ближайшему к Λ. 2.3 Метод обратной итерации Пусть нам дана матрица A. Обозначим соответствующие ее собственному значению λi собственные вектора как xi : Axi = λi xi . Тогда (A − λi E)xi = 0. (2.7) Задача состоит в том, чтобы найти один из собственных векторов xi . Сначала для поиска собственного значения λi необходимо решать уравнение det(A − λi E) = 0. Если мы нашли λi точно, то вычисление xi не составит труда. Однако, если мы получили неточное собственное значение λi ≈ λi , то определитель det(A − λi E) будет отличен от нуля. Так как должна выполняться формула (2.7), то единственным подходящим xi будет xi = 0. Покажем, что даже при неточно вычисленном собственном значении можно вычислить собственный вектор. Фиксируем произвольный вектор b и решим систему: (A − λi E)x = b. (2.8) Утверждается, что решение этой системы будет приближенно равняться искомому собственному вектору: x ≈ xi . Для решения исходной системы построим последовательность векторов по следующему правилу: (2.9) (A − λi E)xk+1 = xk , k = 0, 1, . . . , где x0 = b. При таком задании итерационного процесса достаточно примерно 5 итераций, чтобы с хорошей точностью определить собственный вектор xi , соответствующий собственному значению λi . Покажем, что алгоритм работает правильно. Пусть матрицаX A такова, чтоX существует базис {ξi } из собственных векторов этой матрицы (xi = C · ξi ). Пусть b = βj ξj и x = αj ξj — разложения j j фиксированного вектора b и искомого вектора x по этому базису. Подставим эти разложения в (2.8) X X αj ξj = βj ξj , (A − λi E) X j j j (αj λj − αj λi )ξj = X j βj ξj , 42 Глава 2. ЗАДАЧИ НА СОБСТВЕННЫЕ ЗНАЧЕНИЯ X [αj (λj − λi ) − βj ]ξj = 0. j Только коэффициенты, равные нулю, могут обратить линейную комбинацию базисных векторов в нуль. Следовательно, получаем такое выражение для αj : αj = βj . λj − λ i Отсюда видно, что в разложении x по базису из собственных векторов коэффициент αi при базисном векторе ξi будет велик по сравнению с другими (если λi близко к λi ). Это означает, что каждый шаг итерационного процесса (2.9) приводит нас к вектору, который все больше похож на искомый собственный вектор ξi . Можно говорить, что итерационный процесс сходится к ξi по направлению. Глава 3 Численные методы решения нелинейных уравнений Пусть f (x) — некоторая непрерывная функция, заданная на отрезке [a; b]. Нашей задачей будет поиск корней уравнения f (x) = 0 на этом отрезке. Обычно это проходит в два этапа: сначала проводится локализация корней: выделение небольших отрезков, содержащих только один корень, а потом на этих отрезках проводится уточнение его значения. 3.1 Методы разделения корней Как легко можно убедиться, локализация корней даже для функции одной переменной представляет собой достаточно трудоемкую задачу. Одним из способов является разбиение отрезка [a; b] произвольным образом на N подотрезков [xk ; xk+1 ] : x0 = a < x1 < x2 < . . . < xN = b. Теперь мы считаем значение функции на границах этих маленьких отрезков: f (xk ) = fk . Обратим свое внимание на те из них, для которых выполняется неравенство: fk · fk+1 < 0. (3.1) Оно означает, что значения функции на краях отрезка различны по знаку, а из непрерывности f (x) следует, что она где-то внутри обращается в ноль. Далее повторяем вышеописанную процедуру для всех отрезков, для которых верно (3.1). Упрощенным вариантом предыдущего метода является метод бисекции. Пусть f (a) · f (b) < 0 — — середина отрезка это означает, что внутри [a; b] точно есть корень. Теперь обозначим x0 = a+b 2 [a; b]. Если он не является корнем, то либо f (a) · f (x0 ) < 0, либо f (x0 ) · f (b) < 0. Соответственно, делим пополам [a; x0 ] или [x0 ; b] и так до достижения нужной точности — очевидно, это всегда можно сделать. 3.2 Примеры численных методов Метод простой итерации ∗ Пусть x — корень уравнения f (x) = 0, лежащий на отрезке [a; b]. Зададим τ (x) — некоторую функцию-параметр, не обращающуюся в ноль на [a; b]. Тогда, очевидно, f (x) = 0 ⇐⇒ −τ (x)f (x) = 0 ⇐⇒ x − τ (x)f (x) = x. 43 44 Глава 3. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ УРАВНЕНИЙ Обозначим S(x) = x − τ (x)f (x), тогда получим, что f (x) = 0 ⇐⇒ S(x) = x. (3.2) Теперь будем строить последовательность приближений к x∗ по следующему правилу: xk+1 = S(xk ), k = 0, 1, . . . , задаваясь при этом некоторым начальным приближением x0 . Если предел последовательности {xk } существует: lim xk = x∗ , то x∗ = S(x∗ ), и, согласно (3.2), k→∞ он будет являться корнем исходного уравнения. Очевидно, эта последовательность не всегда будет сходиться, а зависит это от функции S(x), которая в свою очередь определяется параметром τ (x). Из геометрических соображений (нарисовав соответствующие графики) можно сделать предположение, что условие |S 0 (x)| < 1 на некотором отрезке [c; d] внутри [a; b] будет достаточным для существования предела {xk }, если начальное приближение тоже взять из [c; d]. Докажем это позже, а пока рассмотрим один из вариантов данного метода — метод Ньютона. Метод Ньютона1 Итак, пусть x∗ — нуль функции f (x) : f (x∗ ) = 0, (3.3) а f 0 (x) существует, непрерывна и отлична от нуля на [a; b]. Это означает, в частности, что других нулей у f (x) на этом отрезке нет. Теперь перепишем (3.3) следующим образом: f (xk + (x∗ − xk )) = 0. Применим теперь к этому выражению формулу Лагранжа: f (xk ) + f 0 (x)(x∗ − xk ) = 0, x ∈ [a; b]. Для получения формулы итерационного процесса заменим в этом равенстве x на xk , а x∗ — на x . Равенство превратится в уравнение: k+1 f (xk ) + f 0 (xk )(xk+1 − xk ) = 0 Выразим отсюда xk+1 : xk+1 = xk − f (xk ) . f 0 (xk ) (3.4) Проведенное преобразование называется линеаризацией уравнения (3.3). Метод Ньютона называется также методом касательных, так как новое приближение является абсциссой точки пересечения касательной к графику функции f (x), проведенной в точке M (xk , f (xk )), с осью OX. Замечание. Метод Ньютона имеет (когда сходится) квадратичную скорость сходимости:
xk+1 − x∗ = O (xk − x∗ )2 . Это хорошее свойство, однако метод ведь может и не сходиться. Возвращаясь к ранее введенным обозначениям, получим, что S(x) в методе Ньютона считается так: S(x) = x − 1 Разработан И. Ньютоном (1669). f (x) , f 0 (x) 3.2. ПРИМЕРЫ ЧИСЛЕННЫХ МЕТОДОВ то есть τ (x) = 1 f 0 (x) 45 . Мы предположили, что метод простой итерации (т.е. и метод Ньютона) является сходящимся, если |S 0 (x)| < 1. Рассмотрим это неравенство поподробнее: S 0 (x) = 1 − (f 0 (x))2 − f (x)f 00 (x) =⇒ |S 0 (x)| < 1 ⇐⇒ (f 0 (x))2 ⇐⇒ f (x)f 00 (x) < 1. (f 0 (x))2 (3.5) Таким образом, метод Ньютона будет сходится, если неравенство (3.5) будет выполняться на некотором отрезке, содержащем начальное приближение x0 и нужный нам корень x∗ . Часто такое требование приводит к необходимости брать x0 очень близко к x∗ , что не всегда выполнимо. Модифицированный метод Ньютона Если у нас есть проблемы с подсчетом производной на каждом шаге, то можно воспользоваться 1 модифицированным методом Ньютона, где берется τ (x) = f 0 (x 0 ) , то есть xk+1 = xk − f (xk ) , f 0 (x0 ) k = 0, 1, . . . , где x0 — начальное приближение. Однако, в этом методе, как мы можем подсчитать, уже не квадратичная, а линейная скорость сходимости: |xk+1 − x∗ | = O(xk − x∗ ). Замечание. Модифицированный метод Ньютона можно назвать методом одной касательной, так как новые приближения являются абсциссами точек, получающихся при пересечении с осью OX прямых, параллельных касательной к графику f (x) в точке M (x0 , f (x0 )). Метод секущих Когда нет возможности считать производную, просто заменим ее в формуле (3.4) на разностное приближение: f (xk ) − f (xk−1 ) . f 0 (xk ) ≈ xk − xk−1 Тогда получим такую формулу: xk+1 = xk − xk − xk−1 f (xk ), f (xk ) − f (xk−1 ) k = 1, 2, . . . (3.6) — здесь уже требуется задать два начальных приближения: x0 и x1 . Скорость сходимости будет линейной. Этот метод называется методом секущих. Для пояснения названия заметим, что уравнение для секущей, проходящей через точки M 0 (xk−1 , f (xk−1 )) и M 00 (xk , f (xk )) (точки предыдущих приближений), будет выглядеть так: f (xk ) − f (xk−1 ) y − f (xk ) = . x − xk xk − xk−1 Положив y = 0 и x = xk+1 , можно получить формулу (3.6). Это означает, что xk+1 — это абсцисса точки пересечения нашей секущей с осью OX : 46 Глава 3. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ УРАВНЕНИЙ y M 00  M0   xk+1 3.3
xk xk−1 x Сходимость метода простой итерации Как было сказано ранее, в методе простой итерации уравнение для нового приближения xk+1 таково: xk+1 = S(xk ), k = 0, 1, . . . , (3.7) а x0 — заданное начальное приближение. Кроме того, введем такое обозначение для r -окрестности точки a : Ur (a) = {x : |x − a| 6 r}. Теперь мы готовы к формулировке теоремы о сходимости этого метода. Теорема 3.1. Пусть для некоторых r и a функция S(x) удовлетворяет на множестве Ur (a) условию Липшица с константой q ∈ (0; 1) : |S(x0 ) − S(x00 )| 6 q|x0 − x00 | ∀ x0 , x00 ∈ Ur (a), причем |S(a) − a| 6 (1 − q)r. Тогда уравнение x = S(x) имеет на множестве Ur (a) единственное решение, а метод простой итерации (3.7) сходится к этому решению при любом начальном приближении из Ur (a), причем для погрешности на k-й итерации справедлива оценка: |xk − x∗ | 6 q k |x0 − x∗ | — то есть скорость сходимости линейная. Доказательство. Возьмем начальное приближение x0 из множества Ur (a). Индукцией покажем, что остальные xj тоже принадлежат множеству Ur (a). Пусть xj ∈ Ur (a), докажем, что следующий член последовательности xj+1 ∈ Ur (a). |xj+1 − a| = |S(xj ) − a| = |S(xj ) − S(a) + S(a) − a| 6 |S(xj ) − S(a)| + |S(a) − a| из того, что функция S(x) Липшиц-непрерывна и условия теоремы получаем: |xj+1 − a| 6 q|xj − a| + (1 − q)r 6 qr + (1 − q)r = r. 3.3. СХОДИМОСТЬ МЕТОДА ПРОСТОЙ ИТЕРАЦИИ 47 Теперь покажем, что последовательность {xk } имеет предел, являющийся решением уравнения x = S(xk ), причем это решение единственно. Сначала установим сходимость, для этого оценим разность двух соседних элементов: k+1 |xj+1 − xj | = |S(xj ) − S(xj−1 )| 6 q|xj − xj−1 | 6 . . . 6 q j |x1 − x0 |. Покажем, что выполняется критерий Коши сходимости числовой последовательности: k+p |x k −x |= p X k+j x k+j−1 −x
j=1 6 p X q k+j−1 6 p X xk+j − xk+j−1 6 j=1 1 k 1 |x − x | = q |x − x | · p X q j−1 k 1 < q |x − x | · j=1 j=1 ∞ X q j−1 = j=1 qk |x1 − x0 |. 1−q Последнее выражение не зависит от p, и его можно сделать меньше любого ε > 0, и мы можем вычислить, с какого k(ε) будет выполнена эта оценка. Таким образом, числовая последовательность {xk }∞ k=0 сходится при k → ∞ к некоторому ∗ x ∈ Ur (a). Покажем, что этот предел является корнем уравнения x = S(x). Запишем итерационную форму нашего уравнения xk+1 = S(xk ), и перейдем к пределу при k → ∞. Левая часть xk+1 , как было уже показано, сходится к x∗ ∈ Ur (a), а правая часть S(xk ) в силу Липшиц-непрерывности S(x) сходится к S(x∗ ). Таким образом, решение существует. Покажем единственность найденного корня. Пусть существуют два решения: x∗ и x∗ . Рассмотрим модуль их разности: |x∗ − x∗ | = |S(x∗ ) − S(x∗ )| 6 q|x∗ − x∗ |, откуда вытекает, что x∗ = x∗ , иначе возникает противоречие, так как q < 1. В завершение докажем оценку на погрешность |xk − x∗ | = |S(xk−1 ) − S(x∗ )| 6 q|xk−1 − x∗ | 6 · · · 6 q k |x0 − x∗ |. Теперь теорема полностью доказана. Замечание 1. В теореме получена оценка ∀p |xk+p − xk | 6 qk |x1 − x0 |, 1−q qk |S(x0 ) − x0 |. 1−q Потребуем, чтобы xk отличалось от x∗ не более, чем на ε. Так как то есть, переходя к пределу при p → ∞, |x∗ − xk | 6 qk |S(x0 ) − x0 | 6 ε =⇒ |x∗ − xk | 6 ε, 1−q то число итераций, необходимых для достижения такой точности, можно подсчитать так:   (1−q)ε ln |S(x 0 )−x0 | . k(ε) =    ln q   Замечание 2. В условиях теоремы сделано предположение, что S(x) Липшиц-непрерывна: |S(x0 ) − S(x00 )| 6 q|x0 − x00 | ∀ x0 , x00 ∈ Ur (a). Это достаточно слабое ограничение, но его сложно проверять. Тем не менее, если функция дифференцируема, а ее производная ограничена той самой константой q, то условие Липшица будет выполнено: |S(x0 ) − S(x00 )| = {применяя формулу Лагранжа} = |S 0 (ξ)| · |x0 − x00 | 6 q|x0 − x00 | — это и есть условие Липшиц-непрерывности. 48 Глава 3. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ УРАВНЕНИЙ 3.4 Метод Эйткена Пусть метод имеет линейную скорость сходимости, то есть для погрешности выполнена следующая оценка: |xk − x∗ | 6 q k |x0 − x∗ | =⇒ xk − x∗ ≈ aq k (3.8) — мы выразили погрешность через некоторую константу a и параметр q. Допустим, что данное соотношение выполняется точно. Рассмотрим погрешность на трех соседних итерациях. xk − x∗ = aq k ; k+1 x − x∗ = aq k+1 ; k+2 x − x∗ = aq k+2 . Из этих условий легко найти x∗ . Для этого вычтем из второго уравнения первое, а из третьего второе, тогда получим: xk+1 − xk = aq k (q − 1); k+2 x − xk+1 = aq k+1 (q − 1). Вычтем одно уравнение из другого: xk+2 − 2xk+1 + xk = aq k (q − 1)2 . Свяжем это уравнение с предыдущим следующим соотношением: (xk+2 − xk+1 )2 a2 q 2k+2 (q − 1)2 = = aq k+2 = xk+2 − x∗ . k+1 k − 2x +x aq k (q − 1)2 xk+2 Если бы равенство (3.8) выполнялось точно, то можно было бы получить x∗ через три последних приближения: (xk+2 − xk+1 )2 . x∗ = xk+2 − k+2 x − 2xk+1 + xk Тем не менее, выражение, стоящее в правой части, приближает x∗ намного лучше, чем xk+2 . Обозначим его (xk+2 − xk+1 )2 xk+2 = xk+2 − k+2 , x − 2xk+1 + xk и будем считать это равенство алгоритмом построения подправленной итерационной последовательности с элементами xk . Если эту операцию (вычисление подправленных значений) проводить достаточно часто, то новая последовательность из xk сходится к точному решению значительно быстрее, чем исходный метод. Построение подправленных значений последовательности называется методом Эйткена. 3.5 Сходимость метода Ньютона Запишем итерационный процесс: xk+1 = xk − f (xk ) . f 0 (xk ) Он является модификацией метода простой итерации, тогда условием сходимости этого итерационного процесса будет неравенство |S 0 (x)| 6 q < 1, где S(x) = x − (f 0 (x))2 − f (x)f 00 (x) f (x) . S (x) = 1 − , поэтому f 0 (x) (f 0 (x))2 |S 0 (x)| 6 q < 1 ⇐⇒ f (x)f 00 (x) 6 q < 1. (f 0 (x))2 3.5. СХОДИМОСТЬ МЕТОДА НЬЮТОНА 49 Но в силу того, что мы ищем корень уравнения f (x) = 0, существует такая окрестность, где S 0 (x) 6 q < 1, но в общем случае эта область будет мала, то есть нужно подбирать начальное приближение достаточно близко расположенным к корню. Теорема 3.2 (о сходимости метода Ньютона). Пусть x∗ — простой вещественный корень уравнения f (x) = 0, а функция f (x) — дважды дифференцируема в некоторой окрестности Ur (x∗ ), причем первая производная нигде не обращается в нуль. Тогда, следуя обозначениям 0 < m1 = inf x∈Ur (x∗ ) |f 0 (x)|, M2 = |f 00 (x)|, sup x∈Ur (x∗ ) при выборе начального приближения x0 из той же окрестности Ur (x∗ ) такого, что M2 |x0 − x∗ | = q < 1, 2m1 итерационная последовательность xk+1 = xk − f (xk ) , k = 0, 1, . . . f 0 (xk ) будет сходиться к x∗ , причем для погрешности на k-м шаге будет справедлива оценка: |xk − x∗ | 6 q 2 k −1 |x0 − x∗ |. (3.9) Доказательство. Основным этапом нашего доказательства будет получение xk+1 из xk . Получим оценку погрешности на (k + 1)-й итерации из определения итерационного процесса: (xk − x∗ )f 0 (xk ) − f (xk ) f (xk ) ∗ − x = = f 0 (xk ) f 0 (xk ) F (xk ) = {обозначим F(x) = (x − x∗ )f 0 (x) − f(x)} = 0 k . f (x ) xk+1 − x∗ = xk − Преобразуем F (xk ) по формуле Ньютона-Лейбница: k k ∗ Zx k ∗ ∗ F (ξ)dξ = {F(x ) = f(x ) = 0} = F (x ) = F (x ) + x∗ k Zx Zx F (ξ)dξ = x∗ (ξ − x∗ )f 00 (ξ)dξ. x∗ Применив к последнему интегралу формулу среднего значения (ξk ∈ [x∗ ; xk ]), получим k F (xk ) = f 00 (ξk ) Zx (ξ − x∗ )dξ = f 00 (ξk ) (xk − x∗ )2 . 2 x∗ Подставив запись для F (xk ) в выражение для погрешности, получим xk+1 − x∗ = f 00 (ξk ) (xk − x∗ )2 . 2f 0 (xk ) Так как вторая производная по модулю ограничена сверху, а первая — снизу, то из последнего равенства следует, что метод Ньютона имеет квадратичную скорость сходимости. Докажем оценку на погрешность по индукции. 50 Глава 3. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ УРАВНЕНИЙ База индукции. Рассмотрим погрешность при k = 1 : x1 − x∗ = f 00 (ξ0 ) |x0 − x∗ |2 M2 |x0 − x∗ | 0 ∗ 6 {ξ ∈ U (x )} 6 |x − x∗ | = q|x0 − x∗ |. r 2f 0 (x0 ) 2m1 Таким образом, база индукции верна. Предположение индукции. Пусть оценка (3.9) выполняется для некоторого k : |xk − x∗ | 6 q 2 k −1 |x0 − x∗ |. Индуктивный переход. Докажем, что она выполняется для k + 1. Согласно показанному ранее, |xk+1 − x∗ | = |f 00 (ξ k )| |xk − x∗ |2 . 2|f 0 (xk )| ξk ∈ Ur (x∗ ), так как ξk выбиралась в соответствии с теоремой о среднем, и поэтому принадлежит отрезку Ur (x∗ ). Получим верхнюю оценку, используя предположение индукции: |f 00 (ξ k )| M2 (xk − x∗ )2 M2 2k −1 2 |xk − x∗ |2 6 6 (q ) · |x0 − x∗ |2 . k 2|f (x )| 2m1 2m1 Вспомним, что мы обозначали M2 0 |x − x∗ | = q, тогда получим: 2m1 |xk+1 − x∗ | 6 q 2 k+1 −1 |x0 − x∗ |. Откуда следует, что оценка верна, и, следовательно, теорема доказана (q < 1 по предположению, и при k −→ ∞ правая часть стремится к нулю, а это значит, что последовательность сходится к x∗ ). Замечание. M2 |x0 − x∗ | < 1. Но как это проверить, ведь мы не знаем 2m1 ∗ точного решения x ? Можно поступить так. 1. В условии теоремы мы требуем, чтобы Рассмотрим условие M2 |x0 − x∗ | < 1. 2m1 (3.10) Распишем f (x0 ) = f (x0 ) − f (x∗ ) = f 0 (x)(x0 − x∗ ), тогда |x0 − x∗ | = |f (x0 )| |f (x0 )| 6 . |f (x)| m1 Подставим эту запись для |x0 − x∗ | в (3.10), тогда из неравенства M2 |f (x0 )| <1 2m21 следует, что выполняется условие (3.10). Таким образом, зная m1 и M2 , можно подбирать x0 , исходя из этого неравенства (подбираем достаточно малую окрестность на этапе локализации корней, и дальше работаем с ней; если окрестность велика, уточняем расположение корня). 2. В условии требовалось, чтобы x∗ был простым вещественным корнем; если же x∗ — корень кратности p, то метод Ньютона будет иметь квадратичную скорость сходимости, если некоторым образом подправить итерационную последовательность. 3.6. РЕШЕНИЕ СИСТЕМ НЕЛИНЕЙНЫХ УРАВНЕНИЙ 3.6 51 Решение систем нелинейных уравнений Перейдем к поиску численных методов для решения систем нелинейных уравнений. Пусть имеется n уравнений следующего вида: fi (x1 , x2 , . . . , xn ) = 0, i = 1, n. (3.11) Поиск решений данной системы обычно проводится в два этапа: сначала проходит разделение решений (термин имеет тот же смысл, что и локализация корней при решении одиночных уравнений), а затем на полученных участках производится их уточнение. Как мы увидим, для поиска решений систем будут применяться те же методы, что и для поиска корней одиночных уравнений, но с небольшими модификациями. Метод простой итерации Итак, мы решаем систему вида (3.11). Построим итерационный процесс для нахождения точного решения (обозначим его x∗ ) на отрезке [a; b]. Для этого потребуем, чтобы оно являлось решением такой системы уравнений: (3.12) xi = Si (x1 , x2 , . . . , xn ), i = 1, n, где Si — некоторая функция. Будем задавать Si (x) в виде: Si (x1 , x2 , . . . , xn ) = xi − τi (x1 , x2 , . . . , xn )fi (x1 , x2 , . . . , xn ). где τi — функция-параметр, не обращающаяся в нуль в некоторой окрестности x∗ . Легко проверить, что в данном случае x∗ будет решением (3.12). Теперь зададим итерационный процесс следующим образом (xkl — координаты вектора xk , который должен будет сходиться к решению): = Sl (xk1 , xk2 , . . . , xkn ), xk+1 l k = 0, 1, . . . При этом задается вектор x0 — начальное приближение. Итак, мы ожидаем сходимость:  k→∞ −→ x∗ .  xk = (xk1 , xk2 , . . . , xkn )T Для одиночных уравнений достаточным условием сходимости было выполнение неравенства |S 0 (x)| < 1 всюду на рассматриваемом отрезке. Посмотрим, что нам потребуется в данном случае для уменьшения погрешности на каждом шаге: zik+1 = xk+1 − x∗i = Si (xk1 , xk2 , . . . , xkn ) − Si (x∗1 , x∗2 , . . . , x∗n ). i Пусть у функций Si существуют первые производные. Тогда мы можем применить обобщенную формулу Лагранжа и получим: Si (xk1 , xk2 , . . . , xkn ) − Si (x∗1 , x∗2 , . . . , x∗n ) = n X ∂Si (ξ1k , ξ2k , . . . , ξnk ) k (xl − x∗l ). ∂xl l=1 Соберем все производные в одну матрицу: Ak = (akij ) = формула может быть переписана так: z k+1 = Ak z k .  ∂Si (ξ1k , ξ2k , . . . , ξnk ) ∂xj  . Тогда предыдущая 52 Глава 3. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ УРАВНЕНИЙ Из этого следует, что ||z k+1 || = ||Ak z k || 6 ||Ak || · ||z k ||. Таким образом, если для всех k на нашем отрезке [a; b] выполняется неравенство ||Ak || 6 q < 1, (3.13) то последовательность ||z k || будет сходится по норме. Этого нам будет  достаточно.  ∂Si Заметим, что условие (3.13) будет выполнено, если матрица A = max будет иметь норму, [a; b] ∂xj меньшую единицы. Тогда, очевидно, ||Ak || 6 ||A||, и процесс будет сходиться. Такое требование можно удовлетворить, подбирая параметры τi . В частности, связывая их с производными функций fi (x), можно получить метод Ньютона. Метод Ньютона Здесь мы требуем от функций fi (x1 , x2 , . . . , xn ) существование первых производных. Согласно определению x∗ , fi (x∗1 , x∗2 , . . . , x∗n ) = 0, i = 1, n. (3.14) Теперь с помощью преобразований этого тождества получим формулу для итерационного процесса. Для этого зададимся k-м приближением xk = (xk1 , xk2 , . . . , xkn )T и зафиксируем некоторое i ∈ [1; n]. Заметим, что (3.14) можно переписать так: fi (xk1 + (x∗1 − xk1 ), xk2 + (x∗2 − xk2 ), . . . , xkn + (x∗n − xkn )) = 0. Применив обобщенную формулу Лагранжа, получим: fi (xk1 , xk2 , . . . , xkn ) + n X ∂fi (ξ1k , ξ2k , . . . , ξnk ) ∗ (xl − xkl ) = 0. ∂xl l=1 Теперь, заменив ξ k на xk , а x∗l — на новое приближение xlk+1 , получим такую формулу для поиска : n X ∂fi (xk1 , xk2 , . . . , xkn ) k+1 fi (xk1 , xk2 , . . . , xkn ) + (3.15) (xl − xkl ) = 0. ∂xl xk+1 l l=1 Обозначив ∆xkl = xk+1 l − xkl , перепишем ее так: n X ∂fi (xk1 , xk2 , . . . , xkn ) ∆xkl = −fi (xk1 , xk2 , . . . , xkn ). ∂xl l=1 Это система линейных уравнений для поиска ∆xk вида A∆xk = −f. Решая ее, мы находим ∆xk , а затем и xk+1 = xkl + ∆xkl для l = 1, n. l Теорему о сходимости данного метода мы напишем неформально и примем без доказательства. Теорема 3.3. В достаточно малой окрестности искомого корня итерационный процесс по методу Ньютона (задаваемый формулой (3.15)) сходится, если определитель матрицы A не обращается в этой окрестности в нуль. При этом скорость сходимости — квадратичная. Замечание. Данный итерационный процесс сходится быстро: для достижения неравенства (обычный для метода Ньютона критерий завершения ИП) ||xk+1 − xk || < ε при ε = 10−5 достаточно всего 3-5 итераций. Однако он требует большого объема вычислений — ведь на каждом шаге нам приходится решать систему уравнений. 3.6. РЕШЕНИЕ СИСТЕМ НЕЛИНЕЙНЫХ УРАВНЕНИЙ 53 Примеры решения нелинейных уравнений Пример 3.1. (Метод простой итерации) Пусть f (x) = x3 − x − 1. Найдем корень уравнения x3 − x − 1 = 0 на отрезке [−2; 3]. Для начала проведем локализацию корней: введем на отрезке сетку и посчитаем значение функции в ее узлах: −2 −7 xi f (xi ) −1 −1 −1 1 −1 2 5 3 23 Так как f (1)f (2) < 0, а наша функция непрерывна, то на отрезке [1; 2] обязательно есть корень уравнения. Будем искать его методом простой итерации. 1. Как уже говорилось, мы ставим в соответствие уравнению f (x) = 0 уравнение x = S(x), где S(x) наиболее часто берется в виде S(x) = x − τ (x)f (x). Выберем τ (x) = τ > 0. Тогда S(x) = x − τ (x3 − x − 1). Соответственно, итерационный процесс задается так: xk+1 = S(xk ) = xk − τ ((xk )3 − xk − 1). Для его сходимости, согласно замечанию к теореме 3.1, нам было достаточно выполнения неравенства |S 0 (x)| < 1 на [1; 2]. Посмотрим, какие ограничения это условие даст на τ : |S 0 (x)| < 1 ⇐⇒ |1 − τ (3x2 − 1)| < 1 ⇐⇒ −1 < 1 − τ (3x2 − 1) < 1. (3.16) При x ∈ [1; 2] и для любого положительного τ правое неравенство в (3.16) верно всегда. Левое неравенство перепишется так: 2 . τ< 2 3x − 1 При x ∈ [1; 2] знаменатель дроби достигает минимума при x = 2. Отсюда итоговое ограничение на τ таково: 2 2 = . τ< 3 · 22 − 1 11 1 1 3 9 . В этом случае S(x) = x − (x − x − 1), max |S 0 (x)| = . Это число — 11 11 11 x∈[1; 2] знаменатель геометрической прогрессии, характеризующей скорость сходимости процесса. Как видно, оно не очень мало, и сходится все медленно. Насколько медленно, можно понять из таблицы первых приближений: Возьмем τ = k 1 2 3 4 5 — в данном примере x∗ ≈ 1.32472. xk 1.1 1.16991 1.22161 1.25784 1.28218 1.29802 xk+1 = S(xk ) 1.16991 1.22161 1.25784 1.28218 1.29802 1.30812 54 Глава 3. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ УРАВНЕНИЙ 2. Можно по-разному выбирать функцию S(x) : x3 − x − 1 = 0 =⇒ x = — и можно взять S(x) = √ 3 √ 3 x + 1. x + 1. При этом скорость сходимости будет выше, так как S 0 (x) = 1 3(x + 1) 2 3 1 < 2 3(1 + 1) 3 ≈ 0.2 = q. 3. С другой стороны, x3 − x − 1 = 0 ⇐⇒ x = x3 − 1. Однако, если взять S(x) = x3 − 1, то S 0 (x) = 3x2 > 1 на [1; 2], что противоречит достаточному условию сходимости. Пример 3.2. (Метод простой итерации и метод Ньютона) Будем искать корень уравнения x3 − 1 = 0. Зная, что x∗ = 1 подходит, не будем утруждать себя локализацией корней, а просто покажем, как проходят итерационные процессы различных методов. Метод простой итерации. Возьмем отрезок [ 32 , 43 ] и зададим на нем итерационный процесс: xk+1 = S(xk ), где зададим S(x) = x−τ (x3 −1). Рассмотрим, какие ограничения накладываются на τ в этом примере: S 0 (x) = 1 − 3τ x2 . Потребуем выполнение достаточного условия сходимости: |S 0 (x)| < 1 ⇐⇒ −1 < 1 − 3τ x2 < 1. Отсюда ограничение на τ таково: 0<τ < Возьмем τ = 0.25, тогда S(x) = x − процесса, взяв x0 = 1.1 : k 1 2 3 4 5 2 3x2 = x= 34 3 . 8 x3 − 1 . Распишем первые несколько шагов итерационного 4 xk 1.1 1.01725 1.00409 1.00101 1.00025 1.00006 xk − x∗ 0.1 0.01725 0.00409 0.00101 0.00025 0.00006 Видно, что скорость сходимости не очень высокая. Однако на третьем и четвертом шаге уже можно применить коррекцию Эйткена: xk+2 = xk+2 − Тогда получим такую таблицу: (xk+2 − xk+1 )2 . − 2xk+1 + xk xk+2 3.6. РЕШЕНИЕ СИСТЕМ НЕЛИНЕЙНЫХ УРАВНЕНИЙ xk − x∗ 0.1 0.01725 0.00409 0.000068 0.000001 xk 1.1 1.01725 1.00409 1.000068 1.000001 k 1 2 3 4 55 Видно, что точность существенно возросла. Метод Ньютона. Согласно канонической форме метода Ньютона, xk+1 = xk − f (xk ) . f 0 (xk ) Зная, что f (x) = x3 − 1, а f 0 (x) = 3x2 , получим такую формулу: xk+1 = xk − (xk )3 − 1 . 3(xk )2 Таблица приближений такова: k 1 2 3 xk 1.1 1.00882 1.00008 1.000000006 xk − x∗ 0.1 0.00882 0.00008 0.000000006 Здесь уже заметно, что скорость сходимости — квадратичная. Пример 3.3. (Метод Ньютона для систем уравнений) Рассмотрим следующую систему нелинейных уравнений:  F (x, y) = x2 + y 2 − 4 = 0; G(x, y) = xy − 1 = 0. (3.17) На первом этапе (разделение корней) получаем 4 корня. Это можно показать, нарисовав графики соответствующих кривых. 2 1 3 4 Возьмем начальное приближение для поиска одного из корней таким: x0 = 2, y 0 = 0. Решать систему будем методом Ньютона. Линеаризованные уравнения имеют вид: ∂F k k (x , y )∆xk + ∂x ∂G k k (x , y )∆xk + ∂x ∂F k k (x , y )∆y k ∂y ∂G k k (x , y )∆y k ∂y = −F (xk , y k ); k = 0, 1, . . . k k = −G(x , y ), 56 Глава 3. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ УРАВНЕНИЙ Учитывая вид системы (3.17), можно найти частные производные: ∂F = 2x, ∂x ∂G =y, ∂x ∂F = 2y; ∂y ∂G =x. ∂y Подставив их в (3.17), можно найти приращения, решив систему:  2xk ∆xk + 2y k ∆y k = 4 − (xk )2 − (y k )2 ; y k ∆xk + xk ∆y k = 1 − xk y k . и подставить в определение итерационного процесса:  k+1 x = xk + ∆xk ; y k+1 = y k + ∆y k . Получим следующую последовательность итерационных приближений: N 1 xk 2 yk F G −1 2 2 0.5 0.25 3 1.93 0.517 −0.0077 −0.0022 ∆xk 1 − 15 ... ∆y k 0.5 1 60 ... Как видно из приведенной таблицы, процесс очень быстро сходится. Об этом можно судить по величине функций F (x, y) и G(x, y), так как мы решаем уравнения F (x, y) = 0 и G(x, y) = 0. Как же работает метод Ньютона? Дадим геометрическую интерпретацию этого метода. В каждой текущей точке (xk , y k ) к поверхностям z = F (x, y) и z = G(x, y) строятся касательные плоскости. Потом рассматриваем линии пересечения этих плоскостей с плоскостью z = 0 — это две прямые. И в заключение, определяем точку пересечения полученных прямых как новое, k + 1 значение. Глава 4 Интерполяция и приближение функций В этой главе мы будем восстанавливать функцию по значениям в некоторых заданных точках. Итак, пусть на отрезке [a; b] задан набор точек (n + 1 точка). Эти точки называют узлами интерполяции. Занумеруем их в следующем порядке a = x0 < x1 < . . . < xn = b, и пусть в каждой из этих точек известно f (xk ) = fk , k = 0, n. Наша задача заключается в том, чтобы вычислить приближенные (с некоторой точностью) значения функции между значениями в узлах интерполяции. Опираясь на эти сведения, построим на том же отрезке функцию Φ(x), которую назовем интерполянтой — она и будет служить приближением исходной функции. Не любая Φ(x) нам подойдет — потребуем, чтобы она была легко вычислима, и совпадала с исходной функцией в узлах интерполяции: Φ(xk ) = f (xk ), k = 0, n. Иногда, правда, от последнего условия отказываются — тогда говорят о построении наилучшего приближения. Впрочем, вся терминология будет объясняться по ходу дела. Построим интерполянту. Введем базисные функции ϕi (x), i = 0, m — линейно независимые элементы в нашем пространстве функций. Φ(x) будем строить как линейную комбинацию базисных функций ϕi (x) : m X Φ(x) = ai ϕi (x). i=0 В качестве базисных функций мы можем выбрать следующие: 1. Степенные функции: ϕi (x) = xi ; 2. Тригонометрические функции: в случае интерполяции периодических функций с периодом, к πix примеру, 2l лучше взять sin πix l и cos l ; 3. Дробно-полиномиальные функции. Интерполяция алгебраическими многочленами Одним из наших первых требований было совпадение значений интерполянты и исходной функции в узлах сетки. Записав интерполянту через базисные функции (пусть их будет n + 1 — столько же, сколько узлов), получим такие уравнения: n X ai ϕi (xk ) = f (xk ), k = 0, n. i=0 57 (4.1) Глава 4. ИНТЕРПОЛЯЦИЯ И ПРИБЛИЖЕНИЕ ФУНКЦИЙ Относительно ai мы получили СЛАУ с n + 1 неизвестными, так как ϕi (xk ) и f (xk ) заданы. Условием разрешимости в данном случае является то, что определитель системы (4.1) отличен от нуля. При этом решение будет существовать и оно будет единственным. Определитель матрицы системы будет выглядеть так: ϕ0 (x0 ) ϕ1 (x0 ) . . . ϕn (x0 ) ϕ0 (x1 ) ϕ1 (x1 ) . . . ϕn (x1 ) (4.2) . . . . . . . ϕ0 (xn ) ϕ1 (xn ) . . . ϕn (xn ) Будем брать в качестве базисных функций степенные. Тогда (4.2) — это определитель Вандермонда: 1 1 x0 x1 . . 1 xn x20 x21 . . . x2n ... ... . . ... xn0 xn1 = (x1 − x0 )(x2 − x0 ) . . . (xn − x0 )(x2 − x1 )(x3 − x1 ) . . . (xn − xn−1 ) = xnn Y = (xm − xk ) 6= 0, если xi 6= xj ∀ i 6= j. 06k n — сплайнами). Теперь обсудим случай, когда m < n (узлов, в которых известно значение функции, больше, чем неизвестных коэффициентов). В этом случае возникает понятие наилучшего приближения: мы отказываемся от требования совпадения значения интерполянты и табличной функции, а требуем минимизировать некоторый функционал от, к примеру, вектора погрешностей в узлах сетки: ~r = (Φ(x0 ) − f (x0 ), Φ(x1 ) − f (x1 ), . . . , Φ(xn ) − f (xn )) , где x0 , . . . , xn — точки, в которых нам задана функция f (x). В качестве функционала рассмотрим норму этого вектора, которую можно задавать по-разному: ! 21 n X (Φ(xi ) − f (xi ))2 ; ||~r||(1) = i=0 ||~r||(2) = max |Φ(xi ) − f (xi )|. i В каждом случае мы пытаемся подобрать коэффициенты ai в задании Φ(x) так, чтобы минимизировать эту норму. Эту задачу называют поиском наилучшего среднеквадратичного приближения (а метод ее решения — методом наименьших квадратов) или поиском наименьшего равномерного приближения соответственно. Приведем пример решения такой задачи. 68 Глава 4. ИНТЕРПОЛЯЦИЯ И ПРИБЛИЖЕНИЕ ФУНКЦИЙ Наилучшее приближение в гильбертовом пространстве p Пусть H — евклидово пространство функций с нормой ||f || = hf, f i, а ϕi (i = 0, n)— его линейно независимые элементы. Нашей задачей будет поиск наилучшего приближения ci ∈ R, ϕ = c0 ϕ0 + . . . + cn ϕn , для элемента f ∈ H. Оценкой точности будет служить величина погрешности ||f − ϕ||. Определение. Элемент ϕ, e доставляющий минимум этой норме (для которого верно равенство min ||f − ϕ|| = ||f − ϕ||), e называется элементом наилучшего приближения. ϕ Покажем, что такой элемент существует и единствен. В качестве примера можно взять пространство L2 [a; b] (пространство функций, интегрируемых со своими квадратами) и в нем некоторую функцию f. В качестве скалярного произведения берут обычное в L2 : 1 Zb hg, f iL2 = Zb f (x)g(x) dx, соответственно ||f ||L2 = a 2 2 f (x) dx . a Продолжим рассуждения в общем виде. Наша задача — это минимизировать ||f − ϕ||, подобрав соответствующую функцию ϕ. Приведем эту норму (будем работать с ее квадратом) к более удобному виду: * + n n X X 2 ||f − ϕ|| = hf − ϕ, f − ϕi = f − cl ϕl , f − ck ϕk = = hf, f i − n X cl hϕl , f i − n X cl hϕl , f i + = { Введем обозначения = ||f ||2 − 2 l=0 cl fl + ck hϕk , f i + n X n X n n X X cl ck hϕl , ϕk i = l=0 k=0 cl ck hϕl , ϕk i = l=0 k=0 l=0 n X k=0 k=0 l=0 = ||f ||2 − 2 l=0 n X Rb     fl = hϕl , f i =    akl = hϕk , ϕl i = n X n X f (x)ϕl (x) dx; a Rb } = ϕk (x)ϕl (x) dx. a ck cl akl = l=0 k=0   c = (c0 , c1 , . . . , cn )T ; = { Обозначим f = (f0 , f1 , . . . , fn )T ; } =  A = (akl ). = ||f ||2 − 2 c, f + hAc, ci = ||f ||2 + F (c). Соединив первое и последнее равенство, получим: ||f − ϕ||2 = ||f ||2 + F (c). (4.18) Таким образом, задача о минимизации ||f − ϕ|| свелась к задаче минимизации функции F (c) от вектора переменных c. Из определения матрицы A следует, что она симметрична. Покажем, что она положительно определена, то есть ∀c 6= 0 hAc, ci > 0. 4.2. НАИЛУЧШЕЕ ПРИБЛИЖЕНИЕ ТАБЛИЧНОЙ ФУНКЦИИ 69 Если взять в равенстве (4.18) f ≡ 0, то получим, что hAc, ci = ||ϕ||2 > 0. Предположим, что существует вектор y 6= 0 такой, что hAy, yi = 0. Но это будет означать, что ||ϕ|| = 0. Так как ϕ — линейная комбинация линейно независимых элементов ϕi , то это возможно тогда и только тогда, когда эта комбинация тривиальна — то есть yi = 0 для всех i. Отсюда делаем вывод, что y = 0 =⇒ ∀c 6= 0 hAc, ci > 0. Таким образом, матрица A положительно определена. Это позволяет воспользоваться следующей теоремой. Теорема 4.4. Пусть даны матрица A такая, что A = AT > 0, и f — некоторый вектор (соответствующей размерности). Тогда у функции F (c) = hAc, ci − 2 f , c векторного переменного c точка минимума существует и единственна, причем элемент c реализует этот минимум тогда и только тогда, когда он является решением системы: Ac = f . (4.19) Доказательство. Сначала докажем утверждение об эквивалентности. Достаточность. Пусть c — решение системы (4.19). Покажем, что для любого ненулевого вектора v F (c + v) > F (c) : F (c + v) = hA(c + v), c + vi − 2 f , c + v = hAc, ci + hAc, vi + hAv, ci + hAv, vi − 2 f , c − 2 f , v = = {A = AT =⇒ hAc, vi = hAv, ci} = F (c) + hAv, vi + 2 v, Ac − f = F (c) + hAv, vi . Так как A > 0, то hAv, vi > 0. Это означает, что c — точка минимума . Достаточность доказана. Необходимость. Пусть c — точка минимума F (c). Фиксируем произвольный ненулевой вектор y и положим v = λy (λ — параметр). Тогда, согласно проведенным в доказательстве достаточности преобразованиям, F (c + v) = F (c) + hAv, vi + 2 v, Ac − f = F (c) + λ2 hAy, yi + 2λ y, Ac − f . (4.20) Обозначим выражение, стоящее в правой части равенства, за g(λ). Из принятых условий следует, что в точке λ = 0 функция g(λ) достигает минимума. Она, очевидно, дифференцируема, поэтому g 0 (0) = 0. Продифференцировав (4.20) по λ, получим, что
Ac − f , y = 0. 2λ hAy, yi + 2 y, Ac − f λ=0 = 0 ⇐⇒ Из произвольности выбора y следует, что Ac − f = 0. Это означает, что c — решение системы (4.19). Необходимость доказана. Осталось заметить, что из доказанной эквивалентности следует существование и единственность точки минимума функции F (c). Это вытекает из того, что матрица A положительно определена, поэтому система (4.19) имеет единственное решение. Теорема полностью доказана. 70 Глава 4. ИНТЕРПОЛЯЦИЯ И ПРИБЛИЖЕНИЕ ФУНКЦИЙ Алгоритм построения наилучшего приближения Опираясь на доказанную теорему, можно построить алгоритм нахождения наилучшего приближения для функции f ∈ L2 [a; b]. Он будет выглядеть так: 1. Выбираем (n + 1) линейно независимый элемент ϕk , k = 0, n из L2 [a; b]. 2. Строим матрицу A = (akl ) скалярных произведений: Zb ϕk (x)ϕl (x) dx. akl = a 3. Готовим вектор скалярных произведений f = (f0 , f1 , . . . , fn ), где fi находятся так: Zb fi = f (x)ϕi (x) dx, i = 0, n. a c0 , e c1 , . . . , e cn ), решая систему уравнений Ac = f . 4. Ищем вектор коэффициентов c = (e 5. Строим элемент ϕ e: ϕ e=e c0 ϕ0 + e c1 ϕ1 + . . . + e cn ϕn . Он будет являться наилучшим приближением согласно доказанной теореме. Теперь посмотрим, насколько точно ϕ e приближает f, то есть оценим ||f − ϕ||. e Для этого нам понадобится лемма. Лемма. Пусть ϕ e — элемент наилучшего приближения для f. Тогда hf − ϕ, e ϕi e =0 — то есть ϕ e ортогонален (f − ϕ). e Доказательство. Подставим представление ϕ e через ϕk в искомое скалярное произведение: * + n n n n X n X X X X hf − ϕ, e ϕi e = f− e ck ϕk , e cl ϕl = e cl hf, ϕl i − e ck e cl hϕk , ϕl i = k=0 = n X e cl fl − l=0 n n X X l=0 l=0 k=0 l=0 e ck e cl akl = {согласно старым обозначениям } = c, f − hAc, ci = f − Ac, c . k=0 l=0 Согласно построению ϕ, e c находился из условия Ac = f , поэтому получаем, что hf − ϕ, e ϕi e = 0. Лемма доказана. 4.2. НАИЛУЧШЕЕ ПРИБЛИЖЕНИЕ ТАБЛИЧНОЙ ФУНКЦИИ 71 Теперь можно оценить отклонение: ||f − ϕ|| e 2 = hf − ϕ, e f − ϕi e = hf − ϕ, e f i − hf − ϕ, e ϕi e = hf − ϕ, e f i = hf, f i − hϕ, e fi . Согласно лемме, hf, ϕi e = hϕ, e ϕi e , поэтому ||f − ϕ|| e 2 = hf, f i − hϕ, e ϕi e = ||f ||2 − ||ϕ|| e 2. Замечание. Если {ϕk } — ортонормированная система, то есть hϕ ek , ϕ el i = δkl , тогда A = E. Отсюда следует, что ck = hf, ϕk i = fk , и для наилучшего приближения получается простая формула: ϕ e= n X fk ϕk . i=0 В этом случае коэффициенты ck называются коэффициентами Фурье, а построенный элемент ϕ e — многочленом Фурье. Пример 4.3. Построим для функции из предыдущего примера наилучшее приближение. Итак, f (x) задана таблично в точках x0 , x1 = x0 + h, x2 = x0 + 2h. Обозначим F0 = f (x0 ), F1 = f (x1 ), F2 = f (x2 ). Приближать будем снова многочленами: ϕ0 (x) = 1, ϕ1 (x) = x − x1 . Напомним, что, применяя метод наименьших квадратов, мы нашли такое приближение: Φ(x) = F2 − F0 F0 + F1 + F2 + (x − x1 ) 3 2h Для поиска наилучшего приближения для f последовательно пройдем по построенному нами алгоритму: 1. ϕ0 (x) = 1, ϕ1 (x) = x − x1 . xZ1 +h 2. Подсчитаем akl = ϕk (x)ϕl (x) dx : x1 −h a00 = 2h; xZ1 +h a10 (x − x1 )2 (x − x1 ) dx = 2 = a01 = x1 −h xZ1 +h a11 (x − x1 )2 dx = = x1 −h  2h =⇒ A =  (x − x1 )3 3 x1 +h = 0; x1 −h x1 +h = x1 −h h3 2h3 h3 + = . 3 3 3  2h3  . 3 3. На данном этапе возникают сложности, так как, не зная f (x), мы не можем точно вычислить f0 и f1 . Будем вычислять их приближенно: f0 через формулу среднего значения, а f1 — по 72 Глава 4. ИНТЕРПОЛЯЦИЯ И ПРИБЛИЖЕНИЕ ФУНКЦИЙ формуле Симпсона (как известно, она дает маленькую погрешность): xZ1 +h f0 f (x) dx = f (x) · 2h ≈ 2h = F0 + F1 + F2 ; 3 x1 −h xZ1 +h f1 Zb b−a (G0 + 4G1 + G2 )} ≈ f (x)(x − x1 ) dx ≈ { G(x) dx ≈ 6 = a x1 −h ≈ h2 2h [F0 (−h) + 4F1 · 0 + F2 · h] = (F2 − F0 ) 6 3 4. Теперь решаем систему Ac = f . Матрица A — диагональная, поэтому ее решение запишется просто:   F0 + F1 + F2 F0 + F1 + F2     ; c0 = 2h ; c0 =  e  2he 3 3 ⇐⇒   2h3 F2 − F0 h2    e  c1 = . e c1 = (F2 − F0 ). 2h 3 3 5. Вычислив коэффициенты e c0 и e c1 , можем записать построенное приближение: F0 + F1 + F2 F2 − F0 ϕ(x) e = + (x − x1 ). 3 2h Оно совпало с построенной ранее функцией Φ(x). Совпадение это не случайно: мы очень неточно вычислили f0 , хотя могли бы этого избежать. Вычислим f0 , применяя формулу Симпсона: xZ1 +h f (x) dx ≈ h [F0 + 4F1 + F2 ]. 3 x1 −h Заново решив систему Ac = f , получим такие выражения для e c0 и e c1 :  F + 4F + F 1 2   ; c0 =  e 6  F2 − F0  e  c1 = . 2h Эти коэффициенты дадут более точное приближение: F0 + 4F1 + F2 F2 − F0 ϕ(x) e = + (x − x1 ). 6 2h Естественно задаться вопросом: «А насколько оно точнее?». Легко показать, что ϕ e отличается от ϕ e на константу: h2 F0 − 2F1 + F2 h2 F0 − 2F1 + F2 ϕ(x) e − ϕ(x) e = = · = fxx,1 , 6 6 h2 6 где fxx,1 — вторая разностная производная f (x) в точке x1 . Ранее было показано, что такое же значение принимает квадрат нормы погрешности между ϕ e и f — то есть ||ϕ(x) e − f (x)||2 : ||ϕ(x) e − f (x)||2 = 1 (F0 − 2F1 + F2 )2 . 6 Неформально можно сказать, что f −ϕ e≈ϕ e − ϕ. e Столь большие отклонения возникают из-за неточности при вычислении f . Отсюда следует вывод, что алгоритм построения наилучшего приближения слишком зависит от методов приближенных вычислений, и лучше использовать единые формулы, например, для подсчета интегралов — ту же формулу Симпсона. При этом получаются неплохие результаты. Глава 5 Численные методы решения краевых задач Нашей первой задачей будет поиск численных методов решения задачи Коши:  du   = f (t, u(t)), 0 < t < T ; dt   u| = u . t=0 Для начала необходимо построить соответствующую дискретную модель. Для этого разобьем весь отрезок [0; T ] на точки ωτ = {tn = nτ, n = 0, 1, . . . Tτ }, где τ — диаметр дискретной сетки. Обозначим значения искомой функции — un = u(tn ), приближенное решение — yn , и погрешность на n-й итерации как zn = yn − un в узлах сетки (zn , yn , un — сеточные функции). От будущего алгоритма потребуем как можно более точного воспроизведения функции u — для этого нужно, чтобы погрешность zn была мала. Обсудим понятие сходимости приближенного решения к точному. Фиксируем точку tn и построим последовательность сеток ωτ такую, что точка tn является узлом для сеток с номерами m > k, то есть при сгущении сетки только добавляются новые узлы. На каждой из этих сеток строится сеточная функция yn . τ →0 Определение. Сеточная функция yn сходится к решению un в узле tn , если |zn | −→ 0. Определение. Сходимость на отрезке означает сходимость в каждой точке этого отрезка. Определение. Пусть погрешность по порядку роста ведет себя как |zn | = O(τ p ), тогда приближенное решение имеет p-й порядок точности. Если мы будем приближать производную в узлах сетки ее разностным аналогом: yn+1 − yn u0 (tn ) = , τ то исходное уравнение примет такой вид: yn+1 − yn = f (tn , yn ), n = 0, 1, . . . (5.1) τ Получившиеся уравнения относительно yi называются разностной схемой. Из этого аналога нашего непрерывного уравнения можно определить значения yn во всех точках сетки, если стартовать с фигурирующего в условии задачи y0 .  yn+1 = yn + τ f (tn , yn ); n = 0, 1, . . . y0 = u 0 . 73 74 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ Получившаяся схема называется явной, так как (n+1)-е приближение выражается через n-е явно, непосредственно (не требуется решать уравнение). Рассмотрим другой вариант, в котором мы заменяем обыкновенную производную на разностную производную назад. Тогда функция f (t, y) будет браться в точке (tn+1 , yn+1 ) : yn+1 − yn = f (tn+1 , yn+1 ). τ При этом относительно yn+1 возникнет нелинейное уравнение:  yn+1 = τ f (tn+1 , yn+1 ) + yn , y = 0, 1, . . . ; y0 = u 0 . То есть на каждом шаге требуется решать нелинейное, вообще говоря, уравнение, и формула для y соответственно является неявной. Отсюда и название — неявные разностные схемы. un+1 − un + f (tn , un ) называется невязкой (погрешностью Определение. Значение ψn = − τ аппроксимации) расчетной схемы (5.1) в узле tn . Примечание. Вообще говоря, далее под невязкой мы будем понимать разность левой и правой частей расчетной схемы при подстановке в нее точного решения. Вышеприведенная формула выражает невязку для методов Рунге-Кутта с параметром m = 1. τ →0 Определение. Разностная схема аппроксимирует исходное ОДУ в узле tn , если ψn −→ 0 (этот предельный переход связан с последовательностью сеток). Определение. Разностная схема имеет p-й порядок аппроксимации в точке tn , если выполнено равенство: |ψn | = O(τ p ). Данную терминологию будем использовать для сравнения методов. 5.1 Сходимость методов Рунге-Кутта Методами Рунге-Кутта1 называется семейство методов, общий вид разностных схем которых задается так: m X yn+1 − yn = σi Ki (y), (5.2) τ i=1 где величины Ki вычисляются по следующим формулам: K1 (y) = f (tn , yn ); K2 (y) = f (tn + a2 τ, yn + b21 τ K1 (y)); K3 (y) = f (tn + a3 τ, yn + b31 τ K1 (y) + b32 τ K2 (y)); ... m−1 X Km (y) = f (tn + am τ, yn + bmi τ Ki (y)). i=1 Параметры ai , bij , σi выделяют конкретный метод Рунге-Кутта. В этих методах нужно вычислять значения функции f в промежуточных точках сетки. Их количество определяется параметром m, а соответствующие методы называются m-этапными методами. Схемы с m > 5 используются крайне редко (чаще всего используют 4-этапные методы). 1 Основная идея — К. Рунге (1885). 5.1. СХОДИМОСТЬ МЕТОДОВ РУНГЕ-КУТТА 75 Какие ограничения накладываются на параметры методов Рунге-Кутта для того, чтобы обеспечить сходимость? Попробуем ответить на этот вопрос. Потребуем, чтобы разностная схема (5.2) аппроксимировала исходное ОДУ в соответствии с ввеτ →0 денным определением ( |ψn | −→ 0). Невязка в данном случае будет иметь вид: m ψn = − un+1 − un X + σi Ki (u). τ i=1 Разложим un+1 в ряд Тейлора с центром в точке tn : un+1 = un + u0n τ + O(τ 2 ). Рассмотрим выражение для Ki (аналогично разложив его в ряд Тейлора): Ki (u) = f (tn , un ) + O(τ ). Подставим эти формулы в выражение для невязки: ψn = −u0n + f (tn , un ) m X σi + O(τ ) = f (tn , un ) −1 + i=1 m X ! σi + O(τ ). i=1 Очевидно, для того, чтобы разностная схема аппроксимировала исходное ОДУ (с порядком апm X σi = 1; тогда невязка будет равна проксимации p = 1, то есть ψn = O(τ )), достаточно, чтобы i=1 O(τ ). Таким образом, можно ввести первое ограничение на параметры метода (5.2): m X σi = 1. i=1 Теорема 5.1 (О сходимости методов Рунге-Кутта). Пусть метод Рунге-Кутта аппроксимирует исходное уравнение, тогда приближенное решение yn сходится к точному un , и порядок точности приближенного решения совпадает с порядком аппроксимации разностной схемы исходного ОДУ. Доказательство. Будем предполагать, что функция f (t, u) Липшиц-непрерывна по второму аргументу, то есть ∀ u1 , u2 |f (t, u1 ) − f (t, u2 )| 6 L|u1 − u2 |, где L — некоторая константа. Рассмотрим функцию погрешности: zn = yn − un . Выразим из нее yn = zn + un и подставим в левую часть (5.2). Тогда получим: m m X zn+1 − zn un+1 − un X =− + σi Ki (u) + σi (Ki (y) − Ki (u)). τ τ i=1 i=1 m m X un+1 − un X + σi Ki (u) = ψn , и обозначим ψn = σi (Ki (y) − Ki (u)). τ i=1 i=1 Оценим для разных i выражение |Ki (y) − Ki (u)| : Заметим, что − i=1: i=2: |K1 (y) − K1 (u)| = |f (tn , yn ) − f (tn , un )| 6 L|yn − un | = L|zn |. |K2 (y) − K2 (u)| = |f (tn + a2 τ, yn + b21 τ K1 (y)) − f (tn + a2 τ, un + b21 τ K1 (u))| 6 6 L|yn − un + b21 τ (K1 (y) − K1 (u))| 6 L(|zn | + τ |b21 | · |K1 (y) − K1 (u)|) 6 6 L(|zn | + τ |b21 |L|zn |) 6 {обозначим b = max bij } 6 L|zn |(1 + τ bL). i=2, n j=1, m−1 76 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ i=3: Если тоже самое проделать для i = 3, то получим такую оценку: |K3 (y) − K3 (u)| 6 L|zn |(1 + τ bL)2 . Теперь перейдем к общей оценке. Докажем, что |Kl (y) − Kl (u)| 6 L|zn |(1 + τ bL)l−1 , l = 1, m. Пусть эта оценка верна для некоторого i : |Ki (y) − Ki (u)| 6 L|zn |(1 + τ bL)i−1 , докажем ее для i + 1 :  |Ki+1 (y) − Ki+1 (u)| = f tn + ai+1 τ, yn + τ i X  b(i+1)j Kj (y) − j=1  −f tn + ai+1 τ, un + τ i X    X  i b(i+1)j Kj (u) 6 L yn − un + τ b(i+1)j (Kj (y) − Kj (u))  6 j=1 j=1 6 {воспользуемся неравенством bij 6 b и подсчитаем сумму геометрической прогрессии} 6     i i X X j−1 (1 + τ bL) = |Kj (y) − Kj (u)| 6 L |zn | + τ bL|zn | 6 L |yn − un | + τ b j=1 j=1   bL)i = L|zn | 1 + τ bL 1−(1+τ = L|zn |(1 + τ bL)i . 1−(1+τ bL) Таким образом, оценка (5.3) действительно имеет место. Получим теперь оценку на |ψn | : |ψn | 6 m X |σi | · |Ki (y) − Ki (u)| 6 {обозначим σ = max |σi |} 6 i=1, m i=1 6σ m X L|zn |(1 + τ bL)i−1 = σL|zn |(1 + τ bL)m−1 m. i=1 Оценим два последних множителя: m(1 + τ bL)m−1 6 {(1 + y)α 6 eαy } 6 meτ bL(m−1) 6 {τ 6 T} 6 meT bL(m−1) — обозначим за α, тогда |ψn | 6 |zn |σLα. Откуда, учитывая, что zn+1 − zn = ψn + ψn , получаем оценку на погрешность: τ |zn+1 | 6 |zn | + τ |ψn | + τ |ψn | 6 |zn | + τ |ψn | + τ |zn |σLα = = |zn |(1 + τ σLα) + τ |ψn | 6 |zn−1 |(1 + τ σLα)2 + τ |ψn−1 |(1 + τ σLα) + τ |ψn |. Применив эту же операцию n − 1 раз, получим: |zn+1 | 6 |z0 |(1 + τ σLα)n + τ n X |ψj |(1 + τ σLα)n−j . j=0 Если мы обозначим ψ = max |ψj | и учтем, что z0 = y0 − u0 = 0, то получим: j=0, n |zn+1 | 6 ψτ n X j=0 (1 + τ σLα)j 6 ψτ max (1 + τ σLα)j (n + 1). j=0, n (5.3) 5.2. МЕТОДЫ РУНГЕ-КУТТА ВТОРОГО ПОРЯДКА АППРОКСИМАЦИИ 77 В силу того, что τ n = T, |zn+1 | 6 ψT eT σLα , откуда следует, что если наша схема аппроксимирует на всей сетке ОДУ, то есть |ψi | → 0, то имеет место сходимость (|zn | → 0 ). Кроме того, если наша разностная схема аппроксимирует исходное ОДУ с p-м порядком аппроксимации (ψn = O(τ p )), то погрешность имеет соответственно p-й порядок: |zn | = O(τ p ). Таким образом, теорема полностью доказана. Теперь свяжем требование на порядок аппроксимации с количеством промежуточных точек, в которых требуется вычислять значение f (t, y). 5.2 Методы Рунге-Кутта второго порядка аппроксимации Рассмотрим семейство методов Рунге-Кутта при m = 2. Схема для вычисления приближенного значения ( yn+1 ) будет выглядеть так:  yn+1 − yn    τ    = σ1 K1 (y) + σ2 K2 (y); K1 (y) = f (tn , yn ); K2 (y) = f (tn + a2 τ, yn + τ b21 K1 (y)). Найдем, что является достаточным условием для достижения 2-го порядка точности. Для этого рассмотрим выражение для невязки: ψn = − un+1 − un + σ1 K1 (u) + σ2 K2 (u), τ (5.4) и потребуем, чтобы ψn = O(τ 2 ). Для начала распишем дробь в правой части (5.4), применив формулу Тейлора: un+1 − un 1 1 τ2 τ = (u(tn + τ ) − un ) = (un + u0n τ + u00n + O(τ 3 ) − un ) = u0n + u00n + O(τ 2 ). τ τ τ 2 2 Согласно постановке, K1 (u) = f (tn , yn ). Обозначим это число за fn . Используя это обозначение, разложим выражение для K2 (u) по формуле Тейлора для функции двух переменных: K2 (u) = f (tn + a2 τ, un + τ b21 fn ) = fn + ft0 (tn , yn )a2 τ + fu0 (tn , un )b21 fn τ + O(τ 2 ). Подставив это выражение в (5.4), получим такое выражение для невязки: ψn = −u0n − u00n τ + σ1 fn + σ2 (fn + ft0 (tn , yn )a2 τ + fu0 (tn , un )b21 fn τ ) + O(τ 2 ). 2 Так как u — точное решение, то для него справедливы формулы:  0 u = f (t, u); u00 = ft0 + fu0 u0 = ft0 + fu0 f. С их использованием выражение для невязки перепишется так: ψn = fn (−1 + σ1 + σ2 ) + ft0 (tn , yn )(− τ τ + σ2 a2 τ ) + fu0 (tn , yn )fn (− + σ2 b21 τ ) + O(τ 2 ). 2 2 78 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ Нетрудно подобрать такие σ1 , σ2 , a2 , b21 , чтобы первые три слагаемых обратились в ноль — это и даст требуемую оценку для невязки. Коэффициенты σ1 , σ2 , a2 , b21 должны быть такими, что  σ1 + σ2 = 1;      1 σ2 a2 = ; 2      σ b = 1. 2 21 2 (5.5) Из последних двух уравнений следует, что a2 и b21 равны друг другу. Обозначим это число за a. Теперь обозначим σ2 = σ, и тогда из первого уравнения будет следовать, что σ1 = 1−σ. Это позволяет переписать систему (5.5) так:   σ1 = 1 − σ;  σa = 1 . 2 Таким образом, методы с расчетной схемой следующего вида: yn+1 − yn = τ ((1 − σ)f (tn , yn ) + σf (tn + aτ, yn + af (tn , yn )τ )) 1 имеют 2-й порядок аппроксимации, если выполняется условие σa = . Согласно доказанной теореме, 2 построенное решение будет иметь 2-й порядок точности. 1 Можно привести примеры таких методов. При σ = 1, a = получается такая схема: 2 yn+1 − yn = τ f (tn + А при σ = τ τ f (tn , yn ) , yn + ). 2 2 1 , a = 1 — вот такая: 2 yn+1 − yn = τ [f (tn , yn ) + f (tn+1 , yn + τ f (tn , yn )] . 2 Несмотря на бесконечное количество схем заданного порядка точности, нам может не подойти ни одна. Это связано с тем, что методы Рунге-Кутта не являются устойчивыми, и при их использовании накапливается машинная погрешность, в конце вычислений сравнимая с полученными величинами. Подробнее о вычислительной устойчивости речь пойдет в следующих разделах. Методы Рунге-Кутта четвертого порядка точности Приведем без вывода расчетную схему 4-го порядка точности в методе Рунге-Кутта с параметром m=4:  1 yn+1 − yn   = (K1 (y) + 2K2 (y) + 2K3 (y) + K4 (y));   τ 6      K1 (y) = f (tn , yn );    τ τ K2 (y) = f (tn + , yn + K1 (y));  2 2    τ τ    K3 (y) = f (tn + , yn + K2 (y));   2 2    K4 (y) = f (tn + τ, yn + τ K3 (y)). Все методы Рунге-Кутта требуют возможность вычислять значение функции в произвольной точке. Теперь рассмотрим методы, в которых этого делать не надо. 5.3. ОПИСАНИЕ МНОГОШАГОВЫХ МЕТОДОВ 5.3 79 Описание многошаговых методов Определение. Линейным m-шаговым методом называется метод с расчетной схемой следующего вида: a0 yn + a1 yn−1 + a2 yn−2 + . . . + am yn−m = b0 fn + b1 fn−1 + . . . + bm fn−m . (5.6) τ где ai , bj — параметры метода, а yn−k и fn−i означают следующее:  yn−k = y(tn−k ); fn−i = f (tn−i , yn−i ). Таким образом, для реализации m-шагового метода на первом шаге требуется знать значения y0 , y1 , . . . , ym−1 . Значение y0 можно взять равным u0 — начальному условию, а вот для вычисления y1 , . . . , ym−1 применяют методы Рунге-Кутта соответствующего порядка точности. Заметим также, что в методе используются только табличные данные о f (x) — то есть уметь вычислять функцию f на промежуточных точках в общем случае не требуется, а может понадобиться только для получения y1 , . . . , ym−1 . Если в схеме (5.6) коэффициент b0 равен нулю, то в правой части fn не присутствует, и соответствующий метод называется явным (по тем же причинам, что и раньше). Если же b0 6= 0, то метод называется неявным (как ни странно, тоже по тем же причинам, что и раньше); возникает нелинейное уравнение относительно yn , которое в общем случае решается методом Ньютона. Очевидно, что метод не изменится, если выражение (5.6) домножить на какую-нибудь ненулевую m X константу. Поэтому устраним неоднозначность, введя условие нормировки: bi = 1. Покажем, что в i=0 этом случае правая часть уравнения (5.6) будет аппроксимировать правую часть дифференциального уравнения исходной задачи: f (tn , un ) − m X bi f (tn − iτ, u(tn − iτ )) = {разлагая слагаемые в сумме в ряд Тейлора} = i=0 = f (tn , un ) − m X bi [f (tn , un ) + O(τ )]gg = fn (1 − i=0 m X bi ) + O(τ ) = O(τ ) i=0 — это и означает аппроксимацию. Теперь выведем достаточные условия для достижения k-го порядка аппроксимации исходной функции. Для этого рассмотрим выражение для невязки: m X ψn = − m X ai u(tn − iτ ) i=0 + τ m X bi f (tn − iτ, u(tn − iτ )) = − ai un−i i=0 i=0 + τ m X bi f (tn−i , un−i ). i=0 Теперь разложим составляющие равенства в ряд Тейлора: u(tn − iτ ) = k (j) X un j=0 f (tn−i , un−i ) j! (−iτ )j + O(τ k+1 ); = {u0 = f(t, u)} = u0n−i = u0 (tn − iτ ) = k−1 X j=0 (j+1) un j! (−iτ )j + O(τ k ). Подставив эти формулы в выражение для невязки, получим:   m k m k−1 (j) X X X X u(j+1) 1 un n ψn = − ai  (−iτ )j  + bi (−iτ )j + O(τ k ). τ i=0 j! j! i=0 j=0 j=0 80 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ Поменяем порядки суммирования в двойных суммах, при этом из первой суммы вынесем отдельно составляющую при j = 0, а во второй сделаем замену j = j + 1. Тогда выражение для невязки перепишется так: ψn = − m k m k (j) m (j) X X 1X 1 X un X un un ai − ai (−iτ )j + bi (−iτ )j−1 + O(τ k ) = τ i=0 τ j=1 j! i=0 (j − 1)! i=0 j=1 = − m k m (j) X un X un j−1 X ai + τ (−i)j−1 (iai + jbi ) + O(τ k ). τ i=0 j! j=1 i=0 Заметим, что если все суммы подбором коэффициентов ai , bj обратить в нуль, то для невязки будет справедлива оценка: ψn = O(τ k ). Таким образом, достаточным условием k-го порядка аппроксимации будет выполнение системы равенств:  X m   ai = 0;   i=0 (5.7) m X   j−1  (−i) (iai + jbi ) = 0, j = 1, k.  i=0 Рассмотрим отдельно последнее условие при j = 1 : m X i=0 Согласно условию нормировки, m X iai + m X bi = 0. (5.8) i=0 bi = 1, поэтому (5.8) перепишется так: i=0 m X i=0 iai = −1 ⇐⇒ m X iai = −1. i=1 Добавив это уравнение и условие нормировки в систему (5.7), получим окончательный вариант достаточного условия k-го порядка аппроксимации:  m X   bi = 1;     i=0     m   X   ai = 0;    i=0 (5.9) m  X    iai = −1;     i=1     m  X    ij−1 (iai + jbi ) = 0, j = 2, k.  i=1 Мы получили систему из k + 2 линейных уравнений, решив которую, мы получим параметры, определяющие метод k-го порядка аппроксимации. Система содержит 2m + 2 неизвестных. Чтобы она не была переопределенной, потребуем, чтобы k + 2 6 2m + 2. Таким образом, порядок аппроксимации m -шагового линейного метода не может превышать 2m — для неявного метода. Если же метод явный, то одним неизвестным в системе становится меньше, и максимально возможный порядок аппроксимации будет равен 2m − 1. Перейдем к практическим примерам. 5.4. МЕТОДЫ АДАМСА И ГИРА 5.4 81 Методы Адамса и Гира Определение. Методы Адамса2 — семейство m -шаговых линейных методов решения задачи Коши, в которых берется a0 = 1, a1 = −1, a2 = a3 = . . . = am = 0. Таким образом, общая формула для нахождения приближения yn выглядит так: yn − yn−1 = b0 fn + b1 fn−1 + . . . + bm fn−m . τ (5.10) Посмотрим, к чему приведет требование m-го порядка аппроксимации для данного метода. Ранее было показано, что достаточным условием будет выполнение системы равенств:  m X    bi = 1;     i=0    m  X    ai = 0;   i=0 m X     iai = −1;    i=1    m  X    ij−1 (iai + jbi ) = 0,   j = 2, k. i=1 Второе и третье равенство следуют из определения методов Адамса. Подставив значения ai в последнее условие, получим такую систему:  m X    bi = 1;   i=0 (5.11) m X   j−1   j i bi = 1, j = 2, k.  i=1 Эти уравнения на коэффициенты bi — достаточные условия k-го порядка аппроксимации. Их k штук, а должны они определять m + 1 неизвестное. Очевидно, чтобы система была разрешима, необходимо выполнение неравенства: k 6 m + 1. Таким образом, максимально возможный порядок аппроксимации не может превышать m + 1. Если мы потребуем, чтобы k = m + 1, то система (5.11) даст единственное решение — то есть схема, отвечающая максимально возможному порядку аппроксимации, одна. Взглянув на формулу (5.10), легко заметить, что при b0 = 0 она становится явной относительно yn . Очевидно, что в этом случае система (5.11) содержит меньше неизвестных, и разрешима она будет уже при k 6 m. Если же b0 6= 0, то схема является неявной, и для нахождения yn приходится решать такое, в общем случае нелинейное, уравнение: yn − τ b0 f (tn , yn ) = yn−1 + τ m X i=1 В общем случае оно решается методом Ньютона. 2 Впервые предложены Дж. Адамсом (1855). bm fn−i . 82 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ Пример 5.1. Положим b0 = 0, m = 1. Это означает, что схема для вычисления yn будет явная, поэтому максимально возможный порядок аппроксимации будет равен m, то есть 1. Пусть будет так. От системы (5.11) останется только первое уравнение — на b1 , и из него просто получается, что b1 = 1. Общая схема метода будет такова:   yn − yn−1 = f (tn−1 , yn−1 ), n = 1, 2, . . . τ  y0 = u0 — начальное условие. Пример 5.2. Здесь возьмем b0 = 0, m = 3. Максимально возможный порядок аппроксимации будет равен 3 — его и потребуем. При этом система (5.11) для нахождения bi будет записана так:  23   ; b1 =    12    b1 + b2 + b3 = 1;  4 2(b1 + 2b2 + 3b3 ) = 1; ⇐⇒ b2 = − ;   3  3(b1 + 4b2 + 9b3 ) = 1.      b3 = 5 . 12 Общая схема такова: yn − yn−1 23fn−1 − 16fn−2 + 5fn−3 = , n > 3. τ 12 Значение y0 , как и раньше, берется равным u0 , а y1 и y2 обычно ищутся методами Рунге-Кутта. В данном случае они должны быть не менее, чем третьего порядка точности (так как k у нас равно трем). Пример 5.3. Рассмотрим пример неявного метода Адамса при m = 1. Тогда можно взять k, равное 2 — максимально возможному порядку аппроксимации. Точно так же решаем систему:  b0 + b1 = 1; 2b1 = 1.  1   b0 = ; 2 ⇐⇒   b = 1. 1 2 и получаем такую схему:   yn − yn−1 = fn + fn−1 , τ 2  y0 = u 0 . n>1 В данном случае yn приходится находить из, вообще говоря, нелинейного уравнения. По-другому его можно записать так: τ τ yn − f (tn , yn ) = yn−1 + fn−1 . 2 2 Обычно его решают методом Ньютона, где в качестве начального приближения берут yn−1 . Пример 5.4. Последним примером на метод Адамса будет неявный метод, с m = 2 и k = 3. Получаем для нахождения bi такую систему:  5      b0 = 12 ;    b0 + b1 + b2 = 1;  2 2(b1 + 2b2 ) = 1; ⇐⇒ b1 = ;   3  3(b1 + 4b2 ) = 1.      b2 = − 1 . 12 5.4. МЕТОДЫ АДАМСА И ГИРА 83 Отсюда получаем общую схему:  yn − yn−1 5fn + 8fn−1 − fn−2   = , n>2  τ 12  y = u0 ;   0 y1 ищется методом Рунге-Кутта 3-го порядка точности. yn ищется из неявной формулы по методу Ньютона. Методы Гира Определение. Методами Гира называется семейство линейных m -шаговых методов, в которых заранее определяется b0 = 1, b1 = b2 = . . . = bm = 0. Общая расчетная формула будет такова: a0 yn + a1 yn−1 + . . . + am yn−m = fn . τ Заметим, что все эти методы являются неявными, так как yn приходится находить из нелинейного, вообще говоря, уравнения: m X a0 yn − τ f (tn , yn ) = − ai yn−i . i=1 В системе достаточных условий для k-го порядка аппроксимации  m X    bi = 1;     i=0    m  X    ai = 0;   i=0 m X     iai = −1;    i=1    m  X    ij−1 (iai + jbi ) = 0,   j = 2, k. i=1 первое равенство будет выполняться всегда. Подставив в последнее условие значения из определения методов Гира, получим такую систему:  m X   ai = 0;     i=0    m  X iai = −1;   i=1    m  X    ij ai = 0, j = 2, k.  (5.12) i=1 Она содержит m + 1 неизвестное и состоит из k + 1 уравнения. Таким образом, чтобы обеспечить разрешимость системы и требование k-го порядка аппроксимации, приходится ограничивать k числом m. При k = m порядок аппроксимации будет максимально возможным, а схема, им определяемая, тоже будет одна. 84 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ Пример 5.5. В первом простейшем примере берем m = 1, k = 1. Система для нахождения ai будет такова:   a0 + a1 = 0; a0 = 1; ⇐⇒ a1 = −1. a1 = −1. Отсюда получаем простейшую разностную схему: ( y −y n n−1 = f (tn , yn ), τ y0 = u 0 . n > 1; Пример 5.6. Взяв в этом примере m = 2 и k = 2, мы получим такую систему уравнений:   a0 + a1 + a2 = 0; a1 + 2a2 = −1;  a1 + 4a2 = 0.  3  a0 = ;   2  ⇐⇒ a = −2; 1     a = 1. 2 2 Соответствующая схема будет выглядеть так:  3yn − 4yn−1 + yn−2   = τ f (tn , yn ),  2 n>2  y = u0 ;   0 y1 ищется методом Рунге-Кутта 2-го порядка точности. То есть, для определения yn мы имеем такое нелинейное уравнение: 1 3 yn − τ f (tn yn ) = 2yn−1 − yn−2 . 2 2 Пример 5.7. В этом примере мы возьмем m = 3, k = 3. Тогда последнее условие в системе (5.12) разобьется на два уравнения, и мы получим такую систему:  11   a0 = ;    6   a + a + a + a = 0;   1 2 3     a1 = −3; a1 + 2a2 + 3a3 = −1; 3 ⇐⇒ a1 + 4a2 + 9a3 = 0;   a2 = ;     2  a1 + 8a2 + 27a3 = 0.      a3 = − 1 . 3 Этому будет соответствовать такие расчетные формулы:  11yn − 18yn−1 + 9yn−2 − 2yn−3   = τ f (tn , yn ),  6 n>2  y = u0 ;   0 y1 , y2 ищутся методом Рунге-Кутта 3-го порядка точности. Замечание. На практике используются методы Гира вплоть до десятого порядка аппроксимации (и соответствующей точности). Это связано с тем, что эти методы обладают свойством вычислительной устойчивости, о котором мы поговорим в следующем разделе. 5.5. УСТОЙЧИВОСТЬ ЧИСЛЕННЫХ МЕТОДОВ РЕШЕНИЯ ЗАДАЧИ КОШИ 5.5 85 Устойчивость численных методов решения задачи Коши Мы будем рассматривать численные методы для поиска функции, являющейся решением такой задачи Коши:  du  = f (t, u), t > 0; dt  u t=0 = u0 . Для обоснования дальнейших действий сначала проведем теоретические рассуждения. Будем считать, что функция u является решением задачи Коши: du dt    u t=0 = f (t, u), t > 0; (5.13) = u0 . Пусть функция U — решение аналогичной задачи, но с «возмущенными» начальными данными: dU dt    U t=0 = f (t, U ), t > 0; (5.14) = u0 + u0 . Представим ее следующим образом: U = u + u, где u — функция-погрешность. Исследуем u, подставив представление для U в (5.14) и разложив f (t, u) в ряд Тейлора по второй переменной: ( du du + dt dt u t=0 + u t=0 = f (t, u + u) = f (t, u) + fu0 (t, u)u + Rf , t > 0; = u0 + u0 . Используя (5.13), получим: du dt ( u t=0 = fu0 (t, u)u + Rf , t > 0; (5.15) = u0 . Из курса «Дифференциальные уравнения» известно, что если fu0 (t, u) < 0, то функция u(t) монотонно стремится к нулю на бесконечности. В этом случае говорят, что задача (5.13) устойчива по начальным данным. Теперь вернемся к численным методам. Если их применять к возмущенной задаче (при условии fu0 < 0 ), то естественно требовать, чтобы функция u получалась убывающей (или, что тоже самое, функция U приближалась к u, то есть погрешность не накапливалась). После проведенных выкладок понятно, что это то же самое, что требовать получение убывающей функции при решении системы (5.15). Поэтому методы обычно тестируют на такой модельной задаче: ( u du dt t=0 = λu, t > 0; λ = const < 0; = u0 , справедливо ожидая, что если при ее решении будет получаться убывающая на бесконечности функция, то метод, примененный к системе (5.14), даст функцию, не сильно отличающуюся от u. Такие методы называются устойчивыми. Примечание. Мы требуем убывания от функции, если u0 , задающее начальное условие, положительно. Если же оно отрицательно, то логично требовать возрастание, или, точнее, стремление к нулю слева. 86 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ Будем называть метод условно устойчивым, если он устойчив при некоторых значениях своих параметров. Теперь рассмотрим несколько связанных с этой модельной задачей примеров. Пример 5.8. Явный одношаговый метод Адамса первого порядка аппроксимации. Разностная схема имеет вид: yn+1 − yn = f (tn , yn ). τ Правая часть, согласно дифференциальному уравнению, выражается как f (tn , yn ) = λyn , тогда схема будет иметь следующий вид: yn+1 = (1 + τ λ)yn . Как уже говорилось, метод называется устойчивым, если сеточная функция yn не возрастает по n. В нашем примере |yn+1 | = |1 + τ λ| · |yn | не будет возрастать, если |1 + τ λ| 6 1, то есть при −2 6 τ λ 6 0. Правая часть неравенства выполняется всегда, так как τ > 0, а λ < 0. Из этого следует, что явный одношаговый метод Адамса удовлетворяет условию устойчивости лишь при τ 6− 2 2 = . λ |λ| Этот метод условно устойчивый, и это не очень хорошо, так как, чтобы не набирать погрешность, надо учитывать ограничение на τ, где величина λ зависит от поведения функции f (t, u) (λ — разностный аналог производной функции f (t, u)). То есть надо выбирать очень мелкий шаг интегрирования в соответствии с поведением функции f. Пример 5.9. Одношаговый метод Гира первого порядка аппроксимации. yn+1 − yn = f (tn+1 , yn+1 ). τ Согласно модельной задаче, f (tn+1 , yn+1 ) = λyn+1 , откуда yn+1 (1 − τ λ) = yn . Для устойчивости метода мы требуем, чтобы |yn+1 | 6 |yn |. Это неравенство выполнено всегда, так как τ > 0 и λ < 0, то есть сеточная функция не возрастает при любых шагах интегрирования. Мы получили, что данный метод Гира является абсолютно устойчивым. Неоспоримым достоинством этого метода по сравнению с методами Адамса является то, что шаг интегрирования мы выбираем с оглядкой лишь на требуемую точность, а к недостаткам можно отнести то, что на каждом шаге приходится решать нелинейную систему уравнений. При рассмотрении модельной задачи у нас возникают так называемые сеточные уравнения. В общем виде их можно записать так: yn+1 = an yn + an−1 yn−1 + . . . + a0 y0 . В правой части стоит линейная комбинация значений сеточной функции. По сути дела, это некоторое рекуррентное соотношение. В общем виде оно решается с помощью характеристического многочлена, но мы ограничимся случаем, когда решение можно задать в виде yi = q i . На основе такого представления решений аппарат исследования на устойчивость таков. Если удается показать, что |q| 6 1 для любых параметров метода (от них будет зависеть рекуррентное соотношение), то общее решение будет монотонно убывать. Если же существует хотя бы один набор параметров такой, что |q| > 1, то условие устойчивости будет нарушено для данного метода. 5.5. УСТОЙЧИВОСТЬ ЧИСЛЕННЫХ МЕТОДОВ РЕШЕНИЯ ЗАДАЧИ КОШИ 87 Рассмотрим пример. Проверим на устойчивость явные одношаговые методы Адамса. yn+1 − yn = f (tn , yn ). τ Решение сеточного уравнения будем искать в виде yn = q n . Так как схема сводится к виду yn+1 = (1 + τ λ)yn , то q будет равным 1 + τ λ. Для устойчивости метода необходимо, чтобы |1 + τ λ| 6 1, но это условие мы уже получали выше. Таким образом, пока наш метод проверки на устойчивость ничего нового не дал, но мы увидим в полученных формулах закономерности. Если применить этот аппарат к методам Гира, опять получаем знакомую формулу: |q| = 1 . |1 − τ λ| В общем случае для сеточных уравнений, в которых задано более двух узлов, основание решения q может быть, вообще говоря, комплексным числом. Таким образом, проведем некоторое обобщение аппарата исследования устойчивости — будем считать, что в модельной задаче λ ∈ C. Еще раз применим метод Адамса к решению модельной задачи, но уже с комплексным λ. Обозначив в получившихся выкладках τ λ = µ ∈ C, получим: |1 + µ| 6 1. =µ −1 µ <µ Для всех точек, лежащих внутри области, соответствующие методы Адамса устойчивы. Аналогично для методов Гира: |1 − µ| > 1. Как нетрудно заметить, метод устойчив во внешней области. Определение. Назовем областью устойчивости численного метода решения задачи Коши для ОДУ ту область значений µ, в каждой точке которой |q| 6 1. Определение. Численный метод называется А-устойчивым3 , если его область устойчивости содержит отрицательную полуплоскость µ. В соответствии с этим определением, метод Гира — А-устойчив, а метод Адамса не является Аустойчивым. Переход с оси на комплексную полуплоскость привел к тому, что можно показать, что явных Аустойчивых методов не существует, а среди неявных нет А-устойчивых методов выше второго порядка аппроксимации. 3 от абсолютно устойчивый. 88 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ =µ µ <µ 1 Все это вытекает из «плохого» определения устойчивости. Попробуем немного исправить эту ситуацию. Определение. Линейный многошаговый метод называется А( α )-устойчивым, если область его устойчивости содержит угол |arg(−µ)| < α. =µ µ α α <µ −µ = rei(π−ϕ) =⇒ Заметим, что µ = reiϕ =⇒ arg(−µ) = π − ϕ. Можно доказать, что для введенного определения устойчивости выполняется следующее утверждение. Утверждение 5.1. Среди явных линейных m-шаговых методов нет А(α )-устойчивых. Это следует учитывать при решении химических, биологических задач, где часто решение имеет экспоненциальный вид. Среди неявных методов существуют А(α )-устойчивые методы. Рассмотрим пару примеров таких методов. Пример 5.10. Двухшаговый метод Гира второго порядка аппроксимации. 1 3 yn − 2yn−1 + yn−2 = τ f (tn , yn ). 2 2 Пример 5.11. Четырехшаговый метод Гира четвертого порядка аппроксимации. 1 (25yn − 48yn−1 + 36yn−2 − 16yn−3 + 3yn−4 ) = τ f (tn , yn ). 12 Найдем область устойчивости для рассмотренного выше двухшагового метода Гира второго порядка аппроксимации. Как уже делали ранее, зададим f (tn , yn ) = λyn и, обозначим, τ λ = µ, тогда 5.5. УСТОЙЧИВОСТЬ ЧИСЛЕННЫХ МЕТОДОВ РЕШЕНИЯ ЗАДАЧИ КОШИ 89 сеточное уравнение будет выглядеть так: 3 1 yn − 2yn−1 + yn−2 = µyn . 2 2 (5.16) Как мы уже договаривались, нас интересуют только решения вида yn = q n . Подставив это значение в (5.16), получим: 1 3 2 q − 2q + = µq 2 . 2 2 или: 1 3 2 µ = − + 2. 2 q 2q Найдем те µ, для которых |q| < 1. При этом комплексная плоскость разбивается на две области — ту, где метод устойчив, и — где неустойчив. При этом, так как на границе |q| = 1, то q = eiϕ и уравнение границы примет вид: 3 1 3 1 1 − 2e−iϕ + e−i2ϕ = − 2 cos ϕ + 2i sin ϕ + cos 2ϕ − i sin 2ϕ = 2 2 2 2 2 3 1 = − 2 cos ϕ + (2 cos2 ϕ − 1) + i(2 sin ϕ − sin ϕ cos ϕ). 2 2 µ(ϕ) = Сделаем замену переменной x = cos ϕ (x ∈ [−1; 1]), получим: p µ = 1 − 2x + x2 ± i(2 − x) 1 − x2 . =µ |q| < 1 |q| > 1 4 <µ |q| = 1 Вообще то неплохо было бы отметить здесь еще область определения, а то, например, при µ = 23 итерационный процесс неопределен. Проверить область определения мы √ предоставим читателю. Под1±i 1/2 ставим точку µ = − 32 , она лежит во внешней части плоскости, а q = , |q| < 1 и, таким образом, 3 мы получили, что внешняя область является областью устойчивости данного метода Гира. Хотя при реализации неявных методов требуется больше времени на просчет одного шага, но их плюс в том, что мы можем выбирать произвольный шаг — тот, который нам нужен. Численные методы решения задачи Коши для систем ОДУ Исследуем методы решения задачи Коши для систем ОДУ. Будем искать приближенное решение для такой системы: ( ut = f (t, u), t > 0; u|t=0 = u0 , где u = (u1 , u2 , . . . , un ) и f = (f1 , f2 , . . . , fn ). 90 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ Методы Рунге-Кутта и линейные одношаговые Чтобы стало яснее, приведем пример задачи:  ut =    vt = u(0) =    v(0) = методы легко переносятся на системы уравнений. f (t, u, v); g(t, u, v); u0 ; v0 . Покажем, как к этой задаче применяется метод Рунге-Кутта четвертого порядка аппроксимации (и точности соответственно). Обозначим u(tn ) = yn , а v(tn ) = zn , тогда формулы для подсчета следующих сеточных значений будут таковы: 1 yn+1 − yn = (K1 + 2K2 + 2K3 + K4 ); τ 6   zn+1 − zn = 1 (M + 2M + 2M + M ). 1 2 3 4 τ 6    Параметры Ki , Mi , i = 1, 2, 3, 4 вычисляются по следующим формулам: K1 K2 K3 K4 = f (tn , yn , zn ), M1 = g(tn , yn , zn ); = f (tn + τ2 , yn + τ2 K1 , zn + τ2 M1 ), M2 = g(tn + τ2 , yn + τ2 K1 , zn + τ2 M1 ); = f (tn + τ2 , yn + τ2 K2 , zn + τ2 M2 ), M3 = g(tn + τ2 , yn + τ2 K2 , zn + τ2 M2 ); = f (tn + τ, yn + τ K3 , zn + τ M3 ), M4 = g(tn + τ, yn + τ K3 , zn + τ M3 ). Мы долго и упорно говорим, что условно устойчивые методы хуже устойчивых, но не привели не одного примера, объясняющего разницу между ними. Исправим ситуацию. Представим, что мы моделируем два процесса, причем один из них протекает существенно быстрее другого. Система уравнений для такой задачи будет иметь вид:  0 u + a1 u1 = 0, a1 > 0;    10 u2 + a2 u2 = 0, a2 > 0; u1 (0) = u0 ;    u2 (0) = u0 , а точные решения соответственно:  u1 u2 = u0 e−a1 t ; = u0 e−a2 t . Как мы уже говорили, a1  a2 — константы суть характерное время протекания первого и второго процессов. Нас интересует описание системы при достаточно больших t > t∗ . Пусть мы используем для решения этой задачи условно устойчивый метод (например, какой-нибудь из методов Адамса). При этом на шаг τ мы должны сделать ограничение τ λ > −2, где роль параметра λ играют a1 и a2 , то есть λ = −a1 или λ = −a2 . Неравенства должны выполняться одновременно, то есть τ 6 a21 (в силу того, что a1  a2 ). Теперь представим, что нам важно исследовать первый процесс в то время, когда второй уже не влияет на систему. Тем не менее, от поведения второго процесса зависит шаг аппроксимации (он будет уменьшаться), поэтому весьма вероятно, что мы будем выполнять «лишнюю» работу. Отсюда можно сделать вывод, что в случае решения систем ОДУ, описывающих разномасштабные процессы, надо использовать абсолютно устойчивые методы (например, методы Гира). Опишем еще одно свойство систем уравнений. В качестве примера системы возьмем такую: ∂u ∂t    u t=0 = f (t, u), = u0 . 0 < t 6 T; (5.17) 5.6. ИНТЕГРО-ИНТЕРПОЛЯЦИОННЫЙ МЕТОД ПОСТРОЕНИЯ РАЗНОСТНЫХ СХЕМ 91 u u0 | t∗ u2 (t) u1 (t) t Поставим ей в соответствие матрицу из производных (якобиан):   ∂fj (t, u) A(t, u) = ∂ui и обозначим ее собственные значения за λk (t), t ∈ [0; T ]. Определение. Система (5.17) называется жесткой, если выполнены два условия: Re λk (t) < 0 ∀ k, t ∈ [0; T ]; max |Re λk (t)| sup t∈[0; T ] k min |Re λk (t)|  1 — много больше единицы. k Такие системы чаще всего решают неявным абсолютно устойчивым методом. 5.6 Интегро-интерполяционный метод построения разностных схем Перейдем к рассмотрению более сложных краевых задач. Для начала исследуем применение численных методов для решения такой модельной задачи:   (k(x)u0 (x))0 − q(x)u(x) + f (x) = 0, 0 < x < l; −k(0)u0 (0) + βu(0) = µ1 ; (5.18)  u(l) = µ2 . Это краевая задача для обыкновенного дифференциального уравнения второго порядка с заданными функциями k, q, f и неизвестной функцией u. Известно, что если выполняются условия   k(x) > k0 > 0; q(x) > 0;  β > 0, то решение задачи (5.18) существует и единственно. Задача (5.18) содержит уравнение параболического типа. Обычно такие уравнения возникают при исследовании распределения температуры в тонком стержне или в диффузионных процессах. Решение системы типа (5.18) проходит в несколько этапов. Сначала ей сопоставляется дискретная модель, а на ее основе строится разностная схема. Существует несколько методов построения таких схем, и первым мы рассмотрим интегро-интерполяционный метод. Название его происходит от того, что в процессе построения соответствующей разностной схемы мы переходим от интегральных соотношений к интерполяционным уравнениям. Построение разностной схемы 92 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ Перейдем к построению схемы. Для начала введем на отрезке [0; l] равномерную сетку: l } N ωh = {xi = ih, i = 0, N , h = — очевидно, x0 = 0, xN = l. Теперь введем такое обозначение для средних точек между узлами сетки: xi± 21 = xi ± h . 2 Также обозначим ui = u(xi ) — значение искомой функции в узлах сетки, и W (x) = k(x)u0 (x). Применяя эти обозначения, фиксируем произвольное i ∈ [1; N −1] и проинтегрируем первое уравнение системы (5.18) по отрезку [xi− 12 ; xi+ 12 ] : xi+ 1 2 Z ((k(x)u0 (x))0 − q(x)u(x) + f (x)) dx = 0 ⇐⇒ xi− 1 2 xi+ 1 xi+ 1 Z Z 2 Wi+ 21 − Wi− 12 − 2 q(x)u(x) dx + xi− 1 f (x) dx = 0, (5.19) xi− 1 2 2 где Wi± 21 = W (xi± 12 ). xi+ 1 2 Z Первый интеграл можно приблизить значением ui · q(x) dx. Тогда (5.19) можно переписать как xi− 1 2 приближенное равенство: xi+ 1 xi+ 1 Z Z 2 Wi+ 21 − Wi− 12 − ui · 2 f (x) dx ≈ 0. q(x) dx + xi− 1 (5.20) xi− 1 2 2 Перейдем от интегральных выражений к линейным. Для этого введем новые обозначения: xi+ 1 ϕi = 1 h xi+ 1 2 Z f (x) dx; di = xi− 1 q(x) dx. (5.21) xi− 1 2 Теперь заметим, что u0 (x) = 2 Z 1 h 2 W (x) . Проинтегрировав это равенство на отрезке [xi ; xi+1 ], получим: k(x) x Zi+1 ui+1 − ui = W (x) dx. k(x) xi Заменим это равенство приближенным: x Zi+1 ui+1 − ui ≈ Wi+ 21 xi dx . k(x) 5.6. ИНТЕГРО-ИНТЕРПОЛЯЦИОННЫЙ МЕТОД ПОСТРОЕНИЯ РАЗНОСТНЫХ СХЕМ 93 Отсюда следует, что если обозначить  ai+1 1 = h x Zi+1 −1 dx  , k(x) xi то будут справедливы приближенные равенства:  ui+1 − ui  ;  Wi+ 21 ≈ ai+1 h   W 1 ≈ a ui − ui−1 . i i− 2 h (5.22) Воспользовавшись обозначениями (5.21) и (5.22), приближенное равенство (5.20) можно переписать так: ui − ui−1 ui+1 − ui − ai − hdi ui + ϕi h ≈ 0. (5.23) ai+1 h h Обозначим за yi такие числа, которые при подстановке в (5.23) вместо ui дают точное равенство: ai+1 yi+1 − yi yi − yi−1 − ai − hdi yi + ϕi h = 0. h h (5.24) Найденные из таких уравнений значения yi и будут считаться приближениями к ui . Полученное равенство и будет искомой разностной схемой, однако ее можно переписать и в более компактном виде, заметив, что первые две дроби — не что иное, как разностные производные назад. Обозначив их yx,i+1 и yx,i , перепишем (5.24) в следующем виде: ai+1 yx,i+1 − ai yx,i − hdi yi + ϕi h = 0. Разделив это равенство на h, мы сможем объединить две разностные производные во вторую разностную производную вперед: (ayx )i+1 − (ayx )i ≡ (ayx )x,i . h В итоге мы получим такой вид разностной схемы: (ayx )x,i − di yi + ϕi = 0. (5.25) Мы имеем право написать такие равенства для i = 1, N − 1. Их можно объединить в систему линейных (по построению) уравнений относительно yi . Она будет содержать N − 1 уравнение и N + 1 неизвестное. Необходимые для однозначной разрешимости системы 2 уравнения добавим из краевых условий. Мы можем заменить в последнем равенстве в (5.18) u(l) на yN , тогда получим, что yN = µ2 . (5.26) Для получения последнего уравнения мы выполним те же самые действия, что и при выводе равенства (5.19), но интегрирование будем проводить на отрезке [0; h2 ]. Тогда можно получить такое равенство: h h Z2 Z2 W 12 − W0 − u0 q(x) dx + f (x) dx ≈ 0. (5.27) W 1 2 и W0 мы найдем, заменив в приближенных равенствах u1 − u0 ; h W0 ≈ βu0 − µ1 . W 12 ≈ a1 94 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ uk на yk и получив уравнения: y1 − y0 ; h W0 = βy0 − µ1 . W 21 = a1 После этого, воспользовавшись обозначениями h ϕ0 = 1 h 2 h Z2 f (x) dx, d0 = 1 h 2 Z2 q(x) dx, мы приведем (5.27) к такому виду: h h a1 yx,0 − βy0 + µ1 − d0 y0 + ϕ0 = 0 ⇐⇒ 2 2 a1 yx,0 − βy0 = µ, d0 h2 , (5.28) ϕ0 h2 . µ = µ1 + где β = β + Уравнения (5.25), (5.26) и (5.28) представляют собой окончательный вариант разностной схемы, полученной с использованием интегро-интерполяционного метода. Решение разностной схемы Вторым шагом на пути решения краевой задачи численно становится выбор метода решения построенной схемы. В нашем случае мы получили систему линейных уравнений, для которой метод решения и выбирается. Заметим, что (5.25) можно переписать так:   1 yi+1 − yi yi − yi−1 ai+1 − ai − di yi + ϕi = 0. h h h Собрав коэффициенты при yi , получим: Ai yi+1 − Ci yi + Bi yi−1 = −Fi , i = 1, N − 1, где Ai = ai+1 , Bi = ai , Ci = ai + ai+1 + di h2 . Добавив к этим уравнениям уравнения (5.26) и (5.28), получим систему из N + 1 уравнения. Матрица, задающая эту систему уравнений, будет являться трехдиагональной, а такие системы обычно решаются методом прогонки. Он применим, так как выполнены условия: Ai , Bi > 0, Ci > Bi + Ai — они дают существование и единственность yi , отвечающих уравнениям (5.25), (5.26) и (5.28). 5.7 Метод аппроксимации квадратичного функционала Это другой метод построения разностных схем. Будем рассматривать его для задачи, схожей с (5.18), но с более простыми краевыми условиями:  (k(x)u0 (x))0 − q(x)u(x) + f (x) = 0, 0 < x < 1; u(0) = u(1) = 0. Известно, что решение такой задачи эквивалентно поиску u, минимизирующих функционал Z1 J[u] = [k(x)(u0 (x))2 + q(x)u2 (x) − 2f (x)u(x)] dx. 5.8. КОРРЕКТНОСТЬ РАЗНОСТНОЙ СХЕМЫ 95 Задав на отрезке [0; 1] равномерную сетку, разобьем интеграл по всему отрезку на сумму по подотрезкам: N Zxi X   J[u] = k(x)(u0 (x))2 + q(x)u2 (x) − 2f (x)u(x) dx. i=1x i−1 Для упрощения поиска экстремума заменим обыкновенную производную на ее разностный аналог: N Zxi X   2 J[u] ≈ k(x)ux,i + q(x)u2 (x) − 2f (x)u(x) dx. i=1x i−1 1 Воспользовавшись обозначением ai = h Zxi k(x) dx, получим: xi+1  Zxi N  X
2 ai ux,i h + J[u] ≈ q(x)u2 (x) − 2f (x)u(x) dx . i=1 xi−1 Интегралы посчитаем по формуле трапеций, заменив всюду ui на yi — приближенные значения: J[u] ≈ J[y] = Jh (y0 , y1 , . . . , yN ) = = {y0 = yN = 0} = N X 2 h+ ai yx,i  N  X h 2 2 ai yx,i h + (qi yi2 − 2fi yi + qi−1 yi−1 − 2fi−1 yi−1 ) = 2 i=1 N −1 X (qi yi2 − 2fi yi )h. i=1 i=1 Мы свели задачу о поиске элемента, минимизирующего функционал, к поиску чисел yi , доставляющих минимум функции многих переменных — при этом, правда, мы потеряли в точности. Необходимым условием экстремума будет равенство нулю всех частных производных: ∂Jh = 0 ⇐⇒ ∂yi ⇐⇒ ⇐⇒ 1 1 2ai+1 yx,i+1 (− · h) + 2ai yx,i ( · h) + (2qi yi − 2fi )h = 0 ⇐⇒ h h ai+1 yx,i+1 − ai yx,i − qi yi + fi = 0 h (ayx )x,i − qi yi + fi = 0. Последнее уравнение — это уже часть итоговой разностной схемы (осталось добавить только краевые условия). Можно заметить, что она похожа на схему, возникающую в интегро-интерполяционном методе, однако коэффициенты различны, да и свойства схем тоже довольно сильно различаются. 5.8 Корректность разностной схемы Напомним несколько определений. Определение. Разностная схема аппроксимирует исходную дифференциальную задачу в точке xi [на всей сетке], если погрешность аппроксимации в этой точке [соответственно, норма погрешности аппроксимации] стремится к нулю [соответственно, тоже к нулю] с уменьшением h : h i h→0 h→0 ψ(xi ) −−−→ 0 ||ψ||h = ||ψ||C(ωn ) = max |ψi | −→ 0 . i 96 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ Определение. Разностная схема аппроксимирует исходную дифференциальную задачу с p-м порядком аппроксимации, если ψi = O(hp ) в точках xi , или в целом на сетке, если ||ψi ||h = O(hp ). Будем поступать так же, как и в случае выяснения порядка аппроксимации для задачи Коши в разностной задаче аппроксимации. Можно показать, что в целом на сетке схема, построенная интегро-интерполяционным методом, будет иметь второй порядок аппроксимации. Вычисления громоздки и мы их опускаем, заметив, что невязка имеет порядок не хуже O(h2 ) даже в крайних узлах: ψ0 = O(h2 ), ψN = 0. Если мы требуем такой порядок аппроксимации, то можно сэкономить на вычислении параметров, вычисляя ai , di , ϕi по формуле прямоугольников, при этом получим следующие значения параметров:      ai = k xi− 12 ; di = q(xi );   ϕi = f (xi ). Если считать параметры по квадратурной формуле трапеций, то получим следующие выражения:  1 1 1  + ; =    a k(x ) k(x i i−1 i)   


1 di = q xi− 12 + q xi+ 12 ;  2       ϕi = 1 f x 1
+ f x 1

. i− 2 i+ 2 2 Рассмотрим вопрос о сходимости приближенного решения к точному. Как обычно, обозначим zi = yi − ui — погрешность и напомним несколько определений. h→0 Определение. Приближенное решение yi сходится к точному в точке xi , если zi −→ 0. h→0 Определение. Приближенное решение yi сходится к точному на всей сетке, если ||zi ||h −→ 0. Определение. Если величина погрешности zi в каждой точке (или на всей сетке) есть O(hp ), то метод имеет p-й порядок точности. Можно установить, что сеточная норма ||z||C(ωi ) = O(h2 ). Доказательство этого утверждения можно посмотреть в книге [1]. Подставим в расчетную схему yi = ui + zi : (azx )x,i − di zi = −(aux )x,i + di ui − ϕi . Аналогичную операцию проведем для граничных условий:  (azx )x,i − di zi = −ψi ;    −a1 zx,0 + βz0 = −ψ;    zN = 0. Как нетрудно заметить, задача для погрешности имеет ту же структуру, что и исходная разностная схема, с заменой правой части на невязку. После преобразований системы, которые мы снова опускаем, можно получить, что ||z||h 6 M1 ||ψ||h . 5.9. ЯВНАЯ РАЗНОСТНАЯ СХЕМА ДЛЯ УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ 97 Разностная схема имеет 2-й порядок аппроксимации (||ψ||h = O(h2 )), а, следовательно, и второй порядок точности ||z||h = O(h2 ). Подробные указания на то, как это получить, можно найти в [1]. Проделав те же самые действия, можем получить оценку на приближенное решение в сеточной норме (по аналогии — уравнения очень похожи). ||y||h 6 M1 (||ϕ||h + |µ1 | + |µ2 |), (5.29) — такой оценки и следовало ожидать. Теперь мы можем дать несколько определений. Определение. Если для решения разностной задачи выполняется оценка (5.29), то решение называется устойчивым по правой части. Это определение устойчивости разностной задачи является непосредственным следствием общего определения устойчивости. Определение. Задача называется корректно поставленной по Адамару4 , если: 1) решение существует; 2) решение единственно; 3) решение непрерывно зависит от входных данных (устойчиво по правой части). По аналогии запишем определение для разностной схемы. Определение. Разностная схема называется корректной, если: 1) решение существует; 2) решение единственно; 3) решение устойчиво по правой части. Теорема 5.2. Пусть дифференциальная задача корректно поставлена, и разностная схема, соответствующая этой дифференциальной задаче, корректна. Тогда, если разностная схема аппроксимирует исходную задачу, то решение разностной схемы сходится к решению исходной дифференциальной задачи и порядок аппроксимации совпадает с порядком точности. Доказательство. Как уже записывали ранее, оценка на приближенное решение такова ||y||h 6 M1 ||ϕ||h . В силу линейности разностной схемы, оценка на погрешность будет: h→0 ||z||h 6 M1 ||ψ||h −→ 0. Таким образом, так как ||ψ||h = O(hp ), решение разностной схемы сходится и имеет p-й порядок точности. Далее работать с теоремой мы будем по следующему плану. Сначала исследуем разностную схему на аппроксимацию (аппроксимирует ли она исходное ОДУ). Затем проверяем схему на устойчивость, и потом уже можно пользоваться теоремой, что мы и будем с успехом делать. 5.9 Явная разностная схема для уравнения теплопроводности Мы будем рассматривать достаточно простые задачи, в которых решение можно построить и аналитическими методами, но наша задача — изучать численные методы. Рассмотренная ниже техника решения краевых задач легко обобщается на более сложные случаи, не имеющие аналитического решения. 4 Понятие впервые предложено Ж. Адамаром (1923). 98 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ Запишем краевую задачу для  ut (x, t)    u(x, 0) u(0, t)    u(1, t) уравнения теплопроводности (УТ): = = = = uxx (x, t) + f (x, t), u0 (x), µ1 (t), µ2 (t), 0 < x < 1, 0 < t 6 T ; 0 6 x 6 1; 0 6 t 6 T; 0 6 t 6 T. Проведем дискретизацию области изменения независимого переменного:  ωh = xj = jh, j = 0, N , h = N1 ;  T . ωτ = tn = nτ, n = 0, K, τ = K Параметры N и K характеризуют «густоту» сетки. Теперь построим семейство линий j = 0, N ; n = 0, K. x = xj , t = tn , Будем рассматривать точки пересечения этих линий. Разделим узлы на две группы — граничные узлы (в них заданы дополнительные условия — краевые и граничные условия) и внутренние узлы. t T × žŸ ¡ ¢£ ¤¥ ¦§ ¨© ª« ¬­ ®¯ °± ²³ ´µ ‚ƒ „ †‡ ˆ‰ Š‹ Œ Ž ‘ ’“ ”• –— ˜™ š› hi jk lm no pq rs tu vw xy z{ |} ~ € NO PQ RS TU VW XY Z[ \] ^_ `a bc de fg 45 67 89 :; <= >? @A BC DE FG HI JK LM    ! "# $% &' () *+ ,- ./ 01 23
          × τ × × × × × × œ × × × × × × × × × × × h × × × × × × × × × × 1 x После дискретизации строим некоторый аналог исходного уравнения — разностную схему. Сначала введем обозначения: u(xj , tn ) = unj ; f (xj , tn ) = ϕnj , а для частных производных возьмем такие приближения: unj+1 − 2unj + unj−1 ; j , tn ) h2 un+1 − unj j ut |(x , t ) ≈ . j n τ — эти конструкции возникают при применении интегро-интерполяционного метода с формулой прямоугольников. Подставив эти формулы в краевую задачу, получим ее алгебраический аналог:  n+1 n n yj − yjn yj+1 − 2yjn + yj−1    = + ϕnj , j = 1, N − 1, n = 0, K − 1;   τ h2    yj0 = u0 (xj ), j = 0, N ; (5.30)   n  y0 = µ1 (tn ), n = 0, K;      n yN = µ2 (tn ), n = 0, K. uxx |(x ≈ 5.9. ЯВНАЯ РАЗНОСТНАЯ СХЕМА ДЛЯ УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ 99 Будем исследовать класс схем для решения задачи (5.30). Нас будут интересовать следующие вопросы: 1) существование и единственность решения; 2) методы получения решения разностной схемы; 3) как соотносятся разностная схема и исходная дифференциальная задача (т. е. аппроксимация); 4) есть ли сходимость приближенного решения к точному. Будем рассматривать (пытаться разрешить) уравнение относительно yjn+1 . Множество узлов дискретной сетки с одинаковым t = const назовем временным слоем. Первое уравнение из (5.30) можно переписать так:  
2τ τ n n yjn+1 = 1 − 2 yjn + 2 yj+1 + yj−1 + τ ϕnj , j = 1, N − 1, n = 0, K − 1. h h При n = 0 получим: yj1 =  1− 2τ h2  yj0 +
τ yj+1 + yj−1 + τ ϕ0j , 2 h j = 1, N − 1. В правой части все значения известны из краевых условий. Поэтому мы можем получить искомую сеточную функцию на всем первом временном слое. Аналогично можно рассчитать второй и последующие слои. Таким образом, данная разностная схема решается по слоям, и понятно, что решение существует и единственно. Так как все формулы явные, то вся разностная схема является явной разностной схемой. Эта схема построена по четырем узлам дискретной сетки, что было необходимо для представления первой производной по t и второй по x. Четырех узлов оказалось достаточно, и схема получилась простой. Определение. Совокупность узлов дискретной сетки на базе которых построено разностное уравнение, называется шаблоном разностной схемы. Далее мы будем изучать более сложные схемы с большим числом узлов в шаблоне, а пока вернемся к исследованию только что полученного разностного уравнения. Изучим поведение невязки: ψjn = − un+1 − unj unj+1 − 2unj + unj−1 j + + ϕnj . τ h2 Для этого разложим unj±1 и un+1 в ряд Тейлора в точке unj : j 2 3 unj±1 = unj ± unx,j h + unxx,j h2 ± unxxx,j h6 + O(h4 ); un+1 j = unj + unt,j τ + O(τ 2 ). Тогда получим выражение для невязки: ψjn = (−unt,j + unxx,j + fjn ) + ϕnj − fjn + O(τ + h2 ). 100 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ Выражение в скобках равно нулю в силу того, что unj — точное решение уравнения теплопроводности. Если взять функции ϕnj так, что ϕnj = fjn + O(τ + h2 ) (мы допускаем некоторый произвол при выборе ϕnj ), то для невязки будет справедлива оценка: ||ψ||h = O(τ + h2 ) — то есть наша разностная схема аппроксимирует исходное ДУ со вторым порядком аппроксимации по h и с первым порядком аппроксимации по τ. Рассмотрим теперь вопрос о сходимости приближенного решения к точному. Как обычно, выразим искомое решение yjn через точное решение и погрешность: yjn = unj + zjn . Тогда разностная схема примет вид:  n+1 n n zj − zjn zj+1 − 2zjn + zj−1    + ψjn , j = 1, N − 1, n = 1, K − 1; =   τ h2    zj0 = 0, j = 0, N ;    z0n = 0, n = 0, K;      n zN = 0, n = 0, K. Выразим из первого уравнения погрешность на (n + 1)-м временном слое: zjn+1 = (1 − 2τ n τ n n )zj + 2 (zj+1 + zj−1 ) + τ ψjn . 2 h h Получим оценку на погрешность: |zjn+1 | 6 1 − 2 τ τ |z n | + 2 h2 j h n n zj+1 + zj−1
+ τ ψjn 6  1−2 τ  τ + 2 max zjn + τ max ψjn j j h2 h2 — это неравенство выполняется для любого j, поэтому:  τ  τ max zjn+1 6 1 − 2 2 + 2 2 max zjn + τ max ψjn . j j j h h Норма погрешности на n-м временном слое считается так: zjn C(ω ) = max zjn . Тогда на (n+1)-м n j временном слое норма погрешности оценивается как  τ τ  zjn+1 C(ω ) 6 1 − 2 2 + 2 2 zjn C(ω ) + τ ψjn C(ω ) . n n n+1 h h τ Если предположить, что 1 − 2 2 > 0 (накладываем ограничения на шаг), то h zjn+1 C(ωn+1 ) 6 zjn C(ωn ) + τ ψjn C(ωn ) , причем эта оценка выполняется для любого n. Применив оценку рекурсивно n раз, получим: zjn C(ωn ) 6 zj0 C(ω0 +τ ) n−1 X ψjk C(ωk ) k=0 но zj0 C(ω0 ) = 0 согласно постановке задачи, поэтому: zjn C(ωn ) 6τ n−1 X k=0 ψjk C(ωk ) . , 5.9. ЯВНАЯ РАЗНОСТНАЯ СХЕМА ДЛЯ УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ 101 Напомним, что у нас получена следующая оценка на невязку: ψjn = O(τ + h2 ). Обозначим ψjk C(ω ) = Mk (τ + h2 ), max Mk = M , тогда получим: k k=0, n−1 zjn C(ω ) n 2 6 τ (τ + h ) n−1 X Mk 6 τ nM (τ + h2 ). k=0 Напомним, что τ n = tn и не зависит от τ, тогда обозначим τ nM = M — не зависящая от τ константа, и получим: zjn C(ω ) 6 M (τ + h2 ) n — то есть полученная разностная схема имеет первый порядок точности по τ и второй — по h. Структура разностной схемы и задачи для погрешности одинакова, разница только в том, что в разностной схеме правая часть равна f, а в задаче для погрешности — ψ. Значит, по аналогии (то есть n−1 X ϕkj C(ω ) . Это дискретный анапроведя те же действия) получим, что yjn C(ω ) 6 yi0 C(ω0 ) + τ k n k=0 лог принципа максимума для уравнения теплопроводности (подробнее он описан в курсе «Уравнения математической физики»). Он говорит о том, что решение краевой задачи для уравнения теплопроводности устойчиво по начальным данным и по правой части. Но вспомним, что эта оценка получена при ограничении шагов дискретизации hτ2 6 12 , то есть разностная схема, которую мы исследуем, скорее всего, является условно устойчивой. Далее мы рассмотрим некоторый метод исследования разностной схемы на устойчивость, называемый методом гармоник. Сопоставим рассматриваемой разностной схеме однородное разностное уравнение: n n yjn+1 − yjn yj+1 − 2yjn + yj−1 = . τ h2 (5.31) Исследуем полученное уравнение на решения вида: yjn = q n eijhϕ , (5.32) где q, ϕ — некоторые параметры. Подставим такое yjn в (5.31) и сократим: q−1 τ q−1 τ q−1 τ q
1 ihϕ e − 2 + e−ihϕ =⇒ 2 h 1 = (2 cos(hϕ) − 2) =⇒ h2 2 = − 2 · 2 sin2 hϕ =⇒ 2 h τ = 1 − 4 2 sin2 hϕ 2 . h = Из формулы (5.32) нетрудно заметить, что если |q| 6 1, то ограниченность начального условия влечет ограниченность в любой момент времени n. То есть сеточная функция будет устойчива. Если же |q| > 1 при каких-то τ и h, то решения разностной схемы yjn будут расти с ростом n. Отсюда следует необходимое условие сходимости разностной схемы — |q| 6 1. В нашем случае оно перепишется так: 1−4 τ hϕ sin2 6 1. h2 2 102 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ |q| > 1 |q| < 1 |q| = 1 Раскрыв модуль, получим: τ hϕ sin2 6 1. 2 h 2 Правое неравенство выполнено всегда, перепишем второе неравенство: −1 6 1 − 4 hϕ 1 τ sin2 6 , h2 2 2 или 1 τ 6 h2 2 sin2 hϕ 2 . Оно должно выполняться для всех ϕ. Взяв минимум по правой части, получим окончательное ограничение на параметры схемы: τ 1 6 . h2 2 1 , тогда τ должно быть меньше 12 10−4 . Это не очень хорошо. Для примера рассмотрим h = 100 Если верхнюю границу отрезка, на котором мы ищем функцию, взять t∗ = 1, то количество шагов по времени будет N = τ1 ≈ 20000, а это, понятно, немало. Как исправить этот глобальный дефект? Можно ли, меняя шаблон, на котором происходит аппроксимация исходного дифференциального уравнения, менять соответствующий разностный аналог? 5.10 Неявная разностная схема для уравнения теплопроводности Изменим шаблон: Уравнения разностной схемы перепишутся следующим образом:  n+1 n+1 n+1 yj − yjn yj+1 − 2yjn+1 + yj−1    = + ϕn+1 , j = 1, N − 1, n = 0, K − 1;  j 2  τ h    yj0 = u0 (xj ), j = 0, N ;    y0n+1 = µ1 (tn+1 ), n = 0, K − 1;      n+1 yN = µ2 (tn+1 ), n = 0, K − 1. (5.33) Тогда невязка будет иметь вид: ψjn+1 = −
un+1 − unj 1 j n+1 n+1 + 2 un+1 + un+1 . j+1 − 2uj j−1 + ϕj τ h (5.34) 5.10. НЕЯВНАЯ РАЗНОСТНАЯ СХЕМА ДЛЯ УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ 103 n Как и в предыдущем случае, разложим все un+1 j±1 и uj в ряд Тейлора: un+1 j±1 unj h2 h3 ± un+1 + O(h4 ); xxx,j 2 6 2 = un+1 − un+1 j t,j τ + O(τ ). n+1 = un+1 ± un+1 j x,j h + uxx,j Теперь подставим эти разложения в формулу (5.34) и в правой части добавим и вычтем fjn+1 . Тогда невязка посчитается так: n+1 n+1 ψjn+1 = (−un+1 ) + ϕn+1 − fjn+1 + O(τ + h2 ). t,j + uxx,j + fj j Если функцию ϕn+1 взять равной fjn+1 с точностью O(τ + h2 ), то выражение для невязки сильно j сократится и примет вид: ψjn+1 = O(τ + h2 ). Как видно, эта разностная схема имеет первый порядок аппроксимации по τ и второй по h. Для исследования на устойчивость воспользуемся методом гармоник. Сопоставим нашему разностному уравнению однородное уравнение: n+1 n+1 yjn+1 − yjn yj+1 − 2yjn+1 + yj−1 = . τ h2 Проделаем те же действия, что и в предыдущем случае. Подставив в качестве решения yjn = q n eijhϕ и сократив множители, получим такое уравнение относительно параметров этого решения: eihϕ − 2 + e−ihϕ q−1 =q τ h2 4 q−1 hϕ = −q 2 sin2 . τ h 2 =⇒ Выразим отсюда q :   4τ hϕ q 1 + 2 sin2 =1 h 2 =⇒ q= 1 4τ 1 + 2 sin2 h hϕ 2 . Очевидно (так как знаменатель всегда больше или равен единице), что эта неявная разностная схема абсолютно устойчива (устойчива при любых значениях τ и h ). Посмотрим, как можно получить решение разностного уравнения из системы (5.33). Перепишем его так:   τ n+1 2τ τ n+1 y − 1 + yjn+1 + 2 yj+1 = −yjn + τ ϕn+1 . (5.35) j−1 j 2 2 h h h Обозначим для удобства Aj = τ , h2 Bj = τ , h2  Cj = 1+ 2τ h2  , Fjn = yjn − τ ϕn+1 . j Тогда (5.35) будет выглядеть постройнее: n+1 n+1 Aj yj−1 − Cj yjn+1 + Bj yj+1 = −Fjn . Рассмотрим эту систему при n = 0 : 1 1 Aj yj−1 − Cj yj1 + Bj yj+1 = −Fj0 , j = 1, N − 1 1 — это система линейных алгебраических уравнений относительно y 1 = (y01 , y11 , . . . , yN ). Перепишем ее в компактном виде: M y1 = F 0 , 104 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ где матрица F 0 состоит из элементов Fj0 (которые полностью определяются начальными условиями), а M — из коэффициентов Aj , Bj , Cj , расположенных на трех диагоналях (то есть матрица имеет трехдиагональный вид). Следовательно, применим метод прогонки, и мы можем найти сеточную функцию на первом временном слое. Поступая так дальше, мы сможем определить искомую сеточную функцию на всех временных слоях («послойно» применяя метод прогонки). 5.11 Разностная схема с весами для уравнения теплопроводности Рассмотрим теперь не минимально возможный шаблон, а «избыточный». Будем аппроксимировать производные в шести узлах. «Избыточность» схемы скомпенсируем введением некоторого параметра — весового множителя. Соответствующая этому шаблону разностная схема такова:  n+1 yj − yjn      τ      yj0      y0n+1      n+1 yN = σ n+1 n+1 n n yj+1 − 2yjn+1 + yj−1 yj+1 − 2yjn + yj−1 + (1 − σ) + ϕnj , h2 h2 j = 1, N − 1, n = 0, K − 1; (5.36) j = 0, N ; = u (xj ), = µ1 (tn+1 ), n = 0, K − 1; = µ2 (tn+1 ), n = 0, K − 1. — это так называемая разностная схема с весами для уравнения теплопроводности. Список вопросов остается тем же. Рассмотрим вопрос об аппроксимации. Выясним, как ведет себя невязка: n+ 1 def ψj 2 = un+1 − unj j n n − + σun+1 xx,j + (1 − σ)uxx,j + ϕj . τ Как мы уже поступали, разложим функции unj и un+1 в ряд Тейлора: j n+ 12 un+1 j = uj unj = uj n+ 12 τ 1 n+ 1  τ 2 + utt,j 2 + O(τ 3 ); 2 2 2 1 n+ 1  τ 2 n+ 1 τ − ut,j 2 + utt,j 2 + O(τ 3 ). 2 2 2 n+ 12 + ut,j Тогда получим: n+ 12 ψj n+ 12 = −ut,j n n + O(τ 2 ) + σun+1 xx,j + (1 − σ)uxx,j + ϕj . (5.37) 5.11. РАЗНОСТНАЯ СХЕМА С ВЕСАМИ ДЛЯ УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ 105 Теперь в представлении второй разностной производной разложим все вхождения функции в ряд Тейлора с членами до пятого порядка включительно: 1 (u(xj+1 , t) − 2u(xj , t) + u(xj−1 , t)) = h2 h2 h3 h4 h5 = {u(xj±1 , t) = u(xj , t) ± ux (xj , t)h + uxx (xj , t) ± uxxx (xj , t) + uxxxx (xj , t) ± uxxxxx (xj , t) + O(h6 )} = 2 6 4! 5! h2 4 = uxx (xj , t) + uxxxx (xj , t) + O(h ). 12 uxx (xj , t) = Воспользуемся этим разложением для слагаемых в выражении (5.37): h2 n+ 1 + O(h4 ) = {Разложение в ряд Тейлора с центром в точке xj 2 } = 12 τ h2 h2 τ n+ 21 n+ 12 n+ 21 n+ 12 = uxx,j + uxxt,j · + uxxxx,j · · + uxxxxt,j · + O(τ 2 + h4 ); 2 12 12 2 h2 n+ 1 = unxx,j + unxxxx,j · + O(h4 ) = {Разложение в ряд Тейлора с центром в точке xj 2 } = 12 τ h2 h2 τ n+ 21 n+ 21 n+ 12 n+ 12 = uxx,j − uxxt,j · + uxxxx,j − uxxxxt,j · + O(τ 2 + h4 ). · · 2 12 12 2 n+1 n+1 = un+1 uxx,j xx,j + uxxxx,j · n uxx,j Таким образом, выражение для невязки принимает вид: n+ 21 ψj n+ 21 = −ut,j n+ 1 n+ 12 Добавляя и вычитая fj n+ 12 ψj n+ 1 n+ 1 2 + uxx,j2 + uxxxx,j n+ 1 n+ 21 = −ut,j 2 + uxx,j2 + fj h2 1 n+ 12 1 n+ 12 h2 + ϕnj + (σ − )uxxt,j τ + O(τ 2 + h4 ). τ + (σ − )uxxxxt,j 12 2 2 12 , получим эквивалентное выражение: n+ 1 2 + uxxxx,j h2 1 n+ 12 1 n+ 12 h2 n+ 1 + ϕnj − fj 2 + (σ − )uxxt,j τ + O(τ 2 + h4 ). τ + (σ − )uxxxxt,j 12 2 2 12 Согласно уравнению теплопроводности, первые три слагаемых обращаются в нуль: n+ 12 ψj n+ 1 2 = uxxxx,j h2 1 n+ 21 1 n+ 21 h2 n+ 1 + ϕnj − fj 2 + (σ − )uxxt,j τ + O(τ 2 + h4 ). τ + (σ − )uxxxxt,j 12 2 2 12 При σ = 21 схема (5.36) называется симметричной. Тогда в последнем равенстве последние слагаемые обнулятся, и с помощью условия на параметр ϕnj n+ 12 ϕnj = fj + O(τ 2 + h2 ) мы можем достичь такого порядка аппроксимации: n+ 12 ψj = O(τ 2 + h2 ). Теперь вернемся на шаг назад и воспользуемся тем, что ut = uxx + f =⇒ uxxt = uxxxx + fxx . Тогда формула для невязки будет несколько иной:   2 1 h2 n+ 12 1 n+ 12 h2 n+ 1 n+ 1 n+ 1 h ψj 2 = ϕnj − fj 2 − fxx,j2 · + (σ − )τ + uxxt,j + (σ − )uxxxxt,j · τ + O(τ 2 + h4 ). 12 2 12 2 12 1 h2 Взяв σ = − , мы обнулим четвертое слагаемое, а коэффициент при пятом оценится как O(h4 ). 2 12τ Осталось потребовать, чтобы n+ 21 ϕnj = fj n+ 1 + fxx,j2 · h2 + O(τ 2 + h4 ), 12 106 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ тогда порядок аппроксимации будет таков: n+ 12 ψj = O(τ 2 + h4 ). При данном σ выражение (5.36) называется разностной схемой повышенного порядка аппроксимации. При всех остальных σ порядок аппроксимации будет меньше: n+ 12 ψj = O(τ + h2 ) при условии, что n+ 21 ϕnj = fj + O(τ + h2 ). Теперь исследуем схему на устойчивость методом гармоник. Для начала запишем однородное уравнение: n+1 n+1 n n yjn+1 − yjn yj+1 − 2yjn+1 + yj−1 yj+1 − 2yjn + yj−1 =σ + (1 − σ) . τ h2 h2 Подставим в качестве решения yjn = q n eijhϕ . Сокращая степени q, получим: q−1 q 1 = σ 2 (eihϕ − 2 + e−ihϕ ) + (1 − σ) 2 (eihϕ − 2 + e−ihϕ ) ⇐⇒ τ h h q−1 q 1 hϕ hϕ = −4σ 2 sin2 − 4(1 − σ) 2 sin2 . τ h 2 h 2 Отсюда получаем выражение для q : τ τ hϕ hϕ sin2 ) = 1 − (1 − σ)4 2 sin2 =⇒ h2 2 h 2 τ hϕ 1 − (1 − σ)4 2 sin2 h 2 . q= τ 2 hϕ 1 + 4σ 2 sin h 2 q(1 + 4σ Для получения условий на устойчивость мы требуем, чтобы |q| 6 1. В данном случае это дает два неравенства:  τ τ hϕ  2 hϕ  6 1 + 4σ 2 sin2 ;  1 − (1 − σ)4 2 sin h 2 h 2  hϕ τ hϕ τ   6 1 − (1 − σ)4 2 sin2 . −1 − 4σ 2 sin2 h 2 h 2 Первое выполнено всегда, так как τ > 0. Второе перепишется так: −8σ τ hϕ τ hϕ sin2 6 −4 2 sin2 +2 h2 2 h 2 ⇐⇒ σ> 1 − 2 h2 4τ sin2 hϕ 2 — оно должно быть выполнено при любом ϕ. Взяв максимум по правой части, приходим к окончательному условию для σ : 1 h2 . σ> − 2 4τ 1 Значение σ = удовлетворяет этому неравенству. Это означает, что соответствующий метод яв2 ляется абсолютно устойчивым. Последним нашим долгом будет обосновать возможность вычисления приближения по этой схеме. Для этого перепишем (5.36) в таком виде: σ τ n+1 τ τ n+1 n y − (1 + 2σ 2 )yjn+1 + σ 2 yj−1 = −yjn − (1 − σ)τ yxx,j − ϕnj . h2 j+1 h h 5.12. РАЗНОСТНЫЕ СХЕМЫ ДЛЯ УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ ОСОБОГО ТИПА 107 Теперь почти очевидно, что мы можем получить все yjk (j = 0, N , k = 0, K). Действительно, фиксируем n = 0. Тогда мы получим систему линейных уравнений с трехдиагональной матрицей. Правые части уравнений мы можем найти, используя начальные условия (заметим, что ϕnj мы задали при исследовании порядка аппроксимации). После этого применяется метод прогонки, после которого становятся известны все yj1 . Теперь можно увеличить n и снова получить СЛАУ, подставив в правую часть уравнений только что найденные yj1 . Так действуем, пока не найдем все yjk . Мы рассмотрели схему для нахождения численного решения простейшей краевой задачи. Для нее существует более простое аналитическое решение, но в общем случае его может и не быть. В то же время вся методика построения разностных схем и их решение достаточно легко переносятся на более сложные задачи, к которым мы и перейдем. 5.12 Разностные схемы для уравнения теплопроводности особого типа Разностная схема для уравнения теплопроводности с переменными коэффициентами Рассмотрим такую краевую задачу:    ∂u ∂ ∂u   ρ(x, t) = k(x, t) + f (x, t),    ∂t ∂x ∂x u(0, t) = µ1 (t), 0 6 t 6 T ;   u(1, t) = µ2 (t), 0 6 t 6 T ;    u(x, 0) = u0 (x), 0 6 x 6 1. 0 < x < 1, 0 < t 6 T ; Аналитического выражения для решения нет. Тем не менее, известно, что если всюду верны неравенства 0 < c1 6 ρ(x, t); 0 < c2 6 k(x, t), то оно существует и единственно. Теперь будем приближать ut соответствующей разностной производной: ∂u ≈ unt,j , ∂t а производные по x с использованием интегро-интерполяционного метода можно представить следующим образом:    
∂ ∂u 1 uj+1 − uj uj − uj−1 k(x, t) ≈ aux x,j = aj+1 − aj , ∂x ∂x h h h где aj вычисляются по такой формуле:  1 aj =  h Zxj −1 dx   k(x, t) ≈ k(xj− 21 , t). xj−1 Используя шаблон из шести точек, мы построим разностную схему. Выкладки аналогичны предыдущим схемам, поэтому опустим их и приведем только окончательный результат: 
n+1
n yjn+1 − yjn    = σ ayx x,j + (1 − σ) ayx x,j + ϕnj ; j = 1, N − 1, n = 0, K − 1 ρ(x , t) j   τ    yj0 = u0 (xj ), j = 0, N ; (5.38)   n  y0 = µ1 (tn ), n = 1, K;      n yN = µ2 (tn ), n = 1, K. 108 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ Здесь остаются нефиксированными момент времени t в первом уравнении (от него зависят aj и 1 ρ) и параметр метода σ. Взяв σ = и t = tn+ 21 , мы можем получить такую оценку на порядок 2 аппроксимации: n+ 12 ψj = O(τ 2 + h2 ). В противном случае оценка будет похуже: ψjn = O(τ + h2 ). Наконец, общая формула для получения y будет такова: n+1 n+1 Aj yj+1 − Cj yjn+1 + Bj yj−1 = −Fj . При этом правая часть уравнения зависит только от yn . Это означает, что, решая данное уравнение послойно (с n = 0 ) и используя начальные условия, мы можем найти все yjk . Все эти умозаключения абсолютно идентичны тем, что были в конце предыдущего параграфа. Разностная схема для нелинейного уравнения теплопроводности Исследуем случай, когда уравнение теплопроводности в краевой задаче имеет такой вид:   ∂ ∂u ∂u = k(u, x, t) + f (u, x, t). ∂t ∂x ∂x — это означает, что коэффициенты при производных зависят еще и от искомой функции u, причем, вообще говоря, нелинейно. В этом случае рекомендуется использовать неявные разностные схемы, так как они чаще всего абсолютно устойчивы. Приведем пример:
n+1 yjn+1 − yjn = ayx x,j + f (yjn+1 ), τ при этом коэффициенты a зависят еще и от y. Расписав разностные производные, получим такую формулу: " # n+1 n+1 yjn+1 − yjn yj+1 − yjn+1 yjn+1 − yj−1 1 n+1 n+1 = aj+1 (yj ) · − aj (yj ) · + f (yjn+1 ). τ h h h — для каждого слоя это система нелинейных уравнений относительно yjn+1 . Решается она итераци(k) онным методом следующего вида. Если обозначить за yj приближение для yjn+1 , то формула для получения следующего приближения будет такова:   (k+1) (k+1) (k+1) (k+1) − yjn (k) (k) (k) 1 yj+1 − yj yj − yj−1  = aj+1 ( yj ) · − aj ( yj ) · + f ( yj ) τ h h h (k+1) yj — это уже система линейных уравнений с трехдиагональной матрицей. Она, в свою очередь, решается методом прогонки, и мы получаем (k + 1)-е приближение к yjn+1 . Обычно ограничиваются пятью приближениями: (5) yjn+1 = yj . 5.13. РАЗНОСТНАЯ СХЕМА ДЛЯ УРАВНЕНИЯ КОЛЕБАНИЙ 5.13 109 Разностная схема для уравнения колебаний Рассмотрим стандартную краевую задачу на уравнение колебаний:  utt = uxx + f (x, t), 0 < x < 1, 0 < t 6 T ;     0 6 t 6 T;  u(0, t) = µ1 (t), u(1, t) = µ2 (t), 0 6 t 6 T;   u(x, 0) = u0 (x), 0 6 x 6 1;    ut (x, 0) = ψ(x), 0 6 x 6 1. (5.39) Решение такой задачи существует и единственно. Введем дискретную сетку на рассматриваемой области:  ωk = xj = jh, h = N1 , j = 0, N ;  T , n = 0, K . ωτ = tn = nτ, τ = K Заметим, что начальные условия дают нам значения искомой функции на границе прямоугольника. Будем использовать уже привычные нам обозначения: u(xj , tn ) = unj — точное решение в узлах сетки, yjn — искомое приближение. Выпишем аналоги краевых и начальных условий:  n y = µ1 (tn );    0 n yN = µ2 (tn );    0 yj = u0 (xj ). Сопоставим уравнению дискретный аналог в каждом узле. Будем использовать шаблон из пяти точек (это минимально необходимый шаблон для аппроксимации вторых производных по t и x). Приближения построим так:
1 n utt ≈ utt,j = 2 un+1 − 2unj + un−1 ; j j τ
1 uxx ≈ unxx,j = 2 unj+1 − 2unj + unj−1 . h Тогда получим следующий дискретный аналог исходной задачи:  n n ytt,j = yxx,j + ϕnj , j = 1, N − 1, n = 1, K − 1;        y0n = µ1 (tn ), n = 1, K;     n yN = µ2 (tn ), n = 1, K; (5.40)    yj = u (xj ), j = 0, N ;      1    yj − yj = ψ(xj ), j = 0, N . τ Шаблон будет выглядеть следующим образом: Стоит отметить, что для того, чтобы разностная схема была сбалансированной (то есть, чтобы не делать лишних вычислений в одном месте, а потом загрублять их в другом), необходимо, чтобы порядки аппроксимации в уравнении и в краевом условии второго типа были согласованы. Иначе, если использовать низкий (первый) порядок аппроксимации в краевом условии, то вся схема будет аппроксимировать исходную задачу с первым порядком аппроксимации. Постараемся все данные в задаче приблизить со вторым порядком точности. Для этого разложим u1j в ряд Тейлора в точке (xj , 0) с остаточным членом второго порядка малости: u1j − u0j τ = u0t,j + u0tt,j + O(τ 2 ). τ 2 110 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ Постараемся избавиться от u0tt,j . Пусть на границе также выполняется уравнение колебаний, тогда utt = uxx + f (x, 0), и краевое условие примет вид:
τ u1j − u0j = u0t,j + u0xx,j + f (xj , 0) + O(τ 2 ), τ 2 где выражение u0xx,j + f (xj , 0) точно вычислимо из начальных условий. Тогда разностную схему можно переписать в виде:  n n ytt,j + ϕnj , j = 1, N − 1, n = 1, K − 1; = yxx,j        y0n = µ1 (tn ), n = 1, K;     n yN = µ2 (tn ), n = 1, K;    yj0 = u0 (xj ), j = 0, N ;     1     yj − yj = ψ(xj ) + τ u0 (xj ) + f 0
, j = 0, N . j τ 2 xx Порядок аппроксимации краевого условия будет O(τ 2 ). Потребуем в задаче второй порядок аппроксимации и в первом уравнении. Для этого достаточно взять ϕnj = f (xj , tn ) + O(h2 + τ 2 ). Таким образом, разностный аналог исходного ДУ будет его аппроксимировать со вторым порядком по τ и по h. Рассмотрим вопрос о поиске решения. Найдем из построенной разностной схемы выражение для yjn+1 : n yjn+1 = 2yjn − yjn−1 + τ 2 yxx,j + τ 2 ϕnj . (5.41) Покажем,что значение сеточной функции вычислимо на всей сетке. Рассмотрим формулу (5.41) при n = 1 :
τ2 1 1 yj2 = 2yj1 − yj0 + 2 yj+1 − 2yj1 + yj−1 + τ 2 ϕ1j . h Так как yj0 известно из начальных условий, а yj1 можем найти из второго краевого условия, то 2 yj находится явно на всем втором временном слое. Аналогично можно найти yj3 и так далее, то есть определить значение искомой сеточной функции на всей сетке (как и раньше, делая это послойно). Как видно из процесса поиска решения, эта схема — явная, и построенное решение будет единственным. Рассмотрим вопрос устойчивости. Будем использовать привычный нам метод гармоник. Запишем однородное уравнение:
τ2 n n yjn+1 − 2yjn + yjn−1 = 2 yj+1 − 2yjn + yj−1 . h Подставим в него решение вида yjn = q n eijhϕ , тогда получим:   τ2 τ2 2 hϕ ihϕ −ihϕ 2 2 q − 2q + 1 = 2 q(e −2+e ) ⇐⇒ q − 2 − 4 2 sin q + 1 = 0. h h 2 5.14. РАЗНОСТНАЯ АППРОКСИМАЦИЯ ЗАДАЧИ ДИРИХЛЕ 111 Решения этого квадратного уравнения будут такими:   s 2 τ2 hϕ τ2 2 hϕ q1,2 = 1 − 2 2 sin − 1. ± 1 − 2 2 sin2 h 2 h 2 Выделим два случая на дискриминант (подкоренное выражение): 1. D > 0 : Так как свободный член равен единице, то |q1 ||q2 | = 1. В этом случае (можно показать, что корни не могут быть противоположными по знаку) |q1 | или |q2 | больше единицы, поэтому устойчивости не будет при данных значениях параметров. 2. D 6 0 : Аналогично, |q1 ||q2 | = 1. Найдем значение абсолютной величины корня: 2 |q1 | =  τ2 hϕ 1 − 2 2 sin2 h 2 2 τ2 hϕ + 1 − 1 − 2 2 sin2 h 2  2 = 1. Тогда |q1 | = |q2 | = 1, и, следовательно, решение будет устойчивым. Распишем условие неположительности дискриминанта: |1 − 2 hϕ τ2 sin2 | 6 1, h2 2 или, расписав: τ2 hϕ sin2 6 1. h2 2 Правое неравенство выполнено всегда, перепишем левое: −1 6 1 − 2 τ2 hϕ sin2 6 1. 2 h 2 τ2 τ 6 1, что эквивалентно 6 1 (напомним, что в уравнении теплоh2 h 1 τ проводности мы получали ограничение вида 2 6 ). Получается, что построенная схема является h 2 условно устойчивой. Итак, можно выделить несколько особенностей для уравнения колебаний: 1) дополнительное условие на аппроксимацию (из начального условия); 2) решение строится послойно (впрочем, как и для многих других задач); 3) условная устойчивость с простыми ограничениями. В худшем случае получим 5.14 Разностная аппроксимация задачи Дирихле Кратко опишем численное решение задачи Дирихле на уравнение Пуассона. Пусть область, для которой задано уравнение — прямоугольник в E 2 :  ux1 x1 + ux2 x2 = f (x1 , x2 ), 0 < x1 < l1 , 0 < x2 < l2 ;     u(0, x2 ) = µ1 (x2 ), 0 6 x2 6 l2 ;  u(l1 , x2 ) = µ2 (x2 ), 0 6 x2 6 l2 ;   u(x , 0) = µ (x ), 6 x1 6 l1 ;  1 1 1   u(x1 , l2 ) = µ2 (x1 ), 0 6 x1 6 l1 . 112 Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧ Дискретизируем область: ωhx1 ωhx2   l1 = x1i = ih1 , i = 0, N1 , h1 = ; N1   l2 . = x2j = jh2 , j = 0, N2 , h2 = N2 Выглядеть это будет следующим образом (обозначения аналогичны предыдущим рисункам): `2 × × × × × × × × × × × × × ˆ‰ Š‹ Œ Ž ‘ ’“ ”• –— ˜™ š› lm no pq rs tu vw xy z{ |} ~ € NO PQ RS TU VW XY Z[ \] ^_ `a bc de fg 45 67 89 :; <= >? @A BC DE FG HI JK LM    ! "# $% &' () *+ ,- ./ 01 23
          × × × × × × × × × × × × × × × × †‡ jk × × × „ hi × × × ‚ƒ × × × × × × × `1 Запишем дискретные аналоги граничных условий:        y0,j yN1 ,j yi,0 yi,N2 = = = = µ1 (x2j ), µ2 (x2j ), µ1 (x1i ), µ2 (x1i ), j = 0, N2 ; j = 0, N2 ; i = 0, N1 ; i = 0, N1 . (5.42) Производные будем приближать разностными производными со вторым порядком аппроксимации: yx1 x1 ,ij + yx2 x2 ,ij = ϕij , i = 1, N1 − 1, j = 1, N2 − 1, (5.43) где для достижения второго порядка аппроксимации во всем уравнении ϕij считается по формуле: ϕij = f (x1i , x2j ) + O(h21 + h22 ). Найдем формулу для поиска приближенного решения. Явно получить выражение для приближений на следующем слое из (5.43) довольно сложно, поэтому сделаем хитрее. Введем вектор всех неизвестных (заметим, всех, то есть получим огромный вектор), используя самый простой, регулярный способ нумерации: Y = y11 , y21 , . . . , y(N1 −1)1 , y12 , y22 , . . . , y(N1 −1)2 , . . . , y1(N2 −1) , . . . , y(N1 −1)(N2 −1)
T . Аналогично введем вектор правой части: F = ϕ11 , ϕ21 , . . . , ϕ(N1 −1)1 , ϕ12 , ϕ22 , . . . , ϕ(N1 −1)2 , . . . , ϕ1(N2 −1) , . . . , ϕ(N1 −1)(N2 −1)
T . Объединив, используя эти обозначения, уравнения (5.42) и (5.43), получим уравнение: AY = F. Из структуры шаблона следует такой вид матрицы A (она будет иметь пять ненулевых «диагона- 5.14. РАЗНОСТНАЯ АППРОКСИМАЦИЯ ЗАДАЧИ ДИРИХЛЕ 113 лей»; i, j, s — некоторые параметры, зависящие от N1 , N2 ):   a11 a12 ... a1j ...  a21 a22 a23 ... a2j ... 0     0 a32 a a . . . a . . . 0  33 34 3j     . . . . . . . . . . . . . . . . . . . . . . .     0 ... a a a . . . a jj j1  j(j−1) j(j+1)   . . . . . . . . . . . . . . . . . . . . . . . A=     ai1 0 . . . a a a . . . ii i(i−1) i(i+1)     . . . . . . . . . . . . . . . . . . . . . . .     0 . . . a(s−2)(i−2) . . . a a a (s−2)(s−3) (s−2)(s−2) (s−2)(s−1)    0 ... a(s−1)(i−1) 0 ... a(s−1)(s−2) a(s−1)(s−1) a(s−1)s  0 0 ... asi ... as(s−1) ass — так называемая матрица ленточной структуры5 . Мы, кстати, опять получили сеточное уравнение (напомним, подробные указания о том, как их решать, можно найти в [2]). 5 названа в честь ленточного червя (тоже большая гадость). Оглавление 1 Решение систем линейных алгебраических уравнений 1.1 Прямые методы решения СЛАУ. Метод квадратного корня . 1.2 Линейные одношаговые итерационные методы . . . . . . . . 1.3 Сходимость одношаговых стационарных методов . . . . . . . 1.4 Оценка погрешности одношаговых стационарных методов . 1.5 Попеременно-треугольный итерационный метод . . . . . . . 1.6 Чебышевский набор итерационных параметров . . . . . . . . 1.7 Одношаговые итерационные методы вариационного типа . . 1.8 Примеры итерационных методов вариационного типа . . . . 1.9 Двухшаговые итерационные методы вариационного типа . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 4 9 13 17 21 25 29 31 34 2 Задачи на собственные значения 36 2.1 Поиск собственных значений методом вращений . . . . . . . . . . . . . . . . . . . . . . . . 36 2.2 Степенной метод поиска собственных значений . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.3 Метод обратной итерации . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3 Численные методы решения нелинейных 3.1 Методы разделения корней . . . . . . . . 3.2 Примеры численных методов . . . . . . . 3.3 Сходимость метода простой итерации . . 3.4 Метод Эйткена . . . . . . . . . . . . . . . . 3.5 Сходимость метода Ньютона . . . . . . . . 3.6 Решение систем нелинейных уравнений . уравнений . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 43 43 46 48 48 51 4 Интерполяция и приближение функций 57 4.1 Интерполирование кубическими сплайнами . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.2 Наилучшее приближение табличной функции . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5 Численные методы решения краевых задач 5.1 Сходимость методов Рунге-Кутта . . . . . . . . . . . . . . . . . . . 5.2 Методы Рунге-Кутта второго порядка аппроксимации . . . . . . . 5.3 Описание многошаговых методов . . . . . . . . . . . . . . . . . . . 5.4 Методы Адамса и Гира . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Устойчивость численных методов решения задачи Коши . . . . . 5.6 Интегро-интерполяционный метод построения разностных схем . 5.7 Метод аппроксимации квадратичного функционала . . . . . . . . 5.8 Корректность разностной схемы . . . . . . . . . . . . . . . . . . . 5.9 Явная разностная схема для уравнения теплопроводности . . . . 5.10 Неявная разностная схема для уравнения теплопроводности . . . 5.11 Разностная схема с весами для уравнения теплопроводности . . . 5.12 Разностные схемы для уравнения теплопроводности особого типа 5.13 Разностная схема для уравнения колебаний . . . . . . . . . . . . . 5.14 Разностная аппроксимация задачи Дирихле . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . *Литература [1] А. А. Самарский, А. В. Гулин. "Численные методы". [2] А. А. Самарский, Е. С. Николаев. "Методы решения сеточных уравнений". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 74 77 79 81 85 91 94 95 97 102 104 107 109 111