'********************************************************* '* Calculate the coefficients of the characteristic poly-* '* nomial P(l)=A-l*I of a square complex matrix A(i,j) * '* ----------------------------------------------------- * '* Ref.: "Algèbre, Algorithmes et programmes en Pascal * '* by Jean-Louis Jardrin, DUNOD Paris, 1988" [10].* '* ----------------------------------------------------- * '* SAMPLE RUN: * '* (Find the characteristic polynomial of matrix: * '* 1 -1 3 * '* A = -1 i -1-2i * '* i 1 -2+i ) * '* * '* Input matrix size (max 20): 3 * '* * '* Line 1 * '* Column 1 * '* Real part: 1 * '* Imag. part: 0 * '* Column 2 * '* Real part: -1 * '* Imag. part: 0 * '* Column 3 * '* Real part: 3 * '* Imag. part: 0 * '* Line 2 * '* Column 1 * '* Real part: -1 * '* Imag. part: 0 * '* Column 2 * '* Real part: 0 * '* Imag. part: 1 * '* Column 3 * '* Real part: -1 * '* Imag. part: -2 * '* Line 3 * '* Column 1 * '* Real part: 0 * '* Imag. part: 1 * '* Column 2 * '* Real part: 1 * '* Imag. part: 0 * '* Column 3 * '* Real part: -2 * '* Imag. part: 1 * '* * '* Characteristic polynomial P(l): * '* Coefficient for degree 3: * '* -1.000000 + 0.000000 i * '* Coefficient for degree 2: * '* -1.000000 + 2.000000 i * '* Coefficient for degree 1: * '* 3.000000 + 1.000000 i * '* Coefficient for degree 0: * '* 0.000000 + 0.000000 i * '* * '* (The characteristic polynomial of matrix A is: * '* P(l)=-l^3+(-1+2i)l^2+(3+i)l * '* * '* Basic version by J-P Moreau, Paris * '* (www.jpmoreau.fr) * '********************************************************* 'Program CARPOL1 'Note: a complex number is represented by a couple of two real numbers (Ex.:z0 by z01 and z02). DEFDBL A-H, O-Z DEFINT I-N NMAX = 20 NMAXP = 21 DIM A1(NMAX, NMAX), A2(NMAX, NMAX) DIM P1(NMAXP), P2(NMAXP) GOSUB 1000 'Read data GOSUB 2000 'Call routine to calculate coefficients ' print results F\$="###.######" PRINT PRINT " Characteristic polynomial P(l):" FOR k = 1 TO n + 1 PRINT " Coefficient for degree "; n - k + 1; ": " PRINT USING F\$; P1(k); IF P2(k) >= 0 THEN PRINT " +"; ELSE PRINT " -"; END IF PRINT USING F\$; ABS(P2(k)); : PRINT " i" NEXT END 1000 'Read data CLS : PRINT PRINT " Input matrix size (max "; NMAX; "): "; INPUT "", n FOR i = 1 TO n PRINT : PRINT " Line "; i FOR j = 1 TO n PRINT " Column "; j PRINT " real part: "; : INPUT "", A1(i, j) PRINT " imag. part: "; : INPUT "", A2(i, j) NEXT j NEXT i RETURN 1500 'calculate complex trace of matrix C 't01,t02:Complex c01 = 0#: c02 = 0# FOR i = 1 TO n t01 = c01: t02 = c02 'CADD(t0,C(i,i),c0) c01 = t01 + C1(i, i) c02 = t02 + C2(i, i) NEXT RETURN '****************************************************** '* The procedure PCCS calculates the complex coeffi- * '* cients of the characteristic polynomial P(l)=A-l*I * '* of a real square tridiagonal matrix A(i,j). * '* * '* Note: the roots of P(l) are the eigenvalues of * '* matrix A(i,j). * '* -------------------------------------------------- * '* INPUTS: * '* n: zize of matrix (integer) * '* A: complex square matrix (n,n) * '* (see main program). * '* OUTPUT: * '* P: vector of complex coefficients of P(l) * '* * '****************************************************** 2000 'Procedure PCCS ' c0,c1,s,t0,z0,z1: Complex DIM B(NMAX, NMAX), C(NMAX, NMAX) z01 = 0#: z02 = 0# z11 = 1#: z12 = 0# IF n / 2 <> INT(n / 2) THEN 'CPRO(-1.0,z1,P^[1]) P1(1) = -1# * z11: P2(1) = -1# * z12 ELSE P1(1) = z11: P2(1) = z12 END IF FOR l = 1 TO n IF l = 1 THEN FOR i = 1 TO n FOR j = 1 TO n C1(i, j) = A1(i, j) C2(i, j) = A2(i, j) NEXT j NEXT i ELSE FOR i = 1 TO n FOR j = 1 TO n s1 = z01: s2 = z02 FOR k = 1 TO n c01 = s1: c02 = s2 'CMUL(B(i,k),A(k,j),c1) c11 = B1(i, k) * A1(k, j) - B2(i, k) * A2(k, j) c12 = B1(i, k) * A2(k, j) + B2(i, k) * A1(k, j) 'CADD(c0,c1,s) s1 = c01 + c11: s2 = c02 + c12 NEXT k C1(i, j) = s1: C2(i, j) = s2 NEXT j NEXT i END IF GOSUB 1500 'TRC(n,C,c0) 'CPRO(1.0/l,c0,t0) t01 = 1# / l * c01: t02 = 1# / l * c02 'CMUL(t0,P(1),c0) c01 = t01 * P1(1) - t02 * P2(1) c02 = t01 * P2(1) + t02 * P1(1) 'CPRO(-1.0,c0,P(l+1)) P1(l + 1) = -1# * c01: P2(l + 1) = -1# * c02 IF l < n THEN FOR i = 1 TO n FOR j = 1 TO n IF j = i THEN 'CDIF(C(i,j),t0,B(i,j)) B1(i, j) = C1(i, j) - t01 B2(i, j) = C2(i, j) - t02 ELSE B1(i, j) = C1(i, j): B2(i, j) = C2(i, j) END IF NEXT j NEXT i END IF NEXT l RETURN 'end of file carpol1.bas