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

если заменить разложение в ряд Тэйлора f(x)(x – скаляр соответствующим разложением вектор – функции F(x) (x – вектор При этом f'(x заменяется матрицей Якоби Fx, а f''(x матрицей Fxx. Тогда формула метода примет вид

Fx(xk)Δxk+1=- F(xk), xk+1=xk+ Δxk+1.

Алгоритм следует переписать следующим образом

1.Задаемся начальным приближением x0.

2. Вычисляем F(x0 ) и Fx (x0 ).

3.Решаем систему линейных алгебраических уравнений

F(x0 )+Fx (x0 )Δx=0.

4.Находим уточненное приближение x= x0+Δx

5.Проверяем критерий окончания вычислений |x-x0| εзад.

6.Если критерий выполнен, x найденное решение. Иначе x0=x и

переход к п.2.

Доказательство сходимости аналогично скалярному варианту. Метод сходится с квадратичной скоростью. Область сходимости (выбор

начального приближения) определяется условием ||x0- x ||<1/C, где C= c2 , 2c1

||Fx(x)||≥ c1, ||Fxx(x)|| c2, x - точное решение системы.

Критерий останова переписывается следующим образом

||xk-xk-1||<εдоп. .

Возможно использование предложенных ранее модификаций метода.

13.4. Метод простых итераций.

Метод простых итераций легко обобщается для системы уравнений.

Итерационный процесс задается системой разностных уравнений x(k+1)=φ(x(k)) . Условие сходимости принимает вид ||φx ( x )||<1.Критерий окончания итерационного процесса |xn-1-xn| ε(1-q)/q, где q= ||φx ( x )||.

161

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

Проиллюстрируем применение вариантов метода Ньютона – Рафсона для решения систем нелинейных уравнений.

Пример 1.

Ниже приведена программа на языке Matlab, реализующая одну из модификаций алгоритма Ньютона – Рафсона.

1.function nwsys

2.n=2;x=[-.8;8];tol=1e-5;k=0;

3.dx=x;c=.1

4.while ((norm(dx)>tol)&&(k<500))

5.[F1,F2]=mF(x);

6.dx=F2\(-F1);k=k+1;

7.if k>300

a.c=1;

8.end

9.x=x+c*dx;

10.end

11.x,k, [F1,~]=mF(x)

12.function [F1,F2]=mF(x)

a.%F1=[exp(x(1)+x(2))-x(1)^2+x(2)-2;(x(1)+.5)^2+x(2)^2-1];

b.%F2=[x(2)*exp(x(1)+x(2))- 2*x(1),2*(x(1)+.5);x(1)*exp(x(1)+x(2))+1,2*x(2)];

c.F1=[sin(x(1)+x(2))-1.3*x(1)-.1;x(1)^2+x(2)^2-1];

d.F2=[cos(x(1)+x(2))-1.3,2*x(1);cos(x(1)+x(2)),2*x(2)];

e.% F1=[sin(x(1))-x(2)-1.32;cos(x(2))-x(1)+.85];

f.% F2=[cos(x(1)),-1;-1,-sin(x(2))];

Встроке 2 задают размер n системы, начальное приближение и

абсолютную погрешность tol. Выбирая значение переменной "c" в строке 3 c=1 или c<1 можно задать либо основной метод, либо его модификацию,

расширяющую область сходимости. В строке 4 величина "k" ограничивает допустимое число итераций.

В строке 12 задается подфункция, вычисляющая значения компонентов вектор – функции системы (переменная F1) и компонентов матрицы Якоби (переменная F2).

Программа исследования.

1. В строках 2 и 7а задать с=1. Убедится в трудности подбора начального приближения.

162

2. В строках 2 и 7а установить одинаковые значения "c". Варьируя

"c" и "k", отметить факт расширения области сходимости (с=0.1÷1).

3. Установить в строке 2 подобранное в п.2 оптимальное значение

"c", а в строке 7а с=1. Подобрать величину "k" таким образом, чтобы число итераций оказалось минимальным.

Пример 2.

Воспользуемся модификацией метода Ньютона – Рафсона,

использующую конечно – разностную аппроксимацию частных производных в матрице Якоби

f j

(x(k ) )

1

[ f

(x

(k ) ,..., x

(k ) h

(k ) ,..., x (k ) )

x

h (k )

j

1

i

i

n

i

i

f

j

(x (k ) ,..., x (k ) ,..., x

(k ) )];

h (k ) f

(x (k ) ,..., x

(k ) ,..., x (k ) );

1

i

n

i

i

1

i

n

i, j 1, n.

Здесь "k" номер итерации.

При указанном способе выбора приращения hik , модификация называется методом Стеффенсона.

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

1.%Метод Стефенсона

2.function nwsys1

3.n=2;x=[.2;0.8];tol=1e-9;k=0;

4.dx=x;c=.1;

5.while ((norm(dx)>tol)&&(k<500))

6.%Расчет матрицы Якоби

7.for i=1:n %Выбор столбца

i.y=mF(x);x1=x(i);

ii.x(i)=x(i)+y(i);%Расчет приращения очередного аргумента

iii.for j=1:n %Выбор строки

iv.y1=mF(x);%Значение очередной функции при новом аргументе

v.F2(j,i)=(y1(j)-y(j))/y(i);%Значение частной производной

vi.end

vii.x(i)=x1;

8.end

a.%Конец расчета матрицы Якоби

9.F1=mF(x);

10.dx=F2\(-F1);k=k+1;

11.if k>50

a.c=1;

12.end

13.x=x+c*dx;

14.end

163

15.x,k, F1=mF(x)

16.function F1=mF(x)

a.% F1=[exp(x(1)+x(2))-x(1)^2+x(2)-2;(x(1)+.5)^2+x(2)^2-1];

b.% F1=[sin(x(1)+x(2))-1.3*x(1)-.1;x(1)^2+x(2)^2-1];

c.% F1=[sin(x(1))-x(2)-1.32;cos(x(2))-x(1)+.85];

d.% F1=[x(1)^2-x(2)^2+x(3)^2-1.65;

e. % x(1)*sin(x(1))+x(2)*sin(x(2))+x(3)*sin(x(3))-3.1;

f.% x(1)+x(2)+x(3)-sin(x(1)+x(2)+2*x(3))-2.9];

g.% F1=[x(1)^2+x(2)^2+x(3)^2-.3;

h. % cosh(x(1))+cosh(x(2))+cosh(x(3))-3.1;

i.% asin(x(1))+asin(x(2))+asin(x(3))-.3];

j.F1=[sin(x(1)-2*x(2))-x(1)*x(2)+1;x(1)^2-x(2)^2-1];

Воспользовавшись этой программой предлагается провести исследования, аналогичные предложенным в примере1.

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

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

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

Собственным вектором x квадратной матрицы A называется такой ненулевой вектор, воздействие на который матрицы A не меняет его направления Ax=λx.Скаляр λ есть собственное число матрицы А,

соответствующее собственному вектору x. Собственный вектор является решением системы линейных однородных уравнений (A-λE)x=0,

ненулевые решения которого возможны при |A-λE|= χn(λ)=0 . Назовем χn(λ)

характеристическим полиномом. Ограничимся анализом матриц с вещественными коэффициентами A Rn×n . Характеристический полином имеет n корней с учетом их кратности, некоторые из них могут быть комплексными. Такие корни для вещественной матрицы образуют

164

комплексно - сопряженные пары. Каждому собственному числу с алгебраической кратностью α может соответствовать γ собственных векторов, γ α. Назовем γ геометрической кратностью собственного числа.

Если для каждого собственного числа γ=α, говорят о матрице простой структуры. Частными случаями матриц простой структуры являются:

1.Матрицы с попарно различными собственными числами.

2.Симметрические матрицы. У таких матриц дополнительно все собственные числа вещественны.

3.Симметрические положительно определенные матрицы. Для них все собственные числа положительны. Условие положительной определенности (Ax,x)>0 для произвольного вектора x.

Подобным называется преобразование вида B=P-1AP, где P

невырожденная матрица. Подобное преобразование сохраняет собственные значения. Такое преобразование лежит в основе большинства современных методов вычисления собственных чисел. При этом матрицу преобразования P подбирают таким образом, чтобы собственные числа определялись непосредственно по виду матрицы B.

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

Проанализируем обусловленность задачи вычисления собственных чисел. Пусть А* матрица с приближенно заданными элементами aij*≈aij.

Обозначим через λj* собственные числа матрицы А* (j=1,2,…,n).

Теорема1. Пусть А и А* симметрические матрицы, а λj и λj* их собственные значения, упорядоченные по возрастанию. Тогда справедливы следующие формулы для погрешностей

165

max | j * j ||| A A* ||2

1 j n

n

( j * j )2 || A A* ||E

j 1

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

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

Для матриц простой структуры справедлива теорема Теорема2. Пусть P-1AP=D, где D диагональная матрица из

собственных значений А и пусть d=cond2(P)||A-A*||. Тогда каждое собственное значение матрицы А* удалено от некоторого собственного значения матрицы А не более, чем на d.

Перейдем к изложению некоторых методов расчета собственных значений.

14.2.QR алгоритм

Теоретической основой QR – алгоритма служит лемма Шура.

Приведем её формулировку в применении к вещественным матрицам.

Для произвольной квадратной матрицы А существует ортогональная матрица Q, такая, что Q-1AQ есть верхняя блочно треугольная матрица

1.QR алгоритм реализует итерационный процесс нахождения собственных чисел произвольной матрицы. Итерационный шаг состоит из следующих действий:

Ak=QkRr; Ak+1=RkQk; k=0,1,2,…; A0=A – исходная матрица.

Так как Rk=QkTAk, Ak+1=QkTAQk. Таким образом, Ak и Ak+1 подобны 2.Следующая теорема формулирует свойства сходимости.

Пусть все собственные значения различны по абсолютной величине и│λ1│>…>│λn│. Тогда генерируемые QR алгоритмом матрицы сходятся к

166

верхней треугольной. При этом лежащие ниже главной диагонали элементы aij(k) матрицы Ak сходятся к нулю с линейной скоростью,

определяемой отношением собственных значений матрицы A , а именно

││aij(k)│ С│ λi/ λj│k. Заметим, что │ λi/ λj│<1 при i>j.

3.Условия теоремы исключают не только матрицы с кратными

вещественными, но и с комплексными собственными числами. Однако QR

алгоритм работоспособен при более общих предположениях, например,

тогда, когда имеются группы собственных чисел, равных по модулю. В

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

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

2×2.

Если вещественный корень имеет алгебраическую кратность k, то теоретически должен появиться диагональный блок размером k×k. Однако ЭВМ воспринимает кратные λ как различные в пределах вычислительных ошибок. Поэтому на практике, как правило, имеют дело с диагональными блоками размером не более 2×2.

4.Если задача плохо обусловлена QR алгоритм может оказаться неработоспособным.

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

Преобразования осуществляют по следующему алгоритму. Ak+1=Vk AkVkT; k=0,1,…,n-3; A0=A; An-2=H. Vk матрица отражения. Её строят так, чтобы обнулить элементы очередного k-ого столбца, расположенные ниже

главной поддиагонали. Vk=E-2wkwkT; wk=µk(0,0,…0,ak+1,k-s,ak+2,k,…an,k)T.

	n
Здесь первые k компонент нули. s=	ai,k	2	; sign(s)=sign(-ak+1,k);
	i k 1

167

	µk=		1			. Это преобразование численно устойчиво и
			2sk (sk	ak 1,k )

осуществляется за n-2 шагов. В дальнейшем генерируемые QR алгоритмом матрицы остаются матрицами Хессенберга.

Использование матриц Хессенберга снижает число операций на каждой итерации с n3 до n2.

6. Для увеличения скорости сходимости можно использовать сдвиги.

Пусть ηn хорошее приближение к λn. Тогда работают с матрицей H1k=Hk-

ηnE; λ(H1k)= λ(Hk)- ηn. Скорость сходимости определяется соотношением

│ n n │<│ n │, (элемент		a(k ) | n n |k , |λn-ηn| λi-ηn|).
i n	i	ni
				n
			i

По окончании итерации следует осуществить обратный сдвиг – возврат к исходным собственным числам.

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

(алгоритм Фрэнсиса).

7.Критерий останова - малость элементов под главной диагональю: │ai+1,i│<εm(│ai,i│+│ai+1,i+1│); i=1,2,…n-1.

14.3. Метод Якоби для действительных симметричных матриц.

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

168

Назовем итерационным циклом последовательность из n(n 1) 2

шагов (количество циклов k=1,…) , состоящих в переборе элементов alk в

нижней поддиагональной части матрицы, то есть k 1,..., n . l k 1,..., n

На i – ом шаге цикла осуществляется подобное преобразование Ai→Ai+1=TklAiTklT (TklT=Tkl-1). Tkl подбирается так, чтобы обнулить элемент alk матрицы TklAi . Заметим, что при подобном преобразовании симметричность матрицы сохраняется. Поэтому условия обнуления элемента alk совпадают с таковыми для элемента akl . Последнюю задачу мы решали ранее, когда рассматривали QR разложение методом вращения.

Умножение матрицы TklAi на TklT справа может сделать элемент alk

снова ненулевым. Однако доказывается, что Ai →D при k→∞.

Произведение матриц Tkl справа от Ai формируют матрицу, столбцы которой есть собственные векторы. Перед очередным вращением проверяют выбранный элемент alk на малость. Если выполнено условие

100| alk |/|akk| εмаш., полагают элемент малым и вращение не производят.

Критерий окончания нет вращений в цикле, то есть все элементы в поддиагональной части малы.

14.4. Степенной метод.

Возьмем произвольный начальный вектор x0 и построим последовательность векторов, используя формулы

	x(k ) Ax(k 1)
	(k )	(xk , x(k 1) )		.
		(x(k 1)	, x(k 1) )
	1

Справедлива следующая теорема.

Теорема. Пусть А – матрица простой структуры, для которой выполнено условие │λ1│>│λ2│≥│λ3│≥…≥│λn│. Предположим, что в разложении x0 = c1e1+c2e2+…+cnen по базису из собственных векторов c1

169

не равен нулю. Тогда λ1(k)→λ1 при k→∞ и верна оценка относительной погрешности

( (k ) )

| (k )

|

C

k

1

1

2

.

1

| 1 |

1

Доказательство. Вектор x(k) = Akx0. Так как Akei = λikei, то

x(k) = Akx0 = λ1kc1e1+λ2kc2e2+…+λnkcnen .

Положим z(k )

x(k )

и заметим, что z(k) → c1e1

при k → ∞, причем

k

1

n

k

w(k ) z(k ) c1e1

i

ciei .

i 2

1

Так как

i

2

для всех i ≥ 2, то || w(k ) || c

2

k

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

1

1

1

n

i

k

2

k

n

i

k

2

k

|| wk

||

*| ci | *|| ei ||

| ci

| *|| ei || c

.

i 1

1

1

i 2

1

1

n

i

k

c

| ci

| *|| ei || .

2

i 2

Теперь нетрудно установить, что

(k )

(xk , x(k 1) )

(z(k )

, z(k 1) )

(c e ,c e )

при k .

1

1

1

1

(x(k 1) , x(k

(k 1) )

1

1) )

1 (z(k 1) , z

1 (c e ,c e )

1

1

1

1

1

Кроме того

(k )

(z(k )

z(k 1)

, z(k 1) )

(w(k ) w(k 1) , z(k 1) )

1

1

=

.

(z(k 1) , z(k 1) )

(z(k 1) , z(k 1) )

1

Используя неравенство Коши – Буняковского |(x,y)|

||x|| ||y||, имеем

| (w(k ) w(k 1) , z(k 1) ) | || w(k ) w(k 1) || z(k 1) || .

Откуда следует

|

(k )

|

|| w(k ) w(k 1) ||

|| w(k ) || || w(k 1) ||

C

k

1

1

2

.

| |

|| z(k 1) ||

(k 1) ||

|| z

0

1

1

k . Тогда

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

|| w || c

2

k

1

170