/******************************************************** * Calculate the coefficients of the characteristic poly-* * nomial P(l)=A-l*I of a real tridiagonal matrix A(i,j) * * ----------------------------------------------------- * * Ref.: "Algèbre, Algorithmes et programmes en Pascal * * By Jean-Louis Jardrin, DUNOD Paris, 1988" * * [BIBLI 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 ) * * * * C++ version by J-P Moreau, Paris * * (with dynamically allocated vectors). * * (www.jpmoreau.fr) * ********************************************************/ #include #include int k,n; REAL *B,*C,*D,*E,*P; void *vmblock = NULL; /**************************************************** * The procedure 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) * * * * The pointers to a real vector must be declared * * and allocated in the calling program by using the * * functions of module VMBLOCK.CPP. * ****************************************************/ void PCTR(int n, REAL *B, REAL *D, REAL *P) { int i,k; REAL *T; void *vmblock1 = NULL; vmblock1 = vminit(); T = (REAL *) vmalloc(vmblock1, VEKTOR, n, 0); P[2]=D[1]; T[1]=1.0; P[1]=-T[1]; for (k=2; k<=n; k++) { P[k+1]=D[k]*P[k]-B[k-1]*T[k-1]; for (i=k; i>2; i--) { T[i]=P[i]; P[i]=-T[i]+D[k]*P[i-1]-B[k-1]*T[i-2]; } T[2]=P[2]; P[2]=-T[2]+D[k]*P[1]; T[1]=P[1]; P[1]=-T[1]; } free(vmblock1); } void Read_data(int n) { int i; printf("\n Input %d elements of subdiagonal:\n",n-1); for (i=1; i