Вычислительная математика лекции

					2
|| wk || || wk 1	|| c
					1
|| z(k 1) || c*;C		c*
		1			.
2	0	c*
		2

k		2		k 1		2		k		2				2		k
					c			1					c1*			,
		1				1				1				1

З а м е ч а н и е 1. Теорема справедлива для произвольных матриц

при указанном соотношении модулей собственных чисел.

З а м е ч а н и е 2. Если |λ1| > 1, то ||x(k)|| → ∞ при k → ∞ и при вычислениях на ЭВМ возможно переполнение. Если же |λ1| < 1

возможно исчезновение порядка. Для предупреждения этих ситуаций обычно x(k) нормируют так, чтобы ||x(k)|| = 1 для всех k ≥ 1. Возможна

следующая модификация: y(k) =Ax(k-1), λ(k) = (y(k), x(k-1)), x(k) = y(k)/||y(k) ||.

Предполагается, что ||x(0) || = 1. Последовательность y(k) сходится к собственному вектору по направлению.

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

14.5. Метод обратных итераций для вычисления собственных

векторов.

Большинство методов вычисления собственных значений плохо приспособлены для вычисления собственных векторов. Собственные векторы являются ненулевыми решениями вырожденной однородной системы (A – λjE)x = 0. Если вычислено приближенное значение λj,

матрица этой системы оказывается невырожденной и единственное ее

решение равно тривиальному. В методе обратных итераций приближения к собственному вектору определяются последовательным решением систем уравнений (A – λj *E)y(k+1) = x(k) с последующей нормировкой решения x(k+1) = y(k+1)/||y(k+1) ||. В качестве начального приближения x(0) берут произвольный, часто единичный вектор.

Поясним механизм действия метода на примере матрицы простой структуры и в предположении, что λj простое собственное число. Векторы x(0) и y(1) представим в виде линейной комбинации собственных векторов

171

n

x(0) ciei , y(1) i 1

diei .

i 1n

n

Так как ( A *j E) y(1) di (i *j ) , то система уравнений метода

i 1

принимает вид

n

n

di (i *j )ei

ciei

i 1

i 1

Приравнивая коэффициенты при ei , получим

di

ci

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

*

i j

(1)

n

ci

1

n

j

*j

y

ei

(c jej

ciei ).

*

*

*

i 1

j

j

i j

j

i

i

i

Если |λj – λj*| |λi – λj*| для всех i≠j, то второе слагаемое в правой части последней формулы мало по сравнению с первым. Поэтому

y(1)		1	c	e
		*
			j		j
	i	j

оказывается близким по направлению к вектору ej.

Можно показать, что вычисляемый на k – ой итерации вектор

		n			*
			j	j
	x(k ) (k ) (c jej					ciei ) ,
				*
		i j			j
			i

где |β(k)| → |cj-1| при k → ∞. Вектор x(k) сходится к ej по направлению со скоростью геометрической прогрессии, знаменатель которой

		|		* |
	q max		j	j	.
		|		* |
	i j
			i	j

Если |λj – λj*| |λi – λj*| для всех i≠j, то метод обратных итераций сходится очень быстро.

З а м е ч а н и е 1. Так как λj* почти совпадает с λj, то матрица A – λj *E плохо обусловлена и возникают опасения, что большие погрешности при решении систем уравнений повлияют на точность итерационных приближений. К счастью это не так. Возникающая ошибка оказывается почти пропорциональной вектору ej. В данном случае плохая обусловленность не ухудшает, а улучшает ситуацию.

172

З а м е ч а н и е 2. Записав формулу метода в виде

y(k+1) =(A – λj *E)-1 x(k) , замечаем, что метод обратных итераций это степенной метод, примененный к матрице (A – λj *E)-1 .

14.6. Упражнения.

Объектом исследования является метод Якоби для действительных симметрических матриц. Такие матрицы ортогонально подобны диагональной матрице и имеют вещественные собственные значения.

Задача заключается в детализации метода до уровня алгоритма и его программная реализация на языке Matlab.

Как известно, итерационный цикл в методе Якоби состоит в переборе n(n-1)/2 элементов в поддиагональной части матрицы. После выделения очередного элемента aik (i > k) формируют матрицу вращения

Tki, обнуляющую этот элемент при умножении TkiА. Используя матрицу

Tki проводят подобное преобразование TkiАTkiТ. Очередное преобразование не производят, если элемент aik мал. Критерием малости является выполнение условия 100|aik|/|aii| εмаш. Критерием окончания процесса является отсутствие преобразований в цикле.

Алгоритм.

1.Ввод исходных данных: матрица А, заготовка матрицы собственных векторов q1=E, заготовка диагональной матрицы собственных значений D=A, вспомогательная переменная, хранящая число циклов w=0, вспомогательная переменная, подсчитывающая количество отсутствующих преобразований в цикле, s2=0. Критерием окончания процесса является s2≥n(n-1)/2.

2.Пока одновременно выполнены условия s2<n(n-1)/2 и w<wmax,

положить s2=0, w=w+1,проводить итерационный процесс. Процесс заканчивается , если условия перестают выполняться.

В последующих пунктах описываются операции итерационного процесса.

173

3. Осуществить перебор элементов ниже главной диагонали. Для

этого:

3.1Выбрать очередной столбец k=1,..,n-1.

3.2В очередном столбце выбрать очередную строку i=k+1, …,n и

соответствующий элемент aik.

3.2.1. Проверить aik на малость 100|aik|/|aii| εмаш. При выполнении условия малости положить s2=s2+1 и вернутся к п. 3.2.

В противном случае формировать матрицу Tki , осуществлять подобное преобразование, формировать матрицу собственных векторов,

положить s2=0.

c

akk

, s

aik

, t=E, Tki

t.

a2

a2

a2

a2

kk

ik

kk

ik

t

kk

t c, t

ki

s, t

ik

s. D t D tT ,q1 t q1.

ii

Вернуться к п. 3.2.

Ниже приведена соответствующая программа на языке Matlab.

function EiJacob3

a=[1 2 3 4 5;2 8 -7 -2 3 ;3 -7 2 1 5;4 -2 1 7 2;5 3 5 2 0]; %a=[6.25 -1 .5;-1 5 2.12;.5 2.12 3.6];

%a=[1 2 3 4;2 8 -7 -2;3 -7 2 1;4 -2 1 3]; n=length(a);q1=eye(n);s2=0;

w=0;D=a;

while (s2<n*(n-1)/2)&&(w<500) w=w+1;

s2=0;

for k=1:n-1

for i=k+1:n

if 100*(abs(D(i,k)))/abs(D(i,i))<eps s2=s2+1;

else

%Расчет элементов матрицы Tki d=sqrt(D(k,k)^2+D(i,k)^2); c=D(k,k)/d;s=D(i,k)/d;

%Формирование матрицы Tki ( t ) t=eye(n); t(k,k)=c;t(k,i)=s; t(i,i)=c;t(i,k)=-s;

%Расчет собственных чисел (диаг. матрица D) % и собственных векторов (матрица q1')

D=t*D*t';q1=t*q1;

s2=0;

end

end

174

end

end

D,q1',s2,w, [e,L]=eig(a)

Методы минимизации.

15.Введение в минимизацию функций.

Требуется определить минимум вещественной функции n

переменных f(x(1),…,x(n))=f(x), где x n-мерный вектор. Максимум f(x) совпадает с минимумом (-f(x)). Функция f(x(1),…,x(n)) определена на множестве S n-мерного пространства. Если S совпадает со всем n-мерным пространством, задачу минимизации называют безусловной. В противном случае задача имеет ограничения, заданные системой равенств и неравенств. Функцию f(x) часто называют целевой функцией. Если f(x) и

все ограничения линейны, говорят о задаче линейного программирования.

Если хотя бы одна из этих функций не линейна, имеем дело с задачей нелинейного программирования. Задача с ограничениями существенно сложнее безусловной задачи. Ограничимся рассмотрением задачи безусловной минимизации. Обозначим x* точку минимума функции f(x).

Если f(x*)<f(x), когда ||x-x*||< σ, то x* - локальный минимум. Если же f(x*)<f(x), когда x Rn, то x* - глобальный минимум. Будем рассматривать методы поиска локальных минимумов.

15.1.Минимизация функции одной переменной.

15.1.1. Постановка задачи. Определения.

Предположим, что для f(x), определенной на [a,b] имеется единственное значение x=x*, отвечающее минимуму f на[a,b]. При этом для x>x* функция строго возрастает и для x<x* строго убывает. Такая функция называется унимодальной. Минимум может достигаться как внутри промежутка, так и на его краях. В дальнейшем рассматриваются унимодальные функции. Отметим, что используя процедуру локализации,

175

можно выделить промежуток [a,b], внутри которого функция окажется унимодальной.

Если точка минимума x* находится внутри отрезка [a,b],

необходимое условие минимума имеет вид f'(x*)=0. Существование

минимума гарантировано, если одновременно выполнено условие f''(x*)>0.

Обусловленность задачи минимизации.

Пусть значения функции f вычисляются с погрешностью

|f(x)-f*(x)| . Вблизи точки минимума справедливо приближенное равенство f(x)≈f(x*)+(1/2)f''(x*)(x-x*)2.

Оценим минимальное расстояние между точками x и x*, начиная с

которого заведомо будет выполнено неравенство f*(x)>f*(x*). Это расстояние назовем интервалом неопределенности точки x* локального минимума и будем обозначать его ρ.

Имеем f*(x)-f*(x*)=f(x)-f(x*)+[(f*(x)-f(x))-(f*(x*)-f(x*))]≈

≈(1/2)f''(x*)(x-x*)2+[(f*(x)-f(x))-(f*(x*)-f(x*))]≥

≥ (1/2)f''(x*)(x-x*)2-2 .

Следовательно, неравенство f*(x)>f*(x*) выполнено, если

(x-x*)2 4 /f''(x Таким образом, (x-x*)≈2 Δ/ x

Если задача хорошо масштабирована, т е x*~1, |f(x*)|~1, f''(x)~1, то последнее соотношение можно записать в терминах

относительных погрешностей δ(x*)~ δ Отсюда следует, что если δ(f*)~10-m, то δ(x*)~10(-m/2) потеря половины верных значащих цифр в результате

Полученный результат основан на анализе свойств функции f(x). Если для поиска минимума решается нелинейное уравнение f'(x)=0, то радиус неопределенности вычисляется по ранее найденной формуле ρ /|f(2)(x*)|.

15.1.2. Метод деления отрезка пополам.

176

В этом методе используется принцип последовательного сокращения

отрезка локализации, основанный на следующем соображении.

Пусть f унимодальная на отрезке [a,b] функция и a α<γ<β b. Тогда

10)если f(α) f(β), то x* [a,β],

20) если f(α)≥f(β), то x* [α,b],

30) если f(α)≥f(γ) и f(γ) f(β), то x* [α,β].

10. Предположим противное x*>β. Тогда вследствие унимодальности f получим f(α)>f(β), что противоречит условию.

20.Предположим противное x*<α. Тогда вследствие унимодальности f получим f(α)<f(β), что противоречит условию.

30. В силу п.10 имеем x* [a,β], а всилу п.20 x* [α,b]. Следовательно,

x* [α,β].

Дополнительно заметим, если функция унимодальна на отрезке [a,b],

то она унимодальна на любом отрезке [c,d] [a,b].

Обозначим отрезок [a,b] через [a0,b0]. Выберем на отрезке [a0,b0] две симметрично расположенные точки α0=[(a0+b0)/2]-δ, β0= [(a0+b0)/2]+δ.

Здесь δ параметр метода, δ>0. Далее вычисляют f(α0) и f(β0). Сравнение этих значений позволяет сократить отрезок локализации следующим

образом
если f(α0) f(β0), то x	[a1,b1]=[a0,β0];	(	b0 a0	0			);
		1	2		2
если f(α0) >f(β0), то x	[a1,b1]=[α0,b0]. Как и в варианте1					0		;
					1		2

	2δ
a0
α0	β0

b0

Вариант2; ((b0-a0)/2+δ)

Вариант1; ((b0-a0)/2+δ)

177

За очередное приближение к точке минимума в первом случае принимают x1*=α0, во втором x1*=β0. Описанная процедура представляет один шаг итерационной последовательности.

Обозначим через			n=bn-an длину отрезка [an,bn]. Очевидно
n+1= [	n/2]+δ=[(	n-2δ)/2]+2δ. Значения		n образуют цепочку
1=[(	0-2δ)/2]+2δ,	2=[(	1-2δ)/2]+2δ=[(	0-2δ)/22]+2δ,…, n=[( 0-2δ)/2n]+2δ.

Метод половинного деления сходится о скоростью геометрической прогрессии при q=1/2.

n→2δ при n→∞. Если задана погрешность ε определения минимума,

следует положить |xn-x*|<bn-an= n< ε .Причем должно быть ε>2δ.

Рекомендуется задавать δ≤0.2ε.

15.1.3. Метод золотого сечения.

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

Золотым сечением отрезка называют такое его разбиение на две неравные части, что отношение длины всего отрезка к длине его большей части равно отношению длины большей части к длине меньшей части.

Золотое сечение отрезка [a,b] осуществляется каждой из двух симметрично расположенных относительно центра отрезка точек

		2
	a				(b a) 0.382(b a)	,
		3 5

		2
	=a				(b a) 0.618(b a) .
		1 5

178

Вариант 2

0.618			0.382

a	α	β	b
	Вариант 1

Действительно

b a

b

1 5

,

b a

a

1 5

;

b

a

b

a

2

2

Замечательно то, что точка α осуществляет золотое сечение не

только отрезка [a,b], но и отрезка [a,β]

a

a

1 5

.

a

2

Точно также точка β осуществляет золотое сечение не только отрезка

[a,b], но и отрезка [α,b]. Этот факт далее используется.

Очередная (k+1)-я итерация осуществляется аналогично (k+1)-й

итерации метода половинного деления. Новый отрезок локализации определяется по тому же правилу, как и в методе половинного деления.

Если f(αk) f(βk), x	[ak+1,bk+1] =[ak,βk] .
Если f(αk)>f(βk), x	[ak+1,bk+1] =[αk,bk].

Важно то, что какой бы из отрезков [ak,βk] или [αk,bk] не был выбран за очередной отрезок локализации, точка очередного приближения xk+1

совпадает с одной из точек αk или βk (в первом случае xk+1 =αk, а во втором случае xk+1 =βk). Причем, указанные точки уже осуществляют золотое сечение новых отрезков. Требует определения лишь одна дополнительная симметричная точка золотого сечения для нового отрезка. Поэтому

179

достаточно вычислять значение функции лишь в этой одной недостающей точке.

Точка xk+1 отстоит от конца отрезка [ak+1,bk+1] на величину, не

превышающую		2							(bk		ak )					2				k . Таким образом,
		1
					5									1			5
|xk+1-x*|		2				k k 1.Каждая итерация сокращает длину
1				5
																													k
											5 1															2				0 .
отрезка локализации в														раз Поэтому bk-ak= k=
											2															5 1
Таким образом, метод сходится со скоростью геометрической
прогрессии, знаменатель которой																	q				2			0.62.
																							1
																						5
Критерий прекращения итерационного процесса																										k<ε, где
ε – заданная погрешность

15.2.Введение в многомерную минимизацию.

Ограничимся методами безусловной оптимизации.

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

Пусть f(x)=f(x1, x2, …, xn) действительная функция многих переменных, определенная на множестве x Rn.

Точка x* называется точкой глобального минимума, если f(x*)<f(x)

справедливо для всех x Rn. Если неравенство справедливо только в окрестности точки x*, говорят о локальном минимуме. В дальнейшем излагаются методы поиска локальных минимумов.

15.2.2. Поверхности уровня. Градиент и матрица Гессе. Необходимые и достаточные условия локального минимума.

Множество точек, для которых целевая функция принимает постоянные значения f(x)=с, называется поверхностью уровня. В

трехмерном случае для параболоида z=x2+y2 поверхности становятся

линиями уровня – в данном случае окружностями. Для дифференцируемой

180