Компания ALT Linux
Опубликован: 24.03.2015 | Доступ: свободный | Студентов: 550 / 136 | Длительность: 19:00:00
Лекция 8:

Реализация некоторых численных методов

8.4 Итерационные методы

В приближенных или итерационных методах решение системы линейных алгебраических уравнений является пределом итерационной последовательности, получаемой с помощью этих методов. К ним относятся: метод простой итерации, метод Зейделя и др. Итерационные методы выгодны для системы специального вида, со слабо заполненной матрицей очень большого вида порядка 10^3 \dots 10^5.

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

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

При использовании СКМ Maxima вполне обосновано использование и матричной формулировки итерационных методов.

Рассмотрим решение системы Ax = f (A — квадратная матрица, f — вектор правых частей, x — вектор неизвестных). Обозначим A = L + D + U, где L — нижняя треугольная матрица с нулевыми диагональными элементами; D — диагональная матрица; U — верхняя треугольная матрица с нулевыми диагональными элементами.

Для решения этой системы рассмотрим итерационный процесс

{{x}^{i+1}}={{x}^{i}}-{{H}_{i}}(A{{x}^{i}}-f),
где \det ({{H}_{i}})\ne 0, или
{{x}^{i+1}}={{P}_{i}}{{x}^{i}}+{{d}_{i}},\qquad x(0)={{x}_{0}},
где {{P}_{i}}=I-{{P}_{i}}A — оператор i-го шага итерационного процесса;{{d}_{i}}={{P}_{i}}A.

Итерационный процесс сходящийся, если последовательность \left\{ {{x}_{i}}\right\} сходится к решению {{x}^{*}} при любом x_0.

Если матрица не зависит от номера итерации, итерационный процесс называется стационарным:

{{x}^{i+1}}=P{{x}^{i}}+d. ( 8.4)

Необходимым и достаточным условием сходимости стационарного процесса является выполнение условия \rho (P)<1, где \rho (P) — спектральный радиус матрицы P (наибольшее по модулю собственное число матрицы P).

С использованием введённых обозначений метод простой итерации (метод Якоби) даётся формулой:

P=I-{{D}^{-1}}A, \qquad D = diag({a_{ii}}), \qquad H={{D}^{-1}},
а метод Гаусса–Зейделя — формулой:
P =  - {(D + L)^{ - 1}}U,  \qquad H={{(D+L)}^{-1}}.

Рассмотрим поэлементные расчетные соотношения для методов Якоби и Гаусса—Зейделя.

Все элементы главной диагонали матрицы I = D^{-1}A равны нулю, остальные элементы равны -\cfrac{{a}_{ij}}{{a}_{ii}}, \,\,\, i,j =\overrightarrow {1,n}. Свободный член уравнения (8.4) равен \cfrac{{{f}_{i}}}{{{a}_{ii}}.

Таким образом, для метода Якоби итерационный процесс записывается в виде x^{k+1} = Cx^{k} +, где {{C}_{ij}}=-\cfrac{{{a}_{ij}}}{{{a}_{ii}}};\ k=0,1, \dots, \,\, i,j = \overrightarrow {1,n};\ {E_i} = \cfrac{{{f_i}}}{{{a_{ii}}}}.

Для метода Гаусса–Зейделя x^{k+1}=-(D+L)^{-1}Ux^k + (D+L)^{-1}f, или (D + L){x^{k + 1}} =  - U{x^k} + b$, ${x^{k + 1}} = B{x^{k + 1}} + E{x^k} +e, где

B = \left[ {\begin{array}{*{20}{c}}
  0&0&{\dots}&0 \\ 
  {{c_{21}}}&0&{\dots}&0 \\ 
  {{c_{31}}}&{{c_{32}}}&{\dots}&0 \\ 
  {\dots}&{\dots}&{\dots}&{\dots} \\ 
  {{c_{n - 1,1}}}&{{c_{n - 1,2}}}&{\dots}&0 \\ 
  {{c_{n1}}}&{{c_{n2}}}&{\dots}&0 
\end{array}} \right], 
\text{и}\  
E=\left[ \begin{matrix}
   0 & {{c}_{21}} & \dots & \dots & {{c}_{1n}}  \\
   0 & 0 & \dots & \dots & {{c}_{2n}}  \\
   0 & {} & \ddots  & {} & {}  \\
   \vdots  & {} & {} & \ddots  & {{c}_{n-1,n}}  \\
   0 & {} & {} & {} & 0  \\
\end{matrix} \right].

Рассмотрим решение конкретной системы уравнений Ax = b методом Якоби:

\begin{aligned}
 8{x_1} - {x_2} + 2{x_3} &= 8, \\
 {{x}_{1}}+9{{x}_{2}}+3{{x}_{3}}&=18, \\ 
 2{x_1} - 3{x_2} + 10{x_3} &=  - 5.
\end{aligned}

Вычисляем элементы матрицы B и вектора e:

B=\left[ \begin{matrix}
   0 & \cfrac{1}{8} & -\cfrac{2}{8}  \\
   -\cfrac{1}{9} & 0 & -\cfrac{3}{9}  \\
   -0,2 & 0,3 & 0  \\
\end{matrix} \right],
\quad e = \left[ {\begin{array}{*{20}{c}}
  1 \\ 
  2 \\ 
  {0,5} 
\end{array}} \right] .

Вычислим значения x по формуле x^{k+1} = Bx^k + e. Для решения использована следующая последовательность команд Maxima:

  • преобразование заданных матриц
    (%i1)	A:matrix([8,-1,2],[1,9,3],[2,-3,10])$
    	b:matrix([8],[18],[5])$
    	A0:matrix([A[1,1],A[1,1],A[1,1]],
    		[A[2,2],A[2,2],A[2,2]],
    		[A[3,3],A[3,3],A[3,3]])$
    	B:-A/A0+diagmatrix(3,1)$
    	e:b/matrix([A[1,1]],[A[2,2]],[A[3,3]])$ x:e$
  • собственно вычисление решения
    (%i7)	xt:float(B.x+e)$
    	xt:float(B.xt+e)$
    	xt:float(B.xt+e)$
    	xt:float(B.xt+e)$
    	xt:float(B.xt+e)$
    	xt:float(B.xt+e)$
    	xt:float(B.xt+e)$
    	x0:xt$
    	xt:float(B.xt+e)$
    	x1:xt$
    	r:x1-x0$
    	float(r); /* оценка сходимости*/
    	float(A.x1-b); /* оценка невязки*/
(\%o18)\  \begin{pmatrix}-1.2272477756924971\,{10}^{-5}\cr -1.3018148195786949\,{10}^{-4}\cr 6.2047575160040225\,{10}^{-5}\end{pmatrix}\\
(\%o19)\  \begin{pmatrix}2.5427663227794994\,{10}^{-4}\cr 1.738702477211973\,{10}^{-4}\cr 3.6599949035931445\,{10}^{-4}\end{pmatrix}

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

Для решения системы линейных алгебраических уравнений (8.1) итерационным методом её необходимо привести к нормальному виду:

\bar{x}=P\bar{x}+\bar{g} ( 8.5)

Стационарное итерационное правило получаем, если матрица B и вектор \bar{g} не зависят от номера итерации: {{\bar{x}}^{k+1}}=P{{\bar{x}}^{k}}+\bar{g}. Нестационарное итерационное правило получаем если матрица B или вектор \bar{g} изменяются с ростом номера итерации: {{\bar{x}}^{k+1}}={{P}_{k}}{{\bar{x}}^{k}}+{{\bar{g}}^{k}}.

Стационарное итерационное правило обычно называют методом простой итерации. Предел итерационной последовательности является точным решением системы (8.5) или (8.1).

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

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

  • для того чтобы метод простой итерации сходился, достаточно, чтобы какая-либо норма матрицы P была меньше единицы;
  • для того чтобы метод простой итерации сходился, достаточно, чтобы выполнялось одно из следующих условий:
    • \sum\limits_{j=1}^{n}{\left| {{P}_{ij}} \right|}<1,\ i=\overline{1,n};
    • \sum\limits_{i=1}^{n}{\left| {{P}_{ij}} \right|}<1,\ j=\overline{1,n};
    • {{\sum\limits_{i,j=1}^{n}{\left| {{P}_{ij}} \right|}}^{2}}<1.

Для определения скорости сходимости можно воспользоваться следующей теоремой: если какая-либо норма матрицы P, согласованная с данной нормой вектора, меньше единицы, то имеет место следующая оценка погрешности метода простой итерации:

\left\| {{{\bar{x}}}^{*}}-{{{\bar{x}}}^{k}} \right\|<{{\left\| P
\right\|}^{k}}\cdot \left\| {{{\bar{x}}}^{k}} \right\|+\frac{{{\left\| P
\right\|}^{k}}\cdot \left\| {\bar{g}} \right\|}{1-\left\| P \right\|},
где {\bar{x}}^{*} — точное решение системы (8.1).

Другими словами, условие сходимости выполняется, если выполняется условие доминирования диагональных элементов матрицы исходной системы A по строкам или столбцам: \sum\limits_{i,j=1}^{n}{\left| {{a}_{ij}} \right|}<\left| {{a}_{ii}} \right| или \sum\limits_{i,j=1}^{n}{\left| {{a}_{ij}} \right|}\le \left| {{a}_{jj}} \right|.

В этом случае легко можно перейти от системы вида (8.1) к системе (8.5). Для этого разделим i–ое уравнение системы на a_{i,i} и выразим x_i:

{{x}_{i}}=\frac{{{f}_{i}}}{{{a}_{ii}}}-\frac{{{a}_{11}}}{{{a}_{ii}}}{{x}_{1}}
-\ldots -\frac{{{a}_{1n}}}{{{a}_{ii}}}{{x}_{n}},
т.е. для матрицы P будет выполнено одно из условий сходимости, где
P=\left[ \begin{matrix}
   0 & -\cfrac{{{a}_{12}}}{{{a}_{11}}} & \ldots  &
-\cfrac{{{a}_{1n}}}{{{a}_{11}}}
 \\
   -\cfrac{{{a}_{21}}}{{{a}_{22}}} & 0 & \ldots  &
-\cfrac{{{a}_{2n}}}{{{a}_{22}}}
 \\
   \ldots  & \ldots  & \ldots  & \ldots   \\
   -\cfrac{{{a}_{n1}}}{{{a}_{nn}}} & -\cfrac{{{a}_{n2}}}{{{a}_{nn}}} & \ldots  &
0
 \\
\end{matrix} \right].

Пример реализации метода простой итерации средствами Maxima с печатью промежуточных результатов представлен в скрипте ниже:

(%i1)	iterpr(a0,b0,x,n,eps):=block([a,b,x0,i,j,s,sum,p],
	a:copymatrix (a0), b:copymatrix(b0), x0:copymatrix(x),
	sum:1, p:0,
	while sum>eps do (
		sum:0, p:p+1, print("p= ",p," x= ",float(x)),
		for i:1 thru n do (
			s:b[i,1],
			for j:1 thru n do (s:s-a[i,j]*x0[j,1]),
			s:s/a[i,i], x[i,1]:x0[i,1]+s, sum:sum+abs(s)
		),
		x0:copymatrix(x)
	),
	float(x))$