Мой код R BiCG не работает должным образом - PullRequest
Новая
       37

Мой код R BiCG не работает должным образом

0 голосов
/ 28 августа 2018

Я написал код для решения линейного уравнения Ax = b, используя метод BiCG, но решение не сходилось хорошо, и относительный остаток, установленный в условии завершения, не мог быть очищен. Я не знаю, является ли код R неправильным или решение не сходится в принципе, могу ли я получить какое-то мнение? Условие прекращения - когда норма (r, type = "2") / норма (b, type = "2") падает ниже 10 ^ - 12, как написано в перерыве. 

install.packages('SparseM')
library(SparseM)

ganma=0.5

#Create sparse 6, 6 matrix
 i <- c(1,1,rep(2, times=3),rep(3,times=3), rep(4,times=3), rep(5,times=3),6,6) #Line number of nonzero element
 j <- c(1,2,1,2,3,2,3,4,3,4,5,4,5,6,5,6) # Column number of non-zero element
 g <- c(2,1,rep(c(ganma,2,1), time=4),ganma,2) #value of non-zero element

#Preparation of matrix and preparation of b
A <- sparseMatrix(i= i, j= j, x=g,dims = c(6,6))
b= as.vector(A%*%matrix(data=c(1,1,1,1,1,1)))

A = as.matrix(A)

#Solution by BiCG method

x = c(0.8,0.8,0.8,0.8,0.8,0.8)
r = b-A%*%x

rk = r
p = r
pk = rk

xlist = list()
nlist = list()

#repeat part
 for(i in  1:1000000){
  q = A%*%p
  qk = t(Conj(A))%*%pk

  alpha_k = c(as.vector(rk)%*%as.vector(r) / as.vector(pk)%*%as.vector(q))

  x2 = x + alpha_k*p
  r2 = r - alpha_k*q
  rk2 = rk - Conj(alpha_k)*qk

  beta_k =  c(as.vector(rk2)%*%as.vector(r2) / as.vector(rk)%*%as.vector(r))

  p2 = r2 + beta_k*p
  pk2 = rk + Conj(beta_k)*pk


  x = x2
  r = r2
  rk = rk2
  p = p2
  pk = pk2
  xlist = append(xlist,  list(x)) 
  nlist = append(nlist, list(norm(r, type="2")/norm(b, type="2")))

  if(norm(r, type="2")/norm(b, type="2") <= 10**-12) break
}
...