Лабораторная работа №7. Итерационные методы решения систем линейных уравнений. Приложения численных методов решения систем линейных уравнений

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

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

Подробное описание метода Зейделя можно посмотреть в книге Самарский А.А. «Введение в численные методы». Для каждого итерационного метода решение системы линейных уравнений находится в результате многократного повторения одной и той же вычислительной процедуры, которая называется итерационной. Каждое выполнение этой процедуры называется итерацией. В каждой итерации из начального значения вектора решений X в результате выполнения итерационной процедуры получается «новое» приближенное значение вектора решений Z. Если итерационный процесс сходится, то для каждой следующей итерации разность между векторами Х и Z постоянно уменьшается и после достаточно большого количества итераций может стать меньше заданной погрешности определения вектора решений.

Если система N линейных алгебраических уравнений записана в виде:

, i =1, …, N (1)

где – вектора неизвестных и правых частей, – матрица коэффициентов.

Метод Зейделя применяется в одной из двух форм:

, a_ij ≠0, i =1, 2, …, N (2)

или

, a_ij ≠0, i =1, 2, …, N (3)

В расчетных формулах метода Зейделя (вариант (2), и (3)) Z – «новое» значение вектора неизвестных, а Х – начальное значение вектора Х для каждой итерации.

Рассмотрим вычислительную процедуру (2). Распишем формулу (2) для каждой строки i.

Для i =1:

,

Из этого уравнения легко находится первый элемент «нового» вектора решений z ₁.

Для i =2

,

из этого уравнения легко находится второй элемент «нового» вектора решений z₂,

Для i =3

, и т.д.

Для i = k

.

Уравнение для строки с номером k можно записать в виде:

, где k =1, …, N-1.

В этом уравнении первое слагаемое – состоит из известных величин матрицы коэффициентов и уже вычисленных элементов вектора решений z_j. Из него получаем значение k -й компоненты вектора решений – z_k.

Все компоненты «нового» вектора решений получаем в результате выполнения следующего алгоритма:

Цикл по k от 1 до N-1

Конец Цикла

 

Распишем суммы в формуле в теле цикла. В результате получим более детализированный алгоритм:

Цикл по k от 1 до N

S1=0

Цикл по j от 1 до k-1

S1=S1+ a _kj

Конец Цикла

S2=0

Цикл по j от k+1 до N

S2=S2+ a _kj x _j

Конец Цикла

=(b_k-S1-S2)/ a_kk

Конец Цикла

 

Этот алгоритм представляет собой тело одной итерации. После получения «нового» значения вектора решений – Z, вычисляется разность векторов решений до и после итерации – d = |X-Z|. Если эта разность больше величины заданной допустимой погрешности – «новое» значение вектора решений присваивается начальному: X = Z и итерационный процесс продолжается. Начальные значение вектору X задают исходя из смысла задачи. Если d=|X-Z| меньше величины заданной допустимой погрешности – тогда итерационный процесс завершается.

Можно оптимизировать последний алгоритм, объединив циклы для вычислений S1 и S2. Поскольку в цикле для вычисления S1 компоненты получаются с предыдущей итерации, можно заменить на , подобную же замену можно произвести в выражении =(b_k-S1-S2)/ a_kk в теле основного цикла.

Цикл по k от 1 до N

S=0

Цикл по j от 1 до N

Если j≠k

S=S+ a _kj x _j

Конец Если

Конец Цикла

x _k=(b_k-S)/ a _kk

Конец Цикла

 

Теперь можно записать полный алгоритм метода Зейделя с учетом определения максимального количества итераций, требуемой погрешности и т.д.

Алгоритм решения неоднородной системы линейных уравнений методом Зейделя на естественном языке:

n=КоличествоУравнений

ОпределитьМассив a (n,n), b(n), x0(n), x(n)

Max_i = МаксимальноеЗначениеИтераций

ERR=ЗначениеПогрешности

ПолучитьКоэффициенты(@ а, @b, n)

ЗадатьНачальноеЗначение(@X)

 

Цикл по L от 1 до Max_i

Цикл по k от l до N

x0(k)=x(k)

Конец Цикла

 

Цикл по k от 1 до N

S=0

Цикл по j от 1 до N

Если j≠k

S=S+ a (k,j)*x(j)

Конец Если

Конец Цикла

x(k)=(b(k)-S)/ a( k,k)

Конец Цикла

 

DX=0

Цикл по k от 1 до N

DX = DX + (x(k)- Х0 (k))^2

Конец Цикла

 

DX = КвадратныйКорень(DX)

 

Если DX < ERR Тогда

ПечатьРешения (@ x)

ВыходИзЦикла

Иначе

Цикл по k от l до N

x0(k)=x(k)

Конец Цикла

Конец Если

 

Конец Цикла

 

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

Пример реализации алгоритма решения неоднородной системы линейных уравнений методом Зейделя на VFP:

clear

local i, j, k, l, n

local lb, ub, DX

 

 

n = 3

max_i = 100

ErrX = 0.01

 

lb = -0.5

ub = 0.5

dec = 4

 

LOCAL ARRAY a(n,n), b(n), x(n), x0(n)

? "МАТРИЦА КОЕФИЦИЕНТОВ СИСТЕМЫ А(N,N)"

CreateRNDMartix(@a, n, lb, ub, dec, 1)

? "ВЕКТОР ПРАВОЙ ЧАСТИ СИСТЕМЫ B(N)"

CreateRNDVector(@b, n, lb, ub, dec, 1)

? "НАЧАЛЬНЫЕ ЗНАЧЕНИЯ ВЕКТОРА РЕШЕНИЙ X0(N)"

CreateRNDVector(@x0, n, lb, ub, dec, 1)

 

For l = 1 To max_i

 

For k = 1 To n

x(k) = x0(k)

ENDFOR

 

For k = 1 To n

s = 0

For j = 1 To n

If j <> k Then

s = s + a(k, j) * x(j)

EndIf

ENDFOR

x(k) = (b(k) - s) / a(k, k)

ENDFOR

 

DX = 0

For k = 1 To n

DX = DX + (x(k) - x0(k)) ^ 2

ENDFOR

 

DX = Sqrt(DX)

 

? " Итерация l = " + alltrim(STR(l)) + ": DX = " + TRANSFORM(DX,"99.99999999")

If DX < ErrX Then

?

? " РЕШЕНИЕ НАЙДЕНО "

Print_Vector(@x, n, dec)

EXIT For

Else

For k = 1 To n

x0(k) = x(k)

ENDFOR

ENDIF

 

ENDFOR

 

 

procedure CreateRNDMartix

LPARAMETERS a, n, lb, ub, dec, prnt

 

For i = 1 To n

For j = 1 To n

a(i, j) = (ub - lb) * Rand() + lb

a(i, j) = Round(a(i, j), dec)

ENDFOR

ENDFOR

 

If prnt = 1 Then

For i = 1 To n

?

For j = 1 To n

??(a(i, j))

ENDFOR

ENDFOR

 

ENDIF

 

ENDPROC

 

 

procedure CreateRNDVector

LPARAMETERS a, n, lb, ub, dec, prnt

 

For i = 1 To n

a(i) = (ub - lb) * Rand() + lb

a(i) = Round(a(i), dec)

ENDFOR

 

If prnt = 1 Then

For i = 1 To n

? a(i)

ENDFOR

ENDIF

 

ENDPROC

 

procedure Print_Vector

PARAMETERS x, n, dec

local DStr as String

 

DStr = "99."

For i = 1 To dec

DStr = DStr + "9"

ENDFOR

 

For i = 1 To n

? " X(" + ALLTRIM(STR(i)) + ") = " + TRANSFORM(x(i), DStr)

ENDFOR

 

ENDPROC

 

Пример реализации алгоритма решения неоднородной системы линейных уравнений методом Зейделя на VBA:

Sub Zeidel_test()

Dim i, j, k, l, n

Dim lb, ub, DX

 

Dim a(), b(), x(), x0()

 

n = 3

max_i = 100

ErrX = 0.01

 

lb = -0.5

ub = 0.5

dec = 4

 

ReDim a(n, n): ReDim b(n): ReDim x(n): ReDim x0(n)

 

Call CreateRNDMartix(a, n, lb, ub, dec, 1)

Call CreateRNDVector(b, n, lb, ub, dec, 1)

 

Call CreateRNDVector(x0, n, lb, ub, dec, 1)

 

For l = 1 To max_i

 

For k = 1 To n

x(k) = x0(k)

Next k

 

For k = 1 To n

s = 0

For j = 1 To n

If j <> k Then

s = s + a(k, j) * x(j)

End If

Next j

x(k) = (b(k) - s) / a(k, k)

Next k

 

DX = 0

For k = 1 To n

DX = DX + (x(k) - x0(k)) ^ 2

Next k

 

DX = Sqr(DX)

 

Debug.Print " Итерация l = " & l & ": DX = " & DX

If DX < ErrX Then

Debug.Print

Debug.Print " РЕШЕНИЕ НАЙДЕНО "

Call Print_Vector(x(), n, dec)

Exit For

Else

For k = 1 To n

x0(k) = x(k)

Next k

End If

 

Next l

 

 

End Sub

 

Sub CreateRNDMartix(a, n, lb, ub, dec, prnt)

Dim str

 

For i = 1 To n

For j = 1 To n

a(i, j) = (ub - lb) * Rnd() + lb

a(i, j) = Round(a(i, j), dec)

Next j

Next i

 

If prnt = 1 Then

str = ""

For i = 1 To n

For j = 1 To n

str = str & a(i, j) & ": "

Next j

str = str & vbCrLf

Next i

Debug.Print str

End If

 

End Sub

 

Sub CreateRNDVector(a, n, lb, ub, dec, prnt)

For i = 1 To n

a(i) = (ub - lb) * Rnd() + lb

a(i) = Round(a(i), dec)

Next i

 

If prnt = 1 Then

For i = 1 To n

Debug.Print a(i)

Next i

End If

 

End Sub

 

Sub Print_Vector(x(), n, dec)

Dim i

Dim DStr, str

 

DStr = "0."

For i = 1 To dec

DStr = DStr + "0"

Next i

 

str = ""

For i = 1 To n

str = str & "X(" & i & ") = " & Format(x(i), DStr) & vbCrLf

Next i

 

Debug.Print str

 

End Sub

 

Результат работы программ:

МАТРИЦА КОЕФИЦИЕНТОВ СИСТЕМЫ А(N,N)

0,4194: 0,1317: 0,1276:

-0,0715: -0,402: 0,061:

0,1945: 0,4137: 0,3348:

 

ВЕКТОР ПРАВОЙ ЧАСТИ СИСТЕМЫ B(N)

-0,4774

0,0434

0,4162

НАЧАЛЬНЫЕ ЗНАЧЕНИЯ ВЕКТОРА РЕШЕНИЙ X0(N)

-0,0697

0,1779

0,0025

Итерация l = 1: DX = 2,12836329279661

Итерация l = 2: DX = 0,658840965997187

Итерация l = 3: DX = 8,86987750634338E-02

Итерация l = 4: DX = 1,74853648883518E-02

Итерация l = 5: DX = 2,35213495124003E-03

 

РЕШЕНИЕ НАЙДЕНО

X(1) = -1,8062

X(2) = 0,4723

X(3) = 1,7088

 

Замечание: данные примеры реализаций алгоритмов используют системы уравнений, коэффициенты которых сформированы случайным образом, поэтому сходимость алгоритма зависит от «удачно сформированного» вектора начальных значений.

Метод Гаусса-Зейделя

Для решения систем линейных уравнений наряду с прямыми методами решений (метод Гаусса) используются и итерационные методы. Одним из наиболее простых итерационных методов решения систем линейных уравнений является с метод Гаусса-Зейделя.

В соответствии с методом Гаусса-Зейделя исходную систему линейных уравнений

A*X=B, (4)

приводят к виду:

X=B-C*X, (5)

где матрица С=А-Е, Е–единичная матрица.

Аналогичное преобразование используется при нахождении корня нелинейного уравнения методом простой итерации.

Алгоритм метода Гаусса-Зейделя состоит в следующем: на первом этапе задается начальное приближение вектора X (либо его первая компонента), далее по формуле (5) вычисляется новое улучшенное значение X. Итерации повторяются до тех пор, пока не выполнится условие |X_n-1-X_n| < e.

Применение метода Гаусса-Зейделя приводит к решению не для всех матриц А системы (4) как и метод простой итерации. Условие сходимости метода Гаусса-Зейделя состоит в том, что каждый диагональный элемент должен быть по модулю больше либо равен сумме всех недиагональных элементов, причем хотя бы для одного элемента должно выполняться строгое неравенство.

Описание алгоритма решения системы линейных уравнений методом Гаусса-Зейделя:

Задаем точность вычислений tochn и максимальное количество итераций метода Гаусса-Зейделя max_k.

Задаем коэффициенты исходной системы линейных уравнений a (n,n) и правые части b (n).

Задаем начальные значения решений системы x₀ (n).

В соответствии с методом Гаусса-Зейделя вычисляем вспомогательную матрицу c (n,n) =a (n,n) -e (n,n), где e (n,n)=0 при i<>j, и e (n,n)=1 при i=j.

Задаем равными нулю количество итераций метода – k и количество «правильно» найденных значений x (n) – corr.

Задаем условный цикл по k, с условием повторения тела цикла пока количество «правильно» найденных значений x (n) – corr меньшеколичества уравнений n.

Зануляем значений «правильно» найденных значений corr=0 на данной итерации и увеличиваем счетчик цикла k=k+ 1.

В соответствии с методом Гаусса-Зейделя, используя найденные ранее значения вспомогательной матрицы c (n,n), вычисляем произведение матрицы c (n,n) и x_{k- 1}(n).

Вычисляем x_k (1 …n) =b (1 …n)- c (n,n) ×x_{k- 1}(n).

Подсчитываем количество «правильно» (т.е. с заданной точностью) вычисленных значений x_k (1 …n).

Копируем вычисленные значения x_k (1 …n) в x₀ (1 …n) для следующей итерации метода.

Выводим на экран значения номера итерации, количество «правильно» вычисленных x_k (1 …n), а также сами значения x_k (1 …n).

Задаем ветвление: если значение счетчика итераций k равно максимально допущенному количеству max_k задаем значение «правильно» вычисленных x_k (1 …n) – corr равное количеству уравнений системы n.

Конец итераций по k.

Задаем ветвление: если счетчик итераций метода k достиг максимально допустимого количества, то выводим на экран «Решение не найдено за max_k итераций», иначе – вывод на экран сообщения «Решение найдено за k итераций» и вычисленные значения решений системы линейных уравнений - x_k (1 … n).

Конец алгоритма.

 

Пример алгоритма решения системы линейных уравнений методом Гаусса-Зейделя на естественном языке:

n=размерность_системы

dec=количество_знаков_после_запятой

tochn=1/(10^dec-1)

max_k=максимальное_число_итераций

 

Определить Массив a(n,n),c(n,n),b(n),x0(n),x(n)

 

Цикл по i от 1 до n

Цикл по j от 1 до n

a(i,j)=Если(i=j,Случайный(),0.1*Случайный())

c(i,j)=0

Конец Цикла

b(i)=10*Случайный()

Конец Цикла

 

Вызвать Печать_Матрицы(@a,@b,n,dec)

 

Цикл ПО i от 1 до n

x0(i)=b(i)

Конец Цикла

 

Цикл по i от 1 до n

Цикл по j от 1 до n

c(i,i)=a(i,i)-Если(i=j,1,0)

Конец Цикла

Конец Цикла

 

k=0

corr=0

Цикл Пока corr<n

corr=0

k=k+1

Цикл по i от 1 до n

s=0

Цикл по j от 1 до n

s=s+c(i,j)*x0(j)

Конец Цикла

x(i)=b(i)-s

Конец Цикла

 

Цикл по i от 1 до n

Если Модуль(x(i)-x0(i))<tochn Тогда

corr=corr+1

Конец Если

x0(i)=x(i)

Конец Цикла

 

Печать " k= "+k

Печать " correct x "+corr

 

Вызвать Печать_Вектор(@x,n,dec)

 

corr=Если(k>max_k-1,n,corr)

 

Конец Цикла Пока

 

Если k=max_k Тогда

Печать " Решения не найдены за " + max_k+ " итераций "

Иначе

Печать " Необходимая точность достигнута за " + k+ " итераций"

Печать " Решения системы найдены "

 

Вызвать Печать_Вектора(@x,n,dec)

Конец Если

 

Определить Процедуру Печать_Вектора

Параметры x,n,dec

Определить Массив x(n)

 

Цикл по i от 1 до n

Печать x(i)

Конец Цикла

Конец Определения Процедуры

 

Определить Функцию Печать_Матрицы

Параметры a,b,n,dec

Определить Массив a(n,n),b(n)

 

Цикл по i от 1 до n

Цикл по j от 1 до n

Печать a(i,j)

Конец Цикла

Печать ":" + b(i)

Конец Цикла

Возврат

Конец Определения Функции

 

Пример реализации алгоритма решения системы линейных уравнений методом Гаусса-Зейделя на VFP:

CLEAR

 

n=3

dec=3

tochn=1/(10^dec-1)

max_k=500

 

DIMENSION a(n,n),c(n,n),b(n),x0(n),x(n)

 

FOR i=1 TO n

FOR j=1 TO n

a(i,j)=IIF(i=j,RAND(),0.1*RAND())

c(i,j)=0

ENDFOR

b(i)=10*RAND()

ENDFOR

 

Print_Matrix(@a,@b,n,dec)

 

FOR i=1 TO n

x0(i)=b(i)

ENDFOR

 

FOR i=1 TO n

FOR j=1 TO n

c(i,i)=a(i,i)-IIF(i=j,1,0)

ENDFOR

ENDFOR

 

k=0

corr=0

DO WHILE corr<n

corr=0

k=k+1

FOR i=1 TO n

s=0

FOR j=1 TO n

s=s+c(i,j)*x0(j)

ENDFOR

x(i)=b(i)-s

ENDFOR

 

FOR i=1 TO n

IF ABS(x(i)-x0(i))<tochn then

corr=corr+1

ENDIF

x0(i)=x(i)

ENDFOR

 

?

? " k= "+STR(k,5,0),")"

? " correct x "+STR(corr,5,0)

 

Print_Vector(@x,n,dec)

 

corr=IIF(k>max_k-1,n,corr)

 

ENDDO

 

IF k=max_k then

?

? " Решения не найдены за " + STR(max_k,dec,0)+ " итераций "

 

ELSE

?

? " Необходимая точность достигнута за " + STR(k,dec,0)+ " итераций "

? " Решения системы найдены "

 

Print_Vector(@x,n,dec)

ENDIF

 

PROCEDURE Print_Vector

LPARAMETERS x,n,dec

DIMENSION x(n)

?

FOR i=1 TO n

?STR(x(i),dec*2-2,dec)

ENDFOR

ENDPROC

 

FUNCTION Print_Matrix

LPARAMETERS a,b,n,dec

DIMENSION a(n,n),b(n)

 

FOR i=1 TO n

?

FOR j=1 TO n

?? STR(a(i,j),dec*2+2,dec)

ENDFOR

?? ":" + STR(b(i),dec*2+2,dec)

ENDFOR

RETURN

ENDFUNC

В данной реализации значения коэффициентов матрицы системы уравнений a (n,n) и правых частей уравнений b (n) заданы случайным образом, а начальные значения вектора решений равны вектору правой части системы x₀ (n)= b (n).

Пример реализации алгоритма решения системы линейных уравнений методом Гаусса-Зейделя на VBA:

Sub Gauss_Zeidel()

Dim a(), b(), x0(), x()

 

n = 3

dec = 3

tochn = 1 / (10 ^ dec - 1)

max_k = 500

 

ReDim a(n, n): ReDim b(n): ReDim x0(n): ReDim x(n)

 

For i = 1 To n

For j = 1 To n

a(i, j) = IIf(i = j, Rnd(), 0.1 * Rnd())

Next j

b(i) = 10 * Rnd()

Next i

 

Debug.Print "КОЕФИЦИЕНТЫ СИСТЕМЫ А(N,N) И B(N)"

Call Print_Matrix(a(), b(), n, dec)

Debug.Print

 

For i = 1 To n

x0(i) = b(i)

Next i

 

For i = 1 To n

a(i, i) = a(i, i) - 1

Next i

 

k = 0

corr = 0

Do While corr < n

corr = 0

k = k + 1

For i = 1 To n

s = 0

For j = 1 To n

s = s + a(i, j) * x0(j)

Next j

x(i) = b(i) - s

Next i

 

For i = 1 To n

If Abs(x(i) - x0(i)) < tochn Then

corr = corr + 1

End If

x0(i) = x(i)

Next i

 

Debug.Print " Итерация k = " & k & " Найдено " & corr & " из " & n & " решений "

 

'Call Print_Vector(x(), n, dec)

 

corr = IIf(k > max_k - 1, n, corr)

 

Loop

 

If k = max_k Then

Debug.Print

Debug.Print " Решения не найдены " & max_k & " итераций "

 

Else

Debug.Print

Debug.Print " Необходимая точность достигнута за " & k & " итераций "

Debug.Print " Решения системы найдены "

 

Call Print_Vector(x(), n, dec)

End If

 

End Sub

 

 

Sub Print_Vector(x(), n, dec)

Dim i

Dim DStr, str

 

DStr = "0."

For i = 1 To dec

DStr = DStr + "0"

Next i

 

str = ""

For i = 1 To n

str = str & "X(" & i & ") = " & Format(x(i), DStr) & vbCrLf

Next i

 

Debug.Print str

 

End Sub

 

Sub Print_Matrix(a(), b(), n, dec)

Dim i

Dim DStr, str

 

DStr = "0."

For i = 1 To dec

DStr = DStr + "0"

Next i

 

str = ""

For i = 1 To n

str = str & vbCrLf

For j = 1 To n

str = str & Format(a(i, j), DStr) & ": "

Next j

str = str & ": " & Format(b(i), DStr)

Next i

Debug.Print str

End Sub

 

Результат работы программы на VBA:

КОЕФИЦИЕНТЫ СИСТЕМЫ А(N,N) И B(N)

 

0,799: 0,096: 0,094:: 1,143

0,063: 0,599: 0,025:: 8,602

0,044: 0,025: 0,379:: 5,266

 

Итерация k = 1 Найдено 0 из 3 решений

Итерация k = 2 Найдено 0 из 3 решений

Итерация k = 3 Найдено 0 из 3 решений

Итерация k = 4 Найдено 0 из 3 решений

Итерация k = 5 Найдено 0 из 3 решений

Итерация k = 6 Найдено 0 из 3 решений

Итерация k = 7 Найдено 0 из 3 решений

Итерация k = 8 Найдено 0 из 3 решений

Итерация k = 9 Найдено 1 из 3 решений

Итерация k = 10 Найдено 1 из 3 решений

Итерация k = 11 Найдено 1 из 3 решений

Итерация k = 12 Найдено 1 из 3 решений

Итерация k = 13 Найдено 1 из 3 решений

Итерация k = 14 Найдено 1 из 3 решений

Итерация k = 15 Найдено 2 из 3 решений

Итерация k = 16 Найдено 2 из 3 решений

Итерация k = 17 Найдено 2 из 3 решений

Итерация k = 18 Найдено 2 из 3 решений

Итерация k = 19 Найдено 3 из 3 решений

 

Необходимая точность достигнута за 19 итераций

Решения системы найдены

X(1) = -1,800

X(2) = 14,013

X(3) = 13,206