（３．１，３．２はプログラム例なし）

■３．３　反復法による固有値

表3-1 べき乗法のプログラム
import math
def array(N1, N2=0):#配列
    if N2!=0: return [[0 for i in range(N2)] for i in range(N1)]
    return [0 for i in range(N1)]
def transM(A):#転置行列
    N1=len(A); N2=len(A[0]);C=array(N2,N1)
    for i in range(N1):
        for j in range(N2):C[j][i]=A[i][j]
    return C
def mltScalarV(S, B):
    N=len(B); C=array(N)
    for i in range(N): C[i]=S*B[i]
    return C
def mltM(A,B):#行列乗算
    NA1=len(A); NA2=len(A[0]);NB1=len(B); NB2=len(B[0]);NN=NA2
    C=array(NA1,NB2)
    if NN>NB1: NN=NB1
    for i in range(NA1):
        for j in range(NB2):
            C[i][j]=0
            for k in range(NN): C[i][j]+= A[i][k]*B[k][j]
    return C
def PowerIter(A, EPS):#べき乗法による固有値
    N=len(A); X=array(N);XTemp=array(2,N)
    IterMax = 100; K1 = 1;  K2 = 0; XTemp[K2][0] = 1
    for Iter in range(IterMax):
        KK = K2; K2 = K1; K1 = KK ## K1：前回計算，K2;今回
        for i in range(N): # 行列式の乗算
            T = 0
            for j in range(N):
                T += A[i][j] * XTemp[K1][j]
            XTemp[K2][i] = T
        T = 0   #ノルム1=1による方法
        for i in range(N): T += XTemp[K2][i]*XTemp[K2][i]
        lamb = math.sqrt(T)
        if abs(lamb) < EPS : return [] #固有値が0のとき収束しない
        for i in range(N): XTemp[K2][i] /= lamb
        T1 = 0; T2 = 0
        for i in range(N):
            DX = XTemp[K2][i]- XTemp[K1][i]
            T1 +=DX * DX; T2 += XTemp[K2][i] * XTemp[K2][i]
        if abs(T2) < EPS: return [] #固有ベクトルが0のとき収束しない
        E = T1 / T2
        if E<EPS: break
    for i in range(N): X[i] = XTemp[K2][i]
    return [lamb, X]
def printVect(Evect, Mvect, Lvect):
    N=len(Evect)
    print("\n行番号  固有値(λ) 行列*固有ベクトル λ*固有ベクトル")
    print("-"*45)
    for i in range(N):
        print("%5d, %8.5f,    %8.5f,         %8.5f"%
              (i,Evect[i],Mvect[i][0],Lvect[i])) 
A=[[4,2,1,3],[2,3,1,2],[1,1,2,1],[3,2,1,1]]
R=PowerIter(A, 0.00000000001)
if R==[]: print("収束しません")
else:
    Evalue=R[0]; Evect=R[1]
    print("固有値(λ) =",Evalue)
    printVect(Evect,mltM(A,transM([Evect])),
              mltScalarV(Evalue, Evect))

表3-2 逆行列による最小固有値の計算
def setInvTemp(A):#逆行列用作業領域設定
    N1=len(A); N2=N1*2; C=array(N1,N2)
    for i in range(N1):
        for j in range(N1):
            C[i][j]=A[i][j];C[i][j+N1]=0
            if i==j: C[i][j+N1]=1
    return C
def maxDiag(A,k):
    Max = abs(A[k][k]);IR = k #----ピボット入れ替え
    for i in range(k+1,len(A)):
        ab=abs(A[i][k])
        if ab>= Max : Max = ab; IR = i
    return IR
def replaceLine(k1, k2,A):#行入れ替え
    D=A[k1]; A[k1]=A[k2]; A[k2]=D
def invM(MAT):#ピボット選択消去法による逆行列（改良版）
    A=setInvTemp(MAT);N=len(A);N2=N*2;ESP=0.0000001
    for k in range(N):# 前進消去
        k1 = k + 1; IR=maxDiag(A,k); Max=abs(A[IR][k]) 
        if Max<ESP:print("解を求めることができません");return []
        if IR != k :replaceLine(k,IR,A)
        for i in range(k + 1,N):
            ALFA = A[i][k] / A[k][k]
            for j in range(k + 1,N2):A[i][j]-=ALFA*A[k][j]
    for i in range(N-1, -1, -1):# 後退代入
        for k in range(N):
            T = A[i][N+k]
            for j in range(i + 1,N):T -=A[i][j] * A[j][N+k]
            A[i][N+k] = T / A[i][i]
    B=array(N,N)
    for i in range(N):
        for j in range(N): B[i][j]=A[i][j+N]
    return B
def PowerIterInv(A, EPS):#逆行列を使ったべき乗法による最小固有値
    B=invM(A); R= PowerIter(B, EPS); R[0]=1/R[0]
    return R
A=[[4,2,1,3],[2,3,1,2],[1,1,2,1],[3,2,1,1]]
R=PowerIterInv(A, 0.000001)
if R==[]: print("収束しません")
else:
    Evalue=R[0]; Evect=R[1]
    print("固有値(λ) =",Evalue)
    printVect(Evect,mltM(A,transM([Evect])),
              mltScalarV(Evalue, Evect))

表3-3  LR分解による固有値の計算
def array(N1, N2=0):#配列
    if N2!=0: return [[0 for i in range(N2)] for i in range(N1)]
    return [0 for i in range(N1)]
def mltM(A,B):#行列乗算
    NA1=len(A); NA2=len(A[0]);NB1=len(B); NB2=len(B[0]);NN=NA2
    C=array(NA1,NB2)
    if NN>NB1: NN=NB1
    for i in range(NA1):
        for j in range(NB2):
            C[i][j]=0
            for k in range(NN): C[i][j]+= A[i][k]*B[k][j]
    return C
def copyArray(A):#配列複写(表2-8)
    N1=len(A); N2=len(A[0]); C=array(N1,N2)
    for i in range(N1):
        for j in range(N2): C[i][j]=A[i][j]
    return C
def LUdecom(A):#LU分解
    N=len(A)
    L=array(N,N);U=array(N,N); L[0][0] = 1
    for j in range(N): U[0][j] = A[0][j]
    for i in range(1,N):
        U[i][0]=0;L[0][i]=0;L[i][0]=A[i][0]/U[0][0]
    for i in range(1,N):
        L[i][i]= 1; T = A[i][i]
        for k in range(i):T-=L[i][k]*U[k][i]
        U[i][i] = T
        for j in range(i+1,N):
            U[j][i]=0; L[i][j]=0; T=A[j][i]
            for k in range(i): T-=L[j][k]*U[k][i]
            L[j][i] = T/U[i][i]; T=A[i][j]
            for k in range(i): T-=L[i][k]*U[k][j]
            U[i][j] = T
    return [L,U]
def eigenLR(MAT,EPS):
    A=copyArray(MAT);IterMax=100; E2=EPS*EPS; N=len(A)
    for R in range(IterMax):
        LU=LUdecom(A);A=mltM(LU[1],LU[0]); E = 0
        for i in range(1,N):
            for j in range(i-1): E+=A[i][j]*A[i][j]
        if E<E2: break
    return A
A=[[4,2,1,3],[2,3,1,2],[1,1,2,1],[3,2,1,1]]
B=eigenLR(A,0.0000000001)
for i in range(4): print(B[i][i])

表3-4  QR分解による固有値の計算
def QRdecom(A):#QR分解
    N=len(A);X=array(N); Q=array(N,N); R=array(N,N)
    for k in range(N):
        for i in range(N):X[i] = A[i][k]
        for j in range(k-1):
            Rjk = 0
            for j in range(N):Rjk +=A[i][k]* Q[i][j]
            R[j][k] = Rjk; R[k][j] = 0
            for i in range(N): X[i]-= Rjk * Q[i][j]
        T = 0
        for i in range(N):T+= X[i] * X[i]
        R[k][k] = math.sqrt(T)
        for j in range(N): Q[j][k] = X[j] / R[k][k]
    return [Q,R]
def eigenQR(MAT,EPS):
    A=copyArray(MAT);IterMax=100; E2=EPS*EPS; N=len(A)
    for R in range(IterMax):
        QR=LUdecom(A); A=mltM(QR[1],QR[0]); E = 0
        for i in range(1,N):
            for j in range(i-1):E+=A[i][j]*A[i][j]
        if E<E2: break
    return A
A=[[4,2,1,3],[2,3,1,2],[1,1,2,1],[3,2,1,1]]
B=eigenQR(A,0.0000000001)
for i in range(4): print(B[i][i])

■３．４　行列変換法による固有値

表3-5  ヤコビ法による固有値の計算
def printM(title,A,form="%10.5f"):#行列表示
    print(title)
    for i in range(len(A)):
        S="%5d:" % (i+1)
        for j in range(len(A[i])):
            S+=form %A[i][j]
        print(S)
def printV(title,A,form="%10.5f"):#行列表示
    S=title
    for i in range(len(A)):
        S+=form %A[i]
    print(S)
def eigenJacobi(MAT):
    A=copyArray(MAT); N=len(A); D=array(N);W=array(N,N)
    for k in range(N): W[k][k]=1; D[k]=A[k][k]
    MaxIter = 50
    for Iter in range(MaxIter):
        OK=True
        for j in range(N-1):
            jj=j+1
            for k in range(jj,N):
                T = abs(float(D[j])) + abs(float(D[k]))
                if (abs(A[j][k]) + T) != T :   
                    OK = False
                    T =(D[k]-D[j])/(2.0*A[j][k])
                    if T>=0 :T = 1.0/(T+math.sqrt(T*T+1.0))
                    else:    T = 1.0/(T-math.sqrt(T*T+1.0))
                    C = 1.0/math.sqrt(T*T+1); S = T*C; T *= A[j][k]
                    D[j] = D[j]-T; D[k] = D[k]+T; A[j][k] = 0
                    for i in range(j):
                        X = A[i][j]; Y = A[i][k]
                        A[i][j] = X*C - Y*S; A[i][k] = X*S + Y*C
                    for i in range(j + 1, k):
                        X = A[j][i]; Y = A[i][k]
                        A[j][i] = X*C - Y*S; A[i][k] = X*S + Y*C
                    for i in range(k + 1, N):
                        X = A[j][i]; Y = A[k][i]
                        A[j][i] = X*C - Y*S; A[k][i] = X*S + Y*C
                    for i in range(N):
                        X = W[j][i]; Y = W[k][i]
                        W[j][i] = X*C - Y*S; W[k][i] = X*S + Y*C
        if OK: break
    if OK:
        for i in range(N-1):
            k=i
            for j in range(i + 1, N):
                if D[j] > D[k] : k = j
            T = D[i]; D[i] = D[k]; D[k] = T
            T = W[i]; W[i] = W[k]; W[k] = T
        for i in range(N):
            for j in range(N): A[i][j]=W[j][i]
        return [D,A]
    else: print("収束しません"); return[]
def mltEigenV(D,W):#固有値配列×固有ベクトル行列
    N=len(D); C=array(N,N)
    for i in range(N):
        for j in range(N):
            C[j][i]=D[i]*W[j][i]
    return C
A=[[4, 2, 1, 3], [2, 3, 1, 2], [1, 1, 2, 1], [3, 2, 1, 1]]
DW=eigenJacobi(A); D=DW[0]; W=DW[1]
printV("  λ =",D); printM("   W =",W)
printM(" A*W =", mltM(A,W)); printM("λ*W =",mltEigenV(D,W))

表3-6  三重対角化による固有値の計算
copyArrayについては表3-3, 
mltEigenV，printM，printVについては表3-5
を参照してください。）---
def tridiagonalize(MAT,EPS):
    A=copyArray(MAT)
    N=len(A);V=array(N); W=array(N+1);D=array(N);E=array(N)
    for k in range(N-2):
        D[k] = A[k][k]
        for i in range(k+1,N): W[i] = A[i][k]
        S = 0
        for i in range(k+1,N):S += W[i]*W[i]
        if abs(S) > EPS:
            if W[k+1]>=0: S = math.sqrt(S)
            else: S = -math.sqrt(S)
            E[k+1] = -S; W[k+1] += S
            S = 1 / math.sqrt(W[k+1]*S)
            for i in range(k+1, N): W[i] *= S
            for i in range(k+1, N):
                S = 0
                for j in range(k+1, i+1):S += A[i][j]*W[j]
                for j in range(i+1, N): S += A[j][i] * W[j]
                V[i] = S
            S = 0
            for i in range(k+1,N): S +=W[i]*V[i]
            S /= 2
            for i in range(k+1,N): V[i] -= S * W[i]
            for i in range(k+1,N):
                for j in range(k+1,i+1):
                    A[i][j] = A[i][j] - W[i]*V[j]-W[j]*V[i]
            # 固有ベクトルを求めるときWの内容をAに移す
            for i in range(k+1,N):  A[i][k] = W[i]
        else: E[k+1] = 0
    D[N-2]= A[N-2][N-2]; D[N-1] = A[N-1][N-1]#制限値に注意
    E[0] = 0; E[N-1] = A[N-1][N-2]
    #以下，固有ベクトル計算用
    for k in range(N-1, -1,-1):
        if k < N - 1:
            for i in range(k+1,N): W[i] = A[i][k]
            for i in range(k+1,N):
                S = 0
                for j in range(k+1,N): S += A[i][j]*W[j]
                V[i] = S
            for i in range(k+1, N):
                for j in range(k+1,N): A[i][j] -= V[i]*W[j]
        for i in range(N): A[i][k] = 0
        A[k][k]= 1
    return [D,E,A]
def Zero(D,E,i):
    T=0
    if i>0:T=abs(D[i])+abs(D[i-1])
    return (T+abs(E[i]))==T
def eigenHouseholder(MAT,EPS):
    MaxIter = 20; N=len(MAT)
    DEA=tridiagonalize(MAT,EPS); D=DEA[0]; E=DEA[1]; A=DEA[2]
    E[0] = 1.0; H = N-1
    while Zero(D, E, H):H-=1
    while H > 0:
        E[0] = 0.0; L = H - 1
        while not Zero(D, E, L): L -=1
        for Iter in range(MaxIter):
            W = (D[H-1]-D[H])/2; S = math.sqrt(W*W+E[H]*E[H])
            if W < 0 : S = -S
            X = D[L]-D[H]+E[H]*E[H]/(W+S); Y = E[L+1]; Z = 0
            for k in range(L,H):
                if abs(X) >= abs(Y):
                    T = -Y/X; C = 1/math.sqrt(T*T+1); S = T*C
                else:
                    T = -X/Y; S = 1/math.sqrt(T*T+1)
                    if T < 0 : S = -S
                    C = T*S
                W = D[k]-D[k+1]; T = (W*S+2*C*E[k+1])*S
                D[k]-=T;  D[k+1]+= T
                E[k] = C*E[k]-S*Z; E[k+1] = E[k+1]*(C*C-S*S)+W*S*C
                for i in range(N):
                    X = A[k][i]; Y = A[k+1][i]
                    A[k][i] = C*X - S*Y; A[k+1][i] = S*X + C*Y
                if k < N - 2 :
                    X = E[k+1]; Y = -S*E[k+2];Z = Y;E[k+2] = C*E[k+2]
            if Zero(D, E, H):break
        E[0] = 1.0; H -=1
        while Zero(D, E, H): H -= 1
    for k in range(N-1):#以下，固有ベクトル
        j = k
        for i in range(k+1,N):
            if D[i] > D[j] :L = i
        T = D[k]; D[k] = D[j]; D[j] = T
        T = A[k]; A[k] = A[j]; A[j] = T
    return [D,transM(A)]  
A=[[4,2,1,3],[2,3,1,2],[1,1,2,1],[3,2,1,1]]
Z=eigenHouseholder(A,0.0000000001)
D=Z[0]; W=Z[1]
printV("  λ =",D)
printM("   W =",W)
printM(" A*W =",mltM(A,W))
printM("λ*W =",mltEigenV(D,W))