Реферат: То сходимость итерационного процесса может быть значительно увеличена, если использовать следующий алгоритм
§ 3. Методы решения систем линейных и нелинейных уравнений
accel1:=x-f/d1-der2(x)*f*f/(2*d1 3)
end.
Метод Эйткена - Стеффенсона. Если функция F(x) непрерывна и трижды непрерывно дифференцируема на отрезке [a, b], причем |F'(x)| < 1, то сходимость итерационного процесса может быть значительно увеличена, если использовать следующий алгоритм. На каждом шаге вычисления проводятся в три этапа: находим , затем на основе рассчитываем и далее по трем значениям хn ,, определяем приближение по следующей формуле:
. (1.21)
Итерации заканчиваются, когда знаменатель становится близким к нулю. Описанный алгоритм вычисления следующего приближения уже по имеющемуся реализован в про-
цедуре-функции accel2. Значение и тип параметров ясны из текста процедуры-функции accel2, которая, в свою очередь, обращается к функции, возвращающей значение F(x).
Function accel2(x:real) : real;
var f1,f2 : real;
begin
{ ** func вычисляет значение функции в точке x **}
f1:=func(x);
f2:=func(f1);
accel2:=(x*f2-f1*f1)/(f2-2*f1+x)
end; { *** accel2 ***} .
Сравнительные расчеты, выполненные с использованием процедур accel1, accel2 (по формулам 1.20 и 1.21) и newt, реализующей метод Ньютона (1.19), показали, что при использовании "ускоряющих" процедур для некоторых функций удается добиться ускорения сходимости в 3 - 5 раз.
^ § 3. МЕТОДЫ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ И НЕЛИНЕЙНЫХ УРАВНЕНИЙ
Система m линейных алгебраических уравнений с n неизвестными в общем виде может быть записана следующим образом:
(1.22)
или в матричном виде AX = B, где A - прямоугольная матрица размерности m´n, X - вектор n-го порядка, B - вектор m-го порядка. Решением системы (1.22) называется такая упорядоченная совокупность чисел которая обращает все уравнения системы в верные равенства. Две системы называются эквивалентными (равносильными), если множества их решений совпадают.
Система линейных уравнений называется совместной, если она имеет хотя бы одно решение, и несовместной - в противном случае. Совместная система называется определенной, если она имеет единственное решение, и неопределенной - в противном случае. Система является определенной, если rang A = rang B, где матрица B, полученная из матрицы A добавлением столбца свободных членов, называется расширенной.
Если матрица A - квадратная и det A 0, то она называется неособенной (невырожденной), при этом система уравнений, имеющая неособенную матрицу A, совместна и имеет единственное решение.
Eсли уравнения (1.22) являются нелинейными относительно неизвестного вектора Х, то соответствующая система, записанная в векторной форме
(1.23)
называется системой нелинейных уравнений. Она может быть также представлена в координатном виде:
1 < k < n.
Многообразие численных методов решения систем линейных алгебраических уравнений можно разделить на два класса: прямые (или точные) и итерационные (или приближенные) методы. В настоящем параграфе рассматриваются наиболее эффективные алгоритмы, реализующие ряд методов из обоих классов. Однако заметим, что нелинейные системы решают только итерационными методами, один из которых (метод Ньютона) рассматривается в настоящем параграфе.
^ 3.1. МЕТОД ИСКЛЮЧЕНИЯ ГАУССА ДЛЯ СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ
Из курса линейной алгебры [Крылов и др., 1972; Курош, 1962; Фадеев, Фадеева, 1963 и др.] известно, что решение системы линейных уравнений можно просто найти по правилу Крамера - через отношение определителей. Но этот способ не очень удобен для решения систем уравнений с числом неизвестных > 5, т.е. когда найти определитель сложно, а при числе неизвестных > 10 нахождение определителя с достаточно высокой степенью точности становится самостоятельной вычислительной задачей. В этих случаях применяют иные методы решения, среди которых самым распространенным является метод Гаусса.
Запишем систему линейных уравнений (1.22) в виде
(1.24)
Если матрица системы верхняя треугольная, т.е. ее элементы ниже главной диагонали равны нулю, то все хj можно найти последовательно, начиная с хn, по формуле
. (1.25)
При j > k и аjj 0 этот метод дает возможность найти решение системы.
Метод Гаусса для произвольной системы (1.22) основан на приведении матрицы А системы к верхней или нижней треугольной. Для этого вычитают из второго уравнения системы первое, умноженное на такое число, при котором а21= 0, т.е. коэффициент при х1 во второй строке должен быть равен нулю. Затем таким же образом исключают последовательно а31 , а41 , ..., аm1 . После завершения вычислений все элементы первого столбца, кроме а11, будут равны нулю.
Продолжая этот процесс, исключают из матрицы А все коэффициенты аij, лежащие ниже главной диагонали.
Построим общие формулы этого процесса. Пусть исключены коэффициенты из k - 1 столбца. Тогда ниже главной диагонали должны остаться уравнения с ненулевыми элементами:
Умножим k-ю строку на число
и вычтем ее из m-й строки. Первый ненулевой элемент этой строки обратится в нуль, а остальные изменятся по формулам
(1.26)
Произведя вычисления по формулам (1.26) при всех указанных индексах, исключим элементы k-го столбца. Такое исключение назовем циклом, а выполнение всех циклов назовем прямым ходом исключения.
После выполнения всех циклов получим верхнюю треугольную матрицу А системы (1.24), которая легко решается обратным ходом по формулам (1.25). Если при прямом ходе коэффициент аjj оказался равен нулю, то перестановкой строк перемещаем на главную диагональ ненулевой элемент и продолжаем расчет.
Если предположить, что аjj 0, то тогда можно применить подпрограмму gaus, решающую систему из n линейных уравнений. Это одна из подпрограмм, традиционно использующаяся в библиотеках стандартных программ для ЕС-ЭВМ. Ее перевод на язык PASCAL выполнен авторами.
^ Формальные параметры процедуры. Входные: n (тип integer) - порядок системы; а (тип real) - массив коэффициентов системы размером n´n; b (тип real) - массив-строка свободных членов. Выходные: х (тип real) - массив-строка, в который помещается решение системы.
procedure gaus (n:integer; a : mas; b : mas1;
var x : mas1);
type mst = array [1..n] of real;
mss = array [1..n] of mst;
var a1 : mss; b1 : mst;
i, j, m, k : integer; h : real;
begin
for i := 1 to n do
begin
b1[i] := b[i];
for j := 1 to n do
a1 [i,j] := a[i,j];
end;
for i := 1 to n-1 do
for j := i+1 to n do
begin
a1[j,i] :=- a1[j,i] / a1[i,i];
for k := i+1 to n do
a1[j,k] := a1 [j,k] + a1[j,i]*a1[i,k];
b1[j] := b1[j] + b1[i]*a1[j,i];
end;
x[n] := b1[n] / a1[n,n];
for i := n-1 downto 1 do
begin
h := b1[i];
for j := i+1 to n do
h := h - x[j]*a1[i,j];
x[i] := h / a1[i,i];
end;
end.
Если контроль за невырожденностью матрицы А необходим, то можно предложить воспользоваться другой процедурой gaus1.
В отличие от первой процедуры здесь в выходных параметрах есть переменная tol, возвращающая 0 при нормальном завершении работы процедуры, или 1, если на главной диагонали один из элементов равен 0, или 2, если матрица А размерностью больше, чем 50´50.
procedure gaus1 (n:integer; a : mas; b : mas1;
var x : mas1; var tol : integer);
type mst = array [1..50] of real;
mss = array [1..50] of mst;
var a1 : mss; b1 : mst;
i, j, m, k : integer; h : real;
begin
for i := 1 to n do
begin
b1[i] := b[i];
for j := 1 to n do a1 [i,j] := a[i,j];
end;
tol := 0;
if n > 50 then
begin
tol := 2;
exit;
end;
for i := 1 to n-1 do
for j := i+1 to n do
begin
if abs(a1[i,i])< 1.0e-10 then
begin
tol := 1;
exit;
end;
a1[j,i] :=- a1[j,i] / a1[i,i];
for k := i+1 to n do
a1[j,k] := a1 [j,k] + a1[j,i]*a1[i,k];
b1[j] := b1[j] + b1[i]*a1[j,i]
end;
if abs(a1[n,n]) < 1.0e-10 then
begin
tol := 1;
exit;
end;
x[n] := b1[n] / a1[n,n];
for i := n-1 downto 1 do
begin
if abs(a1[i,i]) < 1.0e-10 then
begin
tol := 1;
exit;
end;
h := b1[i];
for j := i+1 to n do h := h - x[j]*a1[i,j];
x[i] := h / a1[i,i];
end;
end.
По поводу приведенных алгоритмов можно сказать, что они не самые эффективные и пригодны для решения систем, имеющих невырожденные матрицы А и В с рангом не более 50, составленные из приблизительно одинаковых коэффициентов. Если между наибольшим и наименьшим из коэффициентов расстояние более чем 104, то эти алгоритмы применять не рекомендуется, так как погрешность вычислений будет такова, что может либо привести к вырожденной матрице А, либо дать в результате вектор Х с большой вычислительной погрешностью.
Для проверки работы процедур решим систему линейных уравнений методом Гаусса с точностью до 0.0001:
Коэффициенты, промежуточные результаты и окончательное решение приводятся в табл. 1.14. Подобные таблицы удобно использовать для ручного счета и проверки хода решения.
Таблица 1.14
Коэффициенты при неизвестных
Cвободные
Cтрочные
x1
x2
x3
x4
члены
суммы
0.68
0.21
0.11
-0.08
0.55
-0.13
0.84
0.15
-0.11
0.27
0.28
0.50
0.88
0.80
-0.06
0.12
2.15
0.44
0.83
1.16
2.85
-0.01
1.44
0.61
1
0.0735
-0.1618
0.1176
3.1618
4.1912
-0.1454
-0.18319
0.15590
0.30398
0.26220
-0.51290
-0.8247
0.0729
-0.1106
-0.22398
-0.48220
1.41290
-0.88015
-0.97847
0.94530
1
-2.09060
5.6719
1.54040
6.1217
-1.47643
-0.18697
4.79139
-0.99480
0.79920
1.17230
4.1136
-0.0095
1
-3.24410
-0.54110
-2.7851
-1.60130
1.07110
-0.5302
1
-0.66890
0.3311
2.8264
-0.3337
-2.7110
-0.6689
: РЕШЕНИЕ
^ 3.2. МЕТОД ГЛАВНЫХ ЭЛЕМЕНТОВ
На практике применяют множество вариантов вычислительных схем метода Гаусса. Например, при приведении матрицы к верхней треугольной форме выбирают наибольший элемент (в строке или столбце), уменьшая вычислительную погрешность за счет деления на не самый маленький элемент. Такая вычислительная схема называется методом Гаусса с выбором ведущего элемента. Если же выбирать при приведении матрицы самый большой (по модулю) элемент из всех оставшихся, то такая схема будет называться методом Гаусса с выбором главного элемента. Последняя схема относится к наиболее популярным. Главное ее отличие от метода Гаусса, рассмотренного в п. 3.1, состоит в том, что при приведении матрицы А к верхней (или нижней) треугольной форме ее строки и столбцы переставляют так, чтобы наибольший из всех оставшихся элементов матрицы стал ведущим, и на него выполняется деление. Если матрица хорошо обусловлена, то в методе Гаусса с выбором главного элемента погрешности округления невелики.
Описанный алгоритм, как и алгоритм метода Гаусса, представленный в п. 3.1, имеет множество программных реализаций и всегда есть в библиотеках стандартных программ. Один из относительно эффективных алгоритмов есть в БСП 1 на ЕС ЭВМ для транслятора FORTRAN-77. Он реализован в виде процедуры sinq. Перевод на язык PASCAL выполнен авторами.
^ Формальные параметры процедуры. Входные: n (тип integer) - целое положительное число, равное порядку исходной системы; aa (тип real) - массив из n´n действительных чисел, содержащий матрицу коэффициентов системы (aa[1] = а11; aa[2] = а21; aa[3] = a31; ...; aa[n]= an1; aa[n +1]= =a12; aa[n + 2] = a22; ...; aa[2*n] = аn2; ...; aa[n*n] = аnn); b (тип real) - массив из n действительных чисел, содержащий столбец свободных членов исходной системы (b[1] = =b1; b[2] = b2; ...; b[n] = bn). Выходные: b (тип real) - массив из n действительных чисел (он же был входным) при выходе из подпрограммы, содержащий решения исходной системы (b[1]= x1; b[2] = x2; ...; b[n] = xn); ks (тип integer) - признак правильности решения или код ошибки: если ks = =0, то в массиве b содержится решение исходной системы; если ks = 1, то исходная система не имеет решения (главный определитель ее равен нулю).
procedure sinq (var a: mas11; var b : mas1;
var ks : integer; n : integer);
var tol, biga, aa, save : real;
jj,jy,ij,it,j,i, i1, imax, k, iqs, ixj, ixjx : integer;
jjx, ny, ia, ib, io, i2, ix, jx, ky : integer;
label Cont10;
begin tol := 0.0;
ks := 0;
jj := -n;
for j := 1 to n do
begin jy := j + 1;
jj := jj + n + 1;
biga := 0.0;
it := jj - j;
for i := j to N do
begin
ij := it + i;
aa := abs (a[ij]);
if (abs(biga)-aa) < 0.0 then
begin
biga := a[ij];
imax := i; end;
end;
if (abs(biga)-tol) <= 0.0 then
begin
ks := 1;
exit;
end;
i1 := j + n*(j-2);
it := imax - j;
for k := j to n do
begin
inc (i1,n);
I2:= i1 + it;
save := a[i1];
a[i1] := a[i2];
a[i2] := save;
a[i1] := a[i1] / biga;
end;
save := b[imax];
b[imax] := b[j];
b[j] := save / biga;
if j=n then goto Cont10;
iqs := n*(j-1);
for ix := jy to n do
begin
ixj := iqs + ix;
it := j - ix;
for jx := jy to n do
begin
ixjx := n*(jx-1) + ix;
jjx := ixjx + it;
a[ixjx] := a[ixjx] - a[ixj]*a[jjx];
end;
b[ix] := b[ix] - b[j]*a[ixj];
end;
end;
Cont10: ny := n - 1;
it := n*n;
i := 1;
j :=1;
for ky := 1 to ny do
begin ia := it - ky;
ib := n - ky;
io := n;
for k := 1 to ky do
begin
b[ib] := b[ib] - a[ia]*b[io];
Dec (ia,n);
Dec (io);
end;
end;
end.
Для проверки процедуры и сравнения с работой процедур в п. 3.1 решалась система уравнений, приведенная в качестве примера в п. 3.1. Вычисления проводились с точностью до 10-5. Результаты работы программы приведены далее.
Исходная матрица А
Матрица В
0.680000 0.050000 -0.110000 0.080000
0.210000 -0.130000 0.270000 -0.800000
-0.110000 -0.840000 0.280000 0.060000
-0.080000 0.150000 -0.500000 -0.120000
2.150000
0.440000
-0.830000
1.160000
Преобразованная матрица А
Матрица В
1.000000 0.210000 -0.110000 -0.080000
0.000000 1.000000 -0.145441 0.155882
0.000000 0.000000 1.000000 0.258130
0.000000 0.000000 0.000000 1.000000
3.161765
0.579636
-2.851572
-0.669070
РЕШЕНИЕ СИСТЕМЫ
2.826351 -0.333733 -2.711759 -0.669070 .
^ 3.3. РЕШЕНИЕ СИСТЕМ НЕЛИНЕЙНЫХ УРАВНЕНИЙ МЕТОДОМ НЬЮТОНА
Очень распространенной является вычислительная задача нахождения некоторых или всех решений системы (1.23) из n нелинейных алгебраических или трансцендентных уравнений с n неизвестными.
Обозначим через Х вектор-столбец (х1, х2, ..., хn)T и запишем систему уравнений в виде формулы (1.23) F(Х) = 0, где F = (f1, f2, ..., fn)T.
Подобные системы уравнений могут возникать непосредственно, например при конструировании физических систем, или опосредованно. Так, к примеру, при решении задачи минимизации некоторой функции G(х) часто необходимо определить те точки, в которых градиент этой функции равен нулю. Полагая F = grad G, получаем нелинейную систему.
Основная идея метода Ньютона состоит в выделении из уравнений линейных частей, которые являются главными при малых приращениях аргументов. Это позволяет свести исходную задачу к решению последовательности линейных систем. При решении единственного уравнения (n=1) получаем рассмотренный в п. 2.3 метод Ньютона (1.19).
Метод Ньютона для n уравнений применим только тогда, когда могут быть вычислены все частные производные функций fi по переменным хi. Пусть J(х) обозначает матрицу Якоби и ее (i, j)-й элемент есть значение производной в точке xi. Как и в одномерном случае (n = 1), метод Ньютона начинается с произвольного Х, обозначенного Х0. Далее F линеаризуют для Х0, разлагая его в ряд Тейлора и учитывая лишь первые члены:
F(Х)= F(Х0) + J(Х 0)(Х - Х0) + ...
Линейное приближение к F около Х0в операторной форме задается
L(Х) = F(Х0 ) + J 0(Х - Х0), где J0 = J(Х0).
Чтобы найти следующее приближение Х к решению системы F(Х) = 0, решают уравнение F(Х0) +J0(Х1 - Х0) = 0. Естественно, решение можно также записать в форме X1 = =X0- (J0) -1F(X0), сильно напоминающей одномерную формулу метода Ньютона.
Однако для большинства систем из n уравнений с n неизвестными вычисление обратной матрицы (J0) -1не является необходимым, а наоборот, предпочтительнее решать линейную систему относительно поправки Х1 - Х0.
В общем случае, имея Хk, можно найти Хk+1 прибавлением к Хk поправки Хk+1 - Хk, полученной из решения линейной системы
F(Хk ) + J k(Хk+1 - Хk) = 0. (1.27)
Сходимость итерационного процесса (1.27) доказывается теоремой Дж.Форсайта и др. (1980), которую сформулируем неформально.
Пусть r - решение системы F(Х) = 0 такое, при котором J(n) не вырождена и вторые частные производные функции F непрерывны вблизи r. Тогда, если Х0 достаточно близко к r, то ньютоновы итерации сходятся. Более того, для еk = хk - r при k® ¥ отношение
ограничено. Cходимость процесса будет 2-го порядка.
Как и в одномерном случае, здесь основная проблема состоит в удачном выборе начального приближения, которое желательно было бы выбрать достаточно близко к предполагаемому корню, чтобы могла начаться быстрая сходимость.
На практике такое приближение достигается или очень большим везением (удачно выбран Х0), или мужеством и настойчивостью исследователя (выполняется очень много итераций до того, как процесс начнет быстро сходиться).
Еще одно замечание по поводу матрицы Якоби. Если вычисление производной функции уже в одномерном случае бывает более сложной задачей, чем отыскание значения самой функции, то для системы n уравнений вычисление F'(Х) во много раз более трудоемко, чем вычисление F(Х) (не забывайте, что это матрицы, а не просто функции !).
Попытки избежать перечисленные трудности превратились в отдельную вычислительную задачу. Прямое обобщение метода секущих оказалось неудовлетворительным, так как приближения к J(Х), получающиеся как n-мерный аналог метода секущей (метод хорд), часто оказываются вырожденными. Популярны так называемые квазиньютоновские методы, которые начинаются с очень трудных начальных итераций, но по мере приближения к решению точность их возрастает. Эти методы широко используют в задачах многомерной оптимизации.
Полная информация по рассмотренным методам имеется в работах Островского (1963) и Ортега, Рейнболда (1975), последнюю можно считать полным справочником по методам решения систем нелинейных уравнений.
В данной работе можно воспользоваться подпрограммой zeroin, подробно описанной в п. 2.6 с небольшими изменениями. Но здесь приводится другая процедура, взятая из БСП для БЭСМ (язык FORTRAN), переведенная на язык Pascal авторами. Эту подпрограмму используют традиционно для решения систем не выше 10-го порядка, она работает достаточно эффективно и быстро.
^ Формальные параметры процедуры. Входные: n (тип integer) - количество уравнений системы (совпадает с количеством неизвестных, причем n < 10); x (тип real) - массив из n действительных чисел, содержащий перед обращением newts начальное приближение решения; funcf - имя внешней подпрограммы, при помощи которой выполняют вычисления текущих значений функции F по заданным значениям, расположенных в элементах массива x, и размещение их в элементах массива y; funcg - имя внешней подпрограммы, предназначенной для вычисления по заданным значениям x из массива {x} элементoв матрицы dF/dх, размещенной в двухмерном массиве a размером [n´n]; eps (тип real) - значение из условия окончания итерационного процесса. Выходные: x (тип real) - массив из n действительных чисел (он же входной) содержит при выходе из newts приближенное значение решения; k (тип integer) - разрешенное количество итераций.
procedure newts (const n: integer;
var x : array of real; eps : real;
var xkit : integer);
var y, x1 : array [0..2] of real;
a : array [0..4] of real;
l, m : array [0..10] of integer;
j, i, nk : integer; s, d, xx,yy : real;
begin
xkit := 0;
repeat
for j := 1 to n do x1[j] := x[j];
xx := x1[1];
yy := x1[2];
FuncF (N, Xx,yy, Y);
FuncG (N, Xx,yy, a);
MinV (A, N, D, L, M);
for j := 1 to n do
begin
x[j] := x1[j];
for i := 1 to n do
begin
NK := J + N*(i-1);
x[j] := x[j] - A[NK]*Y[i];
end;
end;
Inc (xkit);
s := 1.0;
for j := 1 to N do s := s* sqr(X[j]-X1[j]);
s := sqrt(s);
until s > eps;
end.
Внимание: подпрограмма newts содержит обращение к подпрограмме minv, которая входит в БСП БЭСМ. Полный текст подпрограммы приводится ниже (перевод на язык PASCAL выполнен авторами).
procedure MinV (var a : array of real;
n : integer; var d : real;
var L, M : array of integer);
var nk,kk,k,iz,ij,i,ki,ji,jp, jk, ik, kj, jq : integer;
biga, hold : real;
jr : integer; dn : boolean;
begin d := 1.0;
nk := -n;
dn := false;
for k := 1 to n do
begin Inc(nk,n);
l[k] := k;
m[K] := k;
kk := nk + k;
biga := a[kk];
for j := k to n do
begin
iz := n*(j-1);
for i := k to n do
begin
ij := iz + i;
if (abs(biga)-abs(a[ij])) < 0.0 then
begin
biga := a[ij];
l[k] := i;
m[k] := j;
end;
end;
end;
j :=l[k];
if j > k then
begin
ki := k - n;
for i := 1 to n do
begin
Inc(ki,n);
hold := - f[ki];
ji := ki - k + j;
a[ki] := a[ji];
a[ji] := hold;
end;
end;
i := m [k];
if i > k then
begin
jp := N*(i-1);
for j := 1 to n do
begin
jk:= nk + j;
ji:= jp + j;
hold := - a[jk];
a[jk] := a[ji];
a[ji] := hold;
end;
end;
if abs(biga)< 1.0e-10 then
begin
d := 0;
exit;
end;
for i := 1 to n do
begin
if i<> k then
begin
ik := nk + i;
a[ik] := a[ik] / (-biga);
end;
end;
for i := 1 to n do
begin
ik := nk + i;
hold := a[ik];
ij := i - n;
for j := 1 to n do
begin
Inc(ij,n);
if i<>k then
begin
kj := ij - i + k;
a[ij] := hold*a[kl] + a[ij];
end;
end;
end;
kj := k - n;
for j := 1 to n do
begin
Inc (kj, n);
if j <> k then a[kj] := a[kj] / biga;
end;
d := d * biga;
a[kk] := 1.0 / biga;
end;
k := n-1;
repeat
i := l[k];
if i>k then
begin
jq := n*(k-1);
jr := n*(i-1);
for j := 1 to n do
begin
jk := jq + j;
hold := a[jk];
ji := jr + j;
a[jk] := -a[ji];
a[ji] := hold;
end;
end;
j := m[k];
if j>k then
begin ki := k - n;
for i := 1 to n do
begin
ki := ki + n;
hold := a[ki];
ji := ki - k + j;
a[ki] := -a[ji];
a[ji] := hold;
end;
end;
Dec (k);
until k=0;
end.
Для проверки работы процедур решим систему нелинейных уравнений с точностью до 0.001, используя метод Ньютона:
После отделения корней выяснилось, что система имеет решениe по х на интервале 0.4 < х < 0.5, а по y - 0.75 < у < <-0.73. Таким образом, за начальные приближения можно взять: х0 = 0.4; у0 = -0.75. Далее имеем
F(х,у) = sin(2х - у) - 0.4 - 1.2х ;
G(х,у) = 0.8х2 + 1.5у2 -1;
Fx ' = 2 cos (2х - у) - 1.2; Gх' = 1.6х;
Fy' = -cos(2х - у); Gy' = 3у.
Уточнение корней системы проведем методом Ньютона. Взяв за основу предлагаемые процедуры и доопределив процедуры для вычисления F(х), G(х) и F'(х), G'(х) как
Procedure FuncF (n : integer; xx1,xx2 : real;
var y : array of real);
begin
y[1] := xx1 - exp(-xx2);
y[2] := xx2 - exp (xx1);
end.
Procedure FuncG (n : integer; xx1,xx2: real;
var a : array of real);
begin
a[1] := 1.0;
a[2] := - exp(xx1);
a[3] := exp(-xx2);
a[4] := 1.0;
end
получим решение поставленной задачи:
Х1 = 0.273332; Y1 = 1.275009; итерация = 1;
Х2 = 0.273332; Y2 = 1.275009; итерация = 2.
^ 3.4. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ МЕТОДОМ КВАДРАТНОГО КОРНЯ
Метод квадратного корня основан на представлении матрицы А, составленной из коэффициентов системы в форме произведения треугольных матриц, что позволят свести решение заданной системы к последовательному решению двух систем с треугольными матрицами.
Метод квадратного корня применяется для решения системы линейных уравнений, коэффициенты которой образуют эрмитову симметрическую матрицу (эрмитова матрица совпадает с комплексно-сопряженной транспонированной A* = A).
Представим матрицу А в виде А = S*D S, где S - верхняя треугольная матрица с положительными элементами на главной диагонали; S*- транспонированная к матрице S; D- диагональная матрица, на диагонали которой находятся числа (+1) или (— 1).
Если матрицы S и D найдены, то заданная система AХ= = F может быть решена следующим путем:
AX = S* D S X = (S* D) S X = B Y = F, (1.28)
где S*D = В есть нижняя треугольная матрица и Y = S X - вспомогательный вектор.
Таким образом, решение системы АX = F равносильно решению двух треугольных систем ВY = F и SX =Y.
Пусть S = (sik) при i > k и sik 0, sii > 0; S* = (); D = =(dik), d =1; i k. Тогда из сравнения матриц А и S*DS получим
.
Ограничение в сумме получается из учета того факта, что в матрице S ниже главной диагонали элементы обращаются в нуль. Последнее равенство можно переписать несколько иначе, приняв для определенности k = min(k, l):
;
,
откуда окончательно получим формулы для вычисления элементов матриц S и D:
(1.29)
Единственным условием возможности определения sik является skk 0. Для построения матриц полагаем k = 1 и последовательно вычисляем все элементы первой строки s по формуле (1.29); затем полагаем k = 2 - и определяем элементы второй строки и т.д. Когда k = n, тогда найдены все элементы матриц S и D, а следовательно и S*. Затем последовательно выполняем вычисления (1.28):
S*Z = F; D Y = Z; S X = Y
обычным ходом по формулам
(1.30)
при i = 2, 3, ..., n.
Попутно заметим, что определитель матрицы А можно вычислить из выражения
. (1.31)
Этот метод экономичен, требует n3/3 арифметических действий и при больших n ( n > 50) вдвое быстрее метода Гаусса. Если за основу процедуры принять алгоритм П5.6 Дьяконова [1987], то тогда метод квадратного корня может быть реализован c помощью следующей процедуры:
procedure KVK (M, N : integer;
aa: array of real; bb : array of real;
var C: array of real);
type mas1=array [0..4] of real;
mas=array [0..4] of mas1;
var a : mas; j, k, i : integer; s, c1 : real;
begin
i := 0;
for j := 1 to m do
for k := 1 to n do
begin Inc (i); a[j,k] := aa [i]; end;
for j := 1 to N do
begin
for k := j to N do
begin
s := 0.0;
for i := 1 to M do s := s + a[i,j] * a[i,k];
c[k] := s;
end;
c1 := 0.0;
for i := 1 to M do с1 := c1 + a[i,j]*Bb[i];
for i := j to N do a[i,j] := c[i];
c[J] := c1;
end;
a[1,1] := sqrt (a[1,1]);
for j := 2 to N do a[1,j] := a[j,1] / a[1,1];
for i := 2 to n do
begin
s := 0.0;
for k := 1 to i-1 dо s := s + a[k,i]*a[k,i];
a[i,i] := sqrt (a[i,i] - S);
for j := i+1 to N do
begin
s := 0.0;
for k := 1 to i-1 do
S := S + A[K,I]*A[K,J];
A[I,J] := (A[J,I] - S) / A[I,I];
end;
end;
c[1] := c[1] / a[1,1];
for i := 2 to N do
begin s := 0.0;
for k := 1 to i-1do
s := s + a[k,i] * c[k];
c[i] := (c[i] - s) / a[i,i];
end;
c[n] := c[n] / a[n,n];
for i := n-1 downto 1 do
begin s := 0.0;
for k := i+1 to N do s := s + a[i,k] * c[k];
c[i] := (c[i] - s) / a[i,i];
end;
еnd.
^ Формальные параметры процедуры.Входные: m, N (тип integer) - m уравнений и n неизвестных определяют размер матрицы А. Для этой процедуры обязательно выполнение m > n, т.е. количество уравнений должно быть больше количества неизвестных (система переопределена), предельный разрешенный случай m = n; а - матрица, составленная из коэффициентов при неизвестных; В - массив, составленный из столбца свободных членов. Выходные: С - массив, в котором содержится решение системы.
Для проверки правильности работы процедуры решалась система 4´4 линейных уравнений методом квадратных корней с точностью 0.0001:
0.68X1 + 0.05X2 + 0.11X3 + 0.08X4 = 2.15 ;
0.05X1 + 0.13X2 + 0.27X3 + 0.80X4 = 0.44 ;
0.11X1 + 0.27X2 + 0.28X3 + 0.06X4 = 0.83 ;
0.08X1 + 0.80X2 + 0.06X3 + 0.12X4 = 1.16 ,
решение которой, полученное методом квадратных корней, будет следующим
Х 1 = 2.97; Х2 = 1.11; Х3 = 0.74; Х4 = -0.07.
^ 3.5. РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ МЕТОДОМ ИТЕРАЦИЙ
Итерационные методы решения систем линейных уравнений дают возможность вычислить решение системы как предел бесконечной последовательности промежуточных решений. Причем каждое последующее решение^ § 4. АЛГЕБРА МАТРИЦ
Матрицей называется прямоугольная таблица из чисел, содержащая некоторое количество m строк и n столбцов, при этом числа m и n называются порядками матрицы. Если m = n, матрица называется квадратной, а число m = n - ее порядком. Для записи матриц используются следующие обозначения:
, (1.33)
где числа aij называются элементами матрицы A. В случае, если m = n и определитель матрицы detA 0, матрица A называется невырожденной и для нее можно найти обратную матрицу. Обратной по отношению к данной матрице A называется матрица A-1, которая, будучи умноженной как справа, так и слева на A, дает единичную матрицу:
.
Матрица AT, полученная перестановкой строк со столбцами в матрице A, называется транспонированной. Квадратная матрица A называется симметрической, если AT = A, и ортогональной, если AT A = E.
Характеристическим уравнением матрицы A называется матричное уравнение ||A - E|| = 0, в котором i - собственные (или характеристические) числа (или значения) матрицы A. Собственным вектором, отвечающим собственному числу i, называется вектор
,
который удовлетворяет матричному уравнению Ax =i x. Когда говорят о вычислении собственных чисел, то различают полную и частичную проблему собственных чисел. В полной проблеме вычисляют все собственные числа и соответствующие им собственные векторы матриц, а под частичной проблемой обычно понимают задачу нахождения одного или нескольких собственных чисел и соответствующих им собственных векторов.
В настоящем параграфе рассматривается несколько эффективных алгоритмов, реализующих как операции над матрицами, так и задачи решения полной проблемы собственных значений.
^ Задача определения собственных значений и собственных векторов матриц имеет большое значение при решении очень широкого спектра задач. Методы решения указанной задачи по типу примененной вычислительной схемы можно разделить на прямые (точные) и итерационные (приближенные).
В прямых методах по некоторому правилу вычисляются коэффициенты матрицы с заранее известными свойствами, а затем собственные значения находятся как корни характеристического многочлена по какому-либо известному численному методу. После этого определяются собственные векторы, что считается не очень сложной задачей.
В итерационных методах коэффициенты характеристического уравнения не вычисляют, а составляют некоторые итерационные последовательности, позволяющие найти одно или несколько, а иногда и все собственные значения исходной матрицы А. Итерационные методы почти всегда более трудоемкие, однако они надежнее прямых методов, так как менее чувствительны к ошибкам округления.
К прямым можно отнести методы Крылова, Данилевского, Самуэльсона и др., а к итерационным - степенной и метод Якоби, или, как его еще называют, метод вращений. Последний метод считается наиболее эффективным из всех известных [Kpылов и др., 1972].
^ 4.1. ВЫЧИСЛЕНИЕ СОБСТВЕННЫХ ВЕКТОРОВ И СОБСТВЕННЫХ ЗНАЧЕНИЙ МАТРИЦ ПО МЕТОДУ КРЫЛОВА
Пусть А - некоторая квадратная матрица ||аij|| порядка n. Рассмотрим связанное с ней уравнение ||А -Е|| = 0, определитель которого есть алгебраический многочлен степени n от :
||А -Е|| = (-1)n (n - P1n -1Pn -1Pn . (1.34)
Корнями этого многочлена являются собственные значения (собственные числа) матрицы А. Эта система имеет ненулевое решение тогда и только тогда, когда n.
А.Н.Крылов1 предложил эффективный метод построения характеристического многочлена, основанный на применении минимального многочлена матрицы, аннулирующего заданный вектор.
Возьмем произвольный вектор с(0) 0, размерности n и составим последовательность линейно независимых векторов по правилу
с(1) = А с(0); c(2) = A c(1); ...; c(n) = A c(n -1).
Очевидно, что c(n) будет являться линейной комбинацией предыдущих линейно независимых векторов.
Расписав полученную на последнем шаге систему
(1.35)
где сi - соответствующие координаты вектора с(i), а qi - неопределенные коэффициенты, получим неоднородную систему из n алгебраических уравнений. Определитель этой системы
будет отличен от нуля тогда и только тогда, когда с(i) образуют систему линейно независимых векторов [Крылов и др., 1972]. Полученные значения qi будут равны рi - коэффициентам характеристического уравнения матрицы А (1.34).
Для определения qi систему (1.35) решают методом Гаусса (с обязательной проверкой). Если систему не удается привести в ходе решения к треугольному виду, то говорят о зависимости построенной системы векторов, и тогда можно записать только делитель характеристического многочлена матpицы. Сам многочлен решается одним из известных методов решения линейных алгебраических уравнений.
После того как найдены собственные значения i матрицы А, ищут собственные векторы матрицы А, что приводит к решению следующей однородной системы линейных алгебраических уравнений:
(A - E)х = 0. (1.36)
Так как i являются корнями минимального аннулирующего вектора с(0) многочлена (1.34), то решения системы (1.36) представляют в виде линейной комбинации независимых векторов с(i)
x(k) =i1 c(n -1) + i2c(n -2) + ... + in c(0),
где коэффициенты ij удовлетворяют условию Ax(n)=k x(n). Умножая x(n) на А и учитывая, что с(j)=Ас(j-1), получают
i(i1c(n-1)+i2c(n-2)+...+inc(0)) =
= i1 c(n) + i2c(n-1) + ... + in c(1).
Если же учесть, что A)c(n) º 0, то после приведения подобных (с учетом формул 1.35) последнее равенство можно переписать в виде
(qn i1- i in)c(0) + (qn-1 i1 + in - i in-1)c(1) +
+ (qn-2i1+ in-1- iin)c(2) + ... + (q1i1+ i2- ii1)c(n-1) = 0.
1 Библиотека стандартных программ. Поставляется вместе с транслятором языка.
1 Крылов А.Н. О численном решении уравнения высокой степени // ИАН ОМЕН. 1931. № 4. С. 491-539
еще рефераты
Еще работы по разное
Реферат по разное
Секция «Математика и информатика» Руководитель секции
18 Сентября 2013
Реферат по разное
Углеводы. 26
18 Сентября 2013
Реферат по разное
Методи кількісної оцінки економічного ризику
18 Сентября 2013
Реферат по разное
Учебной дисциплины «Численные методы» для направления 010200. 62 «Математика и компьютерные науки»
18 Сентября 2013