２．１　行列の加減算と乗算

■表2-1  Pythonによるリストを配列とみなす
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)]
A=array(3,3); A[1][1]=10
print(A)
 
■表2-2  行列の加算プログラム
def addM(A,B):#行列加算
    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]+B[i][j]
    return C
def addV(A,B):#ベクトル加算
    N1=len(A); C=array(N1)
    for i in range(N1):
        C[i]=A[i]+B[i]
    return C
def addVM(A,B):#ベクトル加算の繰り返しによる行列加算
    N1=len(A); N2=len(A[0]);C=[]
    for i in range(N1):C.append(addV(A[i],B[i]))
    return C
A=[[7,8,9],[6,5,4], [3,2,1]]#
B=[[1,2,3],[1,1,1], [2,2,2]]
C= addM(A,B)     ;print(C)
C=addV(A[0],B[0]);print(C)
C=addVM(A,B)     ;print(C)

■表2-3  行列の減算プログラム
def subM(A,B):#行列減算
    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]-B[i][j]
    return C
def subV(A,B):#ベクトル減算
    N1=len(A); C=array(N1)
    for i in range(N1):
        C[i]=A[i]-B[i]
    return C
def subVM(A,B):#ベクトル減算の繰り返しによる行列減算
    N1=len(A); N2=len(A[0]);C=[]
    for i in range(N1):C.append(subV(A[i],B[i]))
    return C
A=[[7,8,9],[6,5,4], [3,2,1]]
B=[[1,2,3],[1,1,1], [2,2,2]]
C=subM(A,B)      ;print(C)
C=subV(A[0],B[0]);print(C)
C=subVM(A,B)     ;print(C)

■表2-4  行列の乗算プログラム
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 mltV(A,B):#ベクトル乗算
    NA=len(A); NB=len(B);V=0; NN=NA
    if NN>NB:NN=NB
    for i in range(NA): V+=A[i]*B[i]
    return V
def mltCRV(A,B):#列ベクトル*行ベクトル乗算
    NA=len(A); NB=len(B);C=array(NA,NB)
    for i in range(NA):
        for j in range(NB):
            C[i][j]=A[i]*B[j]
    return C
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
A=[[6,5,4],[12,13,10], [18,21,17]]
B=[[2,3,1],[ 3, 4, 4], [ 1, 2, 5]]
C=transM(A); print(C)
C=mltM(A,B)      ;print(C)
V=mltV(A[0],B[0]);print(V)
C=mltCRV(A[0],B[0]);print(C)
C=mltM(transM([A[0]]),[B[0]]);print(C)#上記と同じ意味になる


２．２　ピボット入替え

■表2-5 置換行列
def unitSqMat(N):#単位正方行列
    C=array(N,N)
    for i in range(N):C[i][i]=1
    return C
def replaceMat(N, p,q):#置換行列
    C=unitSqMat(N); D=C[p]; C[p]=C[q]; C[q]=D
    return C
A=[[1,2,3],[4,5,6],[7,8,9]]; print("     A = ",A)
R=replaceMat(3,0,1);print("     R = ",R) #R : 置換行列 
print(" R * R = ", mltM(R,R))
print(" R * A = ", mltM(R,A))
print(" A * R = ", mltM(A,R))

２．３　三角分解(LU分解)とQR分解

■表2-6  Pythonによるクラウト法
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]
A=[[6,5,4],[12,13,10],[18,21,17]]
S=LUdecom(A)
print(S[0]);print(S[1]); print(mltM(S[0],S[1]))

■表2-7  PythonによるQR分解プログラム
import math
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]
A=[[1,-2,1],[-2,2,1],[-2,1,1]]
S=QRdecom(A)
print("    A =",A)
print("    Q =",S[0])
print("    R =",S[1])
print("Q * R =",mltM(S[0],S[1]))

２．４　連立方程式の解法

■表2-8 VBAによるガウスの消去法
def copyArray(A):
    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 GaussianELM(MAT):
    A=copyArray(MAT);N=len(A);N1=N+1;ESP=0.0000001
    for k in range(N):# 前進消去
        if abs(A[k][k]) < ESP:
            print("解を求めることができません")
            return []
        for i in range(k + 1,N):
            ALFA = A[i][k] / A[k][k]
            for j in range(k + 1,N1):A[i][j]-=ALFA*A[k][j]
    B=array(N,1)
    for i in range(N-1, -1, -1):# 後退代入
        T = A[i][N]
        for j in range(i + 1,N):T -=A[i][j] * A[j][N]
        A[i][N] = T / A[i][i]; B[i][0]=A[i][N]
    return B
A=[[1,2,3,4],[-9,-8,-1,0],[-7,-4,4,9]]
B=GaussianELM(A)
print("A=",A)
print("B=",B)
print("A*B=",mltM(A,B))

■表2-9 Pythonによる掃出し法（処理本体のみ：他はガウスの消去法と同じ）
def GaussJordan(MAT):
    A=copyArray(MAT);N=len(A);N1=N+1;ESP=0.0000001
    for k in range(N):# 前進消去
        if abs(A[k][k]) < ESP:
            print("解を求めることができません")
            return []
        k1 = k + 1
        for j in range(k1,N1): A[k][j] /= A[k][k]
        for i in range(N):
            if i != k:
                for j in range(k1,N1):A[i][j]-= A[i][k] * A[k][j]
    B=array(N,1)
    for k in range(N):B[k][0]=A[k][N]
    return B

■表2-10  ピボット選択付きガウスの消去法（処理本体と確認処理のみ）
def GaussianELMPivot(MAT):
    A=copyArray(MAT);N=len(A);N1=N+1;ESP=0.0000001
    for k in range(N):# 前進消去
        k1 = k + 1;Max = abs(A[k][k]);IR = k #----ピボット入れ替え
        if k != N :
            for i in range(k1,N):
                ab=abs(A[i][k])
                if ab>= Max : Max = ab; IR = i
        if Max < ESP :
            print("解を求めることができません")
            return []
        if IR != k :
            D=A[k]; A[k]=A[IR]; A[IR]=D #---ピボット入れ替えここまで
        for i in range(k + 1,N):
            ALFA = A[i][k] / A[k][k]
            for j in range(k + 1,N1):A[i][j]-=ALFA*A[k][j]
    B=array(N,1)
    for i in range(N-1, -1, -1):# 後退代入
        T = A[i][N]
        for j in range(i + 1,N):T -=A[i][j] * A[j][N]
        A[i][N] = T / A[i][i]; B[i][0]=A[i][N]
    return B
A=[[0,5,4,8],[12,0,10,16],[18,21,0,27]]
B=GaussianELMPivot(A)
print("A=",A)
print("B=",B)
print("A*B=",mltM(A,B))

■表2-11  逆行列の計算
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 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 invM(MAT):#逆行列を求める
    A=setInvTemp(MAT);N=len(A);N2=N*2;ESP=0.0000001
    for k in range(N):# 前進消去
        k1 = k + 1;Max = abs(A[k][k]);IR = k #----ピボット入れ替え
        if k != N :
            for i in range(k1,N):
                ab=abs(A[i][k])
                if ab>= Max : Max = ab; IR = i
        if Max < ESP :
            print("解を求めることができません")
            return []
        if IR != k :
            D=A[k]; A[k]=A[IR]; A[IR]=D #---ピボット入れ替えここまで
        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
A=[[0,5,4],[12,0,10],[18,21,0]]
B=invM(A)
printM("A=",A)
printM("B=",B)
printM("A*B=",mltM(A,B))

■表2-12  LU分解を使った連立方程式の解（「LU分解」は表2-6参照）
def simEqLU(A,B):#LU分解による連立方程式の解法
    N=len(B); Y=array(N);X=array(N)
    LU=LUdecom(A);L=LU[0];U=LU[1]#LU分解
    for i in range(N):#下三角行列による変換
        T=B[i]
        for k in range(i): T-=L[i][k]*Y[k]
        Y[i]=T
    #print(Y) #必要なら注釈を外して表示
    for i in range(N-1,-1,-1):#上三角行列による後退代入
        T=Y[i]
        for k in range(i + 1, N):  T -= U[i][k]*X[k]
        X[i] =T/U[i][i]
    return X
def colVector(B):#ベクトルを列ベクトルに変換
    N=len(B); X=array(N,1)
    for i in range(N):X[i][0]=B[i]
    return X
A=[[6,5,4],[12,13,10],[18,21,17]]
B=[8,16,27]
X=simEqLU(A,B)
print("A=",A)
print("B=",B)
print("X=",X)
print("A*X=",mltM(A,colVector(X)))

■表2-13  ピボット選択によるLU分解
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 transMat(A,K):#変換行列
    N=len(A); G=unitSqMat(N); AKK = A[K][K]
    for i in range(K + 1,N): G[i][K] = -A[i][K] / AKK
    return G
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 LUdecompPivot(A):#ピボット選択によるLU分解
    N=len(A); H=unitSqMat(N); multP=unitSqMat(N)
    for i in range(N-1):
        P=replaceMat(N, i,maxDiag(A,i))#置換行列Piの設定
        C=mltM(P,A)# C=Pi*Ai-1; G=transMat(C,i) # 行列Giを求める
        A=mltM(G,C)     # Ai= Gi*C = Gi*Pi*Ai-1
        if i>0: H=mltM(mltM(P,mltM(GB,H)),P)# Hi=Pi*Gi-1*Hi-2*Pi^(-1)
        multP=copyArray(mltM(P,multP)); GB=copyArray(G)
    return [invM(mltM(G,H)),A]
A=[[6,5,4],[12,13,10],[18,21,17]]
R= LUdecompPivot(A)
printM("A=",A)
printM("L=",R[0])
printM("U=",R[1])

■表2-14  ピボット選択によるLU分解の改善
def LUdecompPivot2(A):#ピボット選択によるLU分解(改良版)
    N=len(A); LU=copyArray(A); IPivot=array(N);EPS=0.00000001
    for k in range(N):IPivot[k]=k
    for k in range(N):
        i=k
        for j in range(k + 1,N):
            if abs(LU[IPivot[j]][k])>abs(LU[IPivot[i]][k]):i=j
        j = IPivot[k]; IPivot[k] = IPivot[i];IPivot[i]=j
        pivot=LU[IPivot[k]][k]
        if abs(pivot)< EPS:
            print("%番目の対角項が%です" %(k, pivot))
        for i in range(k+1,N):
            T = LU[IPivot[i]][k]/pivot; LU[IPivot[i]][k] = T
            for j in range(k+1,N):
                LU[IPivot[i]][j]-=T * LU[IPivot[k]][j]
    R=array(N,N)
    for k in range(N): R[k]=LU[IPivot[k]]
    return R
def seperateLU(LU):
    N=len(LU); L=array(N,N);U=array(N,N);
    for i in range(N):
        for j in range(N):
            if j>=i:
                U[i][j]=LU[i][j]
                if j==i:L[i][j]=1
            else: L[i][j]=LU[i][j]
    return [L,U]
A=[[6,5,4],[12,13,10],[18,21,17]]
printM("A=",A);LU=LUdecompPivot2(A);printM("LU=",LU)
LandU=seperateLU(LU)
printM("L=",LandU[0])
printM("U=",LandU[1])

■表2-15  ヤコビの反復法
idef JacobiIter(A, B, iterMax,EPS):
    N=len(A); X=array(N); Y=array(N); EEPS=EPS*EPS
    for i in range(N): X[i] = B[i] / A[i][i]
    for k in range(iterMax):
        for i in range(N):
            T = B[i]
            for j in range(N):
                if i!=j : T-= A[i][j] * X[j]
            Y[i] = T / A[i][i]
        dev = 0
        for i in range(N):
            EI=(Y[i]-X[i])/Y[i];dev = EI*EI;X[i]=Y[i]
        if E<EEPS : break
    return [k+1, X]
A=[[3,1,1],[1,5,2],[1,2,5]]
B=[10,21,30]
R=JacobiIter(A, B, 500, 0.000000000001)
print("繰り返し回数", R[0], "結果", R[1])
C=mltM(A,colVector(R[1]))
print("A*X=", C)

■表2-16  ガウス・ザイデル法
def GaussSeidel(A, B, iterMax,EPS):
    N=len(A); X=array(N); Y=array(N);EEPS=EPS*EPS
    for i in range(N):X[i] = B[i] / A[i][i]
    for k in range(iterMax):
        dev = 0
        for i in range(N):
            T = B[i]
            for j in range(N):
                if i!=j : T-= A[i][j] * X[j]
            Y[i] = T / A[i][i]
            EI = X[i]-Y[i];dev+=EI*EI;X[i]=Y[i]
        if E<EEPS : break
    return [k+1, X]
A=[[3,1,1],[1,5,2],[1,2,5]]
B=[10,21,30]
print("B=",B)
R=GaussSeidel(A, B, 500, 0.000000000001)
print("繰り返し回数", R[0], "結果", R[1])
C=mltM(A,colVector(R[1]))
print("A*X=", C)

■表2-17  SOR法
def SOR(A, B, iterMax, Omega, EPS):
    N=len(A); X=array(N); Y=array(N); EEPS=EPS*EPS
    for i in range(N):X[i] = B[i] / A[i][i]
    for k in range(iterMax):
        dev = 0
        for i in range(N):
            T = B[i]
            for j in range(N):
                if i!=j : T-= A[i][j] * X[j]
            Y[i] = T / A[i][i]; EI = Y[i]-X[i];
            dev+=EI*EI;X[i] += Omega * EI
        if E<EEPS : break
    return [k+1, X]
A=[[0,2,1],[1,3,2],[2,4,1],[1,5,2],[2,6,1],[1,7,0]]
D=[1,2,3,4,5,6]; print("D=",D)
R=SOR(A, B, 500, 1.45, 0.000000000001)
print("繰り返し回数", R[0], "結果", R[1])
C=mltM(A,colVector(R[1]))
print("A*X=", C)

■表2-18  三重対角な行列の連立方程式
def array(N1, N2=0):#配列(表2-1再掲)
    if N2!=0: return [[0 for i in range(N2)] for i in range(N1)]
    return [0 for i in range(N1)]
def printTridiagonal(title, DT,form="%8.4f"):#三重対角行列の表示
    N=len(DT); Nstr=len(form % 0); Bstr=" "*Nstr
    print(title)
    for i in range(N):
	if i==0: S= Bstr+form % DT[i][1]+form % DT[i][2]
        else:
            S=Bstr*i+(form % DT[i][0])+(form % DT[i][1])
            if i<(N-1):S+=form % DT[i][2]
        print("%5d : " % (i+1)+S)

def multTriDiag(A,X):#三重対角行列とベクトルの乗算
    N=len(A); D=array(N); D[0]=A[0][1]*X[0]+A[0][2]*X[1]
    for i in range(1,N-1):
        D[i]=A[i][0]*X[i-1]+A[i][1]*X[i]+A[i][2]*X[i+1]
    D[N-1]=A[N-1][0]*X[N-2]+A[N-1][1]*X[N-1]
    return D

def TriDiag(A, D):#三重対角行列による連立方程式の解法
    N=len(A); X=array(N); CC=array(N)
    btemp = A[0][1]; X[0] = D[0] / btemp
    for j in range(1,N):
        CC[j] = A[j-1][2] / btemp
        btemp = A[j][1] - A[j][0]* CC[j]
        X[j] = (D[j] - A[j][0] * X[j - 1]) / btemp
    for j in range(N - 2,-1,-1): X[j]-=CC[j + 1] * X[j + 1]
    return X