С точки зрения программирования матрица - это двумерный массив чисел. Так как бейсик "МК 85" не поддерживает массивы, организация матрицы сопряжена с некоторыми трудностями: приходится следить за тем, чтобы в участок адресов, отведенный под матрицу, не попадали другие переменные.
Матрица обрабатывается поэлементно. Обход элементов чаще всего организован как двойной вложенный цикл:
FOR I=0 TO A-1
FOR J=0 TO B-1
PRINT "M(";I;",";J;")=";M(I*B+J)
NEXT J
NEXT I
Здесь внутренний FOR изменяет номер столбца J, а внешний - номер строки I.
Транспонирование матрицы сводится к изменению порядка обхода:
FOR J=0 TO B-1
FOR I=0 TO A-1
PRINT "M(";J;",";I;")=";M(I*B+J)
NEXT I
NEXT J
Симметричные матрицы могут храниться в уплотненной форме - в "треугольном" массиве:
FOR I=0 TO A-1
FOR J=I TO A-1
PRINT "M(";I;",";J;")=";M(I*A+J)
NEXT J
NEXT I
Записаны только элементы, находящиеся на главной диагонали и справа от нее. При обращении к матрице должно выполняться отношение I≤J, что можно обеспечить введением новых индексов:
FOR I=0 TO A-1
FOR J=0 TO A-1
K=I:L=J:IF I>J;L=I:K=J
PRINT "M(";I;",";J;")=";M(K*A+L)
NEXT J
NEXT I
Если Вы работаете с разреженной матрицей, нет необходимости выделять память под нулевые элементы - матрицу можно сжать. Вид сжатия определяется характером разреженности. В общем случае, когда в расположении ненулевых элементов нет системы, применима следующая схема кодирования:
FOR I=0 TO A
K=INT N(I+I):L=100*FRAC N(I+I)
PRINT "M(";K;",";L;")=";M(I+I)
NEXT I
Здесь A - число ненулевых элементов матрицы, M(I+I) - значение I-го ненулевого элемента, N(I+I) - его индексы, упакованные в одну ячейку. Номер строки составляет целую часть, а номер столбца - дробную: N(I+I)=K+L/100. Формула записана для случая, когда число столбцов меньше ста. Уплотнение оправдано лишь тогда, когда число нулевых элементов больше, чем ненулевых.
При регулярной разреженности, в частности, для ленточных матриц, рекомендуется использовать частные алгоритмические правила уплотнения, так как они не требуют дополнительных затрат памяти на хранение индексов ненулевых элементов.
Представленные в разделе программы записаны для матриц произвольного размера. Работая с матрицей размера 2 или 3, нет смысла использовать такие длинные программы. Применяйте явное выражение детерминанта. Примеры:
1. Вычисление детерминанта матрицы 3-го порядка.
10 FOR A=1 TO 3:FOR B=1 TO 3:PRINT
A;B;:INPUT A(B+3*A):NEXT B
20 NEXT A:D=E*I*M-E*J*L+F*J*K-F*H*M+G*H*L-G*I*K:PRINT D
Адресация элементов матрицы: при вводе - индексная, при вычислении детерминанта - буквенная. Размещение матрицы в памяти:
жE
F Gц
зH
I Jч
иK
L Mш
2. Решение системы из двух линейных уравнений методом Крамера.
10 INPUT "A11",A,"A12",B,"A21",C,"A22",D,"B1",E,"B2",F
20 G=A*D-B*C:PRINT "X1";(E*D-B*F)/G,"X2";(A*F-E*C)/G
r'=Ar
Программа 16. Матрица вращения в трехмерном пространстве - типичный пример матричного преобразования. Здесь r - исходный радиус-вектор точки; r' - радиус-вектор точки после поворота вокруг начала координат на углы a, b, g; A - матрица вращения:
|
cosa cosb |
-cosa sinb cosg + sina sing |
cosa sinb sing + sina cosg |
|
sinb |
cosb cosg |
-cosb sing |
|
-sina cosb |
sina sinb cosg + cosa sing |
-sina sinb sing + cosa cosg |
10
MODE 4:INPUT "a",A,"b",B,"g",C:F=COS A:I=SIN B:J=COS C
20 K=SIN C:L=SIN A:G=L*K-F*I*J:H=F*I*K+L*J:M=L*I*J+F*K
30 N=F*J-L*I*K:D=COS B:F=F*D:J=J*D:K=-K*D:L=-L*D
40 INPUT "x",A,"y",B,"z",C:FOR D=0 TO 2:O=0:FOR E=0 TO 2
50 O=O+A(E)*F(E+D*3):NEXT E:PRINT CHR (88+D);O:NEXT D:GOTO 40
Размер: 212, O
Пример: a=30; b=45; g=60; x=1; y=0; z=0.
Ответ: x'=0,612372439; y'=0,70710678; z'=-0,353553389.
C=AB
Программа 17. Расчетная формула:
cij=е(k)aikbkj , k=1...r
10 INPUT "m",A,"r",B,"n",C:DEFM
A*B+B*C-18:FOR G=0 TO 1
20 FOR D=0 TO A(G)-1:FOR E=1 TO B(G):PRINT D+1;E;
30 INPUT G(E+D*B(G)+G*A*B):NEXT E:NEXT D:NEXT G:FOR D=1 TO A
40 FOR E=1 TO C:F=0:FOR G=0 TO B-1:F=F+H(G+D*B-B)*G(E+G*C+A*B)
50 NEXT G:PRINT D;E;F:NEXT E:NEXT D
Размер: 177, H(r*(m+n)). Ограничение размерности (n+m)r≤128, в частности, для квадратных матриц: n≤8.
Пример: m=2; r=3; n=1; A=((2, 7, 4), (8, 1, 3)); B=(6, 2, 4) T
Порядок ввода данных: 2; 3; 1; 2; 7; 4; 8; 1; 3; 6; 2; 4.
Ответ: C=(42, 64) T
Детерминант определяется формулой:
det A=е(k) (-1)kХ(l) ai(l,k),j(l,k); l=1...n
где k - номер перестановки,
n - размерность матрицы A,
aij - элементы матрицы.
Обратная матрица A-1 - это матрица, удовлетворяющая соотношению A-1A=E, где E - единичная матрица.
Программа 18. Модифицированный метод Гаусса - Жордана. Выбор ненулевого диагонального элемента обеспечивается перестановкой строк. Результат возвращается в ячейки исходного массива.
10 INPUT E:DEFM
E*E+E-20:FOR A=1 TO E:FOR B=1 TO E
20 PRINT A;B;:INPUT F(A*E+B):NEXT B:NEXT A
30 D=1:FOR A=1 TO E:IF F(A*E+A)¦0;F(A)=A:GOTO 80
40 FOR B=1 TO E:IF F(B*E+A)¦0 THEN 60
50 NEXT B:PRINT 0:END
60 D=-D:FOR C=1 TO E:F=F(A*E+C)
70 F(A*E+C)=F(B*E+C):F(B*E+C)=F:NEXT C:F(A)=B
80 F(A*E+A)=1/F(A*E+A):FOR C=1 TO E:IF C=A THEN 100
90 F(A*E+C)=F(A*E+A)*F(A*E+C)
100 NEXT C:FOR B=1 TO E:IF A=B THEN 130
110 F=F(B*E+A):F(B*E+A)=0:FOR C=1 TO E
120 F(B*E+C)=F(B*E+C)-F*F(A*E+C):NEXT C
130 NEXT B:D=D/F(A*E+A):NEXT A
140 PRINT D:FOR A=1 TO E:FOR B=1 TO E:F=F(A*E+F(B))
150 IF F(B)¦B;F(A*E+F(B))=F(A*E+B):F(A*E+B)=F
160 PRINT A;B;F:NEXT B:NEXT A
Размер: 453, F(n*(n+1)). Допустимая размерность n≤10.
Пример: n=3; A=
ж1 0 1ц
з1 1 0ч
и0 1 1ш
Порядок ввода данных: 3; 1; 0; 1; 1; 1; 0; 0; 1; 1.
Ответ: det A=2; A-1=
ж0.5 0.5 -0.5ц
з-0.5 0.5 0.5ч
и0.5 -0.5 0.5ш
Проверка с использованием свойства обратимости: (A-1)-1 = A.
Повторный запуск программы командой RUN 30 дает исходную матрицу с удвоенной вычислительной погрешностью. Для плохо обусловленных A погрешность может быть значительной.
P(l) = det(A-lE)
Программа 19. Метод Данилевского. Для невырожденной матрицы программа возвращает коэффициенты характеристического многочлена вида
P(l) = ln + е(i)cili-1 , i=1...n.
Вырожденной матрице соответствует несколько независимых многочленов. Тогда программа выводит значения k, по которым можно определить степени многочленов. Степень первого многочлена равна k1, второго - (k2-k1), третьего - (k3-k2) и т.д. Первые k1 выводимых на дисплей коэффициентов относятся к первому многочлену, вторые (k2-k1) - ко второму и т.д. В этом случае коэффициенты выводятся с противоположным знаком.
10 INPUT A:DEFM
A*A+A-20:FOR B=1 TO A:FOR C=1 TO A
20 PRINT B;C;:INPUT F(C+A*B):NEXT C:NEXT B:D=1
30 FOR B=D+1 TO A:E=0:FOR C=B TO A
40 IF E<ABS E(B+C*A);E=ABS E(B+C*A):F=C
50 NEXT C:IF E>0;FOR C=B-1 TO A:E=F(C+A*B):GOTO 70
60 FOR C=D TO B-1:F(C)=E(B+C*A):NEXT C:D=B:PRINT D-1;:GOTO 30
70 F(C+A*B)=F(C+F*A):F(C+F*A)=E:NEXT C:FOR C=D TO A:E=F(B+C*A)
80 F(B+C*A)=F(F+C*A):F(F+C*A)=E:NEXT C:FOR C=D TO A
90 F(C)=E(B+C*A):NEXT C:E=E(B+A*B):FOR C=B-1 TO A
100 F(C+B*A)=F(C+B*A)/E:FOR F=D TO A
110 IF F¦B;F(C+F*A)=F(C+F*A)-F(F)*F(C+B*A)
120 NEXT F:NEXT C:FOR C=D TO A:E=0:FOR F=B TO A
130 E=E+F(F+C*A)*F(F):NEXT F:IF C≤B;IF C¦D;E=E+E(C)
140 F(B+C*A)=E:NEXT C:NEXT B:FOR B=D TO A:F(B)=F(B*A+A):NEXT B
150 PRINT "":FOR B=1 TO A:PRINT B;F(B):NEXT B
Размер: 521, F(n*(n+1)). Ограничение размерности n≤9. Можно модифицировать программу для расчетов с n=10: убрать INPUT A:DEFM A*A+A-20 из строки 10 и строку вывода 150, а перед RUN дать команды A=10 и DEFM 90. После счета коэффициенты можно вызвать из ячеек G...P.
Пример: n=4; A=
ж1 1 4 3ц
з1 1 2 1ч
з0 0 4 3ч
и0 0 2 1ш
Порядок ввода данных:4; 1; 1; 4; 3; 1; 1; 2; 1; 0; 0; 4; 3; 0; 0; 2; 1.
Ответ: k1=2; c1=0; c2=2; c3=2;
c4=5;
что соответствует многочленам: P1(l)
= l2-2l;
P2(l) = l2-5l-2.
Собственные значения матрицы являются корнями характеристического многочлена:
P(li)=0, i=1...n
Для определения собственных значений несимметричных матриц пользуйтесь программой 19, затем 32. Симметричные матрицы приводятся к диагональным, что и реализовано в следующей программе.
Программа 20. Метод вращений с преградами. Матрица приводится к диагональной применением рекуррентной формулы
A0=A; Ak=UkTAk-1Uk;
где Uk - матрица плоского вращения, аннулирующая внедиагональный элемент, по модулю превышающий текущее значение барьера. От цикла к циклу величина барьера уменьшается.
В целях снижения вычислительной трудоемкости программа предусматривает задание погрешности - требуемого числа e верных знаков результата.
10 INPUT "e",B,"n",A:DEFM
A*(A+1)/2-12:E=0:FOR C=1 TO A
20 FOR D=1 TO C:PRINT C;D;:INPUT F:IF C¦D;E=E+2*F*F
30 N(C*(C-1)/2+D)=F:NEXT D:NEXT C:E=SQR E/A:B=E*10--B
40 $=CHR 1:FOR C=2 TO A:FOR D=1 TO C-1:F=N(C*(C-1)/2+D)
50 IF ABS F<E THEN 150
60 $="":G=N(C*(C+1)/2):H=N(D*(D+1)/2):I=(H-G)/2
70 J=1:IF I¦0;J=SGN I*F/SQR (F*F+I*I)
80 J=J/SQR (2*(1+SQR (1-J*J))):I=SQR (1-J*J)
90 FOR K=1 TO A:FOR L=0 TO 1
100 M(L)=C(L)*(C(L)-1)/2+K:IF C(L)<K;M(L)=K*(K-1)/2+C(L)
110 NEXT L:L=N(M)*I-N(N)*J:IF C¦K;N(N)=N(N)*I+N(M)*J
120 N(M)=L:NEXT K:L=2*F*I*J
130 N(C*(C-1)/2+D)=(G-H)*J*I+F*(I*I-J*J)
140 N(D*(D+1)/2)=H*I*I+G*J*J+L:N(C*(C+1)/2)=H*J*J+G*I*I-L
150 NEXT D:NEXT C:IF $="" THEN 40
160 IF E>B;E=E/10:GOTO 40
170 FOR C=1 TO A:PRINT N(C*(C+1)/2):NEXT C
Размер: 570, N(n*(n+1)/2). Ограничение размерности n≤13.
Пример: e=6; n=4; A=
ж2======ц
з8 3====ч
з4 6 1==ч
и5 2 9 0ш
Ответ: l1 =-9,97435311; l2 =18,5917142; l3 =-4,26123764; l4 =1,64387682.
Ax=b
К системам линейных уравнений с квадратной невырожденной матрицей применяют как точные методы, основанные на гауссовском приведении матрицы к треугольному виду, так и приближенные, использующие процесс пошаговой компенсации невязок. Последние хороши тем, что во-первых, не столь подвержены проклятью размерности, а во-вторых, позволяют учитывать разреженность матрицы коэффициентов. В этом разделе представлено по одному методу из каждого класса - метод Якоби и метод Зейделя, программы для которых оказались наиболее компактными. Системы с переопределенной матрицей и системы с вырожденной матрицей предварительно приводятся к нормальной системе (Используйте программу 17). Для переопределенной матрицы нормальная система получается по следующим формулам:
A'x=b'; где A'=ATA; b'=ATb.
Решение x соответствует минимуму скалярного произведения
(Ax-b)(Ax-b).
В случае вырожденной матрицы решается возмущенная система
(ATA+eE)x=ATb.
Оптимальная величина e зависит от конкретного вида матрицы. Полезно решить систему с разными e, чтобы убедиться в устойчивости результата.
Нормализованная система решается с помощью программы 21.
10 INPUT "n",A:DEFM
A*A+A-19:FOR B=1 TO A:FOR C=0 TO A
20 PRINT B;C;:INPUT G(C*A+B):NEXT C:NEXT B
30 FOR B=1 TO A-1:FOR C=B+1 TO A
40 F=G(B*A+B):G=G(B*A+C):IF F=0;IF G=0 THEN 80
50 E=SQR (F*F+G*G):D=G/E:E=F/E
60 FOR F=0 TO A:G=E*G(F*A+B)+D*G(F*A+C)
70 G(F*A+C)=E*G(F*A+C)-D*G(F*A+B):G(F*A+B)=G:NEXT F
80 NEXT C:NEXT B:FOR B=A TO 1 STEP -1
90 D=0:IF B<A;FOR C=B+1 TO A:D=D+G(C)*G(C*A+B):NEXT C
100 G(B)=(G(B)-D)/G(B*A+B):NEXT B
110 FOR B=1 TO A:PRINT B;G(B):NEXT B
Размер: 348,G(n*n+n). Допустимая размерность n≤10.
Пример: n=3; b=(31, 32, 33)T; A=
ж=0=30===1ц
з45 =0 -13ч
и12 13 ==8ш
Порядок ввода данных: 3; 31; 0; 30; 1; 32; 45; 0; -13; 33; 12; 13; 8.
Ответ: x=(1, 1, 1)T
Метод вращений применен также в программах: 81, 146.
Программа 151. Метод Гаусса с выбором максимального элемента.
10 INPUT A:D=A*(A+2):DEFM
D+A-15
20 FOR B=1 TO A:FOR C=1 TO A+1:PRINT B;C;
30 INPUT K(C*A+B):NEXT C:NEXT B
40 FOR B=1 TO A:K(D+B)=B:NEXT B
50 FOR B=1 TO A-1:E=B*A:F=0:FOR H=B+1 TO A
60 I=K(E+K(D+H)):IF ABS I>F;F=ABS I:G=H
70 NEXT H:GOSUB 200
80 J=K(D+G):K(D+G)=K(D+B):K(D+B)=J
90 FOR H=B+1 TO A:G=K(D+H):F=K(E+G)/I:K(E+G)=F
100 FOR K=B+1 TO A+1:C=K*A:K(C+G)=K(C+G)-F*K(C+J)
110 NEXT K:NEXT H:NEXT B
120 FOR B=A TO 1 STEP -1:J=K(D+B):I=K(D-A+J)
130 IF B<A;FOR H=B+1 TO A:I=I-K(H*A+J)*K(H):NEXT H
140 IF I=0 THEN 160
150 F=K(B*A+J):GOSUB 200:I=I/F
160 K(B)=I:NEXT B
170 FOR B=1 TO A:PRINT "x";B,K(B):NEXT B:END
200 IF F¦0;RETURN
210 PRINT "DET=0":END
Размер: 448, L(n(n+3)). Предельная размерность n≤9.
Примечание: Метод Гаусса представлен еще в программе 105.
Программа 22. Метод Зейделя для разреженных матриц. Для обеспечения сходимости следует провести замену переменных так, чтобы для всех строк i=1...n выполнялось:
|aii|>|aij|, j=1...n, j¦i.
Ввод данных. Сначала задается размерность системы n, затем допустимая абсолютная ошибка e. Система считывается построчно. Для i=1...n вводится диагональный элемент матрицы коэффициентов aii, затем свободный член bi, затем ненулевые внедиагональные элементы aij текущей строки. Для ввода aij необходимо сначала указать столбец j, а затем - само значение aij. Для перехода к следующей строке вводите j=0. Начальное приближение x0 задается командой VAC. При необходимости увеличитьточность пользуйтесь командами: B=1E-10, RUN P1, что вызовет выполнение добавочных итераций.
DEFM 3*n+2*r-20
VAC
В программной области P0:
10 INPUT "n",A,"e",B:D=A*3:FOR C=1
TO A:PRINT C
20 INPUT "a",C(C*3),"b",D(C*3)
30 INPUT "j",E:IF E=0;NEXT C:END
40 D=D+2:INPUT "a",E(D):D(D)=C+E/100:GOTO 30
В программной области P1:
10 $=CHR 1:D=A*3:FOR C=3 TO A*3
STEP 3:E=D(C)
20 IF INT F(D)=C/3;D=D+2:E=E-E(D)*E(300*FRAC D(D)):GOTO 20
30 E=E/C(C):IF ABS (E(C)-E)≥B;$=""
40 E(C)=E:NEXT C:IF $="" THEN 10
50 FOR C=1 TO A:PRINT C;E(C*3):NEXT C
Размер: 110+148, F(3*n+2*r). Ограничение размерности при поочередной загрузке областей P0 и P1: 3n=2r≤154, где n - размерность системы, r - количество ненулевых внедиагональных элементов матрицы A.
Пример: e =10-6
м10x1+5x4=15
н10x2+2x3=12
п10x3-3x1=
7==
о10x4+6x2=16
Порядок ввода данных: 4; 1E-6;
10; 15; 4; 5; 0;
10; 12; 3; 2; 0;
10; 7; 1; -3; 0;
10; 16; 2; 6; 0.
Ответ: x=(1, 1, 1, 1) T
Метод Зейделя применен в программах: 80, 101, 102, 145.