МВ (методичка)

Рис. 5 Для сходимости этого метода справедлива теорема 1 на стр.16. Скорость сходимости метода хорд

есть скорость геометрической прогрессии:

|xn− ξ|.≤ q·|x0 − ξ|, q<1

Оценка погрешности приближения:

xn −ξ

≤

f (xn )

.

m

Если хорду заменить секущей, которая

строится по двум последним итерациям, получим сле-

дующую итерационную формулу:

xn+1 = xn −

(xn − xn−1) f (xn )

(4)

f (xn ) − f (xn−1)

Этот метод называется методом подвижных хорд или методом секущих. В отличие от рассмотренных ранее методов, он является двухшаговым. Скорость сходимости этого метода выше, чем у метода хорд, но меньше, чем у метода Ньютона.

3.ЧИСЛЕННОЕ РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ.

3.1. Постановка задачи.

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

a

11

x

1

+ a

12

x

2

+... + a

1n

x

n

= b

1

a21 x1 +a22 x2 +... + a2n xn

= b2

................................................

+an2 x2

+... +ann xn

= bn

,

an1 x1

или в матричной форме Ax = b, где A — матрица коэффициентов, b, x — столбец свободных членов и столбец неизвестных соответственно.

Если матрица А невырожденная, т. е. det(A) ≠ 0 , то система имеет единственное решение. В этом случае решение системы с теоретической точки зрения не представляет труда. Значения неизвестных xi (i = 1, 2, ..., n) могут быть получены, по известным формулам Крамера

xi = det(det(AAi)) ,

где матрица Ai,- получается из матрицы А заменой ее i-го столбца столбцом свободных членов. Но такой способ решения линейной системы с п неизвестными приводит к вычислению п + 1 определителей порядка п, что представляет собой весьма трудоемкую операцию при сколько-нибудь

большом числе п.

Применяемые в настоящее время методы решения линейных систем можно разбить на две группы: точные и приближенные. Точными методами называются такие методы, которые в предположении, что вычисления ведутся точно (без округлений), приводят к точным значениям неизвестных за конечное число шагов. Так как на практике все вычисления ведутся с округлениями, то и значения неизвестных, полученные точным методом, неизбежно будут содержать погрешности. К точным методам относится, например, метод Гаусса.

21

Приближенными методами называются такие методы, которые даже в предположении, что вычисления ведутся без округлений, позволяют получить решение системы (х1, х2, . . . , хп)Т за конечное число шагов лишь с заданной точностью. Точное решение системы в этих случаях может быть получено теоретически как результат бесконечного процесса. К приближенным методам относятся: метод простой итерации, метод Зейделя и др. Каждый из этих методов не всегда является сходящимся в применении к конкретному классу систем линейных уравнений.

3.2.Метод Гаусса.

Наиболее распространенным методом решения систем линейных алгебраических уравнений является метод Гаусса, в основе которого лежит идея последовательного исключения неизвестных. Существуют различные вычислительные схемы, реализующие этот метод. Рассмотрим одну из них.

Для простоты ограничимся рассмотрением системы четырех уравнений с четырьмя неизвестны-

ми

a11x1 +a12 x2 +a13 x3 +a14 x4 = b1a21x1 + a22 x2 +a23 x3 +a24 x4 = b2a31x1 + a32 x2 + a33 x3 +a34 x4 = b3a41x1 + a42 x2 +a43 x3 +a44 x4 = b4

Пусть a11 ≠ 0(ведущий элемент). Разделив первое уравнение системы на a11, получим

x1 + a(1)12x2 + a(1)13x3 + a(1)14x4 = b(1)1

где a(1)1j = a1j /a11, j = 2, 3, 4, 5.

Пользуясь первым уравнением, можно исключить неизвестное х1 из второго, третьего и четвертого уравнений системы. Для этого следует умножить первое уравнение на а21, а31 и а41 и вычесть результаты соответственно из второго, третьего и четвертого уравнений системы.

В результате получим систему трех уравнений:

a(1) x			2	+ a(1) x		3	+ a(1) x		4	= b(1)
	22			23			24			2
a32(1) x2				+ a33(1) x3			+ a34(1) x4			= b3(1)
	(1)	x2		(1)	x3		(1)	x4		(1)
a42				+ a43			+ a44			= b4
где коэффициенты aij(2)										вычисляются по формуле

aij(2) = aij − ai1 a(1)1j, (i = 2, 3, 4; j = 2, 3, 4, 5).

Далее первое уравнение полученной системы делим на a22(1), получим систему двух уравнений с двумя неизвестными и т. д. Таким образом, исходную систему мы привели к эквивалентной системе с треугольной матрицей:

	(1)		x2	(1)	x3		(1)	x4		(1)
x1 + a12				+ a13			+ a14			= b4
		x2 + a23(2) x3					+ a24(2) x4 = b2(2)
					x		+ a(3) x			= b(3)
						3			4
							34			3
								x4		= b4
										(4)

откуда последовательно находим

x3 = b3(3) − a34(3)x4, x2 = b2(2) − a24(2)x4 − a23(2)x3,

22

x1 = b1 (1) − a12 (1)x2 − a13 (1)x3 − a14 (1)x4

Для системы n – го порядка коэффициенты вычисляются по формулам:

aij(k) = aij(k - 1) − aik(k - 1)/akk(k - 1)·akj(k - 1); где

aij(0) = aij;a1j(1) = a1j /a11 ;

bi(k) = bi(k - 1) − aik(k - 1)/akk(k - 1)·bk(k - 1), где

i = 1, . . . ,n; j = i+1, . . . ,n; k = 1, . . . ,n

Итак, решение системы распадается на два этапа:

прямой ход − приведение системы к треугольному виду и обратный ход − определение неизвестных по формулам

n

xn = bn(n) ; xi = bi(n) − ∑aij (n) xj ; i = 1, . . . ,n − 1

j =i +1

Очевидно, рассмотренный метод применим лишь при условии, что все «ведущие элементы» отличны от нуля. Если же какой-либо из них обращается в нуль, то в соответствующей системе достаточно провести перестановку уравнений с тем, чтобы сделать «ведущий элемент» отличным от нуля (разумеется, в предположении, что матрица А невырожденная).

Число N арифметических операций, необходимых для реализации метода Гаусса, примерно пропорционально кубу числа неизвестных.

Пример. Методом Гаусса решить систему:

2x1 + x2 + 4x3 =163x1 + 2x2 + x3 =10x1 +3x2 +3x3 =16

Прямой ход реализуется с помощью преобразований:

2		1	4		16
	3	2	1
A(0) =				; b(0)	= 10	,
	1	3	3
					16

				1		1			2
						2										8
A	(1)	=				1			−5			b	(1)
				0							;			=	−14		;
						2
																8
				0		1			1
						2
					1	1			2
																8
						2
									−10			; b(2) =				− 28
A(2) =					0 1													;
					0 0				26							78

23

	1	1	2
							8
		2
	0		−10		; b(3)		− 28
A(3) =		1				=		;
	0	0	1				3

Обратный ход: решаем систему с треугольной матрицей, начиная с последнего уравнения

		1	1	2
						x				8
			2				1
				−10						− 28
		0	1			x2			=
		0	0	1						3
						x3

x3 = 3;

x2 = −28 + 10x3 = 2; x1 = 8 − 0,5x2 − 2x3 = 1.

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

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

Ах = b	(1)
преобразуется к виду
х = В·х + с	(2)
и ее решение находится как предел последовательности
xn + 1 = В·хп + с	(3)
Всякая система
x = x − H·(Ax − b) (4)
имеет вид (2) и при det(H) ≠ 0	эквивалентна системе (1).

В то же время всякая система (2), эквивалентная (1), записывается в виде (4) с матрицей H =(Е − В)·А-1.

При изучении итерационных процессов нам понадобятся понятия норм векторов и матриц. Если

в пространстве векторов x = (x1, . . . , xm )T введена норма

x

, то согласованной с ней нормой || • || в

пространстве матриц A называют норму

A

=

sup

Ax

(5)

x

x ≠

0

Наиболее употребительны в пространстве векторов следующие нормы:

						m
	x					1 = ∑	x j			(6)
	x					j =1			xj
					∞ = max						(7)
						1 ≤ j ≤ m

24

m	2 = (x, x) ,
x 2 = ∑ x j		(8)
j =1

а согласованными с ними нормами в пространстве матриц являются соответственно нормы

m

A1 = 1max≤ j ≤ m ∑ aij i=1

m

A∞ = 1maxj m ∑ aij

≤ ≤ j=1

A	2	=	max λ		( AT A)
			1 ≤ i ≤ m	i

(9)

(10)

(11),

где λ i(A) − собственные значения матрицы A.

Теорема 1. (о необходимом и достаточном условии сходимости метода простой итерации). Пусть система (2) имеет единственное решение. Итерационный процесс (3) сходится к решению системы

(2) при любом начальном приближении тогда и только тогда, когда все собственные значения матрицы В по модулю меньше 1.

Условия теоремы 1 неудобны для практического использования, поэтому при исследовании сходимости метода простой итерации используют условия теоремы 2.

Теорема 2. (о достаточном условии сходимости метода простой итерации). Если ||В|| < 1, то система уравнений (2) имеет единственное решение и итерационный процесс (3) сходится к решению со скоростью геометрической прогрессии. (для ||•||∞ это условие эквивалентно условию диагонального преобладания матрицы А, см. стр.29).

Теорема3.(о погрешности приближений, полученных методом простой итерации).

Если в итерационном процессе ||В|| < 1, то для согласованной с ней нормы вектора x справедлива следующая оценка погрешности метода простой итерации:

x(k ) − x*

≤

B

x(k ) − x(k −1)

(12)

1

−

B

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

Используя соотношение (3) и свойства нормы можно выписать следующую цепочку неравенств:

║x(k) − x*║ = ║Bx(k – 1) + c − Bx* − c║ = ║Bx(k - 1) − Bx*║ =

= ║Bx(k – 1) + Bx(k) − Bx(k) − Bx*║ ≤

≤║Bx(k – 1) − Bx(k) ║ + ║Bx(k) − Bx*║ ≤

≤║B║·║x(k – 1) + x(k) ║ + ║B║·║x(k) − x*║.

Разрешая неравенство относительно ║x(k) − x*║, получаем

x(k ) − x*

≤

B

x(k ) − x(k −1)

.

1

−

B

На основе этой оценки осуществляется выход из итерационного процесса при вычислении реше-

ния с точностью ε:

25

x(k +1) − x(k )

≤

1−

B

ε

B

3.4. Метод Якоби.

Вернемся к рассмотрению задачи (1). Рассмотрим один из способов ее приведения к виду (2). Если на главной диагонали матрицы А нет нулевых элементов, можно выразить xi из i – го уравнения. Получим систему (1) в виде

x1 = −(a12 x2 + a13 x3 +... + a1n xn −b1) / a11

x2 = −(a21x1

+ a23 x3 +... + a2n xn −b2 ) / a21

(13)

LLLLLLLLLLLLLLLLLL

xn = −(an1x1 + an2 x2 +... + an,n−1xn−1 −bn ) / an1

Итерационный процесс будет иметь вид:

x (k +1)

= −(a

x(k ) + a

x(k ) +... + a

x(k ) −b ) / a

1

12

2 13

3

1n

n

1

11

x (k +1)

= −(a

x(k ) + a

x(k ) +... + a

x −b ) / a

21

2

22 1

23

3

2n n

2

(14)

LLLLLLLLLLLLLLLLLL

x (k +1)

= −(a

x(k ) + a

x(k ) +... + a

x(k )

−b ) / a

n1

n

n1 1

n2 2

n,n−1 n−1

n

Этот метод называется методом Якоби. Можно показать, что достаточным условием сходимости этого метода является диагональное преобладание в матрице А системы (1). Диагональное преобладание в матрице А означает, что

			n
	aii		> ∑		aij			(15)
			j=1
			( j≠i)
Для любого i = 1, 2, …, n
Пример. Исследовать								сходимость метода Якоби для системы трех линейных уравнений с тремя

неизвестными и в случае сходимости получить для нее приближенное решение этим методом. 2x1 − x2 + x3 = −3

3x1 + 5x2 − 2x3 = 1 x1 − 4x2 + 10x3 = 0

Построим итерационный процесс следующим образом:

x1(k + 1)	= −1,5 + 0,5 x2(k) − 0,5 x3(k)
x2(k + 1)	= 0,2 − 0,6 x1(k) + 0,4 x3(k)
x3(k + 1)	= −0,1 x1(k) + 0,4 x2(k) ,

Достаточное условие (15) не выполняется для первой строки А, в этом случае требуется вычис-

лить норму матрицы В.
	0	0.5	−0.5			−1.5
	−0.6	0	0.4			0.2
B =					, c =
	−0.1	0.4	0			0

26

B∞ =1 ; B1 = 0.9 , следовательно, метод сходится.

Точное решение этой системы: x1 = −1,(21); x2 = 1,(16); x3 = 0,(58).

Возьмем x0 = с = (−1,5; 0,2; 0)т , тогда x1(1) = −1,5 + 0,5 x2(0) − 0,5 x3(0) = −1,4; x2(1) = 0,2 − 0,6 x1(0) + 0,4 x3(0) = 1,1; x3(1) = −0,1 x1(0) + 0,4 x2(0) = 0,23;

x1(2) = −1,5 + 0,5 x2(1) − 0,5 x3(1) = − 1,065; x2(2) = 0,2 − 0,6 x1(1) + 0,4 x3(1) = 1,132; x3(2) = −0,1 x1(1) + 0,4 x2(1) = 0,58;

и т.д.

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

Витерационном методе Зейделя, в отличие от метода Якоби, при уточнении i – й компоненты

решения используются уже уточненные значения x1, x2, . . . ,xi - 1. Итерационный процесс имеет вид

x (k +1)

= −(a x(k ) + a x(k ) +... + a

x(k )

−b ) / a

1

12 2

13 3

1n n

1

11

x (k +1)

= −(a

x(k +1)

+ a

x(k ) +... + a

x

−b ) / a

(16)

2

22

1

23

3

2n

n

2

21

LLLLLLLLLLLLLLLLLL

x (k +1)

= −(a

x(k +1)

+ a

x(k +1) +...

+ a

n,n−1

x(k +1)

−b ) / a

n1

n

n1

1

n2 2

n−1

n

Можно представить матрицу А в виде суммы трех матриц

A = L + D + R, где L – левая треугольная матрица, R – правая треугольная матрица, D – диагональная матрица. Тогда для получения процесса (16) можно записать (1) в виде

(L+D)x(k + 1) + Rx(k) = b, отсюда имеем
x(k + 1) = −(L+D)−1 Rx(k) + (L+D)−1 b	(17)

Полученный метод называется методом Зейделя. Очевидно, что метод Зейделя эквивалентен ме-

тоду простой итерации (2), где матрица В = −(L+D)−1 R. Потребуем выполнение условий теоремы 1, получим необходимые и достаточные условия сходимости метода Зейделя: все корни уравнения

det(R + (L + D)λ) = 0

(18)

должны быть по модулю меньше 1. Иногда используются более удобные для проверки достаточные условия.

Пример. Исследовать сходимость метода Зейделя для системы трех линейных уравнений с тремя неизвестными и в случае сходимости получить для нее приближенное решение этим методом.

2x1 − x2 + x3 = −3 3x1 + 5x2 − 2x3 = 1 x1 − 4x2 + 10x3 = 0

Построим итерационный процесс следующим образом: x1(k + 1) = −1,5 + 0,5 x2(k) − 0,5 x3(k)

x2(k + 1) = 0,2 − 0,6 x1(k + 1) + 0.4 x3(k)

27

x3(k + 1) = −0,1 x1(k + 1) + 0,4 x2(k + 1).

Здесь:

0		0	0	2 0		0	0	−1 1
L =	3	0	0	, D = 0	5	0	, R = 0	0	− 2
1		− 4	0	0	0	10	0	0	0
Уравнение (18) имеет вид
2λ		−1	1
3λ	5λ		− 2	= 0,

λ− 4λ 10λ

или :

100λ3 − 3λ2 + 2λ = 0,

λ1 = 0, λ2, 3 – корни квадратного уравнения λ2 – 0.03λ + 0,02 = 0.

λ2, 3 − комплексные числа, по модулю меньше 1 и метод Зейделя сходится.

Возьмем x0 = с = (−1,5; 0,2; 0)т, тогда x1(1) = −1,5 + 0,5 x2(0) − 0,5 x3(0) = −1,4;

x2(1) = 0,2 − 0,6 x1(1) + 0,4 x3(0) = 1,04; x3(1) = −0,1 x1(1) + 0,4 x2(1) = 0,556;

x1(2) = −1,5 + 0,5 x2(1) − 0,5 x3(1) = −1,258; x2(2) = 0,2 − 0,6 x1(2) + 0,4·x3(1) = 1,177; x3(2) = −0,1 x1(2) + 0,4 x2(2) = 0,597; и т.д.

4.ИНТЕРПОЛЯЦИОННЫЙ МНОГОЧЛЕН ЛАГРАНЖА.

4.1. Постановка задачи

Пусть известны значения некоторой функции f в п + 1 различных точках х0, х1, ..., хп , которые обозначим следующим образом: f0, f1, ..., fп .

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

Ln (x) = ∑n	a i x i	(1)
i =0

степени п, который в точках xi принимает заданные значения, т. е.

Ln(xi) = fi , i = 0, 1, ..., п .(2)

Этот многочлен называется интерполяционным. Точки xi называются узлами интерполяции. Для удобства изложения под многочленом степени n мы будем подразумевать многочлен степе-

ни не выше п. Например, если fi = 0, i = 0, 1, . . . , п, то интерполяционный многочлен Ln(x) тождественно равен нулю, но его тоже будем называть интерполяционным многочленом n-й степени.

Приближенное восстановление функции f(x) по формуле Ln(x) называется интерполяцией функции f (с помощью алгебраического многочлена). Если х расположен вне минимального отрезка, содержащего все узлы интерполяции х0, x1 ..., хn , то замену функции f по формуле (2) называют также

28

экстраполяцией. Чтобы построить Ln(x), необходимо найти коэффициентыai , которые являются ре-

шением системы

∑n ai x j i = f ( x j ) , j = 0, 1, . . . , п. (3)

i =0

Мы предполагаем, что все xj различны. Определитель этой системы (определитель Вандермонда) отличен от 0. Следовательно, система (3) имеет решение, и притом единственное. Таким образом, доказано существование и единственность интерполяционного многочлена вида (1).

Непосредственное нахождение коэффициентов ai с помощью решения этой системы уже при сравнительно небольших п, например при п = 20, приводит к катастрофическому искажению результатов вычислительной погрешностью. Получим явное представление интерполяционного многочлена, не прибегая к решению системы (3).

1,	при i = j		.
Пусть δij =	при i ≠	j
0,

Задача интерполирования будет решена, если удастся построить многочлены Фi(x) степени не

выше п такие, что
Фi(xj) = δji	при i, j = 0, 1, . . . , п.
	n
Многочлен	Ln (x) = ∑f (xi )Φi (x)	будет	искомым	многочленом.	Действительно,
	i =0
n	n
Ln (x j ) = ∑ f (xi )Φi (x j ) = ∑ f (xi )δij = f (x j ) .
i=0	i=0

Из условий: Фi(xj) = 0 при i ≠ j следует, что

Фi(x) = А(x − x0) . . .(x − xi - 1)(x − xi + 1). . . (x − xn),

аиз условия Фi(xi) = 1 получаем

А= 1/((xi − x0) . . .(xi − xi - 1)(xi − xi + 1). . . (xi − xn)).

n		− x j
Таким образом, Φi (x) =∏	x		.
	x	− x
i ≠ j i		j

Интерполяционный многочлен (1), записанный в форме

n

n

− x j

Ln (x) = ∑ f (xi )∏

x

(4)

x

− x

i=0

i≠ j i

j

называют интерполяционным многочленом Лагранжа.

Пусть n = 1, тогда L

(x) =

f

x − x1

+

x − x0

f

(5)

1

0

x

− x

x

− x

1

0

1

1

0

и это линейная интерполяция (функцию заменяем прямой, проходящей через точки (х0, f(х0)); (x1, f(х1))). Геометрическая интерпретация линейной интерполяции представлена ниже.

29

(X1 ,f(X1))

Ln(x)

f(x)

(X0,f(X0))

																			Рис.6.
При n = 2
L2 (x) =			(x − x1) (x − x2 )										f0 +	(x − x0 ) (x − x2 )							f1 +
			(x		− x )				(x	0	− x	)		(x	− x )		(x	− x	2	)
		0					1				2			1	0	1					(6)
	(x − x0 ) (x − x1)
+											f2
	(x	− x )			(x	2	− x )
2		0					1

и это квадратичная интерполяция

4.2.Погрешность интерполяции.

Всегда можно написать равенство f(x) = Ln(x) + Rn(x), (7)

где Rn(x) — остаточный член, т. е. погрешность интерполяции. Если относительно функции f(x) ничего не известно, кроме ее значений в узлах интерполяции, то никаких полезных суждений относительно остаточного члена Rn(x) сделать нельзя. Ищем Rn(x) в следующем виде:

Rn(x) = Aωn(x), где ωn(x) = (x − x0)·(x − x1)·. . .·(x − xn).

(8)

Рассмотрим функцию F(z) = f(z) − Ln(z) − Rn(z). A найдем так, чтобы F(x) = 0, т.е.

A = f (x) − Ln (x) , где x – точка, в которой оценивается погрешность.

ω (x)

n

Тогда F(z) обращается в 0 в точках x0, x1,,. . . , xn, x , следовательно, имеет n + 2 нуля. По теореме Ролля между двумя нулями дифференцируемой функции лежит хотя бы один нуль ее производной. Таким образом:

F′(z) обращается в 0 n + 1 раз F″(z) обращается в 0 n раз

. . .

F(n + 1)(z) обращается в 0 один раз. Вычислим последовательно

F′(z) = f′(z) − L′n(z) − A((n + 1)zn + …) F″(z) = f″(z) − L″n(z) − A((n + 1)nzn − 1 + … )

. . .

F(n+1)(z) = f(n + 1)(z) − 0 − A(n + 1)!

Следовательно, существует точка ξ такая, что:

F(n+1)(ξ) = f(n+1)(ξ) − A(n+1)! = 0.

30