ガウス-ザイデル法

• 線型連立方程式を解く
  • 以前はこちらやこちら
  • Gauss-Seidel 法
    • 
    • 
      • , は の対角成分からなる対角行列
  • 対角成分が大きくないといけないらしい

N<-3
A<-matrix(runif(N^2),N,N)
diag(A)<-diag(A)*10    #対角成分が大きくないといけない
x<-runif(N)
b<-runif(N)

D<-matrix(0,N,N)
diag(D)<-diag(A)    #Aの対角成分の行列
B<-A-D      #対角成分以外の要素の行列
diag(D)<-1/diag(D)       #対角行列の逆行列


Niter<-500
for(i in 1:Niter){
x<-D%*%(b-B%*%x)
}

b-A%*%x

• ガウス-ジョルダン法と並べてみる
  • 
  • 拡大係数行列

N<-3
A<-matrix(runif(N^2),N,N)
diag(A)<-diag(A)*10    #対格成分が大きくないといけない
x<-runif(N)
b<-runif(N)

################
#ガウス-ザイデル法#
################
D<-matrix(0,N,N)
diag(D)<-diag(A)    #Aの対角成分の行列
B<-A-D      #対角成分以外の要素の行列
diag(D)<-1/diag(D)       #対角行列の逆行列


Niter<-500
for(i in 1:Niter){
x<-D%*%(b-B%*%x)
}

b-A%*%x
x1<-x
#########################
#ここからがガウス-ジョルダン法#
#########################
C<-cbind(A,b)
N<-nrow(C)
M<-ncol(C)

L<-min(N,M-1)

for(i in 1:L){
C[i,]<-C[i,]/C[i,i]
	for(j in 1:N){
	if(j != i){
	C[j,]<-C[j,]-C[j,i]*C[i,]
	}	
	}
}
C

if(M-N==1){
b-A%*%x
}

x2<-C[,M]

x1
x2
x1-x2