Системы линейных алгебраических уравнений

LA10

Программа LA10 осуществляет разложение матрицы A в произведение двух треугольных матриц L и U по компактной схеме Гаусса с выбором главного элемента по столбцу:
TA = LU,
где L - нижняя треугольная матрица с единичной диагональю, U - верхняя треугольная матрица, T - матрица перестановок.
Матрица TA есть матрица A с переставленными строками. Поскольку в матрице T порядка N в каждой строке и в каждом столбце содержится ровно одна единица и N-1 нулей, то в программе вместо матрицы сохраняется только вектор, содержащий информацию о перестановках строк матрицы A.
Процесс разложения матрицы включает в себя N-1 шаг. На каждом шаге с номером p вычисляются промежуточные величины a*ip:
p-1
a*ip = aip - ∑ likukj    (i = p, p+1, ..., N),
^k=1
где aip, lik, ukp - элементы матриц A, L и U соответственно. Далее находится такой номер строки m, при котором величина |a*ip| максимальна (для i от p до N). Если m≠p, то строки с номерами m и p переставляются. Если a*m,p = 0, то вычисления прекращаются.
Элементы p-го столбца матрицы L и p-ой строки матрицы U вычисляются по формулам:
lip = a*ip / app    (i = p+1, p+2, ..., N),
upp = a*pp,
N-1
upj = apj - ∑ lpkukj    (j = p+1, p+2, ..., N).
^k=1
На последнем, дополнительном шаге, вычисляется элемент uN,N:
N-1
uNN = aNN - ∑ lNkukN.
^k=1
На этом процесс разложения заканчивается.

Вызов программы

CALL LA10(A, AF, T, Det, Error)

Параметры программы

A - входной параметр;
AF, T, Det, Error - выходные параметры.

Real A(1:N, 1:N) - массив, содержащий элементы исходной матрицы;
Real AF(1:N, 1:N) - массив, содержащий элементы треугольных матриц, единичная диагональ нижней треугольной матрицы не хранится. При вызове массивы A и AF можно совмещать в памяти, в этом случае результат работы программы записывается в массив A;
Integer T(1:N) - массив, содержащий информацию о перестановках строк матрицы A;
Real Det - переменная, содержащая значение определителя матрицы A;
Integer Error - индикатор ошибки. Error=0, если ошибок нет и Error=65 или 66, если разложение матрицы не закончено и работа программы завершилась аварийно.

Вернуться к оглавлению    Скачать LA10

LA11

Программа LA11 решает систему линейных уравнений Ax=b с матрицей коэффициентов A, преобразованной при помощи программы LA10. Программа LA10 должна закончиться без ошибок (Error=0).
Решение системы Ax=b разбивается на решение двух треугольных систем Ly=b и Ux=y и выполняется по формулам:
i-1
yi = bi - ∑ lijyj    (i = 1, 2, ..., N)
j=1
для первой системы (прямой ход исключения неизвестных) и
N
xi = (yi - ∑ uijxj)/ujj    (i = N, N-1, ..., 1)
j=i+1
для второй системы (обратный ход исключения неизвестных).
Перед решением первой системы массив свободных членов b переставляется согласно массиву T, который содержит информацию о перестановках матрицы A при работе программы LA10.

Вызов программы

CALL LA11(AF, T, B, X)

Параметры программы

AF, T, B - входные параметры;
X - выходной параметр.

Real AF(1:N, 1:N) - массив, содержащий элементы треугольных матриц L и U;
Integer T(1:N) - массив, содержащий информацию о перестановках строк матрицы A;
Real B(1:N) - массив, содержащий вектор свободных членов;
Real X(1:N) - массив, содержащий вектор решения системы. При вызове массивы B и X можно совмещать в памяти, в этом случае результат работы программы записывается в массив B.

Пример

!Решение системы линейных уравнений
program Test_LA11
use ConvertCyr
use NML
implicit none
integer, parameter:: n=4
real:: A(n,n), AF(n,n), B(n), X(n), D
integer:: T(n), Er
data A/1.0, 1.1, 1.2, 1.4, 1.1, 1.1, 1.2, 1.3, 1.2, 1.2, 1.2, 1.3, 1.4, 1.3, 1.3, 1.3/
data B/11.1, 11.4, 12.1, 13.4/
integer i, j
!begin
  call LA10(A, AF, T, D, Er)
  if (Er>0) then
    print *, WinToDos('Матрица системы особенная. Программа остановлена.')
    stop
  end if
  call LA11(AF, T, B, X)
  print 10
  print 11, ((A(i,j), j=1,n), i, B(i), i=1,n)
  print 12
  print 13, (i, X(i), i=1,n)
  print 14, D
10 format(/'  The Equations Sistem:'/)
11 format('     |',<n>F5.1,' | X',I1,' |',F5.1,' |')
12 format(/'  Solution:'/)
13 format('       X',I1,' = ',F11.8)
14 format(/'  Determinant =',E15.7)
end

 The Equations Sistem:

    |  1.0  1.1  1.2  1.4 | X1 | 11.1 |
    |  1.1  1.1  1.2  1.3 | X2 | 11.4 |
    |  1.2  1.2  1.2  1.3 | X3 | 12.1 |
    |  1.4  1.3  1.3  1.3 | X4 | 13.4 |

 Solution:

      X1 =  3.99987936
      X2 =  3.00012612
      X3 =  2.00011635
      X4 =  0.99988753

 Determinant = -0.9999890E-04

Вернуться к оглавлению    Скачать LA11

LA12

Программа LA12 предназначена для итерационного решения системы линейных уравнений Ax=b с матрицей коэффициентов A, преобразованной при помощи программы LA10. Программа LA10 должна закончиться без ошибок (Error=0).
Уточненное решение x(p+1) получается из приближенного xp следующим образом. Сначала вычисляется вектор невязки
r(p) = b(p) - Ax(p).
Затем, решая систему
Ad(p) = r(p)
с использованием LU разложения матрицы A (прогр. LA10), получаем вектор поправки d(p), который прибавляется к приближенному решению
x(p+1) = x(p) + d(p).
Этот процесс повторяется до получения заданной точности или до достижения заданного количества итераций, начиная с нулевого приближения x(0)=0.
На каждом шаге итерационного процесса, начиная с p>0, проверяется выполнение неравенства
║d(p)║ ≤ ½║d(p-1)║, (1)
где ║d(p)║ = ∑|di(p)| - норма вектора поправок. Если неравенство (1) справедливо, то итерационный процесс продолжается.
Итерационный процесс заканчивается при появлении одной из следующих ситуаций:
   1. Начиная с некоторого p для заданного ε>0 компоненты вектора поправок удовлетворяют неравенству
|di(p)| ≤ ε•|xi(p)|    (i = 1, 2, ..., N),    (2)
в этом случае x=x(p) выдается в качестве решения.
   2. Начиная с некоторого p неравенство (1) перестает выполняться. В этом случае, если выполняется неравенство
║d(p)║ ≤ ε•║x(p)║, (3)
то заданной точности удовлетворяет только норма полученного решения и в качестве решения выдается x=x(p).
   3. Начиная с некоторого p неравенство (1) перестает выполняться и неравенство (3) не выполняется также. Это свидетельствует о том, что итерационный процесс не сходится.
   4. Неравенство (1) выполняется для всех p, но достигнуто максимальное количество итераций. Тогда в качестве решения выдается x=x(p+1).

Вызов программы

CALL LA12(A, AF, T, B, Eps, Max, X, Error)

Параметры программы

A, AF, B, T, Eps, - входные параметры;
Max - входной и выходной параметр;
X, Error - выходные параметры.

Real A(1:N, 1:N) - массив, содержащий элементы исходной матрицы;
Real AF(1:N, 1:N) - массив, содержащий элементы треугольных матриц L и U;
Integer T(1:N) - массив, содержащий информацию о перестановках строк матрицы A;
Real B(1:N) - массив, содержащий вектор свободных членов;
Real Eps - заданная точность решения;
Integer Max - на входе: максимальное количество итераций, которое разрешено выполнить; на выходе: количество итераций, выполненное программой;
Real X(1:N) - массив, содержащий вектор решения системы;
Integer Error - индикатор ошибки;
Error=0, если каждая компонента вектора x удовлетворяет заданной точности;
Error=1, если только норма x удовлетворяет заданной точности;
Error=2, если приближенное решение получено, но норма решения меньше заданной точности;
Error=3, если достигнуто максимальное количество итераций, но норма решения не удовлетворяет заданной точности;
Error=65, если итерационный процесс расходится.

Пример

!Решение системы линейных уравнений с итерационным уточнением
program Test_LA12
use ConvertCyr
use NML
implicit none
integer, parameter:: n=4
real, parameter:: Eps=1.E-7
real:: A(n,n), AF(n,n), B(n), X(n), D
integer:: T(n), Er, Max
data A /1.0, 1.1, 1.2, 1.4, 1.1, 1.1, 1.2, 1.3, 1.2, 1.2, 1.2, 1.3, 1.4, 1.3, 1.3, 1.3/
data B /11.1, 11.4, 12.1, 13.4/
integer i, j
!begin
  call LA10(A, AF, T, D, Er)
  if (Er>0) then
    print *, WinToDos(''Матрица системы особенная. Программа остановлена.')
    stop
  end if
  Max=20
  call LA12(A, AF, T, B, Eps, Max, X, Er)
  print 10
  print 11, ((A(i,j), j=1,n), i, B(i), i=1,n)
  print 111, Er
  print 112, Max
  print 12
  print 13, (i, X(i), i=1,n)
  print 14, D
10 format(/'  The Equations Sistem:'/)
11 format('     |',<n>F5.1,' | X',I1,' |',F5.1,' |')
12 format(/'  Solution:'/)
13 format('       X',I1,' = ',F11.8)
14 format(/'  Determinant =',E15.7)
111 format(/'  Error = ',I2)
112 format(/'  Iterations Quantity = ',I2)
end

 The Equations Sistem:

    |  1.0  1.1  1.2  1.4 | X1 | 11.1 |
    |  1.1  1.1  1.2  1.3 | X2 | 11.4 |
    |  1.2  1.2  1.2  1.3 | X3 | 12.1 |
    |  1.4  1.3  1.3  1.3 | X4 | 13.4 |

 Error =  0

 Iterations Quantity =  1

 Solution:

      X1 =  3.99989152
      X2 =  3.00011444
      X3 =  2.00010252
      X4 =  0.99989986

 Determinant = -0.9999890E-04

Вернуться к оглавлению    Скачать LA12

LA15

Программа осуществляет разложение симметричной матрицы A с ненулевыми главными минорами в произведение диагональной и двух треугольных матриц с единичными диагоналями согласно уравнению
TTAT = LTDL,
где L - верхняя треугольная матрица, LT - транспонированная к ней нижняя треугольная матрица, D - диагональная матрица, T - матрица перестановок.
С целью ограничения роста элементов в процессе преобразования в программе применяется выбор главного элемента на диагонали с одновременной перестановкой строк и столбцов. Такая перестановка сохраняет симметрию всех промежуточных матриц. Матрица перестановок T в программе не вычисляется, а вместо нее информация о перестановках сохраняется в линейном массиве.
Процесс преобразования осуществляется по формулам:
i-1
di = ∑ a2iidk    (i = 1, 2, ..., N),
k=1
i-1
aij = ∑ aikajk/dj    (i = 2, 3, ..., N;   j = 1, 2, ..., N-1).
k=1

Вызов программы

CALL LA15(A, D, Trp, T, Error)

Параметры программы

A - входной и выходной параметр;
Trp - входной параметр;
D, T и Error - выходные параметры.

Real A(1:N, 1:N) - массив, содержащий элементы исходной матрицы. На выходе: ниже главной диагонали располагается матрица L, главная диагональ и элементы выше главной диагонали остаются без изменений;
   Real D(1:N) - массив, который на выходе содержит элементы диагональной матрицы;
   Logical Trp - логическая переменная. Если при вызове программы Trp=true, то программа осуществляет поиск главного элемента с соответствующими перестановками строк и столбцов, если Trp=false, то поиск главного элемента не производится и программа работает без перестановок;
   Integer T(1:N) - массив, который на выходе содержит информацию о произведенных перестановках. Если при вызове программы Trp=false, то массив заполняется последовательно числами от 1 до N.
   Integer Error - индикатор ошибки. Error=0, если ошибок нет и равно 65, если программа завершилась аварийно.

Вернуться к оглавлению    Скачать LA15

LA16

Программа LA16 решает систему линейных уравнений Ax=b с симметричной матрицей коэффициентов A, разложенной в произведение двух треугольных и диагональной матриц при помощи программы LA15. Программа LA15 должна закончиться без ошибок (Error=0).
Решение системы Ax=b разбивается на решение двух треугольных и одной линейной систем LTz=b, Dy=z, Lx=y и осуществляется по формулам:
i-1
zi = bi - ∑ lijzj    (i = 1, 2, ..., N),
j=1

yi = zi/dii    (i = 1, 2, ..., N),
N
xi = yi - ∑ lijyj    (i = N, N-1, ..., 1).
^j=i+1
Перед решением первой системы переставляется массив b, а после решения последней системы переставляется массив x согласно массиву T, который содержит информацию о перестановках строк и столбцов матрицы A при работе программы LA15.

Вызов программы

CALL LA16(A, D, T, B, X)

Параметры программы

A, D, T, B - входные параметры;
X - выходной параметр.

Real A(1:N, 1:N) - массив, содержащий матрицу L (ниже главной диагонали);
Real D(1:N) - массив, содержащий главную диагональ матрицы D;
Integer T(1:N) - массив, содержащий информацию о перестановках строк и столбцов;
Real B(1:N) - массив, содержащий вектор свободных членов;
Real X(1:N) - массив, содержащий вектор решения системы. При вызове массивы B и X можно совмещать в памяти, в этом случае результат работы программы записывается в массив B.

Пример

!Решение системы линейных уравнений с симметричной матрицей
program TestLA16
use ConvertCyr
use NML
implicit none
integer, parameter:: n=9
real:: A(n,n), D(n), B(n), X(n)
integer:: T(n), Er
integer i, j
!begin
  do j=1,n-1
    A(j+1:n:2,j)=1.0
    A(j+2:n:2,j)=0.5
  end do
  forall(i=1:n-1, j=2:n, i<j) A(i,j)=A(j,i)
  forall(i=1:n) A(i,i)=2.0
  A(5,5)=2.2
  do i=1,n
    B(i)=sum(A(i,:))
  end do
  print 10
  print 11, ((A(i,j), j=1,n), i, B(i), i=1,n)
  call LA15(A, D, .true., T, Er)
  if (Er>0) then
    print *, WinToDos('Матрица системы особенная. Программа остановлена.')
    stop
  end if
  call LA16(A, D, T, B, X)
  print 12
  print 13, (i, X(i), i=1,n)
10 format(/'  The Equations Sistem:'/)
11 format('     |',<n>F5.1,' | X',I1,' |',F5.1,' |')
12 format(/'  Solution:'/)
13 format('       X',I1,' = ',F10.7)
end

The Equations Sistem:

   |  2.0  1.0  0.5  1.0  0.5  1.0  0.5  1.0  0.5 | X1 |  8.0 |
   |  1.0  2.0  1.0  0.5  1.0  0.5  1.0  0.5  1.0 | X2 |  8.5 |
   |  0.5  1.0  2.0  1.0  0.5  1.0  0.5  1.0  0.5 | X3 |  8.0 |
   |  1.0  0.5  1.0  2.0  1.0  0.5  1.0  0.5  1.0 | X4 |  8.5 |
   |  0.5  1.0  0.5  1.0  2.2  1.0  0.5  1.0  0.5 | X5 |  8.2 |
   |  1.0  0.5  1.0  0.5  1.0  2.0  1.0  0.5  1.0 | X6 |  8.5 |
   |  0.5  1.0  0.5  1.0  0.5  1.0  2.0  1.0  0.5 | X7 |  8.0 |
   |  1.0  0.5  1.0  0.5  1.0  0.5  1.0  2.0  1.0 | X8 |  8.5 |
   |  0.5  1.0  0.5  1.0  0.5  1.0  0.5  1.0  2.0 | X9 |  8.0 |

Solution:

     X1 =  1.0000013
     X2 =  0.9999984
     X3 =  1.0000014
     X4 =  0.9999987
     X5 =  1.0000011
     X6 =  0.9999989
     X7 =  1.0000004
     X8 =  1.0000001
     X9 =  1.0000004

Вернуться к оглавлению    Скачать LA16

LA17

Программа LA17 предназначена для итерационного решения системы линейных уравнений Ax=b с симметричной матрицей коэффициентов A, преобразованной при помощи программы LA15. Программа LA15 должна закончиться без ошибок (Error=0).
Описание алгоритма программы полностью совпадает с описанием алгоритма программы LA12.

Вызов программы

CALL LA17(A, D, T, B, Eps, Max, X, Error)

Параметры программы

A, D, B, T, Eps - входные параметры;
Max - входной и выходной параметр;
X, Error - выходные параметры.

Real A(1:N, 1:N) - массив, содержащий матрицу L (ниже главной диагонали);
Real D(1:N) - массив, содержащий главную диагональ матрицы D;
Integer T(1:N) - массив, содержащий информацию о перестановках строк и столбцов;
Real B(1:N) - массив, содержащий вектор свободных членов;
Real Eps - заданная точность решения;
Integer Max - на входе: максимальное количество итераций, которое разрешено выполнить; на выходе: количество итераций, выполненное программой;
Real X(1:N) - массив, содержащий вектор решения системы. При вызове массивы B и X можно совмещать в памяти, в этом случае результат работы программы записывается в массив B.
Integer Error - индикатор ошибки;
Error=0, если каждая компонента вектора x удовлетворяет заданной точности;
Error=1, если только норма x удовлетворяет заданной точности;
Error=2, если приближенное решение получено, но норма решения меньше заданной точности;
Error=3, если достигнуто максимальное количество итераций, но норма решения не удовлетворяет заданной точности;
Error=65, если итерационный процесс расходится.

Пример

!Решение системы линейных уравнений с симметричной матрицей
!и с итерационным уточнением решения
program TestLA17
use ConvertCyr
use NML
implicit none
integer, parameter:: n=9
real, parameter:: Eps=1.E-7
real:: A(n,n), D(n), B(n), X(n)
integer:: T(n), Er, Max
integer i, j
!begin
  do j=1,n-1
    A(j+1:n:2,j)=1.0
    A(j+2:n:2,j)=0.5
  end do
  forall(i=1:n-1, j=2:n, i<j) A(i,j)=A(j,i)
  forall(i=1:n) A(i,i)=2.0
  A(5,5)=2.2
  do i=1,n
    B(i)=sum(A(i,:))
  end do
  print 10
  print 11, ((A(i,j), j=1,n), i, B(i), i=1,n)
  call LA15(A, D, .true., T, Er)
  if (Er>0) then
    print *, WinToDos('Матрица системы особенная. Программа остановлена.')
    stop
  end if
  Max=20
  call LA17(A, D, T, B, Eps, Max, X, Er)
  print 111, Er
  print 112, Max
  print 12
  print 13, (i, X(i), i=1,n)
10 format(/'  The Equations Sistem:'/)
11 format('     |',<n>F5.1,' | X',I1,' |',F5.1,' |')
12 format(/'  Solution:'/)
13 format('       X',I1,' = ',F10.7)
111 format(/'  Error = ',I2)
112 format(/'  Iterations quantity = ',I2)
end

  The Equations Sistem:

     |  2.0  1.0  0.5  1.0  0.5  1.0  0.5  1.0  0.5 | X1 |  8.0 |
     |  1.0  2.0  1.0  0.5  1.0  0.5  1.0  0.5  1.0 | X2 |  8.5 |
     |  0.5  1.0  2.0  1.0  0.5  1.0  0.5  1.0  0.5 | X3 |  8.0 |
     |  1.0  0.5  1.0  2.0  1.0  0.5  1.0  0.5  1.0 | X4 |  8.5 |
     |  0.5  1.0  0.5  1.0  2.2  1.0  0.5  1.0  0.5 | X5 |  8.2 |
     |  1.0  0.5  1.0  0.5  1.0  2.0  1.0  0.5  1.0 | X6 |  8.5 |
     |  0.5  1.0  0.5  1.0  0.5  1.0  2.0  1.0  0.5 | X7 |  8.0 |
     |  1.0  0.5  1.0  0.5  1.0  0.5  1.0  2.0  1.0 | X8 |  8.5 |
     |  0.5  1.0  0.5  1.0  0.5  1.0  0.5  1.0  2.0 | X9 |  8.0 |

  Error =  1

  Iterations quantity =  1

  Solution:

       X1 =  0.9999999
       X2 =  0.9999999
       X3 =  1.0000000
       X4 =  1.0000000
       X5 =  0.9999998
       X6 =  1.0000001
       X7 =  1.0000000
       X8 =  0.9999999
       X9 =  1.0000001

Вернуться к оглавлению    Скачать LA17

LA30

Программа решает систему линейных уравнений с трехдиагональной матрицей коэффициентов методом прогонки. Исходная система может быть представлена в виде

где ai - элементы главной диагонали, bi - верхней диагонали, ci - нижней диагонали.
Решение системы проводится в два этапа. В процессе прямого хода прогонки диагональные элементы в каждом уравнении выражаются через внедиагональные и последовательно исключается нижняя диагональ исходной матрицы коэффициентов. Получаем следующую систему уравнений

│ x1 = p2 + q2x2 + r2xn
│ x2 = p3 + q3x3 + r3xn
│ . . . . . . . . (*)
│ xn-2 = pn-1 + qn-1xn-1 + rn-1xn
│ xn-1 = pn + qnxn + rnxn
где
p2 = d1/a1,   p3 = -(p2c2 - d2)/(a2 + q2c2),   . . . ,   pn = -(pn-1cn-1 - dn-1)/(an-1 + qn-1cn-1)
q2 = -b1/a1,   q3 = -(b2)/(a2 + q2c2),   . . . ,   qn = -(bn-1)/(an-1 + qn-1cn-1)
r2 = -c1/a1,   r3 = -(r2c2)/(a2 + q2c2),   . . . ,   rn = -(rn-1cn-1)/(an-1 + qn-1cn-1)
Вычислением коэффициентов pi, qi и ri по этим рекуррентным формулам заканчивается прямой ход прогонки. Для нахождения коэффициентов обратного хода в системе (*) начиная с xn-1 последовательно подставляем нижнее уравнение в верхнее, и тем самым исключаем верхнюю диагональ исходной матрицы коэффициентов. Система уравнений будет преобразована к виду

│ x1 = s1xn + t1
│ x2 = s2xn + t2
│ . . . . . . . . (**)
│ xn-2 = sn-2xn + tn-2
│ xn-1 = sn-1xn + tn-1
где
sn-1 = qn + rn,   sn-2 = qn-1sn-1 + rn-1,   . . . ,   s1 = q2s2 + r2
tn-1 = pn,   tn-2 = qn-1tn-1 + pn-1,   . . . ,   t1 = q2t2 + p2
На этом заканчивается обратный ход прогонки.
Значение переменной xn можно найти из системы двух уравнений

│ xn = dn/an - (bn/an)x1 - (cn/an)xn-1
│ x1 = s1xn + t1

Исключая из них x1, получаем:

xn = -(bnt1 + pncn - dn)/(an + bns1 + cn(qn + rn))

Подставляя найденное решение в систему (**), получаем окончательное решение исходной системы.

Вызов программы

CALL LA30(A, B, C, D, X, Error)

Параметры программы

A, B, C, D - входные параметры;
X, Error - выходные параметры.

Real A(1:N) - массив, содержащий элементы главной диагонали;
Real B(1:N) - массив, содержащий элементы верхней диагонали;
Real C(1:N) - массив, содержащий элементы нижней диагонали;
Real D(1:N) - массив, содержащий вектор свободных членов;
Real X(1:N) - массив, содержащий вектор решения;
Integer Error - индикатор ошибки
Error = 0, если решение системы получено;
Error = 1, если при вызове программы N<3;
Error = 65, если матрица системы вырожденная.

Пример

!Решение системы линейных уравнений с трехдиагональной матрицей
program TestLA30
use ConvertCyr
use NML_LA30
implicit none
integer, parameter:: n=40
real:: A(n), B(n), C(n), D(n), X(n)
integer:: Er
data A/n*2.0/, B/n*1.1/,C/n*1.0/
integer i
!begin
  forall(i=1:n) D(i)=A(i)+B(i)+C(i)
  call LA30(A, B, C, D, X, Er)
  if (Er>=65) then
    print *, WinToDos('Матрица системы особенная. Программа остановлена.')
    stop
  end if
  print 12
  print 13, (i, X(i), i=1,n,4)
  print 13, (n, X(n))
12 format(/'  Solution:'/)
13 format('       X',I2,' = ',F10.7)
end

 Solution:

       X 1 =  1.0000002
       X 5 =  0.9999990
       X 9 =  0.9999995
       X13 =  0.9999995
       X17 =  0.9999990
       X21 =  1.0000001
       X25 =  0.9999995
       X29 =  0.9999988
       X33 =  0.9999998
       X37 =  1.0000007
       X40 =  0.9999996

Вернуться к оглавлению    Скачать LA30

LA31

Программа решает систему линейных уравнений с пятидиагональной матрицей коэффициентов методом прогонки. Исходная система может быть представлена в виде

где ai - элементы главной диагонали, bi и ci - элементы верхних диагоналей, ei и di - элементы нижних диагоналей.
Для нахождения коэффициентов прямого хода диагональные элементы в каждом уравнении выражаются через внедиагональные и последовательно исключаются две нижние диагонали исходной матрицы коэффициентов. Получаем следующую систему уравнений

│ x1 = p3 + q3x2 + r3x3 + s3xn-1 + t3xn
│ x2 = p4 + q4x3 + r4x4 + s4xn-1 + t4xn
│ . . . . . . . . . . . . . . (*)
│ xn-3 = pn-1 + qn-1xn-2 + rn-1xn-1 + sn-1xn-1 + tn-1xn
│ xn-2 = pn + qnxn-1 + rnxn + snxn-1 + tnxn
где
p3 = g1/a1,   p4 = (g2 - D2p3)/A2,   p5 = (g3 - D3p4 - e3p3)/A3,   . . . ,   pn = (gn-2 - Dn-2pn-1 - en-2pn-2)/An-2
q3 = -b1/a1,   q4 = -(b2 + D2r3)/A2,   q5 = -(b3 + D3r4)/A3,   . . . ,   qn = -(bn-2 - Dn-2rn-1)/An-2
r3 = -c1/a1,   r4 = -c2/A2,   r5 = -c3/A3,   . . . ,   rn = -cn-2/An-2
s3 = -e1/a1,   s4 = -(D2s3)/A2,   s5 = -(D3s4 + e3s3)/A3,   . . . ,   sn = -(Dn-2sn-1 + en-2sn-2)/An-2
t3 = -d1/a1,   t4 = -(D2t3 + e2)/A2,   t5 = -(D3p4 + e3t3)/A3,   . . . ,   tn = -(Dn-2tn-1 + en-2tn-2)/An-2
  D2 = d2,   D3 = d3 + e3q3,   . . . ,   Dn-2 = dn-2 + en-2qn-2
  A2 = a2 + D2q3,   A3 = a3 + D3q4 + e3r3,   . . . ,   An-2 = an-2 + Dn-2qn-1 + en-2rn-2

Вычислением коэффициентов pi, qi, ri, si и ti по этим рекуррентным формулам заканчивается прямой ход прогонки. Для нахождения коэффициентов обратного хода в системе (*) начиная с xn-2 последовательно подставляем нижнее уравнение в два верхних и тем самым исключаем две верхние диагонали исходной матрицы коэффициентов. Система уравнений будет преобразована к виду

│ x1 = u1 + v1xn-1 + w1xn
│ x2 = u2 + v2xn-1 + w2xn
│ . . . . . . . . (**)
│ xn-1 = un-1 + vn-1xn-1 + wn-1xn
│ xn-2 = un-2 + vn-2xn-1 + wn-2xn
где
un-2 = pn,   un-3 = pn-1 + qn-1un-2,   . . . ,   u1 = p3 + q3u2 + r3u3
vn-2 = sn + qn,   vn-3 = sn-1 + qn-1vn-2 + rn-1,   . . . ,   v1 = s3 + q3v2 + r3v3
wn-2 = tn + rn,   wn-3 = tn-1 + qn-1wn-2,   . . . ,   w1 = t3 + q3w2 + r3w3

На этом заканчивается обратный ход прогонки.
Значение переменных xn-1 и xn можно найти из совместно решаемых первых двух уравнений системы (**) и уравнений

│ xn-1 = pn+1 + qn+1xn + rn+1x1 + sn+1xn-1 + tn+1xn
│ xn = pn+2 + qn+2x1 + rn+2x2 + sn+2xn-1 + tn+2xn
где
pn+1 = (gn-1 - Dn-1pn - en-1pn-1)/An-1,   qn+1 = -(bn-1 + Dn-1rn)/An-1,
rn+1 = -cn-1/An-1,   sn+1 = -(Dn-1sn + en-1sn-1)/An-1,   tn+1 = -(Dn-1tn + en-1tn-1)/An-1
pn+2 = (gn - Dnpn+1 - enpn)/An,   qn+2 = -(bn + Dnrn+1)/An,
rn+2 = -cn/An,   sn+2 = -(Dnsn+1 + ensn)/An,   tn+2 = -(Dntn+1 + entn)/An
  Dn-1 = dn-1 + en-1qn-1,   An-1 = an-1 + Dn-1qn + en-1rn-1
  Dn = dn + enqn,   An = an + Dnqn+1 + enrn

Исключая переменные x1 и x2, получаем

│ A11 + A12 = B1
│ A21 + A22 = B2
где
A11 = 1 - sn+1 - rn+1v1,   A12 = -(tn+1 + qn+1 + rn+1w1),   B1 = pn+1 + rn+1u1
A21 = -(sn+2 + qn+2v1 + rn+2v2),   A22 = 1 - tn+2 - qn+2w1 - rn+2w2,   B2 = pn+2 + qn+2u1 + rn+2u2

Решая эту систему по методу Крамера, имеем

│ xn-1 = (B1•A22 - B2•A12)/D
│ xn = (-B1•A21 + B2•A11)/D
где D = A11•A22 - A12•A21
Подставляя найденные значения в систему (**), получаем окончательное решение.

Вызов программы

CALL LA31(A, B, C, D, E, G, X, Error)

Параметры программы

A, B, C, D, E, G - входные параметры;
   X, Error - выходные параметры.

   Real A(1:N) - массив, содержащий элементы главной диагонали;
   Real B(1:N), C(1:N) - массивы, содержащие элементы верхних диагоналей;
   Real D(1:N), E(1:N) - массивы, содержащие элементы нижних диагоналей;
   Real G(1:N) - массив, содержащий вектор свободных членов;
   Real X(1:N) - массив, содержащий вектор решения;
   Integer Error - индикатор ошибки
Error = 0, если решение системы получено;
Error = 1 если при вызове программы N<5;
Error = 65, если матрица системы вырожденная.

Пример

!Решение системы линейных уравнений с пятидиагональной матрицей
program TestLA31
use ConvertCyr
use NML
implicit none
integer, parameter:: n=40
real:: A(n), B(n), C(n), D(n), E(n), G(n), X(n)
integer:: Er
data A/n*2.0/, B/n*1.2/, C/n*1.1/, D/n*0.9/, E/n*1.3/
integer i
!begin
  forall(i=1:n) G(i)=A(i)+B(i)+C(i)+D(i)+E(i)
  call LA31(A, B, C, D, E, G, X, Er)
  if (Er>=65) then
    print *, WinToDos('Матрица системы особенная. Программа остановлена.')
    stop
  end if
  print 12
  print 13, (i, X(i), i=1,n,4)
  print 13, (n, X(n))
12 format(/'  Solution:'/)
13 format('       X',I2,' = ',F10.7)
end

 Solution:

       X 1 =  0.9999997
       X 5 =  1.0000002
       X 9 =  1.0000001
       X13 =  1.0000000
       X17 =  0.9999999
       X21 =  0.9999993
       X25 =  0.9999995
       X29 =  0.9999995
       X33 =  0.9999994
       X37 =  0.9999998
       X40 =  1.0000000

Вернуться к оглавлению    Скачать LA31

LA44

Программа LA44 предназначена для итерационного решения системы линейных уравнений Ax=b с симметричной полжительно-определенной матрицей коэффициентов A методом сопряженных градиентов.
Итерационный процесс строится по следующей схеме. Исходя из начального приближения x(0) полагают s(1) = r(0) = Ax(0) - b и последовательно вычисляют рекуррентные соотношения
r(k) = r(k-1) - α(k)As(k), (1)
x(k) = x(k-1) - α(k)s(k), (2)
s(k+1) = r(k) + β(k)s(k), (3)
  k=1, 2, ...,
где
α(k) = (r(k-1), r(k-1))/(As(k), s(k))
β(k) = (r(k), r(k))/(r(k-1), r(k-1)).
Теоретически векторы r^(k) и x^(k) должны удовлетворять условиям
r(k) = Ax(k) - b, (4)
(r(i), r(k)) = 0  для  i ≠ k, (5)
r(n) = 0, x(n) = x∞, (6)
где x∞ - точное решение, а n - размерность матрицы системы. Однако из-за наличия ошибок округления при вычислениях реальный результат отличается от теории и соотношения (4), (5) и (6) не выполняются. Поэтому полученный после n итераций вектор x(n) может значительно отличаться от точного решения. Чтобы получить достаточно близкое приближение к точному решению в общем случае число итераций k должно превосходить размерность системы n. При этом рано или поздно наступит момент, когда невязка r(k) начнет резко уменьшаться и в конечном итоге точность решения будет определяться числом обусловленности матрицы A.
Вычисления прекращаются, как только выясняется, что продолжение итераций не улучшает решения. Критерием окончания процесса служит оценка малости величины ║r(k)║2 = (r(k), r(k)). Поскольку вектор невязок r(k) в алгоритме вычисляется рекуррентно через r(k-1) и s(k), то из-за наличия ошибок округления полученное значение r(k) несколько отличается от величины истинной невязки r*(k), полученной из вектора Ax(k) - b. По этой разнице можно судить о точности вычислений.
В рассматриваемом алгоритме вычисляется величина
σ = ║r*(k) - r(k)║2, (7)
где вектор r(k) определяется из выражения (1), а вектор r*(k) из выражения (4). Величина σ возрастает с ростом k. Как только ║r(k)║2 и σ становятся величинами одного порядка, это означает, что невязки r(k) и r*(k) достигли уровня ошибок округления.
Однако непосредственно использовать неравенство ║r(k)║2 < σ нельзя, т.к. может оказаться, что величина ║r(k)║2 все время чуть-чуть больше σ. Поэтому в алгоритме используется неравенство
║r(k)║2 < ωσ, (8)
в котором ω определяется выражением
ω = exp((k/n)2), (9)
где k - номер итерации, n - порядок системы уравнений.
Постоянно возрастающий коэффициент ω в выражении (8) обеспечивает окончание процесса даже при очень плохо обусловленной матрице A и, кроме того, позволяет завершить вычисления, когда ║r(k)║2 и σ быстро приближаются друг к другу.
Т.к. определение величины σ требует значительного времени из-за необходимости вычислять истинную невязку, то значение σ вычисляется через √-n шагов и остается постоянным до вычисления нового значения σ. Неравенство (8) проверяется на каждом шаге алгоритма.

Вызов программы

CALL LA44(Amv, B, X, Max, Error)

Параметры программы

Amv, B - входные параметры;
X, Max - входные и выходные параметры;
Error - выходной параметр.

Amv - процедура, которая вычисляет произведение Av для произвольного вектора v;
Real B(1:N) - массив, содержащий вектор свободных членов;
Real X(1:N) - массив, который на входе содержит вектор начального приближения, если оно известно. В этом случае параметр Max должен быть больше нуля. Если начальное приближение неизвестно, то параметр Max должен быть отрицательным. На выходе массив содержит вектор решения системы;
Max - параметр, который на входе задает максимально возможное количество итераций. Параметр должен быть больше нуля, если в массиве X содержится начальное приближение и должен быть меньше нуля, если оно не известно. Для определения максимально возможного количества итераций программа использует абсолютное значение параметра. На выходе параметр содержит число итераций, которое потребовалось программе для получения приближения к решению;
Integer Error - индикатор ошибки
Error = 0, если для получения решения потребовалось не более Max итераций;
Error = 1, если для получения решения Max итераций оказалось недостаточно;
Error = 65, если матрица А либо плохо обусловлена и ее наименьшее собственное значение имеет величину порядка ошибки округления, либо не является положительно-определенной. В обоих случаях решение получить нельзя.

Пример

!Решение системы линейных уравнений методом сопряженных градиентов

module KofArray
  integer, parameter:: n=400
  real:: A(n,n)
  !Примененный здесь метод размешения матрицы в памяти
  !допустим только для данной тестовой задачи
end module KofArray

program TestLA44
use ConvertCyr
use KofArray
use NML
implicit none
real:: B(n), X(n), Z(n)
integer:: Max, Er
!external Amv
integer i, j
!begin
  Max=-3*n
  A=0.0
  do i=1,n; A(i,i)=2.0; end do
  do i=1,n-1; A(i,i+1)=0.6; A(i+1,i)=0.5; end do
  do i=1,n-2; A(i,i+2)=0.4; A(i+2,i)=0.3; end do
  do i=1,n-3; A(i,i+3)=0.2; A(i+3,i)=0.1; end do
  do i=1,n-5; A(i,i+5)=0.1; A(i+5,i)=0.2; end do
  do i=1,n; B(i)=sum(A(i,:)); end do
  call LA44(Amv, B, X, Max, Er)
  if (Er>=65) then
    print *, WinToDos('Матрица системы особенная. Программа остановлена.')
    stop
  end if
  do i=1,n
    Z(i)=dot_product(A(i,:),X(:))
  end do
  print 10, Er
  print 11, Max
  print 12
  print 13, (i, X(i), i=1,9)
  print 14, (i, X(i), i=392,400)
  print 15
  print 16, (B(i), Z(i), abs(B(i)-Z(i)), i=1,20)
10 format(/'  Error = ',I2)
11 format(/'  Iterations quantity = ',I3)
12 format(/'  Solution:'/)
13 format('       X',I1,' = ',F10.7)
14 format('     X',I3,' = ',F10.7)
15 format(/10x,'B(i)',12x,'Z(i)',8x,'|B(i)-Z(i)|')
16 format(/3(4x,F12.7))
!end program TestLA44

contains

subroutine  Amv(V,AV)
!use KofArray
real, intent(in):: V(:)
real:: AV(:)
double precision Sum
integer i, j
  do j=1,n
    Sum=0.D0
    do i=1,n
      Sum=Sum+V(i)*dble(A(i,j))
	end do
	AV(j)=sngl(Sum)
  end do
  return
end subroutine Amv

end program TestLA44


  Error =  0

  Iterations quantity = 800

  Solution:

       X1 =  1.0939335
       X2 =  1.0340042
       X3 =  0.9903530
       X4 =  0.9437244
       X5 =  0.9596815
       X6 =  1.0182668
       X7 =  1.0068535
       X8 =  0.9986082
       X9 =  0.9996919
     X392 =  0.9961162
     X393 =  1.0012031
     X394 =  0.9981908
     X395 =  0.9916883
     X396 =  1.0423292
     X397 =  1.0497801
     X398 =  1.0016056
     X399 =  0.9585300
     X400 =  0.9075571

          B(i)            Z(i)        |B(i)-Z(i)|

       3.3000000       3.4949822       0.1949823
       3.8000000       3.8792984       0.0792985
       4.0999999       4.0795093       0.0204906
       4.1999998       4.0866756       0.1133242
       4.1999998       4.1053619       0.0946379
       4.3999996       4.4206786       0.0206790
       4.3999996       4.4113340       0.0113344
       4.3999996       4.4005747       0.0005751
       4.3999996       4.3919854       0.0080142
       4.3999996       4.3952918       0.0047078
       4.3999996       4.4017878       0.0017881
       4.3999996       4.4007087       0.0007091
       4.3999996       4.3998957       0.0001040
       4.3999996       4.4001207       0.0001211
       4.3999996       4.4001508       0.0001512
       4.3999996       4.3997464       0.0002532
       4.3999996       4.3999000       0.0000997
       4.3999996       4.4000793       0.0000796
       4.3999996       4.4001307       0.0001311
       4.3999996       4.4000778       0.0000782

Вернуться к оглавлению    Скачать LA44