Вычислительная математика. Лекции

3.3Метод квадратного корня (метод Холецкого).

Этот метод предназначен для решения систем с симметричной матрицей A, главные миноры которой отличны от 0:

Матрица A представляется в виде произведения трёх матриц

A = S DS;

где S верхняя треугольная матрица (S нижняя треугольная), а матрица D- диагональная матрица, на

равенств s22d2s2j = a2j s12d1s1j находим все элементы второй строчки матрицы S.

Аналогичным способом рассматривая строчки i = 3; 4; :::; n; получаем последовательно числа di, числа sii > 0 и все элементы i-той строчки матрицы S.

Решение СЛАУ. Ax = y: S DSx = y: Обозначим DSx = z: Тогда S z = y система с нижней треугольной матрицей, диагональные элементы которой sii > 0: Решение z этой системы легко получить. Решение x исходной системы найдем как решение системы DSx = z с верхней треугольной матрицей, диагональные элементы которой отличны от нуля.

3.4Метод окаймления построения обратной матрицы.

Метод окаймления состоит в последовательном построении обратных матриц, соответствующих главным минорам матрицы A. Поэтому предполагается, что все главные миноры матрицы отличны от нуля.

Матрица A1 = a11 размерности 1 1 имеет обратную:

18

Из (4) : sT + uT Am1 1 = 0T и sT = uT Am1 1:

Из (1) : Bm 1 + Am1 1vsT = Am1 1 и Bm 1 = Am1 1(Em 1 vsT ):

Покажем, что am uT Am1 1v 6= 0: СЛАУ с матрицей Am имеет единственное решение при любой правой части. Запишем матрицу Am в блочном виде, в правой части возьмем вектор, у которого только

Тогда Am 1x + xmv = 0, x + xmAm1 1v = 0 и uT x + xmuT Am1 1v = 0:

Последнее уравнение системы имеет вид

uT x + amxm = ym:

Сравнивая эти два полученные уравнения, получаем:

xm(am uT Am1 1v) = ym:

Так как координата xm определяется однозначно при любом значении ym, то величина am uT Am1 1v 6=

0:

Итак, продолжая процесс построения матриц Am1; построим и матрицу An 1: Число операций равно n3: Замечание. Пусть в математическую модель кроме параметров x1; x2; :::; xn вводится еще один параметр xn+1; значение которого неизвестно. В случае линейных моделей получаем матрицу An+1: Метод окаймления позволяет получить матрицу An+11 ; исходя из уже построенной матрицы An 1 за (3n2 + 3n 1)

операций.

3.5Метод прогонки решения систем с трехдиагональной

матрицей.

Рассматривается СЛАУ с матрицей размерности (n + 2) (n + 2) :

19

К решению таких систем сводится например задача приближенного решения граничной задачи для линейного дифференциального уравнения второго порядка с переменными коэффициентами:

Для решения систем с трехдиагональной матрицей можно, например, применить метод Гаусса. Учитывая особенности матрицы, можно существенно упростить процесс построения решения.

Коэффициенты i, i i = (1; 2; :::; n + 1) называются прогоночными коэффициентами.

Решение системы xi получим (обратный ход), рассматривая сначала последнее уравнение системы

(i = n + 1): 2xn + xn+1 = 2 вместе с представлением xn = n+1xn+1 + n+1: Отсюда получим

xn+1 = 1+ 2 n+1 , если 1 2 n+1 6= 0

Значения xn; xn 1; :::; x0 вычислим последовательно:

xi = n+1xi+1 + i+1; (i = n; n 1; :::; x1; x0):

20

При вычислении xi мы замечаем, что если значения j ij > 1; то при большой размерности системы погрешности в задании правых частей системы и погрешности значений xi будут возрастать, т.е. процесс вычисления xi будет неустойчивым относительно погрешности вычислений. Поэтому, кроме условий ci iai; 1 2 n+1 6= 0, желательно ввести требование j ij 1:

Теорема. Если jcij jaij + jbij; j 1j 1; j 2j 1; то ci iai 6= 0; j ij 1 для всех i.

Доказательство. Так 1 = 1, то j 1j 1: Далее величины

jci iaij jbij jcij j iaij jbij jaij + jbij j ijjaij jbij = jaij(1 j ij)

и

jci iaij jbij + jaij(1 j ij) jbij > 0:

Тогда

jbij

j i+1j = jcij iai 1;

Итак j 1j 1; j 2j 1; :::; j n+1j 1 и j1 2 n+1j 1 j 2j > 0:

Замечание. Рассмотренный метод прогонки называется методом верхней прогонки. Если исходить из начального представления в виде xi = i 1xi 1 + i 1; то получим метод нижней прогонки. Можно построить и методы встречных прогонок.

Тогда Bk xk+1 xk + Axk = y, и мы получаем итерационный процесс:

xk+1 = xk Bk 1Axk + Bk 1y = (E Bk 1A)xk + Bk 1y = Sxk + Bk 1y:

Матрица S- матрица перехода. Если k E Bk 1A k< 1, то итерационный процесс сходится к решению системы.

3.6Метод простой итерации.

Рассматривается система Ax = y, где матрица E B :

x = Bx + y:

Задавая вектор xo построим вектор x1 = Bx0 + y и далее xk = Bxk 1 + y: Ясно, что верна теорема:

Теорема. Если k B k= q < 1; то существует единственное решение x , последовательность fxkg стремится к x со скоростью геометрической прогрессии для любого начального вектора x0.

Доказательство. Рассмотрим однородную систему x = Bx: Для её решения верно неравенство k x k=k Bx k k B kk x k= q k x k : Тогда x = 0: Следовательно, решение исходного уравнения существует и единственно при любой правой части y.

Так как x = Bx + y и xk = Bxk 1 + y; то

k x xk k=k Bx Bxk 1 k q k x xk 1 k qk k x x0 k :

Замечание. Полученная оценка называется априорной оценкой, так как вектор x неизвестен. Однако, зная xk 1 и xk, легко получить так называемую апостериорную оценку:

k x xk k=k x xk+1 k + k xk+1 xk k

k x xk k + k Bxk Bxk 1 k q k xk xk 1 k;

21

k x xk k 1 q q k xk xk 1 k :

Система x = Bx + y имеет специальный вид. Для СЛАУ с симметричной положительной матрицей можно рассмотреть систему x (E A)x = b: Если собственные числа матрицы A равны 0 < 1 2

::: n; собственные числа матрицы равны 1 i, и норма матрицы E A k E A k2= max j1 ij:

i

Для произвольной матрицы A можно рассмотреть СЛАУ A Ax = A y и построить для этой системы указанный итерационный процесс. Но следует иметь ввиду, что число обусловленности (A A) увеличи-

вается:

(A A) =k A A kk (A A) 1 k= 2(A) > (A):

3.7Метод Зейделя.

Будем предполагать, что диагональные элементы матрицы A отличны от нуля. Матрицу A представим в виде

A = L + D + U;

где матрица D - диагональная матрица, матрицы L и U нижняя и верхняя треугольные матрицы с нулевыми диагоналями.

Метод Зейделя состоит в задании начального вектора x0 и в последовательном вычислении векторов

xk:

(L + D)xk+1 = y Uxk;

Тогда

xk+1 = (L + D) 1Uxk + (L + D) 1y:

При реализации этого метода не требуется строить матрицу (L + D) 1:

Действительно:

a11xk1+1 = (a12xk2 + a13xk3 + ::: + a1nxkn) + y1, получаем xk1+1 a21xk1+1 + a22xk2+1 = (a23xk3 + ::: + a2nxkn) + y2, получаем xk2+1

: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :

an1xk1+1 + an2xk2+1 + ::: + annxkn+1 = yn, получаем xkn+1

и таким образом строим xk+1:

Вариантами метода Зейделя являются методы релаксации

(!L + D) 1xk+1 = (1 !)Dxk !Uxk + !y

с выбором параметра релаксации ! 2 (0; 2): При ! 2 (0; 1) получаем методы нижней релаксации, при

!2 (1; 2) - методы верхней релаксации. Наиболее употребительны методы верхней релаксации. Сходимость метода Зейделя легко доказывается, например в предположении, что матрица A - матрица

PP

22

Глава 4

МЕТОДЫ МИНИМИЗАЦИИ КВАДРАТИЧНОЙ ФУНКЦИИ.

Пусть f(x) непрерывная функция, отображающая множество D Rn в множество R1. Ставятся две задачи:

-найти inf f(x), x 2 D

-найти точку x 2 D, такую что f(x) f(x ), x 2 D.

В этой главе рассматриваются методы решения этих задач на простом примере минимизации квадратичной функции. Одновременно строятся и методы решения систем линейных алгебраических уравнений с положительно определенной матрицей.

4.1Квадратичная функция.

Функция

где точки x(x; :::; xn) 2 Rn, (D = Rn) называется квадратичной функцией.

Мы не будем различать точки x 2 Rn и векторы x 2 Vn с составляющими (x1; x2; :::; xn). Введя в Vn скалярное произведение, получим

f(x) = 12(Ax; x) + (b; x);

где матрица A симметрична. Собственные числа i матрицы A вещественны и 0 i n:

В этой главе рассматриваются матрицы A, имеющие обратные матрицы. Тогда 1 > 0, 1 2 :::

Значения собственных чисел i как правило неизвестны, но на практике известны оценки 1 m > 0,n M. Тогда (k z k=k z k2)

m k x k2 (Ax; x) M k x k2 : : : : : : : : : : : : : : : : : : : : : : : : : : : 10

Значения jf(x)j при неограниченном росте k x k неограниченно возрастают, так как (Ax; x) 1 k x k2, а значения j(b; x)j имеют порядок роста k x k. Таким образом inf f(x) может достигаться на ограниченном множестве.

Рассмотрим

f(x + ) f(x) = 12(A(x + ); x + ) + (b; x + ) 12(Ax; x) (b; x) =

23

Пусть x0- любой вектор. Если Ax0 + b 6= 0, то в точке x0 нет ни максимума, ни минимума: и при малых значениях k x k знак разности

f(x0 + x) f(x0) определяется знаком (Ax0 + b; x), и при замене x на x знак f(x0 + x) f(x0) меняется на противоположный. Поэтому точка x0 не может быть точкой экстремума функции f(x): Таким образом, единственная точка экстремума функции f(x) является единственным решением x системы линейных алгебраических уравнений

Ax = b:

В точке x f(x + x) f(x ) = 12 (A x; x) > 0, следовательно точка x единственная точка минимума функции f(x):

1. Из равенства

k f(x + x) f(x) (Ax + b; x) k= 0(k x k)2

согласно определению производной отображения f в точке x, получаем что вектор Ax + b определяет линейный оператор f0(x) действующий из Vn в R1 по формуле

f0(x) = (Ax + b; ); f0(x) = Ax + b:

2. Дифференцируя функцию f(x) = f(x1; x2; :::; xn) по координатам xi, получаем

и Ax + b = grad f(x) и вектор grad f(x) будем обозначать r(x) :

r(x) = grad f(x) = f0(x) = Ax + b

Вектор grad f(x) = r(x) определяет направление наискорейшего возрастания значений f(x) в окрестности точки x, а вектор r(x) определяет направление наискорейшего убывания этих значений.

Методы минимизации значений функции f(x) в соответствии с поставленными задачами можно также разделить на две группы:

- близость точки x к точке x оценивается величиной k x x k;

- близость точки x к точке x оценивается значениями f(x) f(x ):

В обоих случаях используются градиенты функции f(x); но методы первой группы называются градиентными методами, а методы второй группы - методами спуска.

4.2Градиентные методы.

Пусть выбран начальный вектор x0. Векторы xk строятся последовательно по правилам xk+1 = xk(Axk + b), > 0: Вектор r(xk) = Axk + b обозначим rk: rk = Axk + b:

Возникает вопрос: как выбрать значение ? Оценим

известны только оценки m и M, то в этих оценках следует заменить 1 на m, а n на M.

Таким образом, при 2 (0; 2 ) градиентный метод сходится к точке x со скоростью геометрической

n

прогрессии:

k xk+1 x k qk+1 k x x0 k= Cqk+1; 0 < q < 1

(при = 0 значение q = n 1 ):

n+ 1

24

4.3Метод наискорейшего спуска.

Вэтом методе величина (шаг по направлению антиградиента) в итерационном процессе xk+1 =

xk rk выбирается из условия минимума значения f(xk+1).

Согласно (2): f(xk+1) = f(xk rk) = f(xk) (Axk + b; rk) + 12 (A rk; rk) = = f(xk) (rk; rk) + 12 2(Ar; rk) = '( ):

Значение находим из условия '0( ) = 0 : и оптимальное значение :

= k = (rk; rk) :

(Ark; rk)

Метод наискорейшего спуска строится по правилу

xk+1 = xk krk:

Вычислим (rk+1; rk) = (Axk+1 + b; rk) = (A(xk krk) + (b; rk)) = Axk + bk(Ark; rk) = (rk; rk) k(Ark; rk) = 0:

Таким образом, в методе наискорейшего спуска векторы rk и rk+1 ортогональны: (rk+1; rk) = 0: Построим векторы pk = xk xk 1; где вектор xk построен по методу наискорейшего спуска. Тогда

(rk+1; pk) = (rk+1 krk) = k(rk+1; rk) = 0.

Тогда в методе наискорейшего спуска векторы rk+1 и pk ортогональны:

(rk+1; pk) = 0:

Сходимость метода наискорейшего спуска.

Рассмотрим последовательность векторов fxkg; полученную по методу наискорейшего спуска с начальным вектором x0 :

xk+1 = xk krk:

Ясно, что

f(xk+1) f(x ) f(xk) f(x ):

Построим вектор yk+1 = xk 0rk, где 0 = 1+2 n ; т.е. yk+1 строим согласно градиентному методу с оптимальным значением = 0: Ясно, что

f(xk+1) f(x ) f(yk+1) f(x ):

Значение f(yk+1) f(x ) = 12 (A(yk+1 x ); yk+1 x ); yk+1 x = xk x 0A(xk x ) = (E 0A)(xk x );

A(yk+1 x ) = (E 0A)A(xk x ):

Вектор xk x представим в виде линейной комбинации собственных векторов 'i матрицы A:

f(xk+1) f(x ) q2(f(xk) f(x )) q4(f(xk 1) f(x )) ::: q2(k+1)[f(x0) f(x )]

Тогда

f(xk) f(x ) q2k[f(x0) f(x )];

и последовательность f(xk) сходится к f(x ) со скоростью геометрической прогрессии со знаменателем q2,

где

q = n 1n + 1

Апостериорная оценка f(xk) f(x ).

25

f(xk) f(x ) = f(xk) f(xk+1) + f(xk+1) f(x ) f(xk) f(xk+1)+ +q2(f(xk) f(x ))

и

4.4 Метод сопряженных градиентов.

Пусть x0 начальный вектор. Вектор x1 построим по методу наискорейшего спуска:

x1 = x0 0r0 (если r0 отличен от нулевого вектора, то 0 6= 0), и вычислим векторы r1 = Ax1 + b, p1 = x1 x0. Эти векторы ортогональны. Любой вектор и в частности вектор x x1 можно представить в виде ортогонального разложения x x1 = r1 + p1 + z, где z ? p1, z ? r1.

Метод сопряженных градиентов начинается с минимизации значений функции f(x) на гиперплоскости x = x1 r1 + p1. Значения функции

1

f(x) = f(x1) + (r1; r1 + p1) + 2(A( r1 + p1); r1 + p1) = '( ; ):

Естественно предполагать, что векторы r1 и p1 отличны от нулевого вектора. Ясно, что минимум функции f(x) на гиперплоскости существует и единственен. Условия минимума функции f(x) :

Решение системы существует и единственно: = 1, = 1 и величина 1 > 0: Обозначим векторы

x2 = x1 r1 + p1

r2 = Ax2 + b = Ax1 + b 1Ar1 + 1Ap1 = r1 1Ar1 + Ap1

26

и в силу системы получаем свойства векторов r2 и p1:

Таким образом векторы r1, r1, r2 попарно ортогональны. Для векторов p2 и p1 получаем:

(Ap1; p2) = (Ap2; p1) = (r2; p1) (r1; p1) = (r1; p1) = 0; (Ap2; p2) > 0:

В методе сопряженных градиентов последовательно строятся векторы xk+1 = xk + krk + kpk; k = 1; 2; : : : :

Если окажется, что вектор rk нулевой, то Axk = b и xk = x : В этом случае итерационный процесс заканчивается.

Коэффициенты k и k находятся из условия минимума значений функции f(xk + rk + pk): Из получаемой системы двух линейных уравнений для и получаем единственное решение = k и = k: Ясно, что значение k > 0: Свойства этого решения можно записать в виде: (rk+1; rk) = 0, (rk+1; pk) = 0:

Теорема. Ненулевые векторы r0; r1; r2; :::; rk попарно ортогональны.

Доказательство. Предположим, что теорема верна для r0; r1; r2; :::; rk: (при k = 2 она верна), построим вектор rk+1:

На предыдущих шагах (при i = 1; 2; :::; k 1) ri+1 = ri iAri+ + i(ri ri 1):

Так как значение i 6= 0, то Ari = 1 (ri+1 ri i(ri ri 1)); т.е. вектор Ari есть линейная комбинация

i

векторов (ri 1; ri; ri+1):

Вычислим (rk+1; ; ri) для i = (1; 2; :::; k 2) : (ri; rk+1) = (ri; rk) k(ri; Ark) + k(ri; rk rk 1) =k(ri; Ark) = k(rk; Ari) = k(rk; L(ri 1; ri; ri+1)) =

= 0:

При i = 0 величина (rk+1; r0) = k(r0; Ark) + k(r0; rk rr 1) =

= k(rk; r1 r0 ) = 0 и вектор rk+1 ортогонален всем векторам r0; r1; :::; rk 2:

0

При i = k согласно свойству решения системы уравнений для k; k получаем ортогональность векторов kk+1 и rk.

Остается показать, что (rk+1; rk 1) = 0: Вектор pk = k 1 + k 1pk 1 = = : : : = L(r0; r1; :::; rk 1): Согласно свойству системы для k и k получаем

0= (rk+1;pk ) = (rk+1; L(r0; r1; : : : ; rk 1)) = k 1(rk+1; rk 1)

ивекторы r0; r1; : : : ; rk; rk+1 попарно ортогональны. Теорема доказана.

Ввекторном пространстве Vn существует не более n линейно независимых векторов. Если векторы r0; r1; : : : ; rn 1 ненулевые, то вектор rn- нулевой вектор и вектор xn = x :

Следовательно, метод сопряженных градиентов приводит к вектору x за конечное число итераций, не превосходящее (n 1).

Построенный итерационный процесс позволяет найти решение СЛАУ Ax = b с симметричной положительно определенной матрицей A за конечное число операций и поэтому может быть отнесен к числу точных методов решения таких систем.

Замечание. Для ненулевых векторов pi и pj верно (Api; pj) = 0 при i 6= j, а при i = j (Api; pi) 6= 0:

Действительно (Api; pj) = (ri ri 1; L(r0; r1; : : : ; rj 1)) при i > j. С другой стороны (Api; pj) = (pi; Apj) = (Apj; pi) = 0 при j > i. Значение (Api; pi) k pi k2> 0:

4.5Метод сопряженных направлений.

1.A-ортогональные системы векторов.

Пусть A симметричная положительно определенная матрица размерности n n:

Определение. Система ненулевых векторов fs1; s2; :::; sng векторного пространства Vn называется A- ортогональной, если (Asi; sj) = 0 при i 6= j, а (Asi; si) 6= 0:

Векторы А-ортогональной системы линейно независимы. Действительно, если si = c1sk +c2sm, i 6= k; m, то (Asi; si) = (Asi; c1sk; c2; sm) = 0:

27