(HP71B) Eigenvalues of Symmetric Real Matrix

[Albert Chan]

HP71B, version 2

1. Suggested by Werner, assume user inputed symmetric square matrix. Remove code: A=(A+A^T)/2

2. Bij = tan(θ/2) replace with tan(θ)/2, (|θ|≤pi/4), which use 1 sqrt, instead of 2
This has a slight effect to convergence rate, but likely not noticeable.

For OP example, it converged slightly better.

Cas> float(eigenvals([[0,1,2],[1,2,3],[2,3,4]])) /* for reference */
[0., 6.87298334621, −0.872983346207]

This is 6th iteration, with old Bij = tan(θ/2)

Code:

FNORM(B)= 7.60886994417E-10
Q=
 .859892597285       -.40824829046         .306460526788
 .193822655517        .816496580921        .543843830086
-.472247286237       -.40824829046         .7812271334
L=
-.872983346221       -1.69E-13            -3.663E-12
-1.69203950392E-13    2.E-24               3.73782958716E-12
 1.66E-12             3.73E-12             6.87298334595

This is 6th iteration, with new Bij = tan(θ)/2

Code:

FNORM(B)= 8.7618796429E-11
Q=
 .408248290464        .859892597282        .306460526793
-.816496580928        .193822655521        .543843830107
 .408248290464       -.472247286237        .78122713341
L=
 0                    0                    0
-1.E-15              -.872983346218        1.786E-12
 0                    3.53E-12             6.87298334624

3. Q correction:

I = (I+B)^-1 (I+B)
2 (I+B)^-1 - I = (I+B)^-1 - (I+B)^-1 B = (I+B)^-1 (I-B)

HP71-MathPac-suppl_manual.pdf Wrote:MAT X=INV(A)*B directly computes the solution of the matrix equation A X=B instead of
inverting A then multiplying by B

RHS cost the same as matrix inverse. (although in this case, LHS is cheap to compute too)

4. Combine "MAT L=TRN(Q) @ MAT L=L*A" to MAT L=TRN(Q)*A

10 DESTROY ALL @ OPTION BASE 1 @ INPUT "N=";N
20 DIM A(N,N),B(N,N),L(N,N),Q(N,N),I0(N,N) @ MAT I0=IDN @ MAT Q=I0
30 MAT INPUT A @ MAT L=A
40 DEF FNF(N,D)=N/(D+(1-2*(D<0))*SQRT(N*N+D*D)) ! = tan(atan(N/D)/2)
100 MAT B=ZER
110 FOR I=2 to N @ FOR J=1 to I-1
120 X=L(I,J) @ IF X THEN B(I,J)=FNF(2*X,L(I,I)-L(J,J))*.5
130 NEXT J @ NEXT I
140 MAT L=TRN(B) @ MAT B=B-L @ DISP "FNORM(B)=";FNORM(B)
150 MAT L=I0-B @ MAT B=I0+B @ MAT L=INV(B)*L @ MAT Q=Q*L
160 DISP "Q=" @ MAT DISP Q
170 MAT L=TRN(Q)*A @ MAT L=L*Q
180 DISP "L=" @ MAT DISP L

Comment: this patch make bij closer to tan(θ/2), with little cost, if desired.

Cas> series(tan(x/2),polynom)      → 1/2*x+1/24*x^3+1/240*x^5
Cas> series(tan(x)/2,polynom)      → 1/2*x+1/6*x^3+1/15*x^5

>120 X=L(I,J) @ IF X THEN X=FNF(2*X,L(I,I)-L(J,J))*.5 @ B(I,J)=X-X*X*X