【Julia1.5.2】求解对称正定线性方程组的共轭梯度法
本帖最后由 流星雨 于 2020-11-30 17:13 编辑续watermelon主题：https://www.0xaa55.com/thread-26183-1-1.html中的雅各比迭代法，整理另一种“半迭代半直接”的求解对称正定线性方程组的共轭梯度法。
参考数目：
【1】《数值分析》(李乃成、梅立泉编著，科学出版社，ISBN：978-7-03-032192-3)
采用Julia语言(1.5.2)编制脚本Conjugate Gradient，主要功能均为求解方程组Ax=b。
脚本中的函数Conjugate_Gradient为共轭梯度求解函数，其输入量为矩阵A、尾向量b和迭代偏差ε，输出量为解向量x。
脚本最后提供了以Julia的LinearAlgebra包中附带的QR分解法求得的验证性解向量x2，给出了x2的值、x2与x的误差向量delta、delta的二范数‖delta‖2^2。
结构：
算法原理：
使用说明:
1.安装Julia编译环境JuliaPro-1.5.2-1。可在URL：https://juliacomputing.com/products/juliapro/，下载最新版的JuliaPro-1.5.3-1。
2.安装完毕后用Juliapro打开文件Conjugate gradient.jl，点击▶即可运行。UI右侧Workspace窗口可实时显示各变量值。
3.在脚本文件Conjugate gradient.jl中可以直接更改矩阵A、尾向量b ⃗和迭代偏差ε的值，再次点击▶即可重新计算。
附上参考书【1】给出的使用建议：
1、共轭梯度法只涉及矩阵与向量的乘法和向量内积运算，在理论上对于n阶线性方程组，最多n次迭代即可求得准确解，更像是一种直接解法；
2、对与良态方程，通常所需的迭代次数远小于矩阵阶数n；
3、对于病态问题，迭代次数有所增加，且计算精度可能受限。特别地，当迭代次数=矩阵阶数n时若仍未收敛，则应继续迭代。
源码：
#conjugate gradient Method
using LinearAlgebra
A = [
3 0 1 0 0;
0 5 0 0 0;
1 0 7 0 1;
0 0 0 5 0;
0 0 1 0 3;
]
b = '
ε = 1e-6 #误差
function Conjugate_Gradient(A,b,ε)
N = length(b) #迭代步数
x = ones(length(b),N)*ε #声明临时解向量
r = zeros(length(b),N) #声明临时残向量
d = zeros(length(b),N) #声明临时d向量
α = zeros(N) #声明临时α常数
β = zeros(N) #声明临时β常数
r[:,1] = b - A*x[:,1] #求解r
if r[:,1] == 0#如果初值猜对了，给出x
x[:,N] = x[:,1]
else #如果初值不对，共轭梯度法计算x
d[:,1] = r[:,1] #求解d
global ∑r2 = 1
global k = 1
while ∑r2 >= ε
if k == N #循环次数=矩阵阶数，退出
break
end
#循环次数<矩阵阶数，继续计算
α = r[:,k]'*r[:,k]/(d[:,k]'*A*d[:,k])
x[:,k+1] = x[:,k] + α*d[:,k]
r[:,k+1] = b - A*x[:,k+1]
#求解r的二范数平方
global ∑r2 = 0
for m = 1:length(b)
global ∑r2 = ∑r2 + r^2
end
if (∑r2 <= ε)
break
end
β = r[:,k+1]'*r[:,k+1]/(r[:,k]'*r[:,k])
d[:,k+1] = r[:,k+1] + β*d[:,k]
k = k+1
end
end
return x
end
x = Conjugate_Gradient(A,b,ε)[:,k] #共轭梯度计算结果
x2 = (inv(qr(A).R)*qr(A).Q'*(b)) #QR分解计算结果
delta = zeros(length(b)) #误差向量
Σdelta = 0
for i = 1:length(b) #计算误差向量
delta = (x - x2)/x2
global Σdelta = Σdelta + delta^2
end
Σdelta = Σdelta^0.5 #计算误差二范数
watermelon 发表于 2020-12-26 17:10
Python写的共轭梯度法来啦！
# -*- coding: utf-8 -*-
"""
numpy确实方便！ Python写的共轭梯度法来啦！
# -*- coding: utf-8 -*-
"""
Created on Sat Dec 26 14:54:24 2020
@copyright: Xi'an Jiaotong University - NuTHeL
@author: Zhou Xingguang
"""
import numpy as np
import copy
class ConjugateGradient(object):
def __init__(self, matrix, b, x):
# define the little number
self.__ERR = 1E-6# private parameter
# define the MAX_STEP
self.__MAX_STEP = 1000 # private parameter
self.matrix = np.array(copy.deepcopy(matrix))
self.b = np.array(copy.deepcopy(b)).reshape(np.size(b), 1)
self.x = np.array(copy.deepcopy(x)).reshape(np.size(b), 1)
self.last_x = self.x[:, :]
self.r = np.array(np.zeros((np.size(b), 1), dtype=float))
self.last_r = self.r[:, :]
self.d = np.array(np.zeros((np.size(b), 1), dtype=float))
self.last_d = self.d[:, :]
self.alpha = 0
self.beta = 0
self.time = 0
def CG(self):
self.r = self.b - np.matmul(self.matrix, self.x)
self.last_r = self.r[:, :]
self.d = self.r[:, :]
for time in range(self.__MAX_STEP):
self.alpha = np.linalg.norm(self.r, ord=2)**2 / np.matmul(np.matmul(np.transpose(self.d), self.matrix), self.d)
self.last_x = self.x[:, :]
self.x = self.x + self.alpha * self.d
self.last_r = self.r[:, :]
self.r = self.b - np.matmul(self.matrix, self.x)
self.beta = (np.linalg.norm(self.r, ord=2)**2) / (np.linalg.norm(self.last_r, ord=2)**2)
self.d = self.r + self.beta*self.d
self.alpha = (np.linalg.norm(self.r)**2) / np.matmul(np.matmul(np.transpose(self.d), self.matrix), self.d)
if np.linalg.norm(self.r) <= self.__ERR:
print("iteration has been converged...\n")
self.time = time + 1
return self.x
print("The result is not converged, please recheck the method")
return None
def GenerateExercise(matrix, b):
cnt = np.size(b)
b = -1
b = -1
for i in range(cnt):
if i == 0:
matrix = -2
matrix = 1
elif i == cnt-1:
matrix = -2
matrix = 1
else:
matrix = -2
matrix = 1
matrix = 1
return matrix, b
if __name__ == '__main__':
cnt = 400
matrix = np.zeros((cnt, cnt), dtype=float)
x = np.zeros((cnt, 1), dtype=float)
b = np.zeros((cnt, 1), dtype=float)
matrix, b = GenerateExercise(matrix, b)
CG = ConjugateGradient(matrix, b, x)
result = CG.CG()
print(result)
print(CG.time) 再来一个小弟写的HouseHolder变换的QR分解：
这些数学公式我看了头大。。。
页:
[1]