2 семестр / ПОСОБИЕ_ВычМат

31

на разрешающий элемент ak k. Итак, коэффициенты i уравнения на k шаге преобразовываются по формулам:

			a(k −1)
ai(kj) = ai(kj−1)	–		k j		× ai(kk−1)	,	(*)
			ak(kk−1)
	b(k −1)
bi(k ) = bi(k −1) –		k		× ai(kk−1) ,
	ak(kk−1)

k = 1, 2, ..., n – 1, i = k+1, ..., n, j = k, ..., n.

Например, при k = 1 система (1) преобразовывается к эквивалентной системе (2):

a11x1 + a12 x2 + a13x3... + a1n xn = b1

,ü

(1)

(1)

(1)

(1)

,

ï

,

(2)

a22 x2

+ a23 x3

+ ... + a2n xn

= b2

ý

a(1) x

2

+ a(1) x

+ ... + a(1) x

n

= b(1)

,

ï

32

33

3

3n

n

þ

где ai(1j) = ai j –

a1 j

× ai 1 , bi(1) = bi

–

a1 j

× ai 1 , i = 2, ..., n, j = 1, 2, ..., n.

a11

a11

Таким образом неизвестная x1 исключена из уравнений 2, 3, ..., n. При k = 2 система (2) преобразовывается к эквивалентной системе (3):

a(1)

где a(2) = ai(1j) - 2 j i j a22(1)

a11x1 + a12 x2 + a13x3 + ... + a1n xn = b1

,ü

a(1) x

2

+ a(1) x

+ ... + a(1) x

n

= b(1)

,ï

22

23

3

2n

2

ï

(2)

(2)

xn

(2)

ï

,

(3)

a33 x3

+ ... + a3n

= b3

,ý

............................................ï

a(2) x

+ ... + a(2) x

= b(2)

ï

n

,ï

n3

3

nn

n

þ

× ai(12) , bi(2) = bi(1) -

b2(1)

× ai(12) , i=3, ..., n, j = 2, ..., n.

a22(1)

Таким образом неизвестная x2 исключена из уравнений 3, 4, ..., n. Процедура преобразования повторяется n–1 раз с пересчетом коэффициентов по формулам (*). В итоге на шаге k = n–1 получаем систему (4) треугольного вида:

a11x1 + a12 x2 + a13x3 + ... + a1n xn = b1,ü
a(1) x	2	+ a(1) x	+ ... + a	(1) x		n	= b(1)	,ï
22		23 3		2n			2	ï
		(2)	(2)		xn		(2)	ï	(4)
		a33 x3	+ ... + a3n				= b3	,ý
		............................................ï
			(n−1)				(n−1)	ï
								ï
			ann	xn = bn				.þ

На этом прямой ход заканчивается.

Этап 2 носит название обратного хода. Из последнего уравнения системы (4) следует определение

32

b(n−1)

xn = n − .

ann(n 1)

Вычисление значений неизвестных выполняется по формулам

		1		n
	xk =		× (bk(k −1)	- å ak(kj−1) x j ) , k = n–1, ..., 1.
		akk(k −1)
				j =k +1

2.2.1Алгоритм метода Гаусса

1.Определить значение n – размерность системы линейных уравнений.

2.Определить коэффициенты расширенной матрицы An × (n+1) системы линейных уравнений.

3.Прямой ход. Организовать цикл исключения xk неизвестной: k = 1, 2, …, n–1.

3.1.Организовать цикл, определяющий i уравнение, из которого исключается xk неизвестная: i = k+1, k+2, …, n.

3.1.1.Определить значение вспомогательной переменной d = Aik, равной коэффициенту при неизвестной xk в i уравнении.

3.1.2.Организовать цикл преобразования Аij коэффициентов i уравнения при исключении из него xk неизвестной: j = k, k+1, k+2, …, n+1.

3.1.2.1.Преобразование выполнить по формуле

Aij = Aij – Akj/Akk*d.

3.1.3. Конец цикла по параметру j. 3.2.Конец цикла по параметру i.

4.Конец цикла по параметру k (конец прямого хода).

5.Начало обратного хода. Организовать цикл вычисления значения xk неиз-

вестной по формуле

n

xk = (Ak, n+1 – å Akj x j )/Akk, k = n, n-1, n-2, …, 1.

j=k +1

5.1.Организовать цикл по параметру k: k = n, n–1, n–2, …, 1.

5.1.1.Установить значение xk = Ak, n+1 – равное свободному члену в k уравнении.

5.1.2.Организовываем цикл вычисления разности

n

Ak, n+1 – å Akj x j .

j=k +1

5.1.2.1.Организовать цикл по параметру j = k+1, …, n.

5.1.2.1.1.Вычислить xk = xk – Akj*xj.

5.1.2.2.Конец цикла по параметру j.

5.1.3.Вычислить xk = xk / Akk.

5.2.Конец цикла по параметру k (конец обратного хода).

33

6.Печать вектора корней X.

7.Конец вычислений.

2.2.2 Программа решения системы линейных уравнений методом Гаусса

program pGauss; uses crt;

const

n=3; {определяем размерность системы}

{в массиве А записываем расширенную матрицу системы уравнений} A: array [1..n,1..n+1] of real=((1, 0, -1, 0),(-1, -1, 1, 1),(4, 1, 1, 4));

var

X: array [1..n] of real; i,j,k:byte; c,d:real; begin

clrscr;

{начало прямого хода-приведение матрицы к треугольному виду} for k:=1 to n-1 do {начало цикла исключения k-й неизвестной } begin

for i:=k+1 to n do {начало цикла, определяющего i-е уравнение, из которого исключается k-я неизвестная }

begin

d:=A[i,k]; {коэффициент при k-й неизвестной в i-м уравнении }

for j:=k to n+1 do {цикл преобразования коэффициэнтов при исключении k-й неизвестной}

A[i,j]:=A[i,j]-A[k,j]/A[k,k]*d; end;

end;

{начало обратного хода - вычисление значений неизвестных} for k:=n downto 1 do {цикл вычисленя значения k-й неизвестной} begin

X[k]:=A[k,n+1]; {это свободный член в k-м уравнении}

for j:=n downto k+1 do {цикл вычисления значения k-й неизвестной} X[k]:=X[k]-A[k,j]*X[j];

X[k]:=X[k]/A[k,k];

end;

{конец обратного хода; вывод результатов} for k:=1 to n do

writeln('X[',k,']=',X[k]:12:6);

readkey;

end.

34

2.3 Решение систем линейных уравнений методом простой итерации

Пусть дана система (1) линейных уравнений. Если все диагональные коэффициенты aii ¹ 0, i = 1, 2, ..., n, то систему (1) можно представить в виде

	b	n	aij
xi =	i	- å			x j , i= 1, 2, ..., n.	(5)
	a		a	ii
	ii	j =1
		j ¹i

Введем обозначения:

	éx	ù
	ê 1	ú
x =	êx2	ú	,
	êKú
	ê	ú
	ëxn	û

где
ì		0, если i = j,
ï			aij
aij = í	-			, если i ¹ j,
ï			aii
î

éa a	...a	1n	ù
ê 11 12			ú
a = êa21a22	...a2n		ú	,
ê ..................			ú
ê			ú
ëan1an2 ...ann			û

éb1 ù êb ú b = êêK2 úú , êëbn úû

	bi =	bi	, i = 1, 2, ..., n, j = 1, 2, ..., n,	(*)
		a
		ii

и перепишем систему (5) в виде одного матричного уравнения:
x = ax + b.	(6)

Систему (6) будем решать методом последовательных приближений. Возьмем в качестве начального приближения x(0) вектор b и подставим его в правую часть уравнения (6), получим

x(1) = ax(0) + b.

Продолжая аналогичные вычисления, придем к векторной последовательности приближений:

x(k+1) = ax(k) + b, k = 1, 2, ... .

(7)

Достаточные условия сходимости итераций (7) к решению системы содержит следующая теорема:

Если какая-либо каноническая норма матрицы a меньше единицы, то урав- нение (6) имеет единственное решение x, к которому стремится последователь- ность итераций (7) при любом выборе начального приближения x(0).

Замечание. Для матрицы A = (aij)n × m определены следующие три нормы:

– m-норма – максимальная сумма модулей коэффициентов строк

m

||A||m = max å| aij |;

1 £ i £ n j =1

– l-норма – максимальная сумма модулей коэффициентов столбцов

n

||A||l = max å| aij |;

1 £ j £ m i=1

35

– k-норма – корень квадратный из суммы квадратов всех элементов матрицы

n m

||A||k = å åaij2 .

i=1 j=1

Достаточные условия сходимости процесса итераций для неприведенной системы (1) можно представить в виде:

-по m-норме

m
|aii| > å| aij | (i = 1, 2, ..., n),	(8)

j =1 j ¹i

т. е. модули диагональных коэффициентов для каждого уравнения системы (1) больше суммы модулей всех остальных коэффициентов уравнения;

-по l-норме

n

|ajj| > å| aij | (j = 1, 2, ..., m),

i=1 i¹ j

т. е. модули диагональных коэффициентов для каждого уравнения системы (1) больше суммы модулей всех остальных коэффициентов столбца.

Отклонение приближения x(k) от решения ξ по норме не будет превышать ε,

если

||ξ - x(k)|| ≤

α

||x(k) – x(k-1)|| ≤ ε.

1

-

α

В оценке (9) при вычислении нормы

вектора ||x(k)–x(k-1)|| следует представить

как матрицу Xn × 1. Например, для матрицы

	é 0	0,2	- 0,4ù
α =	ê- 0,25	0	0,25	ú
	ê			ú
	ê	- 0,25	0	ú
	ë- 0,25			û

m – норма равна

(9)

его

		3
||α||m =	max	å| aij | =
1 £ i £ 3		j =1
= max {0 + 0.2 + 0.4,0.25 + 0 + 0.25,	0.25 + 0.25 + 0}= max {0.6,0.5,0.5} = 0.6.

Для вектора

é1.6ù

β = êê- 1úú ê 1 ú ë û

m – норма равна

||β||m = max| βi | = max {1.6,1,1} = 1.6.

i

36

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

||ξ - x(k)|| ≤

α

k

||x(1)

– x(0)|| ≤ ε.

1

−

α

В частности, если выбрать

x(0) = β,

то

и

Следовательно,

Откуда

||α||k+1≤

x(1) = αβ + β

||x(1) – x(0)|| = ||αβ|| ≤ ||α|| ||β||.

||ξ

- x(k)|| ≤

α

k +1

||β|| ≤ ε.

1 −

α

ε

(1–||α||),

(k + 1)ln||α||

≤ ln ε + ln(1–||α||) – ln||β||,

β

k >

ln ε + ln(1 −

α

) − ln

β

– 1.

(10)

α

ln

2.3.1Алгоритм метода итераций

1.Ввести коэффициенты расширенной матрицы An × n+1 системы уравнений

(1)и требование к точности – ε.

2.Проверить выполнение достаточных условий сходимости процесса итераций по m-норме. Если условия (8) не выполняются, конец вычислений.

3.Привести систему (1) к виду (6) по формулам (*).

4.В соответствии с неравенством (10) определить число шагов – k, дающих

наверняка ответ с точностью до ε.

xi( 0 ) = βi, i = 1, 2, ..., n.

5.Организовать цикл вычисления векторной последовательности приближе-

ний (7):

6. xi( j +1 ) = α xi( j ) + β, i = 1, 2, ..., n, j = 0, 1, 2, ... k.

6.Выдать результаты вычислений – вектор x.

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

37

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

program IterSLU; uses crt;

const n = 4; {порядок системы уравнений}

e = 0.0000000001; {точность приближения}

a:array[1..n, 1..n+1] of real = {расширенная матрица неприведенной системы} ((-0.77,-0.04,0.21,-0.18,-1.24), (0.45,-1.23,0.06,0,0.88), (0.26,0.34,-1.11,0,-0.62), (0.05,-0.26,0.34,-1.12,1.17));

var

x, x0: array[1..n] of real; {для корней системы уравнений и приближений} i, j, k: byte;

p, s, am, bm, t: real; begin

clrscr;

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

for i:=1 to n do begin p:=0;

for j:=1 to n do {цикл получения суммы модулей недиагональных элементов} if j <> i then p:=p + abs(a[i, j]);

if abs(a[i, i]) < p then

begin writeln('Итерационный процесс расходится'); readkey; Exit; end; end;

{начать процесс итераций} {привести систему к виду x = ax+b} for i:=1 to n do

begin

for j:=1 to n do

if j <> i then a[i, j] := -a[i, j]/a[i, i]; a[i, n+1] := a[i, n+1]/a[i, i]; a[i, i]:=0; end;

{вычислить m-норму матрицы системы уравнений и столбца свободных членов}

for i:=1 to n do begin s:=0;

for j:=1 to n do s:=s + abs(a[i, j]);

if i = 1 then begin am:=s; bm:=abs(a[1, n+1]); end else

begin if s > am then am:=s; if abs(a[i, n+1]) > bm then bm:=abs(a[i, n+1]); end; end;

{определить число итераций для достижения заданной точности}

38

k:=round((ln(e)+ln(1-am)-ln(bm))/ln(am)); {установить начальные приближения равные bi} x0[i]:=a[i, n+1];

{цикл вычисления последующих приближений} t:=0;

repeat

for i:=1 to n do {цикл вычисления нового значения корня xi}

begin s:=0; for j:=1 to n do s:=s + a[i, j]*x0[j]; x[i]:=a[i, n+1] + s; end; t:=t + 1; for i:=1 to n do x0[i]:=x[i];

until t > k; writeln('Получено решение');

for i:=1 to n do write(' x[', i, ']=', x[i]:12:7); writeln; writeln('Проверка решения:');

for i:=1 to n do begin

s:=0; for j:=1 to n do s:=s + a[i, j]*x[j]; s:=s + a[i, n+1]-x[i]; write(s:12:9); end;

writeln; readkey;{задержка выполнения программы до нажатия клавиши} end.

2.4 Решение систем линейных уравнений средствами Excel

Решение систем линейных уравнений в Excel выполним, используя механизм

– циклическая ссылка.

Пример. Решить систему линейных уравнений

–0.77x1 – 0.04x2		+ 0.21x3	– 0.18x4 = –1.24,	(11)
0.45x1	– 1.23x2	+ 0.06x3	= 0.88,
0.26x1	+ 0.34x2	– 1.11x3	= – 0.62,
0.05x1	– 0.26x2 + 0.34x3 – 1.12x4 = 1.17.

Метод итераций реализуется следующим алгоритмом.

2.В меню Сервис выбрать команду Параметры. В появившемся окне Па- раметры выбрать вкладку Вычисления, на которой снять флажок Итерации.

3.Выполните распределение ячеек и запись расчетных формул:

–диапазон A5:E8 – расширенная матрица неприведенной системы уравне-

ний (11);

–диапазон A10:D13 – матрица α приведенной системы уравнений:

0	=-B5/$A$5	=-C5/$A$5	=-D5/$A$5
=-A6/$B$6	0	=-C6/$B$6	=-D6/$B$6
=-A7/$C$7	=-B7/$C$7	0	=-D7/$C$7
=-A8/$D$8	=-B8/$D$8	=-C8/$D$8	0

–диапазон E10:E13 – вектор β свободных членов приведенной системы уравнений:

39

=E5/$A$5 =E6/$B$6 =E7/$C$7 =E8/$D$8

–F5 – требование к точности – значение ε;

–диапазон F10:F13 – сумма модулей элементов строк матрицы α приведенной системы уравнений:

=ABS(A10)+ABS(B10)+ABS(C10)+ABS(D10) =ABS(A11)+ABS(B11)+ABS(C11)+ABS(D11) =ABS(A12)+ABS(B12)+ABS(C12)+ABS(D12) =ABS(A13)+ABS(B13)+ABS(C13)+ABS(D13)

–G10 – формула вычисления m – нормы матрицы α:

=МАКС(F10:F13)

–формула вычисления m – нормы вектора β свободных членов приведенной системы уравнений:

=МАКС(ABS(E10);ABS(E11);ABS(E12);ABS(E13))

–I10 – условие сходимости итерационного процесса:

=ЕСЛИ(G10>1;"нет";"есть")

–G5 – число итераций для достижения заданной точности ε:

=ОКРУГЛ((LN(F5)+LN(1-G10)-LN(H10))/LN(G10);0)

–G16 – счетчик числа итераций:

=G16+1

–диапазон A18:D18 – начальные приближения значений корней:

=ЕСЛИ($G$16=1;E10;A20) =ЕСЛИ($G$16=1;E11;B20) =ЕСЛИ($G$16=1;E12;C20) =ЕСЛИ($G$16=1;E13;D20)

–диапазон A20:D20 – последующие приближения значений корней:

=СУММПРОИЗВ($A10:$D10;$A$18:$D$18)+$E10 =СУММПРОИЗВ($A11:$D11;$A$18:$D$18)+$E11 =СУММПРОИЗВ($A12:$D12;$A$18:$D$18)+$E12 =СУММПРОИЗВ($A13:$D13;$A$18:$D$18)+$E13

–диапазон H5:H8 – проверка решения:

40

=СУММПРОИЗВ(A5:D5;$A$20:$D$20) =СУММПРОИЗВ(A6:D6;$A$20:$D$20) =СУММПРОИЗВ(A7:D7;$A$20:$D$20) =СУММПРОИЗВ(A8:D8;$A$20:$D$20)

4.Выполните расчет:

-в меню Сервис выбрать команду Параметры. В появившемся окне Па- раметры выбрать вкладку Вычисления, на которой установить флажок Итерации;

вполе Предельное число итераций ввести число, обеспечивающее точность вычислений – значение ячейки G5 и нажать кнопку Ок. Выполняется расчет.

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

Замечания:

1.Если Excel регистрирует циклическую ссылку, следует в появившемся окне нажать кнопку Отмена.

2.Если требуется восстановление исходного состояния для начала процесса итераций, то следует на вкладке Вычисления снять флажок Итерации, а затем выполнить двойной щелчок мышкой по ячейке G16. Ячейка открывается для корректировки, формулу не корректируем, а нажимаем на клавишу Enter. Аналогичные действия выполним для ячеек из диапазона A18:D18.

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

1) В меню Сервис выбираем команду Параметры. В появившемся окне Па- раметры выбираем вкладку Вычисления, на которой устанавливаем флажок Ите- рации, в поле Предельное число итераций вводим число 1 и нажимаем кнопку Ok. Выполняется одна итерация.

2)При каждом нажатии клавиши F9 будет выполняться одна итерация.

2.5Решение систем линейных уравнений средствами MathCAD

Если заданы матрица А и вектор В для системы линейных уравнений в мат-

ричной форме А×Х=В, то вектор решения можно получить из очевидного выражения Х=А-1×В.

В MathCAD введена встроенная функция lsolve(А, В), которая возвращает вектор Х для системы линейных уравнений А×Х=В при заданной матрице коэффициентов А и векторе свободных членов В. Известно, что неоднородная система линейных уравнений совместна, если ранг расширенной матрицы равен рангу матрицы системы, т.е. rank(A)= rank(A׀B). Совместная система имеет единственное решение, если rank(A)= rank(A׀B) = n, где n – размерность матрицы А.