'********************************************************* '* Calculate the coefficients of the characteristic poly-* '* nomial P(l)=A-l*I of a real tridiagonal matrix A(i,j) * '* ----------------------------------------------------- * '* Ref.: "Algebre, Algorithmes et programmes en Pascal * '* By Jean-Louis Jardrin, DUNOD Paris, 1988" [10].* '* ----------------------------------------------------- * '* SAMPLE RUN: * '* (Find the caracteristic polynomial of matrix: * '* 4 1 0 0 0 0 * '* 2 4 1 0 0 0 * '* A = 0 2 4 1 0 0 * '* 0 0 2 4 1 0 * '* 0 0 0 2 4 1 * '* 0 0 0 0 2 4 ) * '* * '* Input matrix size (max 100): 6 * '* Input 5 elements of subdiagonal: * '* Element 1: 2 * '* Element 2: 2 * '* Element 3: 2 * '* Element 4: 2 * '* Element 5: 2 * '* Input 6 elements of main diagonal: * '* Element 1: 4 * '* Element 2: 4 * '* Element 3: 4 * '* Element 4: 4 * '* Element 5: 4 * '* Element 6: 4 * '* Input 5 elements of supradiagonal: * '* Element 1: 1 * '* Element 2: 1 * '* Element 3: 1 * '* Element 4: 1 * '* Element 5: 1 * '* * '* Chracteristic polynomial P(l): * '* Coefficient for degree 6: 1.0000 * '* Coefficient for degree 5: -24.0000 * '* Coefficient for degree 4: 230.0000 * '* Coefficient for degree 3: -1120.0000 * '* Coefficient for degree 2: 2904.0000 * '* Coefficient for degree 1: -3776.0000 * '* Coefficient for degree 0: 1912.0000 * '* * '* (The characteristic polynomial of matrix A is: * '* P(l)=l^6-24*l^5+230*l^4-1120*l^3+2904*l^2-3376*l * '* +1912 ) * '* ----------------------------------------------------- * '* Subroutines used: * '* * '* 1000 ReadData: Read elements of input matrix * '* 2000 PCTR : Calculate the coefficients of the * '* characteristic polynomial P(l)=A-l*I * '* of a real square tridiagonal matrix. * '* * '* Basic version by J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '********************************************************* 'Program CARPOL DEFDBL A-H, O-Z DEFINT I-N DIM B(25), C(25), D(25), E(25), P(25) CLS PRINT INPUT " Input size of matrix: ", n GOSUB 1000 'call Read_data(n,C,D,E) ' prepare B vector FOR k = 1 TO n - 1 B(k) = C(k) * E(k) NEXT k ' call routine to calculate coefficients GOSUB 2000 'call PCTR(n,B,D,P) ' print results PRINT PRINT " Characteristic polynomial P(l):" FOR k = 1 TO n + 1 PRINT " Coefficient for degree "; n - k + 1; ": "; PRINT USING "#######.####"; P(k) NEXT k PRINT END 'of main program 1000 'Subroutine ReadData 'integer n, i PRINT PRINT " Input "; n - 1; " elements of subdiagonal:" FOR i = 1 TO n - 1 PRINT " Element "; i; ": "; : INPUT "", C(i) NEXT PRINT " Input "; n; " elements of main diagonal:" FOR i = 1 TO n PRINT " Element "; i; ": "; : INPUT "", D(i) NEXT PRINT " Input "; n - 1; " elements of supradiagonal:" FOR i = 1 TO n - 1 PRINT " Element "; i; ": "; : INPUT "", E(i) NEXT RETURN '***************************************************** '* The subroutine PCTR calculates the coefficients * '* 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) * '* B: vector of codiagonal terms products * '* (see main program). * '* D: vector of main diagonal terms * '* OUTPUT: * '* P: vector of coefficients of P(l) * '* * '***************************************************** 'Subroutine PCTR(n,B,D,P) 'integer n,i,k, ialloc 2000 DIM T(25) P(2) = D(1): T(1) = 1#: P(1) = -T(1) FOR k = 2 TO n P(k + 1) = D(k) * P(k) - B(k - 1) * T(k - 1) FOR i = k TO 3 STEP -1 T(i) = P(i) P(i) = -T(i) + D(k) * P(i - 1) - B(k - 1) * T(i - 2) NEXT i T(2) = P(2): P(2) = -T(2) + D(k) * P(1) T(1) = P(1): P(1) = -T(1) NEXT k RETURN 'end of file carpol.bas