ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 28.05.2024
Просмотров: 25
Скачиваний: 1
Практична робота № 6
Тема: Розв’язування систем лінійних алгебраїчних рівнянь (СЛАР) ітераційними методами.
Ідея ітераційних методів розв’язування СЛАР полягає в послідовному обчисленні векторів рішень, починаючи з деякого вибраного початкового наближення. Якщо ці вектори утворюють послідовність наближень, що сходиться, то розв’язування знайдено. Ефективність ітераційних методів визначається швидкістю збіжності послідовних наближень до розв’язування.
Алгоритм ітераційних методів складається з послідовного виконання трьох кроків:
1. На першому кроці виконується аналіз СЛАР і перетворення її до виду, зручного для ітерації.
2. На другому кроці будується ітераційний процес і задається початкове наближення.
3. На третьому кроці перевіряється збіжність чергового отриманого наближення до точного розв’язку і оцінка похибки отриманого наближення.
3.1 Метод простої ітерації
Алгоритм методу простої ітерації розглянемо на прикладі.
Приклад 1.
Нехай потрібно знайти розв’язок СЛАР.
4х1 – х2 + х3 = 7
4х1 –8х2 + х3 = -21 (1)
-2х1 + х2 + 5х3 = 15
1. Перетворимо систему до виду, зручного для ітерацій.
х1 = (7+ х2 - х3 )/4
х2 = (21+4х1+ х3)/8 (2)
х3 = (15 + 2х1 + х2 )/5
2. Це дозволить побудувати ітераційний процес
х1(k+1) = (7+ х2(k) - х3(k) )/4
х2(k+1) = (21+4х1(k)+ х3(k))/8 (3)
х3(k+1) = (15 + 2х1(k) + х2(k) )/5
Отже, для отримання нових значень вектора наближень ітерація використовує старі значення, отримані на попередньому кроці.
Покажемо, що якщо почати з точки Р0=[1;2;2], то ітерація сходиться до розв’язку х=[2;4;3].
Підставимо значення вектора Р0=[1;2;2] в праву частину кожного рівняння (3) для отримання нових значень.
Виконаємо обчислення значень в Matlab
>> х1=7/4
х1 =1.7500
>> х2=27/8
х2 = 3.3750
>> х3=15/5
х3 = 3
Нове наближення Р1=[1.7500;3.3750;3.00] ближче до розв’язування X=(2;4;3) ніж Р0. Для отримання нових значень наближення в кожне рівняння системи (3) підставляємо координати Р1 і т.д. до отримання розв’язку із заданою точністю.
Формула (3) задає ітераційний процес обчислення нових координат вектора наближень через старі координати.
3. Для перевірки збіжності чергового отриманого наближення Рк до точного розв’язку X=(2;4;3) після кожної ітерації обчислюємо різницю координат векторів Рк і X
||Pк-X||
Обчислення продовжуються поки ця різниця або не досягне заданої похибки eps, або не перевищить заданого числа ітерацій.
Для реалізації алгоритму цього методу в Matlab створимо два m-файли:
Файл програму iter1.m, яка задає вхідні умови і універсальну файл-функцію jaсobi.m, що реалізує метод простої ітерації.
Файл iter1.m
%Розв'язок СЛАР методом простої ітерації
%А – матриця коефіцієнтів початкового рівняння, B – вектор вільних членів, P – вектор початкових наближень, delta – точність обчислень, max – максимальне число ітерацій.
A=[4 -1 1
4 -8 1
-2 1 5];
B=[7;-21;15];
P=[1;2;2];
delta=0.0001;
max=20;
[X]=jacobi(A,B,P,delta,max)
A*X
%Виклик функції, що реалізує ітераційний алгоритм знаходження коренів СЛАР за методом простої ітерації. Результуючий аргумент X – вектор коренів.
[X]=jacobi(A,B,P,eps,max)
A*X
Файл jaсobi.m
% Розв'язок СЛАР методом простої ітерації
function X=jacobi(A,B,P,delta,max)
N=length(B); % N - кількість рівнянь СЛАР
%цикл по заданій кількості ітерацій
for k=1:max
for j=1:N
%цикл обчислення к-го наближення для всіх невідомих xj
% оператор A(j,[1:j-1,j+1:N] вибирає всі елементи в j-й рядку матриці, окрім елементів в j-м стовпці (тобто А(j,j)). X(j) – елементи вектора-рядка X.
X(j)=(B(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j);
end
%для знаходження різниці координат векторів X і P використовуються функції: abs – модуль числа і norm – обчислення норми Евкліда. X' – транспонований вектор-стовпець.
err=abs(norm(X'-P));
relerr=err/(norm(X)+eps);
P=X';
if (err<delta)|(relerr<delta)
break
end
end
X=X';
Розв’язок:
X =
1.9999
3.9999
2.9999
Цей метод може використовуватися для СЛАР певного типу, а саме, для систем, в яких матриця А – строго діагонально домінуюча.
Це означає виконання умови
|akk|> |akj|, j=1,N, j k
Іншими словами, в кожному рядку матриці величина елемента на головній діагоналі повинна перевищувати суму величин решти елементів в рядку.
Дійсно, легко перевірити, що матриця А нашої СЛАР строго домінуюча.
|4|>|-1|+|1|
|-8| > |-1|+|1|
|5|>|-2|+|1|
Приклад 2
Розглянемо приклад іншої СЛАР. Використаємо нашу універсальну функцію jaсobi.m для обчислення кореня, а для встановлення значень вхідних аргументів створимо нову файл-програму iter2.m
Нехай СЛАР задана матрицею А і вектором B виду:
A=[0 b 0 0.3*g k
10*a 0 -a 0
0.1*a+k 2*b -0.1*b 0 b
0 -b -10*с з 0
-a 0 з -20*g 2*g+a]
B=[ 4*р; -10; 20*р; 4; -4*р]
Встановимо значення змінних a,b,k,g,р:
|
a |
b |
с |
d |
k |
р |
g |
|
0.1 |
-0.1 |
0.2 |
1 |
-1 |
1 |
0.1 |
Обчисливши елементи А і В отримаємо:
= [0 -0.1 0 0.03 -1
1 0.1 0 -0.1 0
-0.99 -0.2 0.01 0 -0.1
0 0.1 -2 0.2 0
-0.1 0 0.2 -2 0.3]
B =[4; -10; 20; 4; -4]
1. Аналіз матриці: елемент a11=0. Переставимо рядки 1 і 2 в матриці і відповідні елементи у векторі B:
= [1 0.1 0 -0.1 0
0 -0.1 0 0.03 -1
0.99 -0.2 0.01 0 -0.1
0 0.1 -2 0.2 0
-0.1 0 0.2 -2 0.3]
B =[ -10; 4; 20; 4; -4]
2. Визначимо, чи є А строго діагонально домінуючою.
|1|>|0.1|+|-0.1|
|-0.1| <|0.03|+|-1|
|0.01| <|-0.99|+| -0.2|+|0.1|
|0.2| <|0.1| +|-2|
|0.3| <|-0.1|+|0.2|+|-2|
Матриця не діагонально домінуюча, отже, метод простої ітерації може не сходитися до розв’язку.
3. Задамо у файлі iter2.m вхідні умови для виклику функції jaсobi.m
%Розв’язок СЛАР методом простої ітерації
a=0.1;
b=-0.1;
с=0.2;
d=1;
k=-1;
р=1;
g=0.1;
A=[10*a 0 -a 0
0 b 0 0.3*g k
0.1*a+k 2*b -0.1*b 0 b
0 -b -10*с с 0
-a 0 с -20*g 2*g+a]
B=[-10; 4*р;20*р;4;-4*р]
delta=0.0001;
max=120;
P=[-10; 4; 20; 4; -4];
[X]=jacobi(A,B,P,delta,max)
Або введемо обчислені елементи для А і В:
% iter2.m Розв’язування СЛАР методом простої ітерації
=[1.0000 0.1000 0 -0.1000 0
0 -0.1000 0 0.0300 -1.0000
-0.9900 -0.2000 0.0100 0 -0.1000
0 0.1000 -2.0000 0.2000 0
-0.1000 0 0.2000 -2.0000 0.3000];
B=[-10;4; 20;4;-4];
delta=0.0001;
max=20;
P=[-10; 4; 20; 4; -4];
[X]=jacobi(A,B,P,delta,max)
4. Виконаємо програму iter2.m
>> iter2.m
В результаті отримаємо значення вектора рішень X, яке значно відхиляється від точного розв’язку.
X =
1.0e+023 *
0.0108
0.3317
0.7633
-2.0562
1.1883
5. Ітерації розв’язування легко отримати, якщо видалити ‘;’ в операторі P=X’
%Розв’язок СЛАР методом простої ітерації
function X=jacobi(A,B,P,delta,max)
N=length(B);
for k=1:max
for j=1:N
X(j)=(B(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j);
end
err=abs(norm(X'-P));
relerr=err/(norm(X)+eps);
P=X'
if (err<delta)|(relerr<delta)
break
end
end
X=X';
5. Перевіримо правильність знайденого розв’язку
ans =
1.0e+023 *
0.2496
-1.2831
-0.1882
-1.9047
4.6205
6. Знайдемо розв’язок прямим методом:
>> X=A\B
X =
1.2209
-109.9381
-7.2698
2.2713
7.0619
Як бачимо, розв’язок, знайдений методом прямої ітерації значно відхиляється від розв’язку, отриманого прямим методом \.
Висновок. У випадку не строго діагонально домінуючої матриці метод прямої ітерації може розходитися.
