Inversarea matricei nesingulare este subprogramul cel mai căutat căci nenumărate aplicații folosesc inversa matricei, de la ecuația de regresie până la rezolvarea de sisteme liniare.
Am tot căutat subprograme pentru inversarea unei matrice și dintre toate am găsit unul interesant în cartea PROBLEME DE ANALIZĂ NUMERICĂ REZOLVATE CU CALCULATORUL avându-i ca autori pe acad. Gheorghe MARINESCU, Irona RIZZOLI, Ileana POPESCU și Cristina ȘTEFAN, carte publicată îmn Editura Academiei RSR, în anul 1987 și care la paginile 171-173 prezintă metoda iterativă HOTELLING-BODEWING de inversare a unei matrice. Se construiește șirul lui NEWTON de torma X(n+1) = X(n) +(I-X(n)*A)*X(n), unde A, I, X, sunt matrice pătrate cu m linii și m coloane, iar cifrele cuprinse între paranteze acolo la indici arată numărul iterației. X(n+1) înseamnă matricea X la iterația n+1, iar X(n) înseamnă matricea X la iterația n. I este matricea unitate, A este matricea ce trebuie inversată. Subprogramul are ca prim pas evaluarea expresiei X1=X0 +(I+X0*A)*X0. În final matricea X1 va conține matricea inversă, iar X0 este matricea de start care se obține prin convenție dintr-o evaluare de forma: X0 = 1./(maxaij) * I. Variabila maxaij este de fapt elementul maxim dintre elementele matricei A. În subprogramul INVMAT, X0 se va referi la matricea de la iterația precedentă și X1 se va referi la matyricea de la iterația următoare.
Se observă că aici trebuie să fie apelate alte subprograme precum:
- subprogram de copiere matrice;
- subprogram de adunare matrice;
- subprogram de scădere matrice;
- subprogram de calcul normă matrice;
- subprogram de înmulțire matrice.
Față de programul dat în cartea pe care am citat-o am făcut o serie de modificări, căci m-am gândit ca utilizatorul să fie scutit în a da el acolo valoarea de start X0, eroarea, criteriu de stopare a procesului iterativ și matricea unitate.
SUBROUTINE START(A,N,X0)
MAXAIJ=A(1,1)
DO 10 I=1,N
DO 10 J=1,N
IF(MAXAIJ.GT.A(I,J)) GO TO 10
MAXAIJ=A(I,J)
10 CONTINUE
DO 20 I=1,N
DO 20 J=1,N
20 X0(I,J)=0.
ALFA=1./MAXAIJ
DO 30 I=1,N
30 X0(I,I)=ALFA
RETURN
END
Această subrutină construiește matricea de start X0 ca fiind matricea unitate înmulțită cu inversul elementului maxim din matricea A.
Pentru scăderea a două matrice pătrate având n linii și n coloane se foolosește subprogramul:
SUBROUTINE ADMAT(A,B,C,N)
MAXAIJ=A(1,1)
DO 10 I=1,N
DO 10 J=1,N
C(I,J)=A(I,J)+B(I,J)
10 CONTINUE
RETURN
END
Pentru scăderea a două matrice pătrate având n linii și n coloane se foolosește subprogramul:
SUBROUTINE ADMAT(A,B,C,N)
MAXAIJ=A(1,1)
DO 10 I=1,N
DO 10 J=1,N
C(I,J)=A(I,J)-B(I,J)
10 CONTINUE
RETURN
END
Pentru copierea matricei A în matricea C, ambele având n linii și n coloane se foolosește subprogramul:
SUBROUTINE SCMAT(A,C,N)
MAXAIJ=A(1,1)
DO 10 I=1,N
DO 10 J=1,N
C(I,J)=A(I,J)
10 CONTINUE
RETURN
END
Generarea unei matrice unitate de ordin n se face cu subprogramul:
SUBROUTINE UNIMT(N,UNIT)
DO 20 I=1,N
DO 10 J=1,N
UNIT(I,J)=0.
10 CONTINUE
UNIT(I,I)=1.
20 CONTINUE
RETURN
END
Norma ANOR a matricei A este calculată cu subprogramul:
SUBROUTINE NMMAT(A,N,ANOR)
ANOR=ABS(A1,1))
DO 20 I=1,N
DO 10 J=1,N
IF(ANOR.GT.ABS(A(I,J)) GO TO 10
ANOR=ABS(A(I,J))
10 CONTINUE
RETURN
END
Pentru înmulțirea a două matrice A și B pătrate având n linii și n coloane, rezultatul fiind matricea C tot cu n linii și n coloane se foolosește subprogramul:
SUBROUTINE PRMAT(A,C,N)
MAXAIJ=A(1,1)
DO 20 I=1,N
DO 20 J=1,N
CIJ=0.
DO 10 K=1,N
CIJ=CIJ+A(I,K)*B(K,J)
10 CONTINUE
C(I,J)=CIJ
20 CONTINUE
RETURN
END
Subprogramul pentru calculul matricei inverse AINV a matricei A cu n linii și n coloane este:
SUBROUTINE INVMAT(A,AINV,N,IK)
ITER=0
Q=.01
EPS=.0001
CALL START(A,N,X0)
10 CALL UNIMT(N,U)
CALL PRMAT(A,X0,X1,N)
CALL SCMAT(UNIT,X1,DIF,N)
ITER=ITER+1
CALL PRMAT(DIF,X0,X1,N)
CALL ADMAT(X0,X1,X1,N)
CALL NMMAT(X0,N,ANOR)
IF(Q**ITER*ANOR/(1.-Q).LE.EPS) GO TO 20
CALL CPMAT(X1,X0,N)
GO TO 10
20 CONTINUE
CALL CPMAT(X1,AINV,N)
RETURN
END
Față de varianta inițială, aici sunt date valorile care dau finețea rezultatului, fără a da flexibilitate utilizatorului să obțină el precizia dorită, din dorința de a nu avea mulți paramertrii în listă.
(27 decembrie 2017)
No comments:
Post a Comment