корнюшин.численные методы
.pdf51
Численные методы решения системы (1) условно разбивают на две группы, выделяя прямые и итерационные методы (имеются, естественно, и гибридные методы). Прямые методы позволяют за конечное число действий получить точное решение системы уравнений, если входная информация (правая часть уравнения f и элементы aij матрицы A) задана точно, и вычисления ведутся без округления. Простейший пример прямого метода – метод прогонки. Конечно, прямые методы также дают решение с определенной точностью, которая зависит от погрешностей округления, т.е. от машины, и от характера вычислительной устойчивости, что зависит от самого метода.
Итерационный метод позволяет найти приближенное решение системы путем построения последовательности приближений (итераций), начиная с некоторого начального приближения. Само приближенное решение является результатом вычислений, полученным после конечного числа итераций.
Выбор того или иного численного метода зависит от многих обстоятельств – от имеющихся программ, от вида матрицы, от типа расчета и др. Поясним слова «тип расчета». Возможны разные постановки задачи:
1)найти решение одной конкретной задачи (1);
2)найти решение нескольких вариантов задачи (1) с одной и той же матрицей A и разными
правыми частями f. Может оказаться, что неоптимальный для одной задачи (1) метод является весьма эффективным для многовариантного расчета.
При многовариантном расчете можно уменьшить среднее число операций для одного варианта, если хранить некоторые величины, а не вычислять их заново для каждого варианта. Это, конечно, зависит от машины, от объема ее оперативной памяти.
Из сказанного ясно, что выбор алгоритма должен зависеть от типа расчета, от объема оперативной памяти машины и, конечно, от порядка системы. Качество алгоритма определяется тем машинным временем, которое требуется для нахождения решения системы (1). Естественно выбирать тот метод, для которого время решения минимально по сравнению с другими методами. Однако время расчета зависит от многих факторов: от числа арифметических и логических действий, которые нужно затратить для получения решения с заданной точностью, от быстродействия и объема оперативной памяти машины, от качества программы. При теоретических оценках качества алгоритмов их сравнение приводится по числу Q(ε) арифметических действий, достаточных для нахождения решения задачи с заданной точностью
ε>0.
3.5.2.Прямые методы
3.5.2.1.Метод Гаусса
Имеется несколько вычислительных вариантов метода Гаусса, основанного на идее последовательного исключения. Процесс решения системы линейных алгебраических уравнений
Ax=f, или
N |
|
∑aij x j = f j , i=1, 2,…, N, |
(1) |
j=1
по методу Гаусса состоит из двух этапов.
Первый этап (прямой ход). Система (1) приводится к треугольному виду
B+x=ϕ, (2)
где x=(x1,…, xN) – неизвестный, ϕ=(ϕ1,…, ϕN) – известный вектор, B+ – верхняя треугольная матрица.
Второй этап (обратный ход). Неизвестные xN, xN-1,…, x1 определяются по формулам (2)
п.3.5.1.
Перейдем к подробному изложению метода. Первый шаг метода Гаусса состоит в исключении из всех уравнений, кроме первого, неизвестного x1. Предположим, что a11 ≠ 0,
разделим первое уравнение (1) (i=1) на a11 и запишем систему (1) в виде |
|
x1+b12x2+…+b1NxN=ϕ1, b1j=a1j/a11, j=2, 3,…, N, ϕ1=f1/a11, |
(3) |
|
|
|
|
|
|
|
|
52 |
|
|
|
|
ai1x1+ai2x2+…+aiNxN=fi, i=2, 3,…, N. |
(4) |
|||||||||
Умножим уравнение (3) на ai1, где i – любое из чисел 2, 3,…, N, и вычтем полученное |
|||||||||||
уравнение из I-го уравнения (4): |
|
|
|
|
|
|
|
|
|||
|
(ai2-ai1b12)x2+…+(aiN-ai1b1N)xN=fi-ai1ϕ1, i=2, 3,…, N. |
||||||||||
Вводя обозначения |
|
|
|
|
|
|
|
|
|
|
|
a(1) |
= a |
ij |
−a |
b |
, f (1) |
= f |
i |
−a |
ϕ |
, i, j=2, 3,…, N, |
(5) |
ij |
|
|
i1 1 j |
i |
|
|
i1 1 |
|
|
||
перепишем полученную систему уравнений (эквивалентную системе (1)) в виде
x1+b12x2+…+b1NxN=ϕ1,
ai(21) x2 +... + aiN(1) xN = fi(1) , i=2, 3,…, N.
Первый столбец матрицы этой системы состоит из нулей, кроме первого элемента, равного единице при i=j=1.
Второй шаг состоит в исключении x2 из системы
a22(1) x2 +... + a2(1N) xN = f2(1) ,
.................................. (6)
aN(1)2 x2 +... + aNN(1) xN = f N(1) .
Для этого разделим первое уравнение на a22(1) :
x2+b23x3+…+b2NxN=ϕ2,
ϕ2 = f2(1) / a22(1) ,b2 j = a2(1j) / a22(1) , j = 3,..., N,
умножим его затем на (- ai(21) ) и сложим с уравнением
ai(21) x2 + ai(31) x3 +... + aiN(1) xN = fi(1) ,i = 3,4,..., N.
В результате получим систему
x2+b23x3+…+b2NxN=ϕ2, (7)
|
ai(32) x3 +... + aiN(2) xN |
= fi(2) , |
|
|
|
||||
a(2) |
= a(1) |
−a(1)b |
, f (2) |
= f (1) |
−a(1)ϕ |
2 |
, |
(8) |
|
ij |
ij |
i2 2 j |
i |
|
i |
i2 |
|
|
|
i = 3,4,..., N.
Для x3, x4,…, xN имеем систему (N-2)-го порядка, аналогичную системе (6) (N-1)-го порядка
для x2, x3,…, xN.
Продолжая рассуждения, после (N-1)-го шага (т.е. после исключения x1, x2,…, xN-1) получим
aNN( N −1) xN = fN( N −1) , или xN=ϕN, ϕN = f N( N −1) / aNN( N −1) . |
(9) |
В итоге приходим к системе (2) с верхней треугольной матрицей:
x1 +b12 x2 +b13 x3 +... +b1N xN =ϕ1 , x2 +b23 x3 +... +b2 N xN =ϕ2 ,
.......................................... (10)
xN −1 +bN −1,N xN =ϕN −1 , xN =ϕN .
Обратный ход метода Гаусса состоит в определении всех xi из системы (10) с верхней треугольной матрицей. Нетрудно показать, что изложенный выше метод Гаусса можно применять в том случае, когда все главные миноры матрицы отличны от нуля.
Подсчитаем число умножений и делений в методе Гаусса. Рассмотрим сначала прямой ход. На первом шаге требуется Q1=N2 делений и умножений, второй шаг требует Q2=(N-1)2 действий и т.д. Всего надо сделать N шагов прямого хода, затратив на это
∑kN=1 (N −k +1)2 = ∑sN=1 s2 = N (N +1)(2N +1) / 6 умножений и делений. Для обратного хода, очевидно, нужно сделать N(N-1)/2 умножений. Таким образом, для решения системы уравнений
(1) требуется Q=N(N2+3N-1)/3 умножений и делений. Примерно столько же потребуется сложений.
Рассмотрим в качестве примера систему трех уравнений
53 |
|
2x1+4x2+3x3=4, |
(11) |
3x1+x2-2x3=-2, |
(12) |
4x1+11x2+7x3=7. |
(13) |
Прямой ход. Первый шаг. Разделим первое уравнение на a11=2: |
|
x1+2x2+1,5x3=2. |
(14) |
Умножим (14) на -3 и сложим с (12), затем умножим (14) на -4 и сложим с (13): |
|
-5x2-6,5x3=-8, |
(15) |
3x2+x3=1. |
(16) |
Мы получили систему второго порядка. |
|
Второй шаг. Разделим (15) на -5: |
|
x2+1,3x3=1,6. |
(17) |
Умножим (17) на -3 и сложим с (16): |
|
-2,9x3=-5,8. |
(18) |
Третий шаг. Делим (18) на -2,9: |
|
x3=2. |
|
В результате получили систему |
|
x1+2x2+1,5x3=2, x2+1,3x3=1,6, x3=2
с верхней треугольной матрицей |
|
|
1 |
2 |
1,5 |
|
1 |
|
0 |
1,3 . |
|
|
0 |
|
0 |
1 |
Обратный ход. Из системы последовательно находим: x3=2; x2=1,6-1,3*2=-1; x1=2-2x2- 1,5*x3=1. Таким образом, решение системы (11) – (13) найдено:
x1=1, x2=-1, x3=2.
3.5.2.2. Метод квадратного корня
Этот метод пригоден для систем |
|
Au=f |
(19) |
с симметричной матрицей A. Матрица A разлагается в произведение |
|
A=SDS*, |
(20) |
где S=(sij) – нижняя треугольная (sij=0 при j>i), D=(dij) – диагональная матрицы (dij=diδij, δij
– символ Кронекера, δij=0, i ≠ j, δii=1); без ограничения общности можно положить sii=1 для всех i=1, 2,…, N. Решение уравнения Au=f сводится к последовательному решению двух систем:
В силу (20) имеем |
Sdy=f, S*u=y. |
|
|
|||||
|
|
|
|
|
|
|
|
|
ai1=d1si1, i ≥1, |
|
(21) |
|
|||||
i−1 |
|
|
|
|
|
|
|
|
aii = di sii + ∑dk |
|
sik |
|
2 ,i = 2,3,..., N, |
(22) |
|||
|
|
|||||||
k =1 |
|
|
|
|
|
|
|
|
j−1 |
_ |
|
|
|
|
|
|
|
aij = d j sij + ∑dk sik |
s jk ,i ≥ j +1, j = 2,3,..., N −1, |
(23) |
||||||
k =1 |
|
|
|
|
|
|
|
|
_ |
|
|
|
|
|
|
_ |
|
где sik – комплексно сопряженное к sik число, т.к. S*=( s ji ) . Положим |
||||||||
^ |
|
|
|
|
|
|
|
|
sik = dk sik |
(24) |
|
||||||
и выберем |
i−1 ^ |
|
|
|
|
|
||
|
|
|
_ |
|
|
|||
d1=a11, di = aii − ∑sik |
sik ,i = 2,3,.., N. |
(25) |
||||||
k =1
54
Тогда условия (22) будут выполнены, а условия (21) – выполнены для i=1. Запишем оставшиеся соотношения следующим образом:
^ |
|
|
|
si1 |
= ai1 ,i = 2,3,..., N, |
||
|
^ |
j−1 ^ |
_ |
aij |
= sij |
+∑sik |
s jk ,i = 3,4,..., N; j = 2,3,...,i −1, |
k =1
^
и из (24) найдем sij = sij / d j . Отсюда следует алгоритм.
Находим сначала d1=a11.
Для i от 2 до N выполняем действия:
→ begin
^ |
^ |
|
1.si1 |
= ai1 , si1 = si1/ d1. |
|
2.if |
i = 2 then |
go to3, else |
for j = 2 to |
i −1 do : |
|
^ |
j−1 |
^ |
_ |
sij = aij −∑sik s jk |
|||
|
k =1 |
|
|
|
i−1 |
^ |
_ |
3.di aii −∑sik sik
k =1
→ end.
^
, sij = sij / d j .
Замечание. Вычисляемые элементы sij и di могут быть размещены на месте aij и aii, если элементы матрицы A не нужно сохранять. Дополнительная память требуется для размещения при
^
фиксированном i чисел sij (j=1, 2,…, i-1); очевидно, нужно не более N-1 ячеек памяти для их
хранения.
Алгоритм SDS*-разложения может оказаться численно неустойчивым даже при выборе главного элемента. Однако в двух случаях можно гарантировать устойчивость алгоритма – для положительно определенных матриц A и матриц с диагональным преобладанием.
Определитель матрицы равен
N |
|
det A = ∏dii sii2 . |
(26) |
i=1
Метод квадратного корня требует порядка N3/3 арифметических действий, т.е. при больших N он вдвое быстрее метода Гаусса и занимает вдвое меньше ячеек памяти. Это обстоятельство объясняется тем, что метод использует информацию о симметрии матрицы.
3.5.2.3. Связь метода Гаусса с разложением матрицы на множители
Пусть дана невырожденная матрица A размера N*N. Представим ее в виде произведения
|
|
|
A=BC, A=(aij), B=(bij), C=(cij), |
|
(27) |
|
|
|
|
||||||||
где B и C – треугольные матрицы вида |
|
|
|
|
|
|
|
|
|
|
|
||||||
|
b11 |
0 |
0 |
... |
0 |
|
|
1 |
c12 |
c13 ... |
c1N |
|
|||||
|
b |
b |
0 |
... |
0 |
|
|
|
0 |
1 |
c |
23 |
... |
c |
|
|
|
|
|
21 |
22 |
|
|
|
|
|
|
|
|
|
|
|
2 N |
||
B = |
|
|
b32 |
b33 |
... |
0 |
|
C = |
|
0 |
0 |
1 ... |
c3N |
|
|||
b31 |
, |
|
, |
||||||||||||||
|
|
.... .... .... ... .... |
|
|
.... .... .... ... .... |
|
|||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
|
|
|
bN 2 |
bN 3 |
... |
|
|
|
|
0 |
0 |
0 ... |
|
1 |
|
||
|
bN1 |
bNN |
|
|
|
|
|||||||||||
т.е. bik=0 при k>i, cik=0 при k<i, cii=1. Из (27) следует, |
что aij |
= ∑kN=1 bik ckj . Преобразуем эту |
|||||||||||||||
сумму двумя способами: |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
55
N |
i−1 |
|
|
N |
|
i−1 |
|
∑bik ckj = ∑bik ckj |
+bii cij + ∑bik ckj = ∑bik ckj +bii cij , |
||||||
k =1 |
k =1 |
|
|
k =i+1 |
|
k =1 |
|
N |
j−1 |
|
|
N |
|
j−1 |
|
∑bik ckj = ∑bik ckj |
+bij c jj + ∑ bik ckj = ∑bik ckj +bij c jj . |
||||||
k =1 |
k =1 |
|
|
k = j+1 |
|
k =1 |
|
Отсюда находим |
i−1 |
|
|
|
|
|
|
|
|
|
|
|
|
||
|
bij = aij −∑bik ckj |
при i ≥ j,b11 = a11 ,c11 =1, |
|||||
|
k =1 |
|
|
|
|
|
|
|
|
1 |
|
|
i−1 |
|
|
|
cij = |
|
aij −∑bik ckj |
при i<j. |
|||
|
|
|
|||||
|
|
bii |
k =1 |
|
|
||
Матрицы B и C найдены.
Решение уравнения Au=Bcu=f сводится к последовательному решению уравнений
Bϕ=f, Cu=ϕ.
Построение матриц B и C и нахождение ϕ=B-1f соответствует прямому ходу, а решение уравнения
Cu=ϕ
соответствует обратному ходу метода Гаусса.
3.5.3.Итерационные методы
3.5.3.1.Метод итераций для решения системы линейных алгебраических уравнений
Пусть дана система линейных алгебраических уравнений
Au=f. (1)
Для ее решения выберем некоторое начальное приближение y0 H и последовательно найдем приближенные решения (итерации) уравнения (1). Значение итерации yk+1 выражается через известные предыдущие итерации yk, yk-1,…. Если при вычислении yk+1 используется только одна предыдущая итерация yk, то итерационный метод называют одношаговым методом; если же yk+1 выражается через две итерации yk и yk-1, то метод называется двухшаговым. Будем считать, что A: H →H – линейный оператор, отображающий конечномерное пространство H (с заданным скалярным произведением) в себя.
Важную роль играет запись итерационных методов в единой (канонической) форме. Любой двухслойный итерационный метод можно записать в следующей канонической форме:
B |
yk +1 − yk |
+ Ayk = f , k = 0,1,... для всех y0 H, |
(2) |
|
|||
|
τk +1 |
|
|
где A: H →H – оператор исходного уравнения (1), B: H →H – линейный оператор, |
|||
имеющий обратный B-1, k – номер итерации, τ1, τ2,…, τk+1,… – итерационные параметры, τk+1>0. Оператор B может, вообще говоря, зависеть от номера k; для простоты изложения мы
предполагаем всюду, что B не зависит от k. |
|
||
Если B=E – единичный оператор, то метод |
|
||
|
yk +1 − yk |
+ Ayk = f , k = 0,1,... для всех y0 H |
(3) |
|
|
||
|
τk +1 |
|
|
называют явным: yk+1 находится по явной формуле |
|
||
|
|
yk+1=yk-τk+1(Ayk-f). |
|
В общем случае, при B ≠ E, метод (2) называется неявным итерационным методом: для |
|||
определения yk+1 надо решить уравнение |
|
||
|
Byk+1=Byk-τk+1(Ayk-f)=Fk, k=0, 1,… |
(4) |
|
Естественно требовать, чтобы объем вычислений для решения системы Byk+1=Fk был меньше, чем объем вычислений для прямого решения системы Au=f.
56
Точность итерационного метода (2) характеризуется величиной погрешности zk=yk-u, т.е. разностью между решением yk уравнения (2) и точным решением u исходной системы линейных алгебраических уравнений. Подстановка yk=zk+u в (2) приводит к однородному уравнению для погрешности:
|
|
|
|
|
|
B |
zk +1 − zk |
|
+ Azk |
= 0, k = 0,1,... , z0=y0-u. |
(5) |
|||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|||||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
τk +1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||
Говорят, что итерационный метод сходится в HD, если |
|
|
||||||||||||||||||||||||||||||
lim |
|
|
|
z |
k |
|
|
|
|
D |
= 0, где |
z |
= |
|
|
(Dz, z), D = D* > 0, D : H → H. |
||||||||||||||||
|
|
|
|
|
|
|||||||||||||||||||||||||||
k→∞ |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
D |
|
|
|
|
|
|
|
|
|
|
|
|||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||||
Обычно задают некоторую относительную погрешность ε>0, с которой надо найти |
||||||||||||||||||||||||||||||||
приближенное решение yk, и прекращают вычисления, как только выполняется условие |
||||||||||||||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
yn −u |
|
|
|
D ≤ε |
|
|
|
y0 −u |
|
|
|
D . |
(6) |
|
|||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||
Если n=n(ε) – наименьшее из чисел, для которых (6) выполняется, то общее число арифметических действий, затрачиваемых для нахождения приближенного решения уравнения (1), равно Qn(ε)=n(ε)q0, где q0 – число действий, затрачиваемых для нахождения одной итерации,
т.е. для решения уравнения (4). Задача состоит в минимизации Qn(ε) путем выбора B и параметров
{τk }.
3.5.3.2. Метод простой итерации
Для решения системы уравнений (1) может быть использован метод простой итерации
|
yk(i+)1 |
|
|
|
N |
|
|
|
|
− f (i) |
|
|
|
|
|
||||
|
= yk(i) −τ |
∑aij yk( j) |
,i =1,2,..., N, |
(7) |
|||||||||||||||
|
|
|
|
|
j=1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
где τ>0 - итерационный параметр. Запишем (7) в операторной форме: |
|||||||||||||||||||
|
yk +1 − yk |
+ Ayk |
= f , k = 0,1,..., |
для всех y0 H. |
(8) |
||||||||||||||
|
|
||||||||||||||||||
|
τ |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Сравнение (8) с (3) показывает, что метод простой итерации задается явной одношаговой |
|||||||||||||||||||
схемой с постоянным параметром τk |
≡τ. |
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||
Существуют и другие варианты метода простой итерации, например, такой: |
|||||||||||||||||||
|
|
|
|
(i) |
1 |
|
N |
|
|
( j) |
|
|
|
|
|||||
|
|
|
|
|
|
|
|
|
|
(i) |
|
|
|||||||
|
|
|
|
yk +1 = |
|
|
∑aij yk |
− f |
. |
|
|||||||||
|
|
|
a |
|
|
||||||||||||||
|
|
|
|
|
|
|
ii |
|
j =1 |
|
|
|
|
|
|
||||
Подставляя сюда |
|
|
|
|
|
|
|
j≠i |
|
|
|
|
|
|
|||||
|
|
N |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
N |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
∑aij yk( j) =∑aij yk( j) −aii yk(i) = ( Ayk )(i) −(Dyk )(i) , |
||||||||||||||||||
|
j=1 |
|
|
j=1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
j ≠i |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
где D=(aiiδij) – диагональная матрица, получаем |
|
|
|
|
|||||||||||||||
|
|
|
|
(i) |
(i) |
|
|
1 |
|
|
N |
|
( j) |
|
|
|
|||
|
|
|
|
|
|
|
|
|
|
|
(i) |
|
|||||||
|
|
|
yk +1 = yk |
− |
|
|
|
|
|
∑aij yk |
− f |
, |
|
||||||
|
|
|
|
a |
|
|
|
||||||||||||
|
|
|
|
|
|
|
|
|
|
ii |
|
j =1 |
|
|
|
|
|
||
или, в каноническом виде, |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
|
|
yk +1 − yk |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
|
|
D |
|
+ Ayk |
= f , k = 0,1,...,τ =1. |
|
|||||||||||||
|
|
|
|
||||||||||||||||
|
|
|
|
τ |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Хотя формально эта схема является неявной (B=D ≠ E), однако D=(aiiδij) – диагональная матрица и потому yk+1 определяется по явным формулам.
3.5.3.3. Метод Зейделя
57
Широкое распространение на практике (в особенности в тех случаях, когда информации о матрице A недостаточно) получил итерационный метод Зейделя в одной из двух форм:
i
∑a y( +j) +
ij k 1
j=1
N |
|
∑aij yk( j) = f (i) , aii ≠ 0,i =1,2,..., N, |
(9) |
j=i+1
i−1 |
N |
|
∑aij yk( j) +∑aij yk( +j)1 = f (i) ,i =1,2,..., N. |
(10) |
|
j=1 |
j=i |
|
Из обеих формул компоненты вектора yk+1 находятся последовательно. Так, из (9) последовательно определяем yk(1+)1 , yk(2+)1 ,..., yk( N+1) :
|
|
|
|
|
yk +1 |
= a |
|
f |
|
|
N |
, |
||
|
|
|
|
|
|
(1) |
−∑a1 j yk |
|||||||
|
|
|
|
|
|
(1) |
|
1 |
|
|
|
( j) |
|
|
|
|
|
|
|
|
|
|
11 |
|
|
j=2 |
|
||
(i) |
|
1 |
|
|
(i) |
|
N |
|
|
( j) |
i−1 |
|
||
|
|
|
|
|
|
( j) |
||||||||
yk +1 |
= |
|
|
f |
|
− ∑ aij yk |
|
−∑aij yk +1 |
,i = 2,..., N. |
|||||
a |
|
|
|
|||||||||||
|
|
|
|
|
|
j |
= + |
|
|
|
= |
|
||
|
|
|
ii |
|
|
i |
1 |
|
|
|
j 1 |
|
||
Пользуясь (10), находим последовательно для i=N, N-1,…, 1
|
|
|
|
yk( N+1) |
|
(i) |
|
1 |
|
|
(i) |
yk +1 |
= |
|
f |
|
|
|
|
||||
|
|
aii |
|
|
|
|
|
|
1 |
|
( N ) |
N −1 |
( j) |
|
|
|
|
|
|||||
= |
|
|
|
f |
|
−∑aNj yk |
, |
|
|
a |
|
|
|||||
|
|
|
|
|
= |
|
|
|
|
|
|
NN |
|
j 1 |
|
|
|
|
i−1 |
|
|
|
N |
|
|
|
−∑aij yk( j) − ∑aij yk( |
+j)1 ,i = N −1,...,1. |
|||||||
|
j=1 |
|
|
j=i+1 |
|
|
||
Запишем этот метод в матричной (операторной) форме. Для этого представим матрицу A в
виде суммы
A=A–+D+A+,
где D=(aiiδij) – диагональная матрица размера N*N, A–=( aij− ) – нижняя треугольная
(поддиагональная) матрица с нулями на главной диагонали, aij− =0 при i ≥ j , aij− =aij при j<i, A+=( aij+ ) – верхняя треугольная (наддиагональная) матрица с нулями на главной диагонали, aij+ =0
при j ≤ i, aij+ =aij при j>i. Из определения A–, D, A+, следует
|
i−1 |
(Dy)(i) = aii y(i) , (A− y)(i) = ∑aij y( j) , |
|
|
j=1 |
N |
N |
( A+ y)(i) = ∑aij y( j) , ((A+ + D) y)(i) = ∑aij y( j) . |
|
j=i+1 |
j=i |
Поэтому уравнение (10) можно записать−в виде |
|
|
((A++D)yk+1)(i)+ (A yk)(i)=f(i), i=1, 2,…, N, |
||
или, в векторной форме, |
− |
|
(A++D)yk+1+ A yk=f. |
|
|
После очевидных преобразований |
− |
|
− |
|
|
(A++D)yk+1+ A yk=(A++D)(yk+1-yk)+ (A |
+(A++D))yk=(A++D)(yk+1-yk)+Ayk |
|
запишем метод Зейделя (10) в каноническом виде |
|
|
(A++D)(yk+1-yk)+Ayk=f, k=0, 1, 2,… |
(11) |
|
Сравнивая (11) с (2), видим, что метод Зейделя (10) соответствует
B=D+A+, τ ≡1,
Т.е. схема (11) является неявной. Однако, т.к. B=D+A+ – треугольная матрица, то итерация yk+1 находится по явным формулам. Аналогично записывается и другой вариант метода Зейделя:
(A − +D)(yk+1-yk)+Ayk=f, k=0, 1, 2,…, |
(12) |
когда B=D+ A − – нижняя треугольная матрица. Можно показать, что метод Зейделя сходится, если A –симметричная положительно определенная матрица.
58
3.5.3.4. Метод релаксации
Чтобы ускорить итерационный процесс, запишем метод Зейделя в виде итерационной схемы (2), введя параметр ω, так что
(D +ωA− ) |
yk +1 − yk |
|
+ Ayk = f , k=0, 1,…, для всех y0 H. (13) |
|
|
|||||||||||||||||||||||||||||||||
|
|
|
|
|
|
|||||||||||||||||||||||||||||||||
|
|
|
|
|
ω |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
Сравнивая (13) с (2), видим, что |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
B=D+ω A–, τ=ω. |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||
Преобразуем уравнение (13) к расчетному виду. Учитывая, что |
|
|
|
|
|
|
|
|||||||||||||||||||||||||||||||
|
− |
|
y |
k +1 |
− y |
k |
|
|
|
|
|
− |
|
|
1 |
|
|
|
|
|
|
|
|
− |
|
D |
|
|
||||||||||
(D +ωA |
|
|
) |
|
|
|
|
|
+ A = A |
|
|
+ |
|
|
D y |
|
+ A |
− A |
|
− |
|
y |
|
= |
||||||||||||||
|
|
|
|
|
ω |
|
|
|
|
|
ω |
|
|
|
|
|
||||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
k |
|
|
|
|
|
|
|
|
|
|
|
k +1 |
|
|
|
|
|
ω |
k |
|
|||||||
|
|
|
|
|
− |
|
|
|
1 |
|
|
|
|
|
|
+ |
|
|
|
|
|
1 |
|
|
|
|
|
|
|
|
|
|
|
|
||||
|
= A |
|
|
+ |
|
|
D y |
k +1 |
+ |
A |
|
|
+ 1− |
|
|
D y |
k |
, |
|
|
|
|
|
|
|
|||||||||||||
|
|
|
ω |
|
|
ω |
|
|
|
|
|
|
|
|||||||||||||||||||||||||
имеем |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
− |
|
|
1 |
|
|
|
|
|
|
+ |
|
|
|
|
|
|
1 |
|
|
|
|
|
|
|
|
|
|
|
|
|||||||
A |
|
+ |
|
|
|
D |
yk +1 |
+ |
A |
|
+ |
1 |
− |
|
|
D yk = f . |
|
|
|
|
|
|
||||||||||||||||
|
ω |
|
|
|
|
|
|
|
|
|
||||||||||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
ω |
|
|
|
|
|
|
|
|
|
|
|
||||||||||
Отсюда находим
|
ω |
i−1 |
N |
|
yk(i+)1 = yk(i) + |
|
f (i) −∑aij yk( +j)1 −∑aij yk( j) |
||
|
||||
|
aii |
j=1 |
j=i |
|
,i =1,2,..., N,
т.е. при ω=1 получили формулу метода Зейделя.
Скорость сходимости метода релаксации зависит от параметра ω. Для сходимости метода требуется, чтобы 0<ω<2. При 0<ω<1 итерационный процесс (13) называют методом нижней релаксации, а при 1<ω<2 – методом верхней релаксации.
59
60
Модуль 4. Линейное программирование
4.6. Основы линейного программирования
Во многих областях знания часто возникают задачи оптимизационного типа, когда из многих вариантов возможного использования технических средств или других ресурсов желательно выбрать наилучший в некотором смысле. Такие задачи обычно выражают стремление достигнуть максимума или минимума некоторого, наперед обусловленного показателя, характеризующего в желательном смысле качество выполняемой работы. Например, задача может выражать стремление достигнуть максимума производительной способности, минимума себестоимости производства и т.д. Постановке и методам решения такого рода задач и посвящена настоящая глава.
4.6.1.Основы математического программирования
4.6.1.1.Оптимальное планирование и линейное программирование
Теоретической основой оптимизационных задач является относительно молодая область науки, называемая математическим программированием или оптимальным планированием. Она изучает методы отыскания таких количественных значений некоторых параметров (так называемых параметров управления), которые обеспечивают достижение оптимума выбранного показателя качества. Эти значения параметров управления составляют оптимальный план.
Разумеется, оптимизационная задача по существу возникает из-за ограниченности ресурсов, так что оптимальный план разыскивается среди тех планов, которые не нарушают ограничений, наложенных наличными ресурсами на параметры управления. Несмотря на наличие многих ограничивающих условий, почти всегда существует большое число допустимых вариантов решения задачи, т.е. различных возможных планов. Естественно, что конечная цель должна состоять не просто в построении какого-нибудь допустимого (совместного с ограничениями) плана, а в выборе из всех возможных вариантов оптимального (в некотором смысле) плана.
При ближайшем рассмотрении оказалось, что во многих оптимизационных задачах
показатель качества, который также называют целевой функцией или критерием оптимальности,
достаточно хорошо выражается линейной зависимостью от параметров управления. Кроме того, ограничения, наложенные запасами ресурсов, также выражаются линейными зависимостями от параметров управления. Оптимизационная задача при таких линейных условиях называется задачей линейного программирования, а методы решения этой задачи составляют предмет математической дисциплины, которую называют линейным программированием.
Принято выделять пять этапов постановки и решения задач математического программирования.
Первый этап – выбор параметров управления. Известно, что управление и планирование – сложные процессы, на их развитие оказывает влияние множество факторов, связанных между собой. Задача исследователя состоит в том, чтобы выделить основные, определяющие факторы и второстепенные. Конечно, в математическую модель не следует стремиться вводить все факторы, т.к. она станет слишком громоздкой, и исследование будет затруднительным, а подчас и невозможным. Важно выбрать только те параметры, которые существенны для конечной цели планирования, управления. Что касается прочих, малосущественных параметров, то достаточно иметь общее представление об их возможном влиянии.
Второй этап – выбор показателя качества (целевой функции). При выборе целевой функции следует исходить из конкретной обстановки, в которой рассматривается данная задача. То, что выгодно оптимизировать в одних условиях, может быть невыгодно в других. Например, обстановка может подсказать, что оптимизировать следует объем выполняемой работы (выпускаемой продукции, перевозок и др.). Но в других условиях может потребоваться
