13.
Методы отыскания решений систем нелинейных
уравнений.
13.
1. Постановка задачи.
Требуется вычислить корни системы из n уравнений с n неизвестными
1 1
2 2
1 2
1 2
( ,
,...,
)
0
( ,
,...,
)
0
( ,
,...,
)
0
n
n
n
n
f x x
x
f x x
x
f x x
x
Матричная запись F(x)=0, x=(x
1
, x
2
,…,x n
)
T
, F=(f
1
, f
2
,…,f n
)
T
Ограничимся отысканием вещественных корней.
Понятие сходимости легко обобщается для системы, если рассматривать сходимость последовательности векторов по норме
||x k
||
→||
x
||, если k
→∞. Аналогичным образом обобщаются формулы, вводящие понятие скорости и порядка сходимости. Рассуждения об обусловленности скалярной задачи распространяются на систему, если воспользоваться разложением вектор – функции F(x) в ряд Тэйлора и заменить |f'(x)| на норму матрицы Якоби ||F
x
(x)||.
160
13.2. Разложение в ряд Тэйлора вектор - функции F=(f1, f2,…,fn)T. Ряд Тэйлора для скалярной функции f(x) нескольких переменных
(x=(x
1
, x
2
,…,x n
)
T
) с остаточным членом имеет следующий вид f(x)=f(x
1
,x
2
,…,x n
)=f(x
0
) +
0 0
0 1
2 1
2
( )
( )
( )
[
]
nnf xf xf xxxxxxx
+
2 2
2 2
1 1
2 1
2 1
1 2
1 1
( )
( )
( )
[(
(
)
)
2
nnf uf uf uxx xx xxx xx x
…+
2 2
2 2
1 2
1
( )
( )
( )
(
(
)
)
kkknkkknf uf uf uxxxxxx xxx x
…+
2 2
2 1
2 1
( )
( )
(
(
) )]
nnnnf uf uxxxx xx
x=x
0
+Δx,
1 1
0
iinnxuxxx
, 0<θ
i
<1.
∆x i
=x i
-x
0
В матричной форме это можно записать следующим образом f(x)=f(x
0
)+f x
(x
0
)Δx+(1/2)f xx
(u)Δx
2
. Вектор u лежит в окрестности точки x
0
, f x
=( f x1
, f x2
,…, f xn
) – вектор – строка размером (1
×n), f xx
= (f x
)
x
=((f x1
)
x
, (f x2
)
x
,…,(f xn
)
x
) – вектор – строка размером (1
×n
2
),
1 2
2 2
всего n компонентов,
nx xxxxxx
Δx i
– i-ая компонента вектора Δx.
Если речь идет о вектор – функции, то
F(x)=F(x
0
)+F
x
(x
0
)Δx+(1/2)F
xx
(u) Δx
2
, где F
x
(x
0
) – матрица Якоби размера
(n
×n), F
xx
(u) – матрица размера (n
×n
2
). Заметим, что || Δx
2
||
∞
= || Δx ||
2
∞
13.3. Метод Ньютона. Все рассуждения и доказательства, сделанные при
обосновании метода для скалярного варианта, переносятся на систему уравнений,
161 если заменить разложение в ряд Тэйлора f(x)(x – скаляр соответствующим разложением вектор – функции F(x) (x – вектор
При этом f'(x заменяется матрицей Якоби F
x
, а f''(x матрицей F
xx
Тогда формула метода примет вид
F
x
(x k
)Δx k+1
=- F(x k
), x k+1
=x k
+ Δx k+1
Алгоритм следует переписать следующим образом
1.Задаемся начальным приближением x
0 2. Вычисляем F(x
0
) и F
x
(x
0
).
3.Решаем систему линейных алгебраических уравнений
F(x
0
)+F
x
(x
0
)Δx=0.
4.Находим уточненное приближение x= x
0
+Δx
5. Проверяем критерий окончания вычислений |x-x
0
|
ε
зад
6. Если критерий выполнен, x найденное решение. Иначе x
0
=x и переход к п.2.
Доказательство сходимости аналогично скалярному варианту. Метод сходится с квадратичной скоростью. Область сходимости (выбор начального приближения) определяется условием ||x
0
-
x
||<1/C, где C=
2 1
2
c
c
,
||F
x
(x)||
≥ c
1
, ||F
xx
(x)||
c
2
,
x
- точное решение системы.
Критерий останова переписывается следующим образом
||x k
-x k-1
||<ε
доп.
Возможно использование предложенных ранее модификаций метода.
13.4. Метод простых итераций.
Метод простых итераций легко обобщается для системы уравнений.
Итерационный процесс задается системой разностных уравнений x
(k+1)
=φ(x
(k)
) . Условие сходимости принимает вид ||φ
x
(
x
)||<1.Критерий окончания итерационного процесса |x n-1
-x n
|
ε(1-q)/q, где q= ||φ
x
(
x
)||.
162
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. Убедится в трудности подбора начального приближения.
163 2. В строках 2 и 7а установить одинаковые значения "c". Варьируя "c" и "k", отметить факт расширения области сходимости (с=0.1
÷1).
3. Установить в строке 2 подобранное в п.2 оптимальное значение "c", а в строке 7а с=1. Подобрать величину "k" таким образом, чтобы число итераций оказалось минимальным.
Пример 2.
Воспользуемся модификацией метода Ньютона – Рафсона, использующую конечно – разностную аппроксимацию частных производных в матрице Якоби
( )
( )
( )
( )
( )
1
( )
( )
( )
( )
( )
( )
( )
( )
1 1
1
(
)
[
(
,...,
,...,
)
(
,...,
,...,
)];
(
,...,
,...,
);
,
1, .
j
k
k
k
k
k
j
i
i
n
k
i
i
k
k
k
k
k
k
k
j
i
n
i
i
i
n
f
x
f x
x
h
x
x
h
f x
x
x
h
f x
x
x
i j
n
Здесь "k" номер итерации.
При указанном способе выбора приращения h i
k
, модификация называется методом Стеффенсона.
Ниже приведена соответствующая программа на языке 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
164 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
R
n×n
. Характеристический полином имеет n
корней с учетом их кратности, некоторые из них могут быть комплексными. Такие корни для вещественной матрицы образуют
165 комплексно - сопряженные пары. Каждому собственному числу с алгебраической кратностью α может соответствовать γ собственных векторов, γ
α. Назовем γ геометрической кратностью собственного числа.
Если для каждого собственного числа γ
=α, говорят о матрице простой структуры. Частными случаями матриц простой структуры являются:
1. Матрицы с попарно различными собственными числами.
2. Симметрические матрицы. У таких матриц дополнительно все собственные числа вещественны.
3. Симметрические положительно определенные матрицы. Для них все собственные числа положительны. Условие положительной определенности (Ax,x)>0 для произвольного вектора x.
Подобным называется преобразование вида B=P
-1
AP, где P невырожденная матрица. Подобное преобразование сохраняет собственные значения. Такое преобразование лежит в основе большинства современных методов вычисления собственных чисел. При этом матрицу преобразования P подбирают таким образом, чтобы собственные числа определялись непосредственно по виду матрицы B.
Матрица простой структуры подобна диагональной. Произвольная матрица может быть приведена подобным преобразованием по крайней мере к треугольной или блочно треугольной. В последнем случае собственные числа совпадают с собственными числами диагональных блоков.
Проанализируем обусловленность задачи вычисления собственных чисел. Пусть А
*
матрица с приближенно заданными элементами a ij
*
≈a ij
Обозначим через λ
j
*
собственные числа матрицы А
*
(j=1,2,…,n).
Теорема1. Пусть А и А
*
симметрические матрицы, а λ
j
и λ
j
*
их
собственные значения, упорядоченные по возрастанию. Тогда справедливы
следующие формулы для погрешностей
166
*
*
2 1
*
2
*
1
max |
| ||
||
(
)
||
||
jjj nnjjEjAAAA
Теорема1 означает, что задача вычисления собственных значений симметрической матрицы очень хорошо обусловлена.
Встречаются несимметричные матрицы с чрезвычайно плохо обусловленной задачей вычисления собственных значений.
Для матриц простой структуры справедлива теорема
Теорема2.
Пусть P-1AP=D, где D диагональная матрица из собственных значений А и пусть d=cond2(P)||A-A*||. Тогда каждое собственное значение матрицы А* удалено от некоторого собственного значения матрицы А не более, чем на d. Перейдем к изложению некоторых методов расчета собственных значений.
14.2.QR алгоритм Теоретической основой QR – алгоритма служит лемма Шура.
Приведем её формулировку в применении к вещественным матрицам.
Для произвольной квадратной матрицы А существует ортогональная матрица Q, такая, что Q-1AQ есть верхняя блочно треугольная матрица 1.QR алгоритм реализует итерационный процесс нахождения собственных чисел произвольной матрицы. Итерационный шаг состоит из следующих действий:
A
k
=Q
k
R
r
; A
k+1
=R
k
Q
k
; k=0,1,2,…; A
0
=A – исходная матрица.
Так как R
k
=Qk
T
A
k
, A
k+1
=Q
k
T
AQ
k
. Таким образом, A
k и A
k+1
подобны
2.Следующая теорема формулирует свойства сходимости.
Пусть все собственные значения различны по абсолютной величине и│λ
1
│>…>│λ
n
│. Тогда генерируемые QR алгоритмом матрицы сходятся к
167 верхней треугольной. При этом лежащие ниже главной диагонали элементы a ij
(k)
матрицы A
k сходятся к нулю с линейной скоростью, определяемой отношением собственных значений матрицы A , а именно
│ │a ij
(k)
│
С│ λ
i
/ λ
j
│
k
. Заметим, что │ λ
i
/ λ
j
│<1 при i>j.
3. Условия теоремы исключают не только матрицы с кратными вещественными, но и с комплексными собственными числами. Однако QR алгоритм работоспособен при более общих предположениях, например, тогда, когда имеются группы собственных чисел, равных по модулю. В этом
случае формулировка сходимости усложняется, но основная идея остается прежней. Теперь процесс сходится к блочно треугольной матрице.
Паре комплексно сопряженных чисел соответствует диагональный блок
2
×2.
Если вещественный корень имеет алгебраическую кратность k, то теоретически должен появиться диагональный блок размером k
×k. Однако
ЭВМ воспринимает кратные λ как различные в пределах вычислительных ошибок. Поэтому на практике, как правило, имеют дело с диагональными блоками размером не более 2
×2.
4.Если задача плохо обусловлена QR алгоритм может оказаться неработоспособным.
5.Для уменьшения числа операций в процессе QR разложения матрица предварительно приводится к форме Хессенберга. Матрица
Хессенберга подобна исходной и имеет нули ниже главной поддиагонали.
Преобразования осуществляют по следующему алгоритму. A
k+1
=V
k
A
k
V
k
T
; k=0,1,…,n-3; A
0
=A; A
n-2
=H. V
k матрица отражения. Её строят так, чтобы обнулить элементы очередного k-ого столбца, расположенные ниже главной поддиагонали. V
k
=E-2w k
w k
T
; w k
=µ
k
(0,0,…0,a k+1,k
-s,a k+2,k
,…a n,k
)
T
Здесь первые k компонент нули. s=
2
,
1
ni ki ka
; sign(s)=sign(-a k+1,k
);
168
µ
k
=
1,
1 2 (
)
kkkks sa
. Это преобразование численно устойчиво и осуществляется за n-2 шагов. В дальнейшем генерируемые QR алгоритмом матрицы остаются матрицами Хессенберга.
Использование матриц Хессенберга снижает число операций на каждой итерации с n
3
до n
2 6. Для увеличения скорости сходимости можно использовать сдвиги.
Пусть η
n хорошее приближение к λ
n
. Тогда работают с матрицей H1
k
=H
k
-
η
n
E; λ(H1
k
)= λ(H
k
)- η
n
. Скорость сходимости определяется соотношением
│
nnin
│<│
ni
│, (элемент
( )
|
|
kknnniina
, |λ
n
-η
n
|
≪ λ
i
-η
n
|).
По окончании итерации следует осуществить обратный сдвиг – возврат к исходным собственным числам.
Итерационный процесс может формировать диагональный блок размером 2
×2. Такому блоку соответствует пара комплексно сопряженных собственных чисел. В этом случае следует осуществить двойной сдвиг. Обратный двойной сдвиг должен возвратить вещественные числа в элементах матрицы. Это, однако, может не случиться вследствие вычислительных ошибок. Проблема решается применением специального алгоритма, исключающего использование комплексной арифметики
(алгоритм Фрэнсиса).
7.Критерий останова - малость элементов под главной диагональю:
│a i+1,i
│<ε
m
(│a i,i
│+│a i+1,i+1
│); i=1,2,…n-1.
Скачать 3.86 Mb.