双共轭梯度预处理(BICGSTAB)(python,数值积分)

第十五课 双共轭梯度处理

采取预处理的原因在上一篇中已经介绍预处理共轭梯度
此篇预处理在第十三课稳定双共轭梯度的基础上
算例采用之前的
详细计算过程可见高斯消元解方程

左手边预处理

#线性联立方程的双共轭梯度法(左手边预处理)
import numpy as np
import math
import B
n=3
converged=np.array([False])
precon=np.zeros((n,1))
v=np.zeros((n,1))
r=np.zeros((n,1))
r0_hat=np.zeros((n,1))
p=np.zeros((n,1))
s=np.zeros((n,1))
t=np.zeros((n,1))
#xnew=np.zeros((n,1))
a=np.array([[10,1,-5],[-20,3,20],[5,3,5]],dtype=np.float)
b=np.array([[1],[2],[6]],dtype=np.float)
x=np.array([[1],[1],[1]],dtype=np.float)
tol=1.0e-5
limit=100
print('系数矩阵')
print(a[:])
print('右手边向量',b[:,0])
print('初始猜测值',x[:,0])
#简单主对角线预处理
for i in range(1,n+1):
    precon[i-1,0]=1.0/a[i-1,i-1]
for i in range(1,n+1):
    a[i-1,:]=a[i-1,:]*precon[i-1,0]
#对左手边应用预处理
b[:]=b[:]*precon[:]
r[:]=b[:]-np.dot(a,x)
r0_hat[:]=r[:]
rho0=1.0
alpha=1.0
w=1.0
v[:]=1.0
p[:]=1.0
rho1=np.dot(np.transpose(r0_hat),r)
print('前几次迭代值')
iters=0
while(True):
    iters=iters+1
    converged[:]=(np.linalg.norm(r)<tol*np.linalg.norm(b))
    if converged==True or iters==limit:
        break
    beta=(rho1/rho0)*(alpha/w)
    p[:,0]=r[:,0]+beta*(p[:,0]-w*v[:,0])
    v[:]=np.dot(a,p)
    alpha=rho1/np.dot(np.transpose(r0_hat),v)
    s[:,0]=r[:,0]-alpha*v[:,0]
    t[:]=np.dot(a,s)
    w=np.dot(np.transpose(t),s)/np.dot(np.transpose(t),t)
    rho0=rho1
    rho1=-w*np.dot(np.transpose(r0_hat),t)
    x[:,0]=x[:,0]+alpha*p[:,0]+w*s[:,0]
    r[:,0]=s[:,0]-w*t[:,0]
    if iters<5:
        print(x[:,0])
print('到收敛需要迭代次数',iters)
print('解向量',x[:,0])   

终端输出结果
程序结果与计算结果一致

右手边预处理

#线性联立方程的双共轭梯度法(左手边预处理)
import numpy as np
import math
import B
n=3
converged=np.array([False])
precon=np.zeros((n,1))
v=np.zeros((n,1))
r=np.zeros((n,1))
r0_hat=np.zeros((n,1))
p=np.zeros((n,1))
p1=np.zeros((n,1))
s=np.zeros((n,1))
s1=np.zeros((n,1))
t=np.zeros((n,1))
#xnew=np.zeros((n,1))
a=np.array([[10,1,-5],[-20,3,20],[5,3,5]],dtype=np.float)
b=np.array([[1],[2],[6]],dtype=np.float)
x=np.array([[1],[1],[1]],dtype=np.float)
tol=1.0e-5
limit=100
print('系数矩阵')
print(a[:])
print('右手边向量',b[:,0])
print('初始猜测值',x[:,0])
#简单主对角线预处理
for i in range(1,n+1):
    precon[i-1,0]=1.0/a[i-1,i-1]
for i in range(1,n+1):
    a[i-1,:]=a[i-1,:]*precon[i-1,0]
#对右手边应用预处理
b[:]=b[:]*precon[:]
r[:]=b[:]-np.dot(a,x)
r0_hat[:]=r[:]
x[:]=x[:]/precon[:]
rho0=1.0
alpha=1.0
w=1.0
v[:]=1.0
p[:]=1.0
rho1=np.dot(np.transpose(r0_hat),r)
print('前几次迭代值')
iters=0
while(True):
    iters=iters+1
    converged[:]=(np.linalg.norm(r)<tol*np.linalg.norm(b))
    if converged==True or iters==limit:
        break
    beta=(rho1/rho0)*(alpha/w)
    p[:,0]=r[:,0]+beta*(p[:,0]-w*v[:,0])
    p1[:,0]=p[:,0]*precon[:,0]
    v[:]=np.dot(a,p1)
    alpha=rho1/np.dot(np.transpose(r0_hat),v)
    s[:,0]=r[:,0]-alpha*v[:,0]
    s1[:,0]=s[:,0]*precon[:,0]
    t[:]=np.dot(a,s1)
    w=np.dot(np.transpose(t),s)/np.dot(np.transpose(t),t)
    rho0=rho1
    rho1=-w*np.dot(np.transpose(r0_hat),t)
    x[:,0]=x[:,0]+alpha*p[:,0]+w*s[:,0]
    r[:,0]=s[:,0]-w*t[:,0]
    if iters<5:
        print(x[:,0]*precon[:,0])
x[:,0]=x[:,0]*precon[:,0]
print('到收敛需要迭代次数',iters)
print('解向量',x[:,0])   

终端输出结果

程序结果与计算结果一致