�ݺ�ߣ

Методы вычислений
Метод Гаусса

Кафедра теоретической механики
Юдинцев В. В.

Самарский государственный аэрокосмический университет
им. академика С. П. Королёва
(национальный исследовательский университет)
yudintsev@termech.ru

6 марта 2013 г.

Кафедра ТМ (СГАУ) Методы вычислений 6 марта 2013 г. 1 / 34

Содержание

1 Методы решения систем линейных уравнений

2 Треугольные системы
Обратная подстановка
Прямая подстановка

3 ��-разложение
Метода Гаусса
Матричное описание LU разложения

4 Метод Гаусса с выбором ведущего элемента

5 Оценка погрешности
Плохо обусловленные системы
Нормы матриц и векторов

6 Задачи


Методы решения систем линейных уравнений

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


 a11 x1 + a12 x2 + a13 x3 + a14 x4 + a15 x5 + . . . + a1n xn
 = y1
 a21 x1 + a22 x2 + a23 x3 + a24 x4 + a25 x5 + . . . + a2n xn

= y1
. .
 .
 . .
.

an1 x1 + an2 x2 + an3 x3 + an4 x4 + an5 x5 + . . . + ann xn = yn


Матричная форма:
Ax = B
     
a11 a12 a13 . . . a1n x1 y1
a21 a22 a23 . . . a2n  x2  y2 
A= . . . . , x =  . , y =  . 
     
 .. .
. .
. ... . .  .
. .
.
an1 an2 an3 . . . ann xn yn


Методы решения систем линейных уравнений

Методы решения СЛАУ

Точные методы Приближенные методы
Метод Гаусса Метод простых итераций
Метод прогонки Метод Зейделя
LDL-разложение Метод градиентного спуска
Метод квадратного корня Метод сопряженных
градиентов


Треугольные системы

Треугольные системы Обратная подстановка

СЛАУ с верхней треугольной матрицей

aij = 0 для i > j
     
a11 a12 ... a1n x1 y1
 0 a22 ... a2n   x2   y2 
     
. . . . . .
 ... ...  ·  ...  =  ... 
    
 0 . . . an−1,n−1 an−1,n  xn−1  yn−1 
0 ... 0 ann xn yn




     
a11 a12 ... a1n x1 y1
 0 a22 ... a2n   x2   y2 
     
. . . . . .
 ... ...  ·  ...  =  ... 
    
 0 . . . an−1,n−1 an−1,n  xn−1  yn−1 
0 ... 0 ann xn yn

Решение
xn = yn /ann ,
xn−1 = (yn−1 − xn an−1,n )/an−1,n−1 ,
xn−2 = (yn−2 − xn−1 an−2,n−1 − xn an−2,n )/an−2,n−2 ,
...
n
xi = yi − j=i+1 aij xj /aii



Алгоритм

x(n)=y(n)/A(n,n);
for i=n-1:-1:1
x(i)=(y(i)-A(i,i+1:n)*x(i+1:n))/A(i,i);
end


Треугольные системы Прямая подстановка

Нижняя треугольная матрица

Aij = 0 для j > i
     
a11 0 ... 0 x1 y1
 a21 a22 ... 0     
  x2   y2 

... · =
... ... . . .  . . . . . .
an,1 . . . an,n−1 ann xn yn




     
a11 0 ... 0 x1 y1
 a21 a22 ... 0   x2   y2 
 · = 
... ... ... . . .  . . . . . .
an,1 . . . an,n−1 ann xn yn

Решение
x1 = y1 /a11 ,
x2 = (y2 − x1 a21 )/a22 ,
x3 = (y3 − x1 a31 − x2 a32 )/a33 ,
...
i−1
xi = yi − j=1 aij xj /aii



Алгоритм

x(1)=y(1)/A(1,1);
for i=2:n
x(i)=(y(i)-A(i,1:i-1)*x(1:i-1))/A(i,i);
end


��-разложение Метода Гаусса

Идея метода Гаусса

Приведение СЛАУ Ax = b к эквивалентной треугольной системе

3x1 +5x2 = 9
6x1 +7x2 = 4

После умножения первой строки на 2 и вычитания её из второй
строки:
3x1 +5x2 = 9
−3x2 = −14
3 5 1 0 3 5
=
6 7 2 1 0 −3


��-разложение Матричное описание LU разложения

Матрицы преобразования Гаусса
Для n = 2 определим множитель Гаусса τ :
1 0 a1 1 0 a1 a
= = 1 (1)
− a2
a1 1 a2 −τ 1 a2 0
В общем случае вектор множителей Гаусса определяется как:
τ T = (0, . . . , 0, τk+1 , . . . , τn ), τi = ai /ak , (ak = 0) , i = k + 1, . . . n (2)
k

Матрица преобразования Гаусса Mk = E − τ eT k (3)
    
1 ... 0 0 ... 0 a1 a1
. . . ... ... . . . . . . . . . 
  . . .  . . .
   
0 ... 1 0 . . . 0   ak   ak 
Mk a =   =  (4)
0
 . . . −τk+1 1 . . . 0  ak+1   0 
   
. . . ... ... . . . . . . . . .  . . .  . . .
0 . . . −τn 0 ... 1 an 0


Вектор множителей Гаусса
Алгоритм

Вектор множителей Гаусса
function t = gauss(a)
n=length(a);
t=a(2:n)/a(1);
Преобразование Гаусса

Mk C = (E − τ eT )C = C − τ (eT C)
k k (5)

function C = gauss_app(C,t)
n=size(C,1); % Количество строк матрицы C
C(2:n,:)=C(2:n,:)-t*C(1,:);



Приведение к треугольному виду

Mn−1 . . . M1 A = U (6)
 
1 4 7
A = 2 5 8 
3 6 10
   
1 0 0 1 0 0
M1 = −2 1 0 , M2 = 0 1 0
−3 0 1 0 −2 1
   
1 4 7 1 4 7
M1 A = 0 −3 −6  , M2 M1 A = 0 −3 −6
0 −6 −11 0 0 1



Приведение к треугольному виду
Алгоритм

1 на k шаге есть матрица A(k−1) = Mk−1 . . . M1 A верхняя
треугольная с 1 по k − 1 столбец.
2 множители Гаусса для Mk определяются по матрице-столбцу
(k−1)
A(k−1) (k + 1 : n, k), если ведущий элемент akk =0
k=1;
while A(k,k)~=0 && k<=n-1
t=gauss(A(k:n,k));
A(k:n,:)=gauss_app(A(k:n,:),t);
k=k+1;
end



LU разложение

Mn−1 Mn−2 . . . M1 A = U (7)
A = M−1 M−1 . . . M−1 U = LU
1 2 n−1 (8)

Mk = E − τ (k) eT ,
k M−1 = E + τ (k) eT
k k (9)




Теорема
Для матрицы A ∈ Rn×n существует ��-разложение, если

det(A(1 : k, 1 : k)) = 0 для k = 1 : n − 1

Если ��-разложение существует и A не вырождена, тогда
��-разложение единственно и

det A = u11 u22 . . . unn



Исходная система:
Ax = f
��-разложение:
    
1 4 7 1 0 0 1 4 7
A = 2 5 8  = 2 1 0 0 −3 −6
3 6 10 3 2 1 0 0 1
L Ux = f
y

Решение нижней треугольной системы
Ly = f
Решение верхней треугольной системы
Ux = y

Метод Гаусса с выбором ведущего элемента

Несостоятельность метода

Для матрицы
0 1
A=
1 0
��-разложение не существует, т.к. главная подматрица вырождена.

−10−7 x1 + x2 = 1
x1 + 2x2 = 4

Исключая x1 из первого уравнения и подставляя во второе
x2 = (107 + 4)/(107 + 2). C точностью до 7 знач. цифр
x1 = 0.000000, x2 = 1.000000.
Исключая x1 из второго уравнения и подставляя в первое
x2 = (1 + 4 · 10−7 )/(1 + 2 · 10−7 ). С точностью до 7 знач. цифр
x1 = 1.000000, x2 = 2.000000.



Выбор ведущего элемента

1 Выбор главного элемента в столбце: для k = 2 поиск max |aik |
i
2 Выбор главного элемента в строке: для k = 2 поиск max |aki |
i
   
a11 a12 a13 a14 a11 a12 a13 a14
 0 a22 a23 a24   0 a22 a23 a24 
A(1) =
 0
 A(1) = 
a32 a33 a34   0 a32 a33 a34 
0 a42 a43 a44 0 a42 a43 a44

3 Выбор главного элемента в строке и столбце
 
a11 a12 a13 a14
 0 a22 a23 a24 
A(1) = 
 0

a32 a33 a34 
0 a42 a43 a44



Перестановочная матрица
Перестановка строк

Перестановочная матрица – это матрица, отличающаяся от единичной
лишь перестановками строк:
 
0 0 0 1
1 0 0 0
P= 0 0 1 0
 (10)
0 1 0 0

Матрица взаимных перестановок – единичная матрица с
переставленными двумя строками
 
0 0 0 1
0 1 0 0
P= 0
 (11)
0 1 0
1 0 0 0



     
3 17 10 0 0 1 6 18 −12
A = 2 4 −2  → P1 = 0 1 0 → P1 A = 2 4 −2  →
6 18 −12 1 0 0 3 17 10
   
1 0 0 6 18 −12
→ M1 = −1/3 1 0 → M1 P1 A = 0 −2 2 →
−1/2 0 1 0 8 16
   
1 0 0 1 0 0
→ P2 = 0 0 1 → M2 = 0 1 0 →
0 1 0 0 1/4 1
 
6 18 −12
→ M2 P2 M1 P1 A = 0 8 16 
0 0 6


Оценка погрешности

Оценка погрешности Плохо обусловленные системы

Плохо обусловленные задачи

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

x + 10y = 11
100x + 1001y = 1101
Решение x = 1, y = 1

x + 10y = 11.01
100x + 1001y = 1101

Решение x = 11.01, y = 0.00


Оценка погрешности Нормы матриц и векторов

Нормы векторов

Для векторного n−мерного линейного нормированного пространства
1-норма
n
d 1 = |ui |
i=1

2-норма (евклидова)

n 1/2
2
d 2 = |ui |
i=1

∞-норма
d ∞ = max |ui |
i≤i≤n



Нормы матриц
Для векторного n−мерного линейного нормированного пространства
Норма матрицы (подчиненная норме вектора)
Au
A = sup = max Au (12)
u =0 u u =1

Свойства нормы

A+B ≤ A + B (13)
λA = |λ| B (14)
AB ≤ A B (15)
A = 0, тогда и только тогда, когда A = 0 (16)

Норма матрицы A согласована с нормой вектора u, если

Au ≤ A u (17)



Нормы матрицы, согласованные с нормами векторов

n
A 1 = max |aij | (18)
1≤j≤n
i=1
n
A ∞ = max |aij | (19)
1≤i≤n
j=1

A 2 = max λi (AT · A) (20)
1≤i≤n

Для нормы (19):
Au ∞ = max | j aij uj | ≤ max j |aij ||uj | ≤ (max j |aij |)max |uj | =
i i i j
Au ∞
(max j |aij |) u ∞ ⇒ u ∞ ≤ max j |aij |
i i



Число обусловленности матрицы, обусловленность СЛАУ
Пусть правая часть f и невырожденная матрица коэффициентов A
СЛАУ
Ax = f (21)
получили приращения ∆f, ∆A

(A + ∆A)(x + ∆x) = f + ∆f (22)

Тогда, если выполняются условия:
∆A
A = 0, µ < 1, µ = A A−1
A
оценка относительной погрешности решения определяется:

∆x µ ∆f ∆A
≤ + (23)
x 1 − µ ∆A
A
f A



Доказательство

∆x = A−1 (∆f − ∆Ax − ∆A∆x) ⇒
∆x ≤ A−1 ∆f + A−1 ∆A x + A−1 ∆A ∆x
∆A ∆A
∆x ≤ A−1 ∆f
f f + A−1 A A x + A−1 A A ∆x
µ = A−1 A
f
т.к. f = Ax ≤ A x ⇒ A ≤ x
∆A ∆f f ∆A
∆x 1−µ A ≤µ f A +µ A x ≤
∆f ∆A ∆f ∆A
≤µ f x +µ A x =µ f + A x

∆x µ ∆f ∆A
≤ +
x 1 − µ ∆A
A
f A



Число обусловленности матрицы

µ = A−1 A (24)
Число обусловленности характеризует чувствительность решения
СЛАУ к погрешности исходных данных (∆A, ∆f)
µ ≈ 1 . . . 10
µ > 102 – плохо обусловленные системы
Для
x + 10y = 11
100x + 1001y = 1101

1001 −10
A 1 = 1101, A−1 = , A−1 = 1011
1 4

µ = A−1 A > 106


Задачи

Задание 4

1 Напишите программу решения СЛУ методом Гаусса с частичным
(полным* +5 баллов) выбором ведущего элемента.
2 Оцените относительную погрешность решения если известно, что
относительная погрешность каждого элемента правой части не
превышает 10%.
3 Напишите программу разложения матрицы A на нижнюю L и
верхнюю U треугольные: A = LU
4 Напишите программу вычисления определителя матрицы A


Задачи

Список использованных источников

1 Петров И. Б., Лобанов А. И. Лекции по вычислительной
математике: Учебное пособие - М.: Интернет-Университет
Информационных технологий; БИНОМ. Лаборатория знаний,
2006.
2 Дж. Голуб, Ч. Ван Лоун Матричные вычисления: Пер. с англ. -
М.: Мир, 1999.


�ݺ�ߣ

Численные методы решения СЛАУ. Метод Гаусса.

Recommended

More Related Content

What's hot (7)

Similar to Численные методы решения СЛАУ. Метод Гаусса. (17)

More from Theoretical mechanics department (20)

Численные методы решения СЛАУ. Метод Гаусса.