1
ГЛАВА 1.Численные методы решения систем
алгебраических линейных уравнений.
Классификация методов:
1.Прямые (точные);
2.Итерационные (приближенные).
§ 1. Вычислительные схемы метода Гаусса.
Метод Гаусса является примером вычислительного метода .
Пусть задана система линейных алгебраических уравнений: А*x=f (1) ,
А=
1..n; ji, ),(
ij
а
f=
)(
i
f
.
Будем рассматривать систему (1) в виде:
nnnnnnn
nn
fxaxaxaxa
fxaxaxaxa
...
.
.
.
...
322211
11312212111
(2)
Выберем из системы какое-либо уравнение, а в нем какою-нибудь неизвестную ,
коэффициент которой не равен 0. Не уменьшая общности , можно считать , что это первое
уравнение и неизвестное
1
x
, т.е. предполагаем
0
11
a
. Разделим первое уравнение на
11
a
,
которое будем называть ведущим:
11
1
1
11
1
113122121
,
...
a
f
g
a
a
b
gxbxbxbx
ij
j
nn
(3)
Умножим (3) последовательно на
13121
,...,,
n
aaa
и вычтем последовательно из (2) , тем
самым мы исключим неизвестную
1
x
из второго и т.д. уравнений. Преобразованные
уравнения будут иметь вид:
11
1
1i1ij
1
ij
1
1
3
1
22
1
2
1
2
1
23
1
232
1
22
* ,*а а а где
...
.
.
.
...
gafb
fxaxaxa
fxaxaxa
iij
n
nnnnn
nn
(4)
Систему можно рассматривать , как систему с n-1 неизвестными
n
xxx ,...,,
32
и с ней
поступим аналогичным образом как с системой (2). Продолжая этот процесс ,
предположим ,что он возможен до n –го шага . Тогда на n-том шаге получим
2
min
m
ii
mj
m
im
m
ij
m
ij
m
nn
m
mnm
m
nm
m
mn
m
nmm
m
mn
mnmnnmmm
gaff
baaa
fxaxa
fxaxa
gxbxbx
1
1
11
11111
11
,
,...
.
.
...
(5)
Предположим , что n-ый шаг –это последний шаг преобразования.
1.
Если m=n , то соединив все первые уравнения до n-го шага включительно , получим
систему :
nn
nnnnn
nn
gx
gxbx
gxbxbx
.
....
...
111
112121
(6)
Из последнего уравнения найдем х:
nn
nnnnn
nn
xbxbgx
gxbx
gx
121211
111
...
...
(7)
Процесс нахождения
n
x
(n=1..j) называется обратным ходом метода Гаусса. Процесс
приведения (2)
(7) – прямым ходом метода Гаусса.
2.
Пусть m<n:
n – последний возможный дополнительный шаг , т.е. на следующем шаге мы не можем
найти ведущего элемента отличного от 0
m
n
m
m
mnmnnmmm
f
f
gxbxbx
0
...
0
...
1
11
Если
1...nmj ,0
j
f
, то на n-том шаге исходная система может быть записана в виде:
mnmnm
nn
nn
gxbx
gxbx
gxbxbx
...
...
222
112121
(10)
Это означает, что неизвестные
m
xxx ,, ,
21
можно было выразить через
nmm
xxx ,, ,
21
,
а тогда (1) имеет множество решений.
Если m<n и хотя бы один из коэффициентов
j
f
1
, j=m+1...n ,то (1) не имеет решений.
На практике часто осуществляется контроль вычислений. Для этого на ряду с системой
А*x= f (1) параллельно решают А*у=d (11) ,
1
1,
j
ijjj
miafd
, т.е. в каждой строке
3
суммируют свободный член и коэффициенты матрицы А данной строки. Если при
решении (1) и (11) окажется ,что
1
jj
xy
,то система решена правильно.
ЗАМЕЧАНИЕ 1:
в методе Гаусса жестко закреплены все арифметические действия и внесение ошибки
приводит к непредсказуемому результату. Число всех умножений и делений равно :
)16(
3
2
nn
n
S
n
.
ЗАМЕЧАНИЕ 2:
схема метода Гаусса с выбором максимального элемента по строке может оказаться , что
на первом шаге коэффициент
0 но ,0
11
а
, в этом случае находят
njaа
ijj
1,max
0
1
. Объявляют
0
1 j
а
- ведущей и исключают неизвестные
0
j
x
из
второго , ... ,n-го уравнения. Аналогично поступают и с другими неизвестными.
3.
Схема с выбором максимального элемента по столбцу :
1
0
i
а
- ведущий ,но
niaа
ii
1,max
11
0
и исключают неизвестную
1
x
из 1-го , ...,
1,1
00
ii
,..., n-го уравнений.
4.
Схема с выбором максимального элемента по всей матрице:
00
ji
a
- ведущий, но
njiaа
ijji
1,,max
00
.В этом случае неизвестную
0
j
x
исключают из
1-го ,...,
1,1
00
ii
,...,n уравнений системы.
5.
Схема Жордана или схема оптимального исключения:
1.
x
x
x
xx
xx
xx
xx
xx
x1
0
0
1
2.
xx
xx
xx
x
x
xx
xx
10
01
00
10
1
..................
n.
1x0
1x
0x1
000
10
1
x
x
xx
xxx
Таким образом, в методе Гаусса исходная матрица приводится к треугольному виду, а
матрица Жордана – к диагональному. Число операций приблизительно в 2 раза меньше
числа операций в матрице Жордана.
6.
Некоторые другие точные методы решения системы линейных алгебраических уравнений.
а) метод ортогонализации
4
В этом методе исключают решения , определяемые как ортогональные к
подпространству, натянутому на векторы строки матрицы исходной системы и имеет
последнюю координату 1.
б) метод окаймления
Основан на том, что система более низкого порядка уже имеет решение и с помощью
формул решение системы более высокого порядка сводится к решению системы более
низкого порядка.
в) метод отражения
Основан на разложении матрицы А на произведение матрицы унитарной на правую
треугольную. Причем , здесь унитарная матрица образуется как произведение нескольких
квадратных матриц, так называемых матриц отражения. В результате решения исходной
системы сводится к последовательному решению систем с простыми матрицами.
г) метод сопряженных градиентов
Содержит в себе как точную, так и итерационную процедуры. Решение – вектор
последовательных приближений к решению системы, обращает в ноль один из
ортогональных векторов (невязок системы).
§ 2 Метод квадратного корня.
Дана система Аx=f (1) с матрицей А, удовлетворяющей условиям:
1. detA=0;
2. матрица А – эрмитова:
)(
* Т
АААА
Если матрица А не является эрмитовой, то без предварительного приведения матрицы А
к виду
fААхА
**
, то метод квадратного корня применять нельзя.
Однако, такое преобразование влечет за собой дополнительное большое количество
операций, тогда метод квадратного корня применять не стоит.
Схема метода квадратного корня построен на идеи построения матрицы А в виде
произведения треугольной и диагональных матриц, а именно, нужно найти такую
треугольную матрицу S
nn
n
n
S
S
SSS
S
000
00
2
12111
и диагональной матрицы D
nn
d
d
D
0
0
11
, причем
1
ij
d
, чтобы имело место равенство
DSSA
*
(2)
*
S
- нижняя треугольная матрица, D – диагональная, S- верхняя треугольная.
Формулы, позволяют записать матрицу А в виде (2) найдутся в книге.
Предположим, что мы получили представление А в виде (2), тогда решение системы (1)
осуществляется по следующей схеме:
)(,
*
ij
вВDSВ
Sx=y
B- известная матрица, У- неизвестный элемент, тогда
fВyDSхSAх
*
(3)
Для определения вектора У имеем систему (3) с нижней треугольной матрицей и решение
этой системы равносильно обратному ходу метода Гаусса. Определить из (3) вектор У,
5
решая систему : Sx=y (4) c верхней треугольной матрицей. Таким образом, основное в
методе квадратного корня – это по известным формулам вычисления коэффициентов
матриц S и D, а далее решение двух систем с треугольной матрицей. По объему операций
метод квадратного корня приблизительно в два раза эффективнее метода Гаусса. Также
положительным моментом является, если матрица А является разряженной, то S- также
будет разряженной. Недостатками метода квадратного корня является то, что он
применяется только для определения типа матрицы может исказить результат. Т.е. метод
в отличие от итерационного не является самоисправляющимся.
§ 3 Итерационные методы . Некоторые сведения
о нормах векторов матриц. Метод простой итерации.
Нормой вектора х называется сопоставляемое этому вектору действительное число :
1.
00
0,0
xх
;
2 .
xx **
3.
yxyx
Вводить норму можно различными способами:
1. кубическая
),,,(,
21 n
n
xxxxRx
2.
n
i
i
xx
1
11
3. сферическая или евклидова норма вектора
2
1
2
111
)(
i
xx
На ряду с покоординатной сходимостью вектора
*
xx
k
мы будем пользоваться
сходимостью , основанной на следующей теореме :
Для того , чтобы последовательность векторов
k
x
сходилась к
*
x
необходимо и
достаточно ,чтобы
kxx
k
,0
*
.
Опр: нормой матрицы А назовем сопоставляемое этой матрице А число
А
,
удовлетворяющее условиям:
1.
00,0,0 AА
;
2.
A*A*
;
3.
BABA
;
4.
BABA **
.
Будем говорить, что
A
согласованна с данной нормой вектора , если для любой
квадратной матрицы А , размерность которой равна порядку матрицы , выполняется:
xAxA **
.
Среди всех норм матриц , согласованных с данной нормой , выберем наименьшую. Для
этой цели за норму матрицы А выберем максимальную из норм векторов Ах в
предположении , что вектор х пробегает множество всех векторов нормы которых равны
1,т.е.
1,max xAxA
.
6
Для каждой матрицы А , в силу непрерывности норм матриц , этот максимум достигается ,
т. е. всегда найдется вектор
0
x
, такой что
1,
**
xAxA
. Введенную таким образом
норму будем называть подчиненной к данной норме вектора.
Примеры подчиненных норм:
1.
njaA
x
n
i
ij
1,max
:
1
1
1
2.
niaAx
n
j
ij
1,max :
1
1111
3.
значение есобственно оемаксимальн тт.е,0...
0T )( :
321
*
1
111111
n
AATAx
Если А- симметрическая положительно определенная, то
njAAA
j
T
1),(maxА тт,0
111
.
ЛЕММА: модуль каждого собственного значения матрицы не превосходит любой из ее
норм.
Док-во: обозначим
)(A
i
- собственные значения, отвечающие собственному вектору
матрицы А, тогда для любых
i
справедливо равенство
iii
xAx
Aвекторйсобственныx
xAxxAAxкТ
xxAx
i
iiii
iiiii
i
i
,
тт, ..
Лемма доказана.
Рассмотрим с.л.а.у. Аx=f (1),
при условии, что определитель не равен нолю и матрица А – положительно
определена(2). Это значит, что квадратичная форма матрицы больше ноля. Запишем (1) в
каноническом виде :
)(xx
(3)
Для этого (1) можно записать :
fAx
xx
(4)
Для нахождения решения (1) строим итерационный процесс по правилу:
любоеxxHx
kk
01
,)(
(5)
Метод (5) – метод итерации для решения системы (1);
k
x
- итерационная
последовательность;
)(
H
- матрица перехода от
k
x
до
1k
x
;
- стационарный параметр.
В нашем случае метод (5) – одношаговый стационарный метод, т.к.
)(
H
- матрица,
которая не зависит от номера итерации к, если бы
)(
H
=
)(
к
H
,
)(
H
=
)(
кк
H
, то в этом
случае метод называется нестационарным, а матрица перехода – ?нестационарной?.
Достоинства : простота, самоисправляемость.
число
оепроизвольнAEH
fxAEx
,)(
)(
7
Обозначим
*
хх
кк
(6),
к
- погрешность итерационного метода (5). Нам нужно
установить эту погрешность при
к
к
,0
.
kk
лк
H
fxHx
fхНх
)(
)(
)(
1
**
1
(7)
0
0
2
12
01
))((
))(()(
)(
k
k
H
HH
H
(8)
Т.к.
0
х
- произвольный начальный вектор, то можно считать, что
0
- произвольный
вектор, тогда можно положить
0
=
i
z
,
i
z
-собственный вектор матрицы А:
iii
zА
, тогда из (8)
i
k
ii
k
zzH )1())((
к
(9)
Пусть
0k
k
x
- последовательность, построенная по итерационному методу (5),
-
стационарный параметр,
0
х
- произвольное начальное приближение. Тогда для
сходимости последовательности
(1) решение - ,,
**
хкхх
к
, достаточным является
условие:
1)( АЕ
)(
При этом скорость сходимости линейная (метод сходится со скоростью геометрической
прогрессии). Если
Т
АА
- симметричная и сходимость имеет место для любых
0
х
, то
условие
)(
будет также необходимым.
Докажем это утверждение:
необходимость) из (8) следует, что
кН
к
к
к
,0)(
00
достаточность) пусть А – симметрическая и
любая ,0
0
к
.Нужно показать, что
)(
- выполняется, предположим, что
)(
- не выполняется, т.е.
1
, т.к. сходимость
имеет с любым начальным приближением, то в качестве вектора
0
возьмем вектор
p
z
,
тогда
pip
Нz
11max)( ,
0
.
В равенстве (9) вместо i положим р:
p
к
iк
z
1
, т.к. в силу предположения, что
)(
- не выполняется, то при
к
,
01
к
i
(не стремится), что противоречит тому, что метод сходится для любого
начального приближения, следовательно
)(
- выполняется.
Замечание(об оптимизации по
метода итерации 5 ):
мы установили, что сходимость итерационного процесса зависит от величины
, из (11)
следует, что скорость сходимости существенно зависит от максимального по модулю
собственного значения матрицы Н(
). Если
n
...
1
- собственные значения матрицы А, то
ii
H
1max))((max
.
Если собственные значения матрицы А могут быть как положительными, так
отрицательными, то из рис. 1 следует, что
11max
i
и поэтому итерационный
процесс является расходящимся.
8
Обратимся к часто встречающейся ситуации, когда
0
i
,
i
неизвестно:
Mm
i
0
.Скорость сходимости итерационного процесса можно охарактеризовать
величиной
1max)(
.
Рассмотрим задачу минимизации
)(
за счет выбора
: для нахождения минимума
)(
рассмотрим рисунок:
Из этого видно, что минимум
)(
:
mM
mM
Mm
Mm
оптопт
,
2
)1(1
00
Из полученной матрицы следует, что если матрица А имеет большой разброс собственных
значений , т.е.
опт
близко к 1, в силу этого естественная скорость сходимости будет
очень медленной и для получения решения даже с небольшой точностью потребуется
большое число итераций.
ОБЩИЙ НЕЯВНЫЙ МЕТОД ПРОСТОЙ ИТЕРАЦИИ ДЛЯ РЕШЕНИЯ
С.Л.А.У.
В приложении возникает необходимость решения систем высокого порядка, матрицы
которых обладают рядом специфических свойств:
симметричность,
положительно определена,
ленточность и т.д.
Такие системы в частности возникают при использовании метода сеток при решении задач
мат. физики.Наиболее интересные результаты решения таких задач получил Самарский.
Пусть задана система
(1) Ах=f,
njiаА
ij
,1,,
- матрица симметрическая и положительно
определенная, в том смысле, что
),(),(,0, xxxAxААА
Т
Для системы (1) запишем итерационный метод:
параметрfAx
xx
В
kk
, n1,ji, ),(b B,
ij0
1
(2)
Обычно матрица В с определителем не равным нолю выбирается таким образом, чтобы
легко вычислялась матрица
1
В
.В методе (2)к=0,1,2,.....
Перепишем метод (2) в виде:
,...2,1,0,)(
11
1
kfBxABEx
kk
(3)
(3)- общий неявный метод простой итерации для с.л.а.у. В (3)
АВЕН
1
)(
играет роль
матрицы перехода от итерации к к+1. Если
)(
Н
не зависит от
, то процесс,
определяемый (3) , называется стационарным итерационным методом. Если
)(
Н
=
)(
к
Н
,
)(
Н
=
)(
к
Н
,то итерационный процесс называется нестационарным. Мы будем
рассматривать стационарные методы.
Перейдем к изучению свойств последовательности
0k
k
x
и установим при каких условиях
она сходится к
fАх
1*
.
Лемма: пусть выполняются условия:
9
1.
Т
СС
;
2.
0, ВВВ
Т
;
3. ВС=СВ – перестановочны, С – n x n – матрица, тогда для положительно определенной
матрицы С, С >0 необходимо и достаточно выполнение ВС >0.
Док-во: необходимость:
пусть выполнено С>0, покажем, что условие 1. выполняется . Т.к. матрица
0, ВВВ
Т
, то
ее можно представить в виде :
0,
2
1
2
1
2
1
ВВВС
и она симметричная, тогда
2
1
2
1
2
1
2
1
)(
ВСВВВВ
, тогда матрицы ВС и
2
1
2
1
СВВ
- подобны. Т.к. преобразование подобия
матрицы не изменяет ее характеристического многочлена матрицы, то будет верно
соотношение:
)()(
2
1
2
1
СВВВС
,
)(ВС
- собственные значения матрицы ВС;
)(
2
1
2
1
СВВ
- собственные значения матрицы
2
1
2
1
СВВ
.
Покажем, что
2
1
2
1
СВВ
>0:
0))(,(),(
2
1
2
1
2
1
2
1
хВхСВххСВВ
Т
,т.к. С>0 , а это значит, что (Сy,y)>0,
0y
.
0)(
2
1
2
1
2
1
2
1
Т
СВВСВВ
, т.к. матрица положительно определена и симметричная , то все ее
собственные значения положительны
)(
2
1
2
1
СВВ
>0,
)(ВС
>0, поскольку
ВССВВСВС
ТТТ
)(
и все ее собственные значения положительны, то и сама
матрица положительна.
Достаточность:
Пусть выполнены условия 1,2,3 и
. Покажем, что С>0 :
собственные значения
)()( )(
2
1
2
1
СВВВСиВС
, следовательно квадратичная форма
0),(),())(,(),(
2
1
2
1
2
1
2
1
2
1
2
1
yСyхВхСВхВхСВххСВВ
Т
, то из последнего условия будет,
что С>0. Что и требовалось доказать.
Условия, при которых итерационная последовательность
кх
к
,
, построенная по
методу 3 при любом начальном приближении
0k
k
x
→
fАх
1*
сводятся к
теореме Самарского: Если выполнены условия
1.
0, ААА
Т
;
2. .
0В
, то для сходимости последовательности построенной по (3) для любого
начального приближения
,..2,1,0,
0
nRх
n
достаточно выполнения неравенств:
2В>
А (
);
А>0 (
);
3.
Т
ВВ
;
4.
ВААВ
, то (
) и (
) ( условия Самарского) являются также необходимыми для
сходимости последовательности
*
хх
к
(к решению системы (1) при любых
0
х
).
Докажем достаточность:
*
xx
kk
- погрешность приближенного решения
k
x
,
10
fAx
xx
В
fAx
xx
В
k
kk
*
**
1
Вычитая почленно из первого уравнения второе для погрешности
k
x
(6) **
(5) ттогда , (4) 0
1
1
к
1
1
1
кк
кк
кk
kk
ВА
ВАAВ
Если бы выполнялось
кк
кк
,0что , (6) из тт,,0
к1
. Или же
*
хх
к
. Для изучения поведения разности
кк
1
: умножим справа (4) на
(7) 0),(2),2( получим )(2
1111
кк
к
кккккк
АВ
Для получения оценок симметризуем в равенстве (7)
),(),(
),(),(),(),()),((
)),(( ),,)2(( , 0
(8) 0)),((),)2((
),*
2
(2),
2
(2),2(
*
22
11
111111
11
11
кк
11
11
111111
11
кккк
кккккккккккк
ккккк
кккк
к
кккк
кккк
кккккккккккк
кккк
к
АА
ААААА
ААВ
ААВ
ААВ
Т.к.
0АА
Т
равенство (8), которое называют основным энергетическим тождеством
Самарского, примет вид :
(9) 0)),(1)((),)2((
1
11
кккк
кккк
АААВ
Т.к.
0 АА
Т
и В>0 и выполнены условия
и
, то в равенстве (9) первое слагаемое
0
к
0),(),(
11
кккк
АА
(10) . Это неравенство в теории разностных схем
называют основным энергетическим неравенством. Из него следует, что
0
),(
к
кк
А
-
невозрастающая, кроме того, эта последовательность ограничена снизу, т.к.
0),(
кк
А
,
а такая последовательность является сходящейся и
0)),)(,((
11
lim
кккк
k
АА
,
поэтому из (9) следует, что
0
lim
к
k
. Обозначим G=
АВ
2
и
кк
k
z
1
. Тогда , т.к.
G>0, то существует
0
- константа такая, что
(11)
),(
8
),(
2
2
1
1
22
kkk
kk
k
kkkkkkk
z
zGzzzzGz
Учитывая, что
k
k
,0
из (11) следует, что
к
кк
,0
1
и следовательно в
силу (6):
к
к
,0
.
11
Докажем необходимость:
Пусть выполнены условия (1)-(4) теоремы и
*
хх
к
. Нужно полагать, что выполнены
условия
и
.
Докажем от противного: предположим, что
не выполняется , т.е. имеет место
кккккк
АВВАВЕАВЕ
АВ
)2()2()(
02
111
1
(12)
Обозначим через
z,
соответственно собственные значения и отвечающий им
собственный вектор матрицы
)2(
1
АВВ
т.е.
zzАВВ
))2((
1
.
В силу выполнения условий (1)-(4) матрица
)2(
1
АВВ
- симметрическая, кроме того
0)2(,0
1
АВВ
(по предположению), то согласно лемме 1. среди собственных
значений существует хотя бы одно меньше либо равно ноля, обозначим
0
s
,
отвечающий ему собственный вектор
0
х ,0
s
z
-любое, то в качестве
ss
zxxzxx
*
00
*
0
,
, тогда для итерационной последовательности:
11
0,,01
)1(,
....1,0,)2(
1
10
1
1
k
s
sj
k
sк
s
k
sкs
кк
zkz
zz
кАВВЕ
Это противоречит условию сходимости
*
хх
к
. Необходимое условие
:
...1,0,)(
,0
1
1
кАВЕ
А
кк
(13)
Обозначим через
,z
собственный вектор и значения матрицы
АВ
1
, согласно лемме 1.
существует хотя бы одно
ppp
zz ,,0
0
- собственный вектор, отвечающий
p
. Для
погрешности
k
- получим:
0,01
1
kz
p
k
pк
, а это противоречит тому, что
*
хх
к
при любом начальном приближении.
Теорема Самарского устанавливает тот факт сходимости последовательности
0к
к
х
к
fАх
1*
, но ничего не говорит о скорости такой сходимости. Для установления
скорости такой сходимости сформулируем лемму:
Пусть
0
Т
АА
,
0
Т
ВВ
и выполняются условия
,
, сформулированные в теореме
Самарского, тогда
),(),(
11 кккк
ВВ
(14).
Док-во:
Запишем условия
,
Самарского 2В>
А (
),
А>0 (
).
Погрешность
к
у нас вычисляется последовательно по рекуррентной формуле:
0
1
k
kk
AВ
. Запишем это равенство в виде
0))()((
2
1
11
1
kkкk
kk
AВ
и умножим его слева на величину:, тогда
12
получим
(15) 0)),(()),((
2
)),)(
2
((
0)),((
2
)),((
2
)),((
2
)),((
2
)),(())((
111111
111111
111111
кккккккккккк
кккккккккккк
кккккккккккк
ВAAB
AAA
AВВ
В (15) первое слагаемое
0)),)(
2
((
11
кккк
AB
в силу условия (
), а второе
0)),((
2
11
кккк
A
в силу условия (
).
0)),((
11
кккк
В
. Распишем подробнее:
),(),(0),(),(
(16) ),(),(
),( тт,0 тт.к,0),(),(),(),(
1111
11
11111
кккккккк
кккк
к
Т
к
Т
кккккккк
ВВВВ
ВВ
ВВВВВВВВ
Конец док-во.
Если ввести следующее обозначение:
),,(
2
кк
В
к
В
то условие (14) можно записать
0,
1
Т
В
к
В
к
ВВ
(17)
Теорема 2 (об оценке
В
к
):
Пусть выполнены условия :
1.
0
Т
АА
,
0
Т
ВВ
,
2. выполняются условия (
) и (
) теоремы Самарского,
3. А, В,
и некоторое число
согласованы в том смысле, что выполнено неравенство
)(
11
ВАВ
,тогда для погрешности
*
хх
кк
справедливы оценки
10 ,
0
В
к
В
к
(18).
Док-во:
Положим
00
,10,
k
k
k
и подставим данное выражение в равенство:
(19) 0А
1
-АА ,
0)
1
(
:виде в перепишем равенство последнее
0
00
0
1
1
1
11
1
1
k
kk
k
kk
k
kkkk
k
kk
k
k
k
k
k
k
k
kk
В
ВВВ
ВAB
AB
ABAB
AB
к оценке
В
к
можно применить лемму 2, которая приводит к соотношению
В
к
В
к
1
. Из последнего неравенства имеем
13
В
к
к
В
к
В
к
к
В
к
к
к
к
ВВ
к
ВВ
к
В
к
00
0
01
1
но,
Теорема доказана.
Из теоремы 2 следует, что
0
к
к
В
к
и
также, как
к
. Т.е. тем быстрее, чем
1
. Представляет интерес для заданной матрицы А определить свойства матрицы В и
выбрать
так, чтобы значение
было как можно меньше.
Пусть между А и В не существует соответствия
ВАВ
21
(20), В>0,
0
Т
АА
,
0
12
.Обычно
21
,
- постоянные, которые не заданы, а заданы интервалы, в которых
они могут находятся:
222111
0 ,0 MmMm
. Сопоставим
)(
и (20),
получим уравнение для определения параметров
и
.
1212
12
2
1
2
,
1
1
Обозначим
12
/
. Для
можем получить оценки
0
2
1
2
1
m
M
M
m
и тогда
1
1
и это значение будет минимальным, если
2
12
опт
12
12
2
1
Mm
Mm
Mm
m
M
опт
Таким образом можно сказать, что итерационный процесс общего неявного метода
простой итерации имеет наиболее быструю сходимость
12
опт
2
Mm
и этот процесс
можно записать в виде :
f
Mm
хА
Mm
ВВх
кк
1212
1
2
)
2
(
(21).
Последовательность
0к
к
х
, построенная по (21) , стремится к системе
fАх
1*
со
скоростью геометрической прогрессии, знаменатель которой
12
12
Mm
Mm
g
опт
.
Рассмотрим частный случай общего неявного метода простой итерации:
матрицу В нужно выбирать таким образом, чтобы она конкретно отображала свойства и
была связана с матрицей А. Рассмотрим некоторые возможные способы выбора В:
1. метод простой итерации
iinn
aaaadiagB ),,,,(
2211
элемент главной диагонали матрицы А.
1
, тогда формула (3) примет вид :
,..1,0,)(
11
1
kfBxABEx
kk
(22)
Метод (22) будет сходится, если выполняется одно из условий :
14
1)1(max
1)1(max
1
11
1
1
1
1
ij
n
i
ii
ij
ij
n
j
ii
ij
a
a
ABEg
a
a
ABEg
ji
ji
ij
,0
,1
- символ Кронекера.
Таким образом, если в матрице А имеется диагональное преобладание по строкам,
либо столбцам, то итерационный процесс (22) стремится к системе Ax=f с любого
начального приближения. Отметим, что для сходимости (22) достаточно выполнение
условия 2В>A,
0
Т
АА
, это следует из теоремы Самарского.
2. Метод Зейделя
Пусть
0
Т
АА
, тогда она может быть представлена в виде:
T
LDLA
nn
11
11
21
a0
0a
D,
0
00
000
nmn
aa
a
L
В=D+L,
1
, тогда общий неявный метод простой итерации в данном случае примет
вид:
,..1,0,)()(
1
kfxALDxLD
kk
- метод Зейделя для системы Аx=f. Метод
Зейделя сходится для любой с.а.л.у., матрица которой удовлетворяет условию
0
Т
АА
. Метод Зейделя является частным случаем метода релаксации.
3.Метод релаксации.
,,LDВ
параметр. Учитывая это, общий неявный метод простой
итерации в данном случае примет вид:
fAx
xx
LD
k
kk
1
)(
(24) – метод релаксации для системы Ax=f , при этом
если
10
, то метод нижней релаксации,
1
- полной релаксации,
21
-
метод верхней релаксации.
Запишем (24) в виде
))(()()(
,)()(
1
1
1
ALDEHH
fLDxHx
kk
(25)
В (25) имеется возможность выбора параметра
, этот выбор следует производить
таким образом, чтобы
)(
Н
была по возможности < и чтобы выполнилось
1)(
Н
.
Для сходимости метода релаксации (24) достаточно выполнение условий:
1.
0
Т
АА
,
2.
20
.
Покажем, что если выполнятся эти условия, то будут выполнены (
) и (
)
Самарского. Имеем
)0,
АА
-выполняется. Покажем, что выполнено
условие (
): нужно показать, что
ALDВ
)(22
.
Имеем,
),(),)2((),)((),)(2( xAxxDxxxLDLxxLD
T
.
Таким образом, условие (
) Самарского выполняются.
Как уже говорилось, при
1
метод полной релаксации совпадает с методом
Зейделя.
15
§ 4 Метод скорейшего спуска.
Этот метод предназначен для решения с.л.а.у. Ax=f (1) с вещественной,
симметрической, положительно определенной матрицей А.
В методе скорейшего спуска, а также в методе сопряженных градиентов отыскание
решения системы (1) связано с задачей минимизации следующего функционала
F(x)=(Ax,x)-2(f,x) (2), который является квадратичной функцией отпеременных
.,,,
21 n
ххх
Это объясняется тем, что решение системы (1)
fAx
1*
достигает
минимум функционала (2) на множество векторов из вещественного векторного
пространства.
Док-во:
0)),((),(2),(),(2),(
)),(2,(),(),()()(
*******
****
xxxxAxAxxAxxAxxAx
xfxAxxfAxAxxFxF
(3)
при этом равенство в (3) возможно, когда
*
xx
=0,
*
xx
.
Таким образом, задача нахождения решения системы (1) сводится к задаче отыскания
вектора х, доставляющего минимум функционалов F(x).
Лекция 5.
Метод имеет следующую вычислительную структуру. Исходя из
0
x
к решению
*
x
системы
fAx
вычисляется вектор
00
Axfr
и число
00
00
0
,
,
Arr
rr
.
Следующее приближение
1
x
определяется по формуле:
0001
rxx
Вектор
2
x
вычисляем из условия минимума функции
11
rxF
, где
00011
ArrAxfr
.
Исходя из условия минимума функции
11
rxF
мы найдем
1
по формуле
11
11
1
,
,
Arr
rr
и
1112
rxx
. Далее процесс построения последовательных приближений
осуществляется по формулам:
111
kkkkk
ArrAxfr
(4)
kk
kk
k
Arr
rr
,
,
(5)
kkkk
rxx
1
(6)
Отметим, что
k
r
особенно при большом порядке матрицы
A
удобно вычислять по
формуле:
111
kkkk
Arrr
. Чтобы из-за ошибок округления таким образом вычисляемые
k
r
через
несколько шагов не стали сильно отличаться от истинных невязок
k
Axf
, их надо время
от времени вычислять по формуле
kk
Axfr
. В методе скорейшего спуска
ортогонализация векторов невязок системы
k
r
не производится. Формулы (4)-(6) в другой
трактовке с использованием формулы (2) для функционала
)(xF
можно записать в
следующем виде:
)(
1
xgradFxx
kkk
, где
k
определяется из условия минимума функционала по
:
))(())((
min
kkkk
xgradFxFxgradFxF
.
16
Исследуем свойства последовательности векторов
,...,
10
xx
построенных по (4)-(6).
Лемма 1.
Если
nia
i
,1,
некоторые положительные числа, а
i
числа, удовлетворяющие
неравенству
niMm
i
,1,0
, то справедливо:
2
2
1
11
4
1
)(
1
M
m
m
M
a
aa
n
i
i
n
i
i
i
n
i
ii
(7)
Док-во:
Вводим обозначения
mM
a
a
i
i
n
i
i
i
i
,
1
. Тогда (7) запишем в виде
2
11
4
1
M
m
m
M
n
i
i
i
i
n
i
i
(8)
Отметим, что
1
1
n
i
i
и
ni
m
M
M
m
i
,1,
. Так как среднее арифметическое
меньше или равно среднему геометрическому, то
2
11
1
4
1
i
ii
n
i
i
i
i
n
i
i
(9)
Функция
1
на отрезке
m
M
M
m
,
принимает максимальное значение на
концах, то есть при
M
m
или
m
M
и это наибольшее значение равно
M
m
m
M
.
Таким образом
ni
M
m
m
M
i
i
,1,
1
(10)
Из (9) в силу (10) получим
2
2
1
2
2
111
4
1
4
1
4
1
M
m
m
M
M
m
m
M
M
m
m
M
n
i
i
n
i
i
n
i
i
i
i
n
i
i
Введем в рассмотрение понятие: функция ошибок, определив ее формулу
,AxG
(11)
xx
*
вектор ошибок.
Справедлива
Лемма 2
Последовательность значений функции ошибок
),...(),...,(),(
10 k
xGxGxG
, где
k
x
определяется по формулам (4)-(6) и (11) стремится к нулю, при
k
.
Док-во:
В силу (4)-(6) и (11) имеет место равенство:
17
kk
kk
kk
Arr
rr
xGxG
,
,
2
1
и
kkk
rrAxG ,)(
1
. Значит
k
kk
kk
k
kk
k
q
Arr
rr
xG
xGxG
xG
1
,
,
1
2
1
(12)
kkkk
kk
k
rArArr
rr
q
1
2
,,
,
. Оценим снизу величину
k
q
. Пусть
n
,...,
1
- собственные значения
матрицы
A
, удовлетворяющие условию
n
...
21
и
n
uuu ,...,,
21
собственные вектора матрицы
A
, ортогональные
друг к другу и нормированные так, что
1,
ii
uu
,
ni ,1
. По условию
0A
все
собственные значения
0
i
. Предположим, что
1
m
и
M
n
. Разложим вектор
k
r
по
собственным векторам матрицы
A
.
nnk
ucucucr ...
2211
(13)
Так как под
k
r
мы понимаем ненулевой вектор невязок системы, то в разложении (13) не
все
0
i
c
. Из (13) имеем
nnnk
ucucucAr
...
222111
nnnk
ucucucrA
1
2
1
221
1
11
1
...
следовательно
n
i
ikk
crr
1
2
,
;
n
i
iikk
cArr
1
2
,
;
n
i
iikk
crAr
1
211
,
;
Отметим, что вектора
i
u
ортогональные друг к другу.
jiuu
ji
,0,
.
Тогда для
k
q
получим:
n
i
n
i
iiii
n
i
i
k
cc
c
q
1 1
21
2
1
2
)(
. Используя лемму 1 для
k
q
получим
оценку снизу:
2
M
m
m
M
u
q
k
k
k
k
q
xG
xG
1
1
и
1
2
1
mM
mM
xG
xG
k
k
0
)1(2
1
42
1
xG
mM
mM
xG
mM
mM
xG
mM
mM
xG
k
kkk
(14)
Так как
1
mM
mM
, то
kxG ,0)(
.
Теорема
Последовательные приближения
k
xxx ,...,,
10
, построенные по методу скорейшего спуска
(4)-(6) сходятся к решению системы
fAx
со скоростью геометрической прогрессии.
Док-во:
Из леммы 2
kxxxxAxG
kk
,0)),(()(
**
. Так как
0A
, то это означает, что
*
xx
k
,
k
. Определим скорость сходимости.
2
),()(
kkk
mAxG
(15)
kk
kk
Arr
rr
,
,
2
18
Из (15) с учетом (14) получаем
k
k
k
mM
mM
m
xG
m
xG
)()(
0
. Последнее
неравенство означает, что
k
k
,0
со скоростью геометрической прогрессии со
знаменателем
mM
mM
q
.
Замечание
Величина знаменателя
q
зависит от наибольшей и наименьшей границ спектра матрицы
A
. Если
A
плохо обусловлена, то есть
M
и
m
сильно отличаются друг от друга, то
величина
q
будет близка к единице, в следствие чего метод будет медленно сходящимся.
Отметим два свойства приближений метода скорейшего спуска:
1. Невязки двух последовательных приближений ортогональны.
0,
1
kk
rr
2. Каждое последующее приближение не дальше к точному решению, чем
предыдущее.
**
1
xxxx
kk
Метод скорейшего спуска может быть применен к системам с несимметрической
матрицей. После умножения системы
fAx
на
T
A
слева. В отличие от метода
квадратного корня, где такое преобразование неэффективно, в вычислительной системе
метода скорейшего спуска матрица
AA
T
может не вычисляться. Действительно
рассмотрим систему
fAAxA
TT
. В качестве вектора невязок мы должны взять вектор:
kk
Ar
, где
kk
Axfr
. Тогда расчетные формулы примут вид:
kkkkkkk
Arxxx
1
, где
),(
,
),(
,
kk
kk
kk
T
kk
k
AA
rArA
rr
. Метод скорейшего
спуска легко реализуется на ЭВМ и сложности могут возникать только от большого
порядка матрицы
A
.
§5. О погрешности приближенного решения систем линейных алгебраических
уравнений, об обусловленности систем и матриц.
Как уже отмечалось в некоторых методах численного решения систем линейных
алгебраических уравнений о точности полученного приближения решения чаще всего
судят по векторам невязок системы. Однако если для одного класса матриц малость
вектора невязок системы в некоторой метрике означает и малость компонент вектора
погрешностей, то для другого класса такой связи может и не быть. Убедимся в этом.
Рассмотрим
fAx
(1)
*
x
- точное решение,
y
- некоторое приближенное решение. Рассмотрим вектора
yx
*
и
Ayfr
(2)
Их будем называть соответственно вектором погрешности и вектором невязок. Пусть
матрица
A
системы (1) имеет хотя бы одно очень малое по модулю собственное значение
. И пусть
z
- собственный вектор, который соответствует
. Тогда
zfAzAxzxA
**
. Компоненты вектора
zx
*
могут отличаться весьма сильно
от компонент вектора
*
x
хотя бы в силу малости
. Компоненты вектора
zf
будут
мало отличаться от компонент вектора
f
. В связи с этим необходимо ввести такие
соотношения между векторами
и
r
, которые позволяли бы по величине
r
судить более
точно о величине вектора
. В практике вычислений большое значение имеют не нормы
19
векторов
,
r
, а отношения
f
r
x
,
*
, которые являются в некотором смысле
относительными погрешностями. Для количественной характеристики этих соотношений,
а также векторов
и
r
, введено понятие обусловленности систем и матриц. Введем в
рассмотрение величину:
),(sup
*
f
r
x
r
(3)
Если
мало, то из (3)
r
f
x
*
и малость нормы вектора невязок
r
означает
малость нормы погрешности
. В этом случае говорят, что (1) хорошо обусловлена. Если
велико, то малость нормы вектора невязок
r
не означает малости нормы
. В этом
случае говорят, что (1) плохо обусловлена. Число
называется мерой обусловленности
системы (1). По аналогии введем понятие обусловленности матрицы. Из определения
вектора невязок и
, а также из определения нормы матрицы имеем
1
1*
supsupsup
A
r
rA
r
yx
r
rrr
(4)
Учитывая (3) и (4) получаем:
1
*
A
x
f
(5)
Будем рассматривать систему (1) при всевозможных значениях правых частей
f
. Тогда
решением этой совокупности будет некоторое множество
X
векторов
*
x
, отвечающих
соответствующим значениям
f
при одной и той же матрице
A
. Рассмотрим поведение
величины
, определяемое по (5) при всевозможных
Xx
*
, а именно рассмотрим
величину
Xx
*
sup
и обозначим
1
*
1*
*
sup
AA
x
AAx
v
Xx
(6)
v
- число обусловленности матрицы
A
. Из (6) следует, что если
A
близка к особенной,
то
v
будет для такой матрицы велико. В этом случае говорят, что матрица
A
плохо
обусловлена. Если
v
мало, то соответствующую матрицу
A
называют хорошо
обусловленной. Как правило система с хорошо (плохо) обусловленной матрицей будет
хорошо (плохо) обусловлена. Значения
v
зависят от того, каким образом мы определим
норму матрицы
A
, однако неравенство
1v
справедливо при любом выборе норм
матриц. Покажем это:
Имеем
1
1
maxmaxmaxmax
1
i
iii
AAv
. Здесь
i
- собственное
значение
A
.
i
- собственное значение
1
A
. Получим оценку для погрешности
. Такая
погрешность в сильной степени зависит от того, как изменяется решение (1) при малых
изменениях ее коэффициентов и свободных членов, а это означает, что оценка
зависит
от меры и числа обусловленности матрицы системы. Рассмотрим наряду с (1) систему
gBy
(7)
где
B
-заданная матрица,
g
- заданный вектор. Предположим, что
B
и
g
связаны
соответственно с
A
и
f
следующим соотношением
20
ACAB
,
1 qC
;
fg
,
(8)
Отметим, что во многих вычислительных алгоритмах решение системы (1) удовлетворяет
системам вида (7), для которых матрица
C
и вектор
могут быть вычислены реально
или по крайней мере для них известны оценки. Получим оценку погрешности
yx
*
.
Из (7) и (8) получим
fAycE )(
или
............
2222
ccEfccfccEfccE
.
Учитывая, что
Ayfr
из последней формулы получим, что
)......(
22
ccEfccr
Этот вектор можно рассматривать как невязки приближенного решения системы (1). В
силу определения чисел
и
v
имеем
fqq
q
f
ccEfcc
f
r
x
yx
x
1
1
1
......
22
*
*
*
fqq
q
v
1
1
1
(9)
Из (9)
относительная погрешность
*
x
тем меньше, чем меньше
v
. Малость такой
погрешности также сильно зависит от того, как далеко уклоняется матрица
C
и вектор
g
соответственно от матрицы
A
и вектора
f
, то есть зависит от малости величин
q
и
. Из
(9) можно получить оценку для нормы погрешности
. Обозначим
fqq
q
ql
1
1
1
,
. Учитывая, что
yxyyxyx ***
из (9) получим
,1
,
qvl
qvl
y
(10)
В предположении, что
0,1
qvl
. Таким образом мы установили зависимость вектора
погрешности от решения
y
и от величин
и
v
, характеризующих удаленность матрицы
C
и вектора
g
от матрицы
A
и вектора
f
. Большое число методов решения (1) на
преобразовании
A
к некоторому каноническому виду, например, диагональному,
треугольному, ленточному и так далее. Чаще всего такие преобразования выполняются
путем умножения матрицы
A
слева на некоторую невырожденную матрицу М. Выясним,
какой класс матриц М при указанном преобразовании не меняет
v
. То есть нужно
определить класс матриц М, для которых имеет место равенство
AvMAv
(11)
Для любой невырожденной матрицы М справедливы неравенства
AvMvMAv )(
(12)
Так как
AMMA
и
1111
MAMA
. С другой стороны
)()(
1
)(
AvMv
MAv
(13)
Так как
)()()()()(
1
MAvMvMAvMvAv
. Аналогично получим, что
)(
)(
1
)( Mv
Av
MAv
(14)
Из (12)-(14)
21
)()()()(
)(
1
);(
)(
1
max AvMvMAvAv
Mv
Mv
Av
(15)
Из (15) следует, что если при некоторой заданной М правая граница неравенства
достигается для какой-либо матрицы
A
и при этом
)(Mv
велико, то
)(MAv
может стать
достаточно большим. Если положить
1)( Mv
, то из (15) сразу будет следовать (11),
значит все невырожденные матрицы,
v
которых равно единице, не меняют числа
обусловленности матрицы
A
, то есть если
)()(1)( MAvAvMv
. Отметим, что таким
свойством в случае третьей нормы обладают унитарные и ортогональные матрицы.
§6. Проблема собственных значений. Краткий обзор методов численного решения
проблемы собственных значений.
Пусть задана
)(
ij
aA
,
nji ,1,
. Собственным значением матрицы
A
называется число
, такое что система линейных алгебраических уравнений:
xAx
(1)
имеет ненулевое число решений
0x
. Вектор
x
, отвечающий числу
, - собственный
вектор. Рассмотрим
EA
det
. Система (1) имеет нулевое решение
0det EA
(2)
(2) представляет собой уравнение n-ной степени относительно
, со старшим
коэффициентом
n
1
, который можно записать в виде:
0...1
1
1
n
nn
n
pp
(3)
n
nn
pp
...
1
1
- характеристический многочлен матрицы
A
. Само уравнение (3) –
характеристическое, собственное, вековое уравнение матрицы
A
. Корни уравнения (3)
n
,...,
1
- собственные значения матрицы
A
. При решении различных задач возникают
разные требования о собственных значениях и собственных векторах матрицы и это
порождает многообразие проблем и методов решения этой задачи:
1. Для решения ряда задач механики, физики, химии требуется вычисление всех
i
, а иногда и всех собственных векторов. Это задача – полная проблема
собственных значений.
2. В ряде случаев требуется найти лишь минимальное или максимальное по
модулю собственное значение. Возникает в физике. Здесь приходится решать
задачи, эквивалентные задаче отыскания собственных значений матриц
размерности порядка
63
1010
. В таких задачах при малых размерностях
матриц используются итерационные методы, при больших – вероятностные.
3. При исследовании колебательных процессов иногда требуются определить два
максимальные по модулю собственные значения матрицы. Причем меньшее из
них обычно достаточно определить с меньшей точностью.
4. В тех же задачах иногда требуется отыскание собственного значения матрицы,
близкого к заданному значению
0
или же отыскание расстояния от заданного
0
до спектра матрицы.
Формально рассуждая можно было бы сказать, что задачи 2,3,4 являются частным
случаем общей проблемы собственных значений и достаточно ограничиться набором
методов для решения полной проблемы собственных значений. Однако такой подход
приведет к неоправданно большому объему вычислений. Решение задач 2,3,4 обычно
сводят к отысканию максимального по модулю собственного значения некоторой
матрицы В, полученной из
A
, такой что это собственное значение соответствует
отыскиваемому значению матрицы
A
. Рассмотрим случай, когда все
i
R
. Если
требуется определить максимальное или минимальное значение матрицы
A
, то следует
взять матрицу
CEAAg )(
при достаточно больших положительных (отрицательных)
22
значениях
C
максимальному собственному значению
CEA
соответствует
максимальное (минимальное) значение матрицы
CEA
. Для задачи 4. При некотором
0
рассматривают матрицу
2
0
EAcE
. Тогда максимальному по модулю собственному
значению построенной матрицы соответствует искомое значение матрицы
A
. В ряде
случаев в качестве
)(Ag
выбирают:
1
0
)(
EAAg
. При этом
)(Ag
не выписывают в явном виде, а необходимые по ходу
вычисления вектора определяют как решения системы уравнений
yxEA
0
.
Нахождение максимального по модулю собственного значения и отвечающего ему
собственного вектора матрицы
A
.
Предположим, что
A
обладает
i
e
- полная система собственных векторов.
iii
eAe
,
ni ,1
.
ijji
ee
,
.
n
...
321
. Выберем произвольный вектор
0
x
и будем строить последовательность векторов
mm
Axx
1
,
,...2,1,0m
Представим
0
x
в виде:
i
n
i
i
eCx
1
0
. Тогда
i
m
i
n
i
ii
m
n
i
im
eCeACx
11
. Так как у нас
1
равняется максимальному по модулю значению, то можно записать
)(
2111
m
m
m
oeCx
.
Рассмотрим скалярное произведение
m
m
m
m
m
m
mm
oCoeCoeCxx
1
2
2
1
2
121112111
1))(),((,
m
m
mm
oCxx
1
2
2
1
2
111
1,
. Положим:
m
m
m
m
m
mm
mm
m
o
oC
oC
xx
xx
1
2
21
1
2
2
1
2
1
1
2
1
2
1
2
11
1
1
1
1
,
,
. Кроме того имеем
m
m
m
oCx
1
2
3
11
1
.
m
m
i
m
i
m
m
m
oC
eC
x
x
e
1
2
311
1
1
1
,exp
1
1
2
41
1
2
3
11
m
m
m
m
i
m
i
i
oie
o
C
eC
где
m
m
C
11
arg
. Таким образом для
достаточно больших
m
мы можем получить с заданной точностью максимальное по
модулю
m
1
и соответствующий ему вектор
m
e
1
.
Метод Крылова.
Для иллюстрации основной идеи метода Крылов вводит в рассмотрение каноническую
систему обыкновенных дифференциальных уравнений первого порядка с постоянными
коэффициентами.
nn
yayayay
12121111
...'
23
nn
yayayay
22221212
...'
………………………………..
nnnnnn
yayayay ...'
2211
Aa
ij
Характеристическое уравнение этой системы имеет вид:
0...1
1
1
n
nn
n
pp
То есть совпадает с характеристическим уравнением
матрицы
A
. Корни этого уравнения являются собственными значениями матрицы
A
.
Если заданную систему обыкновенных дифференциальных уравнений первого порядка
удается свести к одному дифференциальному уравнению порядка
n
с постоянными
коэффициентами:
ypypypypy
nn
nnn
'...
1
2
2
1
1
По виду этого уравнения можно выписать его характеристическое уравнение:
0...
1
1
n
nn
pp
Крылов указал также на возможность алгебраической интерпретации этой идеи.
Рассмотрим любой вектор
0
c
, согласованный по размерности с размерностью матрицы
A
.
По этому вектору будем строить последовательность векторов
;
01
Acc
;
0
2
12
cAAcc
…
0
; cAc
m
m
до тех пор пока не встретим первый вектор, например,
m
c
, который будет являться
линейной комбинацией предыдущих линейно независимых векторов
110
,...,,
m
ccc
, то есть
пока не будет выполняться равенство
02211
... cqcqcqc
mmmm
. Если
nm
, то
i
q
будут коэффициентами многочлена
степени
m
, являющегося делителем характеристического многочлена степени
n
, то есть
решив уравнение
0...
1
1
m
mm
pp
мы найдем
m
собственных значений матрицы
A
.
Метод Данилевского.
Этот метод построен на известном из линейной алгебры факте о том, что преобразование
подобия
ASS
1
не изменяет характеристического значения матрицы
A
.
nnnnnn
nn
aaaa
aaaa
A
121
1111211
...
.....
...
EASEASSSASSEASS
''
11
Поэтому удачно подобрав преобразование подобия можно надеяться получить матрицу,
собственный многочлен которой легко выписывается по виду этой матрицы. Данилевский
предложил преобразованием подобия приводить исходную матрицу к виду Фробениуса:
01...0
....
0...01
...
21 n
ppp
Тогда характеристический многочлен матрицы
совпадает
с характеристическим многочленом матрицы
A
.
0...
1
1
n
nn
pp
. Преобразование подобия
S
осуществляем по шагам:
121
...MMMS
nn
, где каждая последующая матрица начиная с
1n
M
приводит к
виду Фробениуса только одну строку. Начиная с последней нужно приводить к виду
Фробениуса последнюю строку матрицы
A
, то есть в последней строке должны быть все
24
нули и 1 на главной диагонали. Для этой цели
A
умножаем на
1
1
n
M
и
1n
M
, где
10...00
1
...
.....
.....
00...01
111
2
1
1
1
nn
nn
nnnn
n
nn
n
n
a
a
aa
a
a
a
M
Вид Фробениуса это последняя строка
01...00
.
10...00
...
.....
.....
00...01
121
1
1
nnnnnn
n
aaaa
M
таким образом преобразование
1
1
1
1
AAMM
nn
.
01...00
...
.....
.....
...
1
)1(
11
)1(
22
)1(
11
)1(
1
1
11
1
12
1
11
1
1
nnnnnn
nn
aaaa
aaaa
A
Далее матрица
2
1
1
2
nn
MAM
и до тех пор пока полностью не будет приведена к виду
Фробениуса. После нахождения характеристического многочлена и вычисления его
корней метод Данилевского позволяет уменьшить объем вычислений при определении
собственного вектора. Пусть мы нашли собственное значение
. Тогда собственный
вектор матрицы Фробениуса имеет вид
1
...
2
1
n
n
y
, а собственный вектор матрицы
A
вычисляется по правилу
yMMMSyx
nn 121
...
. Здесь рекомендуется не перемножать
матрицы, а последовательно вычислять произведение матрицы на вектор, то есть
начинаем вычислять справа.
Треугольный степенной метод.
Разрабатывался Бауэром и предназначается для вычисления итерационным путем всех
собственных значений матрицы
A
. При этом предполагается, что для собственных
значений матрицы
A
справедливо распределение
n
...
21
. Пусть
njicc
ij
,1.),(
0
0
- некоторая матрица задаваемая вычислителем. В основе
вычислительной схемы метода лежит последовательное вычисление матриц
1n
M
1n
M
25
1......
.....
0...1
0......1
0...001
)(
2
)(
1
)(
32
)(
31
)(
21
k
n
k
n
kk
k
k
CC
CC
C
C
и вычисление матриц
)(
)(
2
)(
22
)(
1
)(
12
)(
11
......00
.....
.....
......0
......
k
nn
k
n
k
k
n
kk
k
r
rr
rrr
R
Последовательное вычисление матриц по правилу:
1.
110
RCAC
, где
A
заданная матрица,
0
C
- некоторая невырожденная матрица, заданная
вычислителем,
11
, RC
- искомые матрицы. Характер вычислений при определении
коэффициентов матриц
11
, RC
аналогичен обратному ходу метода Гаусса, то есть не
вызывает затруднений.
2.
221
RCAC
Из этого соотношения определяем
22
, RC
…
к.
kkk
RCAC
1
Оказывается, что
i
k
ii
k
r
)(
lim
, то есть диагональные элементы матрицы
R
стремятся к
собственным значениям матрицы
A
. Поэтому при достаточно больших
k
можно
предположить
nir
k
iii
,1,
)(
. Тем самым можем вычислить все собственные значения
матрицы
A
.
Метод вращений.
Пусть
A
эрмитова матрица
*)( AA
. Основой для построения метода служит известная
теорема из алгебры, утверждающая, что если
A
эрмитова матрица, то существует
унитарная матрица
V
(
*
1
VV
), то преобразование подобия с этой матрицей приводит
A
к диагональному виду
AVV
1
(4)
-диагональная матрица, на диагонали которой стоят собственные значения матрицы
A
.
Так как для унитарной матрицы справедливо
*
1
VV
, то последнее равенство можно
записать в виде
AVV *
(5)
Эта формула не может быть использована для прямого вычисления элементов матриц
V
и
, так как она представляет систему
2
n
уравнений с
nn
2
неизвестными. Здесь
2
n
элементов матрицы
V
и
n
элементов матрицы
. Однако можно трактовать задачу
приведения матрицы
A
к диагональному виду, как приближенную задачу в следующем
списке. Предположим, что некоторым преобразованием вида (5) матрица
A
приведется к
некоторой
~
21
2
~
221
112
~
1
~
...
....
...
...
nnn
n
n
26
Предположим также, что недиагональные элементы матрицы
~
таковы, что элементами
2
2
11
n
k
k
можно пренебречь по сравнению с
~
1
2
2,1
22
n
kk
k
~
2
… …
2
1
1
n
k
nkn
~
n
Записав полученное преобразование по правилу
~~
* VAV
(6)
мы можем утверждать, что данным преобразованием мы приводим исходную матрицу к
виду, где на диагонали стоят элементы, совпадающие с заданной точностью с
собственными значениями матрицы
A
. Пусть матрица
T
AA
,
Ra
ij
. Для такой
матрицы метод вращений заключается в построении последовательности матриц
;...;;
210
AAAA
каждая последующая получается из предыдущей при помощи
элементарного шага, состоящего из преобразования подобия путем умножения
предыдущей матрицы на некоторую матрицу, называемую матрицей вращений, имеющей
вид
10
cossin
1
1
sincos
.
001
ij
V
Итерационный процесс продолжается до тех пор, пока вне диагонали элементы матриц
k
A
станут таковыми, что сумма их квадратов в строке будет меньше модуля диагонального
элемента той же строки. Кроме перечисленных выше методов существует большое
количество методов для нахождения собственных значений векторов матриц
специфического вида.
§7. Решение систем нелинейных уравнений. Постановка задачи. Метод итерации.
Задача решения уравнений в достаточно общем виде может быть сформулирована так:
Пусть заданы
X
и
Y
, элементы, которых будем обозначать
yx,
. Природа этих
элементов может быть любой. Это могут быть числа, функции, линии, поверхности и так
далее. Говорят, что в некотором
X
задан оператор
A
со значениями во множестве
Y
,
если для каждого элемента из
X
существует некоторый
YAxy
. Этот оператор можно
рассматривать как преобразование
X
в
Y
.
x
- оригинал,
y
- изображение. Допустим, что
во множестве
Y
взят некоторый
0
y
. Нужно найти такие
Xx
, которые являются для
0
y
оригиналами. Записать такую задачу можно в виде уравнения
0
yxA
(1)
Особое значение для нас будут иметь уравнения, в которых неизвестными величинами
будут числа. Такие уравнения являются частным случаем операторных уравнений, когда
X
и
Y
- численные пространства конечных размерностей. В этом случае
0xf
,
nn
RRf :
(2)
27
0),...,(
...
0),...,(
0),...,(
1
12
11
nn
n
n
xxf
xxf
xxf
(3)
системы вида (3) имеют очень важное значение в науке и ее приложениях. Много задач
сводятся к системам вида (3). Одним из наиболее распространенных методов решения
системы (3) является метод итерации или метод повторных подстановок, который
является общим методом решения нелинейных уравнений и применим к весьма широкому
классу операторных уравнений. Общность метода, его простота и удобная и удобная
реализация сделали метод итерации одним из основных методов вычислительной
математики. Приведем (2) к виду
xx
(4)
Множество значений, которое может принимать
x
обозначим
X
. В приложениях
X
чаще
всего конечный или бесконечный отрезок числовой прямой. Множество значений,
которое может принимать
x
, когда
Xx
, обозначим
Y
. Функцию
x
можно
рассматривать как оператор, преобразующий
X
в
Y
. Уравнение (4) означает, что нужно
найти такие
Xx
, которые при преобразовании множества
X
оператором
переходят
сами в себя, иначе говоря, во множестве
X
нужно найти точки, остающиеся
неподвижными при преобразовании множества
X
оператором
. В методе итерации
последовательные приближения вычисляются по правилу:
kk
xx
1
,...2,1,0, k
(5)
0
x
- начальное приближение.
В случае, когда
X
есть конечный или бесконечный отрезок числовой прямой, задача
решения (1) имеет простой геометрический смысл:
В системе координат проводим
xy
и график
xy
. Выбираем
0
x
.
Рис.1 Рис.2
Если
1)(' x
, то итерационный процесс
расходящийся, как на рис.3
Рис.3
28
Возвращаемся к уравнению
0)( xf
(2)
Преобразовать (2) в (4) наиболее распространенный прием. Записать (2) в виде
)()( xfxAxx
(6)
nijxaxA
ij
,1)),(()(
выбирается так, чтобы она не обращалась в ноль для
Xx
. Вторая
задача, стоящая перед вычислителем: удачный выбор
0
x
.
0
x
надо выбирать достаточно
близко от
*x
. Третья задача:
x
должна быть построена таким образом, чтобы
выполнялись условия Липшица:
10; qxyqyx
Четвертая задача:
0
x
и
q
, и область, в которой ищется решение должны быть согласованы
между собой.
Теорема о сходимости метода итерации.
Пусть
0
x
- начальное приближение к корню
*x
уравнения
0)( xf
.
0
}{
kk
x
итерационная
последовательность, построенная по правилу:
;
1 kk
xx
)()( xfxAxx
Пусть выполняются условия:
1.
)(x
i
ni ,1
определены и непрерывны по всем
j
x
,
nj ,1
в области
}{),(
00
rxxRxrxSS
n
.
2.
x
являются отображением сжатия в
S
, то есть для
Szx ,
выполняется
zxqzx
,
10 q
3. Функция
x
, числа
rq,
и
0
x
согласованы так, что выполняется неравенство:
r
q
xx
1
00
.
Тогда
)
Все
rxSx
k
,
0
.
)
*lim xx
k
k
,
;** xx
),(*
0
rxSx
.
)
Справедлива оценка скорости сходимости
k
k
q
q
xx
xx
1
*
00
.
Док-во:
Докажем
)
методом математической индукции.
rxSx ,
00
очевидно. Предположим, что
rxSxx
k
,,...,
01
. Покажем, что
rxSx
k
,
01
.
011101
... xxxxxxxx
kkkkk
01111
...)( xxqxxqxxxx
k
kkkkkk
(7)
00001
mxxxx
. Применяя оценку (7) ко всем звеньям в неравенстве
rm
q
mmqmq
kk
000
1
0
1
1
...
rxSx
k
,
01
.
)
Для доказательства сходимости
0
}{
kk
x
применим признак Больцано-Коши.
000
1
11
1
...... m
q
q
mqmqxxxxxx
k
kpk
kkpkpkkpk
(8)
29
Правая часть (8) не зависит от
p
и так как
10 q
, то
0
начиная с некоторого
значения
k
правая часть станет меньше
при любом
p
. Поэтому существует
*lim}{lim
0
xxx
k
k
kk
k
. И так как все
rxSx
k
,
0
, то
*x
также будет принадлежать
rxS ,
0
. Чтобы показать, что
*x
является корнем уравнения
xx
надо в левой и
правой части равенства
kk
xx
1
перейти к пределу. Так как
*
1
xx
k
и
*xx
k
, то
** xx
.
)
Мы получили (8)
0
1
m
q
q
xx
k
kpk
(8)
Если устремить
p
, то так как
*xx
pk
, в пределе будем иметь
0
1
* m
q
q
xx
k
k
.
Теорема
Во всяком множестве точек, где для
x
,
yx,
справедливо неравенство:
yxyx
.
Уравнение
xx
имеет не более одного решения.
Док-во:
Допустим противоположное. Будем считать, что это уравнение имеет два различных
решения
yx,
yx
.
yxyxyx
в силу условия теоремы. Получили
противоречие.
Таким образом в сфере
rxS ,
0
при выполнении условий теоремы 1 корень
*x
единственный.
Замечание 1 В методе итерации речь ведется только об одном корне, для другого корня
нужно искать новую область, строить новую
x
, получать другие ограничения.
Замечание 2 Итерационный процесс сходится со скоростью геометрической прогрессии
и, чтобы скорость была удовлетворительной,
x
должна быть выбрана так, чтобы
10 q
.
Замечание 3 В случае первой нормы для того, чтобы
rxS ,
0
- отображение
x
было
отображением сжатия должно выполняться условие
1maxmax
1
1
n
j
j
i
ni
Sx
x
x
.
§8. Об ускорении сходимости простой одношаговой итерации.
Проблема ускорения сходимости является одной из общих задач вычислительной
математики и может возникать всюду, где для разыскиваемого элемента строится
последовательность вычисляемых элементов в той или иной мере приближающихся к
нему. Рассмотрим задачу. Метод простой одношаговой итерации для одного уравнения с
численной неизвестной величиной. Для сходящихся
...2,1,0; nS
n
SS
n
, для которых
погрешность приближения
n
=
SS
n
изменяется по закону, близкому к геометрической
прогрессии
,
1 nn
q
10 q
. Простейшее правило улучшения сходимости было
предложено Эйткеном. Указанный закон изменения
n
весьма прост и по нему легко
найти представление любой
n
S
, для которой он выполняется. Так как в приближенном
равенстве
nn
q
1
ошибка должна быть малой величиной порядка
nn
q
, где
0
n
при
n
. Поэтому точное равенство может быть следующим:
rxS ,
0
30
nnn
q
1
1
.
Если применить указанное равенство несколько раз, начиная с
n
:
)1(...11...111
0210212
2
11
nn
n
nnnnnn
qqq
. Тогда для
n
S
можно получить равенство:
021
1...11
nn
n
n
AqSS
(1)
0
A
,
10;0;0 qn
n
.
За каноническую форму таких равенств естественно принять последовательность, в
которой все
0
n
n
n
AqSS
(2)
Эта последовательность зависит от трех величин
qAS ,,
. Они могут быть найдены по трем
любым значениям
n
S
. Для нас представляет интерес предельное значение
S
. Для этого
достаточно из (2) исключить
qA,
. Выберем три последовательных значения
11
,,
nnn
SSS
.
1
1
n
n
AqSS
n
n
AqSS
1
1
n
n
AqSS
Разделив второе на первое и третье на второе:
SS
SS
q
n
n
1
SS
SS
q
n
n
1
Приравнивая и выражая
S
:
11
2
11
2
nnn
nnn
SSS
SSS
S
(3)
Если (3) применить не к канонической, а к любой близкой к ней, то мы получим не
обязательное значение
S
, а некоторое значение тем более близкое к
S
, чем больше
n
.
Обозначим его
;
2
11
2
11
nnn
nnn
n
SSS
SSS
,...2,1n
(4)
(4) – преобразование Эйткена заданной
n
S
в новой
n
. Оно применимо ко всякой
n
S
, для
которой
02
11
nnn
SSS
для любого
n
или хотя бы для некоторых
n
, начиная с
N
.
Можно ожидать, что
n
будет сходиться к
S
быстрее, чем
n
S
. Возвратимся к методу
простой итерации для одного уравнения
;
1 nn
xx
,...2,1n
(5)
2
1
*'***
nnnn
oxxxx
, то есть можно положить, что
nnn
x
'
1
(6)
и следовательно для итерационной последовательности (5) закон изменения
n
при
условии, что
1'
n
xq
имеет как раз такую форму, для которой для улучшения
сходимости применимо правило Эйткена. В этом она называется правилом Стефенсона.
Лекция 9.
Полученная оценка не может быть лучше, она может быть заменена более наглядной, но
несколько грубой ошибкой:
12
1
2
2
1
n
hxx
n
n
, (7).
формула (7) будет доказана для случая одного нелинейного уравнения.
31
Теорема 2. Пусть для оператора
xf
и начального элемента
0
x
выполнены условия (1)-
(5) теоремы 1.
h
n
t
211
, тогда уравнение
0xf
имеет единственное
решение при
2
1
h
в области
txx
0
.
Так как условия теоремы 1 выполнены, то уравнение
0xf
имеет решение
x
в
области
ttxx
0
, которая составляет часть области предусмотренной теоремой 2.
Нужно показать, что всякое другое решение будет совпадать с
x
.
Рассмотрим случай, когда
h
<
2
1
. Допустим, что существует решение
x
в области
txx
0
. Можно положить, что
00
tttxx
, (8), где
0
<
<
1
. Введём
следующее обозначение:
xfГxxF
0
. Отметим свойства оператора
xF
, которыми
мы будем пользоваться:
xxF
0
0
xF
xfГxF
0
,
выполним оценку
2
2
2
0
2
0000
10
2
000
10
00000001
2
1
2
1
sup
2
1
sup
2
1
tBkxxBk
xxxxxfГxxxxxF
xxxFxFxFxFxFxfГxxxx
Рассмотрим выражение
0
22
111
2
1
2
1
tpB
t
BB
t
B
tkBtBk
Обозначим
ttpt
B
0
1
и
tp
B
t
B
tk
1
2
1
2
.
Заменим
B
на
0
1
tp
, получаем:
1
0
0
000
0
00
0
2
1
1
2
1
tt
tp
tp
tttptttptp
tp
tpttptp
tp
tBk
1
2
1
ttxx
, (8).
Сравнение полученной оценки с оценкой
00
ttxx
говорит о том, что она
получается из оценки (8) заменой
на
2
,
0
t
и
0
x
соответственно на
1
t
и
1
x
,
следовательно, n-кратное применение этой оценки приводит к неравенству:
n
n
n
ttxx
2
<
t
n2
, (9), где
<
1
и не зависит от n.
Поэтому,
0
n
n
xx
, то есть последовательность
xx
n
, но по теореме 1
последовательность
xx
n
, следовательно,
xx
.
32
При
2
1
h
число
может быть равно
1
, но в этом случае корень
t
совпадает с
t
.
Так как последовательность
tt
n
, то из (9) следует, что
0
n
xx
, а из этого следует,
что
x
совпадает с
x
.
§10 Метод Ньютона в случае одного нелинейного численного уравнения.
Метод Ньютона для систем уравнений.
Пусть дано
0xf
,
1
Rx
,
xf
- достаточно гладкая функция. Обозначим через
x
точное решение уравнения
0xf
. Предположим, что каким-либо образом указано для
x
начальное приближение
0
x
. Последовательные приближения в методе Ньютона
строятся по правилу:
n
n
nn
xf
xf
xx
1
, (1),
0
x
,
0n
,
1
,
2
, …
Условиями возможности построения итерационного процесса (1), являются следующие
требования:
1) все
n
x
принадлежат области определения функции
xf
;
2)
0
n
xf
,
0n
,
1
,
2
, …;
Метод Ньютона совпадает с простым одношаговым итерационным процессом для
уравнения
xx
, где
xf
xf
xx
. Метод Ньютона часто называют методом
касательной.
Выясним, как ведут себя последовательные приближения
n
x
вблизи точного решения
x
. Обозначим
nn
xx
. Из равенства
xx
вычтем почленно уравнение (1),
получим:
n
nnn
n
n
nn
xf
xfxf
xf
xf
1
.
Заметим, что
nnnnnnnn
xfxfxfxfxf
2
2
1
0
.
Заменим
nnn
xfxf
на
nnn
xf
2
2
1
, тогда получим:
22
2
1
2
1
2
1
nn
n
nnn
n
o
xf
xf
xf
xf
.
Погрешность в методе Ньютона на каждом шаге изменяется по квадратичному закону.
Теорема (о сходимости метода Ньютона в случае одного уравнения). Пусть выполнены
следующие условия:
1) функция
xf
определена и дважды непрерывно-дифференцируема на отрезке
0
xx
<
, (
), при этом
kxf
,
x
;
2)
0
0
xf
,
B
xf
0
1
;
3)
0
0
xf
xf
;
4)
2
1
kBh
;
33
5)
h
h211
;
тогда верна оценка
12
1
2
2
1
n
hxx
n
n
, (2).
B
t
B
kttP
1
2
1
2
,
1
2
1
2
ht
, где
t
и
B
tP
.
Для нахождения корня уравнения
0
будем строить последовательные приближения
n
nn 1
,
h
h
n
211
,
0
0
.
Заметим, что
0
111
nnnn
, тогда
2
1
2
1111
2
1
2
1
nnnnnnnnn
h
, а
1
nn
h
,
тогда
n
nn
n
n
nn
h
h
12
2
1
1
, (3).
Так как
0
0
,
1
1
, то при
1n
из (3) получаем:
2
1
12
12
h
h
, при
2
1
h
находим
2
:
2
3
1212
.
Если
2n
:
2
12
23
12
h
h
, то при
2
1
1
,
2
3
2
,
2
1
h
будет
4
1
23
, откуда
4
7
4
1
2
3
2323
и так далее. Продолжая этот процесс, получим:
n
n
1
22
, (4).
Исходя из определения
n
и учитывая, что
0
, получим:
111
11
1
1
1
nnn
nn
n
nn
, (5).
По формуле Тейлора
nnnnn
2
1111
2
1
2
1
2
1
2
1
nnn
h
, (*).
Подставим (*) в равенство (5), получим:
2
1
1
2
1
1
1
n
n
n
h
h
, (6).
На основании (4)
n
n
n
h
1
2
1
2
222
1
11
, (7).
Поэтому
2
1
2
2
n
n
n
h
, (8).
Применяя последовательно неравенство (8) для
1n
:
2
1
1
2
h
и учитывая
2
211
h
h
, получим
h2
1
.
Для
2n
:
3
2
12
2
2
1
hh
.
Продолжая этот процесс, получим:
12
1
2
2
1
n
h
n
n
.
34
Учитывая, что
12
1
2
2
1
n
httxx
n
nnn
n
.
Метод Ньютона – это метод линеаризации, позволяющий сводить решение
нелинейных задач к последовательному решению ряда линейных задач. Он является
одним из старейших вычислительных методов решения нелинейных уравнений.
Рассмотрим несколько модификаций этого метода. Все его изменения связаны либо с
увеличением скорости сходимости, либо с упрощением и уменьшением объёма
вычислений, так как метод Ньютона есть частный случай общего метода простой
итерации.
nn
xx
1
,
xf
xf
xx
, то для него справедливо правило увеличения скорости
сходимости Эйткена. Так как основной объём вычислений в методе Ньютона связан с
нахождением и обращением
n
xf
, то модификации связанные с уменьшением объёма
вычислений основаны на упрощении вычисления
n
xf
.
Метод секущих.
В методе секущих
n
xf
заменяется на
1
1
nn
nn
n
xx
xfxf
xf
, тогда
1
11
1
1
1
nn
nnnn
nn
nnn
nn
xfxf
xfxxfx
xfxf
xfxx
xx
, (9).
Равенство (9) имеет простой геометрический смысл:
nnn
xfxM ,
,
111
,
nnn
xfxM
,
через эти точки проведём прямую линию и уравнение прямой, проходящей через эти две
точки, примет следующий вид:
nn
n
nn
n
xx
xx
xfxf
xfy
11
. Правило (9) означает, что за
1n
x
принята абсцисса точки пересечения секущей с осью
ox
.
По сравнению с основным методом Ньютона, этот метод имеет некоторые
особенности: он является двух шаговым, на каждой итерации для организации
вычислений нам нужно значение функции в двух точках, во-вторых, для организации
вычислений нужно два начальных приближения
0
x
,
1
x
корня
x
уравнения
0xf
.
Итерационный процесс (9) может быть осуществлён, если, во-первых, все
n
x
принадлежат
области определения функции
xf
, во-вторых, на каждой итерации выполняется
неравенство:
0
1
nn
xfxf
,
1n
,
2
, …
Для погрешности
nn
xx
можно получить
11
2
nnn
xf
xf
.
В методе секущих погрешность изменяется по закону близкому к квадратичному.
Метод Ньютона с постоянной касательной.
В этой модификации метод Ньютона освобождается от вычисления производных
n
xf
на каждом шаге и пользуется только лишь производной
0
xf
. Итерационный
процесс строится по правилу:
0
1
xf
xf
xx
n
nn
, (10),
,1n
,2
…
35
Для погрешности
nn
xx
можем получить
0
1
1
xf
xf
nn
.
Этот закон близок к геометрической прогрессии со знаменателем
0
1
xf
xf
q
, и если
точка
0
x
близка к решению
x
, то знаменатель
q
близок к нулю, то есть скорость
сходимости достаточно высока. Однако сравнение с основным методом Ньютона говорит
о том, что этот метод сходится гораздо медленнее основного метода Ньютона.
Можно построить также циклический метод Ньютона, в котором производная остаётся
постоянной
s
шагов.
Демпфирированный метод Ньютона.
Рассмотрим этот метод на примере решения систем нелинейных уравнений
0xФ
,
Ф
:
nn
RR
. Итерационный процесс строится по правилу:
kkkkk
xФx
x
Ф
xx
1
,
10
k
,
k
выбираем из условия, чтобы функция
xФxf
i
n
i
i
2
1
2
,
0
i
,
1k
xf
<
k
xf
.
Преимущество этого метода заключается в том, что он сходится к решению с любого
начального приближения, выбранного на промежутке, где существует и при том
единственный корень системы
0xФ
.
Замечание 1. Основными условиями в теореме о сходимости метода Ньютона являются
неравенства
2
1
kBh
и
h
h211
. Выполнение этих неравенств можно
добиться за счёт удачного выбора начального приближения
0
x
. Выбор
0
x
может
осуществляться:
1) вероятностными методами;
2) градиентными;
3) методами оптимизации;
Замечание 2. Ещё раз о целях модификации:
1) упрощение вычислений;
2) уменьшение объёма вычислений;
3) ослабление условий на
0
x
и на функцию
xf
, при которых имеет место
сходимость;
Метод Ньютона для систем.
Рассмотрим систему нелинейных уравнений
0xf
, (11),
f
:
nn
RR
. Запишем её
подробнее:
0,...,
.........................
0,...,
1
11
nn
n
xxf
xxf
, (12)
Метод Ньютона для системы (12) является частным случаем для операторных уравнений
и обобщённым случаем для одного уравнения.
Пусть нам известно
0
x
- начальное приближение к решению системы (12). Итерационный
процесс строится следующим образом: каждое последовательное приближение к решению
36
системы (12) находится из системы линейных алгебраических уравнений составленных по
предыдущему приближению, а именно, решения находятся как решения системы
)()()1(
1
)(
k
i
k
j
k
j
n
j
j
k
i
xfxx
x
xf
, (13),
ni ,1
Формулировка теоремы о сходимости метода (13) зависит от нормы матрицы, в
которой будем рассматривать систему. Сразу рассмотрим случай кубической,
октаэдральной, сферической нормы.
Теорема. Пусть выполнены условия:
nii
xxfxf ,...,
1
,
ni ,1
определены и дважды непрерывно дифференцируемы в
области
)(
:
)0(
ii
xx
,
ni ,1
n
i
ii
xx
1
)0(
n
i
ii
xx
1
2
2
)0(
При этом для вторых производных в области
)(
справедливы неравенства:
a)
K
xx
xf
n
kj
kj
i
1,
2
,
ni ,1
;
b)
i
kj
i
L
xx
xf
2
,
KLLL
n
...
21
,
nkj ,1,
;
c)
2
1,,
2
2
K
xx
xf
n
kji
kj
i
;
Значения
)0()0(
1
,...,
n
xx
являются начальным приближением к решению системы (12)
и для справедливы неравенства:
a)
)0(
xf
i
,
ni ,1
;
b)
n
i
i
xf
1
)0(
;
c)
n
i
i
xf
1
2
2
)0(
;
Матрица Якоби
xf
имеет в точке
)0(
x
определитель
)0(
xfD
отличный от
нуля и если
jk
D
есть алгебраическое дополнение к элементу
n
j
x
xf
)0(
, то:
a)
BD
D
n
j
jk
1
1
,
nk ,1
;
b)
BD
D
n
k
jk
1
1
,
nj ,1
;
c)
BD
D
n
kj
jk
2
1
1,
2
1
;
Для чисел
,,BK
выполняется условие
2
1
2
KBh
;
37
Для числа
выполняется условие
B
h
h211
, тогда:
a) система
0xf
имеет в области
)(
решение
n
xxx ,...,
1
;
b) последовательность Ньютоновских
приближений (13) построена, принадлежит области
)(
и сходится к решению
x
;
c) скорость сходимости оценивается
неравенством:
)()( kk
ii
ttxx
,
ni ,1
;
n
i
kk
ii
ttxx
1
)()(
;
)(
2
1
1
2
)( k
n
i
k
ii
ttxx
;
где
B
h
h
t
211
- меньший корень уравнения
0
1
2
1
2
t
B
Kt
и
)(k
t
последовательность Ньютоновских приближений к корню
t
, построенная при начальном
приближении
0
)0(
t
. Эта теорема есть частный случай теоремы о сходимости метода
Ньютона для операторных уравнений.
§11 Приближённые методы оптимизации.
Градиентные методы.
Рассмотрим задачу минимизации функции
xf
,
n
Rx
, для этой задачи будем
использовать итерационный процесс вида:
kkkk
pxx
1
, (1),
0k
,
1
,
2
…
здесь вектор
k
p
- вектор, определяющий направление движения из точки
k
x
,
k
-
числовой множитель, величина которого определяет длину шага вдоль выбранного
направления движения.
Процесс (1) будет определён, если будут указаны способы построения вектора
k
p
и
вычисления величины
k
на каждом шаге итерации.
Для того чтобы приблизиться к точке
x
, где возможно выполняются с определённой
степенью точности необходимые условия экстремума функции
xf
, естественно
двигаться из точки
k
x
в направлении убывания функции.
Если точка
k
x
не является стационарной точкой или точкой минимума, то существует
бесконечно много векторов
p
определяющих направление спуска, причём каждый из них,
когда
xf
дифференцируемая функция, определяется условием:
0, pxf
k
. Это
вытекает из следующих рассуждений: предположим, что
pxx
k
, тогда
ppxfpxfxfxf
kckk
,
2
,
2
kkkc
xxxx
,
10
.
38
Если
0, pxf
k
, то, во всяком случае, для достаточно малых значений
будет
выполнено неравенство:
k
xfxf
. Проще всего выбирать положительное направление
k
p
, удовлетворяющее условию
0, pxf
k
.
kk
xfp
, где
k
p
- антиградиент
функции
xf
.
Определение. Итерационный процесс
kkkk
xfxx
1
, (2),
,0
k
,0k
…,
называется методом наискорейшего спуска или градиентным методом.
Приступим к изучению свойств алгоритма (2). Способ выбора скалярного множителя
k
:
выбираем
одно и то же на всех итерациях и определяем точку
kk
xfxx
;
вычисляем значение функции
xf
в этой точке;
производим проверку неравенства
kkk
pxfxfxf ,
(3).
10
;
если неравенство (3) выполняется, то значение
выбираем в качестве искомого
k
;
если неравенство (3) не выполняется, то производим дробление
, то есть
умножаем на некоторую константу
10
до тех пор, пока (3) не окажется
справедливым;
Обоснование такого способа выбора
приводит к теореме 1.
Теорема 1. Если функция
xf
ограничена снизу, её градиент удовлетворяет условию
Липшица
yxRyfxf
, (4),
n
Ryx ,
, а выбор
k
осуществляется описанным
выше способом, то для итерационного процесса (2) будет справедлива сходимость:
0
k
xf
при
k
какова бы ни была начальная точка
0
x
.
Теорема 2. Пусть
xf
дважды непрерывно дифференцируемая функция, её матрица
вторых производных удовлетворяет условию
22
, yMyyxfym
, (5),
0m
, а
последовательность
k
x
строится по правилу (2), где
k
выбирается согласно неравенству
(3). Тогда для любого начального приближения
0
x
будет справедливо
xx
k
,
xfxf
k
,
k
, где
x
- точка минимума функции
xf
и для скорости
сходимости справедливы оценки
xfxfqxfxf
k
k 0
и
2
k
k
cqxx
, где с –
некоторая постоянная,
c
,
10 q
.
Рассмотрим градиентный метод. Мы рассмотрели выбор параметра
k
в
итерационном процессе (2) связанный с проверкой неравенства (3), этот способ выбора
длины шага не является единственно возможным.
Рассмотрим другой способ выбора длины
k
. Каждый из таких способов определяет
различные варианты градиентного метода. В процессе доказательства теорем 1 и 2 (мы не
доказывали) установлено, что неравенство (3) будет заведомо выполнено, если
R
1
,
M
12
. Поэтому, если известны постоянные
R
и
M
, характеризующие
минимизируемую функцию
xf
, то в методе (2) можно выбрать
k
, где
R
1
0
либо
M
12
0
. Таким вариантом градиентного метода удаётся
уточнить величину знаменателя
q
в оценках скорости сходимости.
39
Теорема 3. Если функция
xf
удовлетворяет условиям теоремы 2 и для итерационного
процесса (2)
k
, где
M
2
0
, то для скорости сходимости итерационного
процесса (2) справедлива оценка
xxqxx
k
k 0
, где
Mmq
1,1max
,
причём
mM
mM
q
min
достигается при
mm
2
.
Рассмотрим ещё один способ выбора длины шага
. Будем выбирать
k
из условия
минимума функции вдоль направления движения, то есть выбранное значение
k
должно
удовлетворять условию
kkkkk
xfxfxfxf
min
. При таком варианте
градиентного метода все теоремы остаются справедливыми и итерационный процесс
сходится в точке минимума, начиная с любого начального приближения со скоростью
сходимости геометрической прогрессии со знаменателем
mM
mM
q
. Такой вариант
градиентного метода называется методом наискорейшего спуска.
Рассмотрим другие методы.
Пусть
xF
- произвольная симметрическая матрица, удовлетворяющая условию
22
, yRyyxFy
,
0
,
,x
n
Ry
. Будем выбирать вектор
p
(вектор,
определяющий направление движения):
xfxFp
.
0,,
xfpxfxFxfxpxf
,
0
xf
, то есть
x
не является
стационарной точкой. Исходя из этого для
xf
можно построить минимизирующий эту
функцию итерационный процесс
01
xfxFxx
kkkk
,
0
k
,
0k
,
1
, …
В дальнейшем мы будем рассматривать итерационный процесс вида:
kkkkk
xfxFx
1
1
,
0
k
,
0k
,
1
, …
kk
xFxF
1
. Для такой матрицы будет справедливо:
2
1
1
2
1
, yMyyxFym
, где
0
2
1
R
p
m
,
p
M
1
1
,
xfxFxp
1
- вектор,
задающий направление спуска функции
xf
.
Отметим, что для рассмотренных итерационных процессов остаются справедливыми
свойства, сформулированные в теореме 2.
Мы рассмотрели несколько вариантов градиентных методов. Они отличаются друг от
друга способами выбора величины
k
и построением вектора
xp
. Все они имеют
приблизительную скорость сходимости и используются в зависимости от того, что
известно (какие известны постоянные характеризующие заданную функцию
xf
).
В частности
m
и
M
можно выбирать как наименьшее и наибольшее собственное
значение матрицы вторых производных функции
xf
. И если эта матрица хорошо
обусловлена, то есть величины
m
и
M
мало отличаются друг от друга, то постоянная
q
будет близка к нулю, и итерационный процесс будет быстро сходящимся. Если матрица
плохо обусловлена, то
1q
, и итерационный процесс будет медленно сходящимся.
Метод Ньютона.
В градиентных методах для определения направления движения используется лишь
линейный член разложения функции
xf
в ряд Тейлора, то есть используется лишь
наиболее грубая аппроксимация минимизирующих функций. Пусть
xf
(минимум,
который требуется определить) является строго выпуклой
40
yfxfyxf
11
,
10
. Рассмотрим функцию
yxyxyfyxyfyfx
,
2
1
,
. Эта функция представляет собой
квадратичную аппроксимацию функции
xf
. В силу строгой выпуклости функции
xf
,
x
тоже строго выпуклая. Поэтому минимум достигается в единственной точке, причём
вектор минимизирующий функцию
x
определён по формуле
yfyfp
1
.
Направление, заданное вектором
p
, является направлением спуска функции
xf
, так как
0,,
ppyfpyf
. В силу строгой выпуклости функции
xf
, квадратичная
функция
x
в малой окрестности точки
y
аппроксимирует заданную функцию гораздо,
точнее, по сравнению с линейной функцией. Поэтому, если точка
y
находится в
достаточно малой окрестности точки минимума функции
xxf
, то движение из точки
y
в направлении вектора
p
позволяет достичь более существенного убывания и
получить более точное решение функции, чем направление вектора
yfp
-
направление антиградиента.
Поэтому итерационный процесс для построения последовательностей приближённых
функций
xf
kkkkk
xfxfxx
1
1
, (6),
k
0
,
0k
,
1
, …
будет более эффективным, чем градиентный метод. Метод (6) называется методом
Ньютона с регулировкой шага или обобщённым методом Ньютона. Обычный метод
Ньютона соответствует случаю, когда
1
k
. Наиболее часто построение приближённых
последовательностей по методу (6) осуществляется в два этапа:
1) решается система линейных алгебраических уравнений и определяется вектор
k
p
kkk
xfpxf
;
2) последовательность приближений строим по правилу
kkkk
pxx
1
;
Рассмотрим способ выбора параметра
k
:
полагаем
1
и вычисляем
kk
pxx
;
вычисляем значение функции в точке
x
kk
pxfxf
;
производим проверку неравенства
kkk
pfxfxf ,
, (7),
2
1
0
;
Если неравенство (7) выполняется, то значение
принимают в качестве исходного, если
неравенство (7) не выполняется, то производим дробление
, последовательно умножая
на постоянную
10
до тех пор, пока неравенство (7) не будет справедливым.
k
выбирается из условия минимума функции вдоль выбранного направления движения
kkkkkkk
xfxfxfxfxfxf
11
min
.
Метод Ньютона, как это следует из формулы (6), может быть использован лишь для
минимизации таких функций, у которых существует обратимая матрица вторых
производных, причём эта матрица должна быть ограничена снизу, поэтому мы будем
предполагать, что
xf
удовлетворяет условию:
22
, yMyyxfym
, (8),
0m
,
x
,
n
Ry
.
Справедлива теорема 1:
Теорема 1. Если для минимизации функции
xf
, удовлетворяющей условию (8),
используется метод (6), где
k
выбираются из неравенства (7), то независимо от выбора
начального приближения
0
x
, последовательность
k
x
сходится к точке минимума со
41
сверх линейной скоростью сходимости
lNlN
cxx
, где
,N
c
- постоянные,
1
iN
,
i
,
li 0
и
0
j
при
j
.
Предположим, что матрица вторых производных помимо условия (8) также удовлетворяет
условию Липшица:
yxRyfxf
, (9), где
R
- постоянная,
x
,
n
Ry
, тогда
справедлива теорема 2:
Теорема 2. Если функция
xf
такова, что выполняются условия (8) и (9), то
итерационная последовательность, построенная по правилу (6), в которой
k
выбираются
из неравенства (7) независимо от выбора начальной точки, сходится к решению с
квадратичной скоростью сходимости, то есть, справедлива оценка:
2
1
xx
m
R
xx
kk
.
Если сравнивать метод Ньютона и градиентные методы применительно к решению задач
минимизации выпуклых функций, то окажется, что метод Ньютона обеспечивает гораздо
более высокую скорость сходимости последовательных приближений к решению задачи.
Количество вычислений в методе Ньютона на одном шаге гораздо большее, чем в
градиентных методах. Это происходит за счёт необходимости вычисления и обращения
матриц вторых производных. Однако, если нам необходимо получить решение с
достаточно большой степенью точности, то в методе Ньютона времени потребуется
гораздо меньше, чем в градиентных методах. Вследствие чего метод Ньютона оказывается
гораздо более эффективным.
Метод двойственных направлений.
Построим алгоритм, который не содержит вычислений матриц вторых производных,
сохраняя, тем не менее, скорость сходимости метода Ньютона. Между первыми и
вторыми производными может быть установлена связь с помощью формулы Тейлора для
операторов.
xxyxyxfxfyf ,
, (1),
xyxxy
,
.
Равенство (1) показывает, что если вычислить производные функции
xf
в
произвольных, но достаточно близких друг к другу точках
1
x
,
2
x
, …,
1n
x
и определить
квадратную матрицу
A
по правилу:
iiii
xxAxfxf
11
, (2),
ni ,1
, то матрица
A
будет близка к матрице вторых производных, вычисленных в любой из точек
i
x
в
предположении, что векторы
ii
xx
1
линейно независимы. Поэтому, если
последовательность приближений
k
x
сходится к минимуму функции
xf
, то в
достаточно малой окрестности минимума точек
k
x
,
1k
x
, …,
nk
x
они будут близки друг к
другу. Определив матрицу
k
A
с помощью системы уравнений
11
ikikikikk
xfxfxxA
, (3), можно построить
1k
приближение точки
минимума по формуле:
kkkkk
xfAxx
1
1
, (4),
0k
,
1
, …,
0
k
.
Если матрица
k
A
окажется достаточно близкой к матрице
k
xf
, то направление
вектора
kkk
xfAp
1
будет близко к направлению
kkk
xfxfp
1
, то есть
будет направлением убывания функции
xf
.
Приступим к обоснованию методов вида (4). Пусть функция
xf
имеет непрерывные
первую и вторую производные и имеется некоторая последовательность
k
x
минимизирующая функцию
xf
. Она может быть построена одним из градиентных
42
методов. Этой последовательности ставим в соответствие последовательность
k
y
,
определив её по формуле:
kkk
rxy
, (5), где векторы
k
r
удовлетворяют следующим
условиям:
если
k
определитель матрицы, столбцами которой являются векторы
k
k
r
r
, …,
1
1
nk
nk
r
r
, то при
1 nk
должно выполняться неравенство
k
, где
-
любое достаточно малое число;
0
k
r
при
k
;
В остальном выбор векторов
k
r
произволен.
Лемма. Пусть
k
x
- ограниченная последовательность,
0
1
kk
xx
при
k
и
1 nk
, матрица
K
A
определяется системой уравнений
ikikk
erA
, (6),
1,0 ni
,
ikikik
xfyfe
, где
k
r
и
k
y
- элементы
последовательности (5), тогда
0lim
kk
k
xfA
.
Лекция 13.
Теорема. Если
)( xf
дважды непрерывно дифференцируемая функция, удовлетворяющая
условию
,0,),)((
22
myMyyxfym
матрица
k
A
,при
1 nk
определяется из
системы уравнений (6) и удовлетворяет условию
,0))(),((
1
kkk
xfxfA
k
определяется
из условия,
),),(()()(
kkk
pxfxfxf
то независимо от выбора начальной точки
0
x
для последовательности
,..2.1,0,0),(
1
1
kxfAxx
kkkkkk
справедливы
утверждения:
,..2,1,0),()(
1
kxfxf
kk
и
,0
k
k
xx
причем скорость сходимости
сверхлинейна
,..
lNNlN
Cxx
где
NC,
–некоторые константы,
l
lN
,1
и
.0
i
i
Требования, которым должны удовлетворять векторы
,
k
r
используемые для определения
матрицы
k
A
не являются жесткими и предоставляют большую свободу в выборе этих
векторов. В качестве векторов
k
r
можно взять векторы, направленные по координатным
осям, также в качестве векторов
k
r
можно взять
).(
1
1 kkkkkk
xfAxxr
Допустимы и
другие способы построения векторов
.
k
r
Трудоемкость вычисления вектора
k
p
в значительной степени определяет трудоемкость
метода двойственных направлений. Для построения этих векторов, т.е. построения
матрицы
,
1
k
A
используется понятие двойственности, а именно: говорят, что базис
43
1
,..,
nkk
ss
–двойственен к базису
,,..,
1nkk
ll
если выполняется условие
.1),(,,0),(
jiji
lsjils
Используя это понятие, процесс определения матрицы
1
k
A
превращается в
последовательное вычисление некоторых рекуррентных соотношений. Методы
двойственных направлений позволяют решать задачу минимизации строго выпуклой
функции, не зависимо от выбора начального приближения, причем скорость сходимости к
решению сверхлинейна.
Эти методы по скорости сходимости совместимы с методами Ньютона. Однако, для их
реализации в памяти машины необходимо держать фактически две
nn
матрицы.
Глава Решение задачи Коши для обыкновенных дифференциальных уравнений.
§1. Постановка задачи. Метод Эйлера.
Будем рассматривать задачу Коши для систем обыкновенных уравнений
00
0
)(
,),,(
yxy
Xxxyxfу
(1)
Систему (1) запишем покоординатно:
niyxy
niyyxf
dx
dy
ii
ni
i
,1,)(
,1),,..,,(
00
1
(2)
Из теории дифференциальных уравнений известны условия, гарантирующие
существование и единственность задачи Коши.
Предположим, что функции
niyyxf
ni
,1),,..,,(
1
непрерывны по всем аргументам в
замкнутой области
},1,,{
00
nibyyXxxD
ii
Из непрерывности функций
i
f
в замкнутой области следует их ограниченность, т.е.
существует константа
0M
такая, что всюду в области
D
выполняется неравенство
.,1,),..,,(
1
niMyyxf
ni
Предположим, кроме того, что в области
D
функции
i
f
удовлетворяет условию Липшица
по аргументам
n
yy ,..,
1
,
)...(),..,,(),..,,(
1111 nnnini
yyyyLyyxfyyxf
для
точек
),..,,(
1 n
yyx
и
).,..,,(
1 n
yyx
Если выполнены сформулированные выше требования, то в области
D
существует
единственное решение задачи Коши (1)
).(),..,(
11
xyyxyy
nn
В дальнейшем будем предполагать, что решение задачи (1) существует, единственно и
обладает необходимыми свойствами гладкости.
Заменим непрерывную область
Xx ,
0
разностной (или дискретной) областью.
Рассмотрим сетку
},1,0,,,{
101
NNhkkkh
xXxNkxxxx
h
– сетка,
k
x
– узел сетки,
h
– шаг сетки. В данном случае сетка равномерная.
Если
kkk
hxx
1
и
,
ik
hh
хотя бы для одной пары чисел
k
и
,i
то сетка называется не
равномерной.
Непрерывную модель
),( yxfy
заменим дискретной моделью, которая строится
следующим образом: решение задачи рассматривается в узлах сетки, а именно,
вычисляются значения
).(),..,(),(
10 N
xyxyxy
Внутри дискретной модели будем работать только со значениями сеточной задачи, а
именно, будем рассматривать величины
N
zzz ,..,,
10
– решения дискретной модели.
44
Задача состоит в следующем: построить сеточную модель таким образом, чтобы
выполнялось условие
),())(,()(
)(
2
)()()()(
.,0),(
)(
2
1
hrxyxhfxy
xy
h
xyhxyxyxy
Nkxgz
i
kkkiki
hkkikihkiki
kk
(3)
где
).(
2
)(
2
)(
hk
i
k
xy
h
hr
При
,)(,0
kki
zxyh
а
.)(
kk
zxy
Из (3) получим
nizxhfzz
kki
i
k
i
k
,1),,(
1
или в векторном виде:
.,1,0),,(
001
yzNkzxhfzz
kkkk
(4)
(4) – метод Эйлера решения задачи Коши для обыкновенных дифференциальных
уравнений.
Определим погрешность
,,0),( Nkxyz
kkk
нам нужно оценить
k
и установить,
стремится ли
.0
0
h
k
Из (4) имеем
kkkkk
zxhfzz
),(
1
(5)
здесь
k
– погрешность округления.
Будем предполагать, что
1,0, Nk
k
(6)
Из (3) следует
),())(,()()(
1
hrxyxhfxyxy
kkkkk
(7)
где
)(hr
k
– погрешность аппроксимации метода Эйлера.
2
)( Mhhr
k
(8)
где
M
– константа не зависящая от
.h
Из (5) и (7) получим
0
)),(()))(,()(,(
0
1
hrxyxfxyxLf
kkkkkkkk
(9)
Предположим, что функции
),( yxf
– удовлетворяют условию Липшица с постоянной
,B
а именно:
.)()(,(),(
kkkkkk
BxyzBxyxfzxf
Из (9) получим
kkk
hB
1
(10)
где
Из (10) , при
0k
получим
.)1()1(
2
1
hB
hB
При
kkk
k
(,
1
Из (11) так как
(12) – выполняется при обозначим так как правая часть неравенства (12) не зависит от то
из (12) следует
если погрешность округления то тогда из (13) следует, что метод Эйлера имеет первый
порядок сходимости по
.h
Замечание1. Если
),( yxf
константа, то
0
k
и
.0)( hr
k
В этом случае метод Эйлера
будет точным.
Замечание2. Если отрезок
],[
0
Xx
большой и константа
,1B
то чтобы
0
k
нужно брать очень маленьким шаг
,h
но в этом случае из-за большого числа шагов,
увеличивается погрешность округления.
45
Замечание3. Обозначим
,)(
h
Mhh
если
– фиксированная величина, то
оптимальный шаг для минимизации этой величины нужно брать по формуле
.
M
h
опт
В этом случае
.2
M
опт
Замечание4. Модификация метода Эйлера.
.
)),,(),((
2
00
11
)1(
1
yz
zxfzxf
h
zz
kkkkk
i
k
(14)
(14) – метод Эйлера-Мултона или неявный метод Эйлера.
Нам нужно вычислить
1k
z
и
1k
z
есть в правой части. Для практической реализации
используют итерационную обработку этого метода, то есть
.1,0...,2,1,0
,,0
),,(),((
2
)0(
10
)(
11
)1(
1
Nki
zzz
zxfzxf
h
zz
kk
i
kkkkk
i
k
Метод Эйлера является одношаговым методом решения задачи Коши. Для повышения
порядка точности можно использовать многошаговые методы.
14__________________________________________________________________________
§2 Методы типа Рунге –Кутта.
Рассмотрим дифференциальное уравнение
Xxxyxfy
0
),,(
(1)
и начальное условие
.)(
00
yxy
(2)
Здесь
f
– непрерывно дифференцируемая функция
.:
],[:
0
nn
nn
RRy
RRXxf
Предположим, что решение задачи (1),(2) существует и единственно, и обладает
необходимыми условиями гладкости.
Общая характеристика методов Рунге – Кутта:
1. Методы не чувствительны к сетке. Сетка может быть равноотстоящей и не
равноотстоящей.
2. Точность методов может быть любой.
3. Методы универсальные, гибкие, адаптируемые к любому классу задач.
4. Методы одношаговые.
Будем искать значение приближения решения задачи (1),(2) лишь в фиксированных
точках
Nix
i
,0,
этого отрезка.
Выберем равноотстоящую сетку:
].[,,0,
0
0
h
xx
NNiihxx
i
Связь между двумя соседними значениями функции
)(xy
дает следующее очевидное
равенство:
1
1
.)()(
n
n
x
x
nn
tdtyxyxy
(3)
Поскольку при построении одинаковых методов используется информация о решении
лишь на отрезке длинной в один шаг, то в формуле (3) можно опустить индекс,
46
означающий номер шага и переписать формулу (3) в соответствующем виде:
.)()(
hx
x
tdtyxyhxy
Обозначим
yxyhxy )()(
и, учитывая (1), используя замену
,hxt
последнее
равенство можно переписать в виде:
1
0
.))(,()(
dhxyhxfhdttyy
hx
x
(4)
Если значение
)(xy
известно, то чтобы найти
),( hxy
нам нужно найти поправку
.y
Для построения интересующей нас формулы введем три набора параметров:
r
,..,,
32
)(
121
3231
21
,..,,
,
rrrr
)(
.,..,,
21 r
ppp
)( p
С помощью параметров
)(
и
)(
составим величины:
)....,(
),(
),(
1111
12122
1
rrrrrr
kkyhxhfk
kyhxhfk
yxhfk
Если параметры
)(
и
)(
выбраны, то значения
r
kkk ,..,,
21
вычисляются
последовательно.
Заметим, что величины
)...,(
1111
iiiiii
kkyhxhfk
вообще говоря не равны
величинам
.,1)),(,( rihxyhxhf
ii
Однако при удачном выборе
)(
их можно
трактовать как умноженные на
h
приближенные значения интегрируемой функции
f(x+
).(, hxyh
Поэтому можно надеяться, что с помощью параметров
p
нам удастся создать такую
комбинацию величин, которая будет являться квадратурной суммой и позволит найти
приближенное значение интеграла (4):
r
i
ii
kpy
1
.
(5)
Отсюда становится понятен смысл введенных
).(),(),( p
Величины
i
r
i
ir
kpyh
1
)(
представляют собой погрешность приближенного
равенства (5).Если правая часть уравнения (1) является достаточно гладкой функцией, то
)(h
r
будет иметь достаточно большое число производных.
Запишем разложение
)( h
r
в ряд Маклорена:
).(
!
)0(
!
)(
)1(
1
)(
0
h
s
h
j
h
h
s
s
i
s
j
j
r
(6)
Если удастся выбрать
)(),(),( p
так, что
,,0,0)0(
)(
sj
j
r
а
,0)(
)1(
h
s
r
то
погрешность в формуле (5) будет величиной порядка
.
1s
h
),(
)!1(
)(
)1(
1
h
s
h
h
s
s
r
(7)
47
s
порядок или степень точности данного метода типа Рунге-Кутта. Такое определение
точности связано с тем, что если на каждом шаге погрешность расчетной формулы
данного метода типа Рунге-Кутта имеет порядок
,
1s
h
а функция
),( yx
является гладкой
в окрестности решения, то погрешность данного метода (часть погрешности
приближенного решения, определяемая неточностью расчетной формулы) на конечном
отрезке будет величиной порядка
.h
Для построения по методу Рунге-Кутта при данном
r
одношаговых правил возможно более высокого порядка точности
s
выражают величины
ss
hh ,
1
Требуют, чтобы для любой гладкой
),( yxf
обращалось в нуль возможно большее число
этих величин. Иными словами,
)(),(),( p
выбираются исходя из требования, чтобы
разложение
...)(
!3
)(
!2
)(
!1
)()(
32
xy
h
xy
h
xy
h
xyhxyy
(8)
и разложение линейной комбинации
r
i
ii
kp
1
совпадали для
),( yxf
до членов с возможно
более высокими степенями
.h
Записать в общем виде систему уравнений для определения
)(),(),( p
крайне
затруднительно. Поэтому рассмотрим лишь несколько примеров построения одношаговых
правил по методу Рунге-Кутта.
Метод первого порядка точности.
При
1r
приближенное значение формулы (5) имеет вид:
).,(
11`1
yxhfpkpy
Функция может быть записана так:
)()(
),()()(
),()()()(
1
11
11
hxyh
yxfphxyh
yxhfpxyhxyh
)()0(
1
xy
– не зависит от введенного параметра и естественно не может обратиться в
нуль за счет его выбора.
Величина
).,()1(),()()0(
111
yxfpyxfpxy
И эта величина для
),( yxf
обращается в нуль при
.1
1
p
Тогда (5) принимает в случае
1r
вид:
).,( yxhfy
Погрешность метода из формулы (7) в этом случае имеет вид:
).(
!2
)(
2
1
h
h
h
Т. о. построенный одношаговый метод Рунге-Кутта точностью совпадает с методом
Эйлера.
Метод второго порядка точности.
При
2r
формула (5) принимает вид:
).,(),(
1212212211
kyhxhfpyxhfpkpkpy
Нам нужно определить параметры Для этого разложим по степеням
Из (8) будем иметь:
...))(2(
6
)(
2
2
32
fffffffff
h
yff
h
hfy
yxyyyxyxxxyx
(9)
...)(
!1
1212221
kfhf
h
hphfphfpy
yx
(10)
Требуем совпадения равенств (9) и (10) при производных функции
),( yxf
до членов с
возможно более высокими степенями
.h
В результате сравнения в этих разложениях
коэффициентов при получим равенство:
.1
21
pp
48
При
,
2
1
:
22
pf
x
При
.
2
1
:
212
pff
y
Из вида разложений (9) и (10) следует, что для
2r
нельзя при
),( yxf
добиться за счет
выбора введенных четырех параметров совпадения четырех коэффициентов со степенью
.
3
h
Восемь уравнений и четыре неизвестных. Следовательно, правило Рунге-Кутта имеет
при
2r
второй порядок точности. Таких правил существует бесконечное множество.
Например, положив
.1,1,
2
1
2
1
21212
pp
И тогда
).,(),,();()(
2
1
121
3
21
kyhxhfkyxhfkhkky
И такие правила – это аналог квадратурной формулы трапеции для приближенного
вычисления интеграла (3). Положив
).
2
1
,
2
1
(),,(
);(
2
1
,
2
1
,01
121
3
221112
kyhxhfkyxhfk
hkypp
В этом случае правило Рунге-Кутта – аналог квадратурной формулы средних
прямоугольников для приближенных вычислений интеграла (3).
Методы третьего порядка точности.
332211
kpkpkpy
и одна из формул имеет вид:
).2,(),
2
1
,
2
1
(),,(
),()4(
6
1
213121
4
321
kkyhxhfkkyhxhfkyxhfk
hkkky
Методы четвертого порядка точности.
Наиболее употребительным методом решения задачи Коши является следующее правило:
).,(
),
2
1
,
2
(
),
2
1
,
2
(
),,(
),()22(
6
1
34
23
12
1
5
4321
kyhxhfk
ky
h
xhfk
ky
h
xhfk
yxhfk
hkkkky
В литературе этот метод называется методом Рунге-Кутта без ссылок на порядок и
точность.
§3 Сходимость и оценка погрешности одношаговых методов.
Можно считать, что малость ошибки при одном шаге вычислений вообще говоря не
гарантирует малость ошибки при счете на большом числе шагов. При переходе от шага к
шагу, ошибки, допущенные на каждом шаге вычислений, накапливаются и могут быстро
расти. Поведение накопившейся погрешности при этом зависит и от особенностей
решаемой задачи, и от свойств избранного вычислительного метода, а также от
погрешности задания начальных данных и неточностей выполнения вычислений. Оценим
поведение этой погрешности в случае одношагового метода.
Рассмотрим задачу Коши:
49
.)(
,),,(
00
0
yxy
Xxxyxfy
(1)
Для приближенного решения в точках
n
xxx ,..,,
21
решение
)(xy
этой задачи при
использовании одношагового метода можно определить по формуле:
).,,(
1 nnfnn
yxhyy
(2)
Будем предполагать существование в плоскости
xy
выпуклой в направлении оси
y
области
,D
содержащей внутри себя
grad
точного и приближенного решений исходной
задачи и такой, что функция
),( yxf
имеет в
D
непрерывную производную по
y
такую,
что
.
),(
F
y
yxf
При решении реальной задачи вычисления по формуле (2) выполняются,
как правило, не точно (из-за ошибок округления) и величины
n
yyy
~
,..,
~
,
~
21
не будут
удовлетворять разностному уравнению (2), а будут связаны соотношением
nnnfnn
yxhyy
)
~
,,(
~~
1
(3)
Здесь
).(
~
00
xyy
Величина
,
n
определяемая формулой (3) – погрешность округления на
n
–ом шаге.
Лекция 15.
Обозначим через
)(
~
xy
n
решение дифференциального уравнения, удовлетворяющего
начальному условию
,
~
)(
~
nnn
yxу
а через
n
r
погрешность формулы (2) на
n
ом шаге.
Величина
n
r
характеризует локальную погрешность избранного метода или ошибку
приближенного решения на одном шаге без учета погрешности округления и может быть
введена посредством равенства:
.))(
~
,,()(
~
)(
~
1 nnnnfnnnn
rxyxnhxyxy
(4)
Так как
,
~
)(
~
nnn
yxy
то из последнего равенства с учетом (3) получим:
.
~
)
~
,,(
~
)(
~
11 nnnnnnfnnn
ryryxhhyxy
(5)
Поэтому, локальная ошибка приближенного решения с учетом погрешности округления
может быть задана по формуле:
.
~
)(
~
11 nnnnn
ryxy
(6)
Оценив сумму, стоящую в правой части равенства (6), мы можем сделать заключение о
малости ошибки на одном шаге вычислений. Чтобы иметь возможность судить о величине
погрешности при счете на большое число шагов, нужно поучить удобное представление
для разности
.
~
)(
nnn
yxy
Эта разность – погрешность приближенного решения. Она,
очевидно, должна выражаться через погрешность формулы, погрешность округления и
погрешность задания начального условия –
.
~
)(
000
yxy
Для величины
)).(
~
)(
~
()(
~
)(
~
)(
1
10 n
n
i
ininnnnn
xyxyxyxyyxy
(7)
Правая часть равенства (7) выражается через разности различных дифференциальных
уравнений, для оценки этой разности воспользуемся Леммой.
Лемма. Если
)(),(
21
xYxY
решения дифференциального уравнения
),,( yxfy
то
справедливо равенство:
.10
,))()(()()(
))()(()(,(
2121
121
dyxYxYxYx
y
f
YYYY
(8)
Запишем равенства:
)).(,()(
)),(,()(
22
11
xYxfxY
xYxfxY
50
и применим теорему Лагранжа о конечных
приращениях.
.10)),()(()(,())()((
12121
xYxYxYx
y
f
xYxY
Проинтегрировав последнее равенство на отрезке
,
мы получим формулу (8).
На основе доказанной Леммы мы можем утверждать:
если
,0),( yxf
y
то графики решения дифференциального уравнения (1) удаляются друг
от друга;
если
,0),( yxf
y
то графики решения сближаются.
На основании равенства (8) можно записать следующее представление для разностей,
входящих в формулу (7).
,))(
~
(()(
~
)(
0
0
0000
n
x
x
dy
y
f
nn
xyxyxyxy
(9)
.))(
~
)(
~
()(
~
)(
~
11
n
x
i
x
i
dy
y
f
iiiinni
xyxyxyxy
(10)
Здесь
.
))(
~
)(
~
()(
~
,(
,
))(
~
)(()(,(
11000
y
xyxyxyxf
y
f
y
xyxyxyxf
x
f
iiiii
Учитывая, что
.)(
~
)(
~
,)(
~
)(
111
0000
iiiiii
rxyxy
xyxy
Формулу (7) для погрешности приближенного решения можно записать в виде:
.)(
1
110
0
0
n
i
dy
y
f
ii
dy
y
f
n
n
x
i
x
i
n
x
x
r
(11)
Из (11) мы получим зависимость характера накопления погрешности при счете на
n
шагов в зависимости от особенностей исходной задачи.
В случае, когда в области
D
частная производная
0),( yxf
y
влияние погрешности
,
0
начального приближения и каждой из ошибок
,,
11 ii
r
допущенных на предыдущих
шагах возрастает.
Если
,0),( yxf
y
то влияние каждой из этих ошибок наоборот ослабевает.
Если предположить выполнение неравенств:
,,,1,,,
0
0
n
xx
NNkrr
kk
и всюду в области
D
выполняется
,F
y
f
то на основании равенства (11) получим
.)(
1
0
n
i
hFhF
т
inn
r
(12)
Последнюю оценку можно записать в виде:
.
1
1
)(
Fh
hF
hF
n
n
n
r
(13)
При возрастании поведение правой части (13) существенно зависит от
Если
,0F
то с ростом
n
слагаемое
hF
n
уменьшается, а второе слагаемое правой части
(13) остается ограниченным даже при интегрировании на сколь угодно большое число
шагов.
Приведем более грубую, но более наглядную оценку погрешности и выясним поведение
этой погрешности при неограниченном уменьшении шага
h
. Если учесть, что
,
Lnh
hF
in
где
),0,max( FL
то неравенство (12) можно записать в виде:
.))((
Lnh
n
nr
(14)
51
или
,))((
)(
0
0
xxL
n
n
n
h
xx
hr
учитывая, что
,
00
xXxx
n
последняя оценка
примет вид:
.))(((
)(
01
0
xXL
xX
h
r
(15)
На основании последней оценки можно утверждать, что, если выполняются условия:
1)
,0
0
n
2)
,0
2
0
h
h
3)
.0
0
n
n
Таким образом, приближенное решение задачи (1), полученное одношаговым методом
вида (2) равномерно на отрезке
Xx ,
0
сходится к точному решению этой задачи. В
реальных ЭВМ величина
обычно фиксируется, поэтому отношение
h
возрастает с
уменьшением
,h
полученная оценка (5) вообще говоря, является завышенной. Ошибки
округления могут иметь разные знаки и частично компенсировать друг друга. Однако
практическое вычисление показывает, если при значительном уменьшении шага
,h
точность вычислений не уменьшать, то погрешность приближенного решения начинает
расти. Поэтому всегда нужно согласовывать выбор расчетной формулы и величину шага с
требуемой точностью вычислений.
