技术宅的结界

 找回密码
 立即注册→加入我们

QQ登录

只需一步,快速开始

搜索
热搜: 下载 VB C 实现 编写
查看: 1254|回复: 4
收起左侧

【Julia1.5.2】求解对称正定线性方程组的共轭梯度法

[复制链接]

4

主题

7

帖子

392

积分

用户组: 中·技术宅

UID
6285
精华
4
威望
43 点
宅币
208 个
贡献
71 次
宅之契约
0 份
在线时间
4 小时
注册时间
2020-10-2
发表于 2020-11-30 16:41:32 | 显示全部楼层 |阅读模式

欢迎访问技术宅的结界,请注册或者登录吧。

您需要 登录 才可以下载或查看,没有帐号?立即注册→加入我们

x
本帖最后由 流星雨 于 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。

结构:
}K7HIHMQ12C~QN0%B~HR.png

算法原理:
1.png 2.png 3.png 4.png

使用说明:
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 = [1 2 3 4 5]'
ε = 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[0]
    if r[:,1] == 0#如果初值猜对了,给出x
        x[:,N] = x[:,1]
    else #如果初值不对,共轭梯度法计算x
        d[:,1] = r[:,1] #求解d[0]
        global ∑r2 = 1
        global k = 1
        while ∑r2 >= ε
            if k == N #循环次数=矩阵阶数,退出
                break
            end
             #循环次数<矩阵阶数,继续计算
            α[k] = r[:,k]'*r[:,k]/(d[:,k]'*A*d[:,k])
            x[:,k+1] = x[:,k] + α[k]*d[:,k]
            r[:,k+1] = b - A*x[:,k+1]
            #求解r的二范数平方
            global ∑r2 = 0
            for m = 1:length(b)
                global ∑r2 = ∑r2 + r[m,k+1]^2
            end
            if (∑r2 <= ε)
                break
            end
            β[k] = r[:,k+1]'*r[:,k+1]/(r[:,k]'*r[:,k])
            d[:,k+1] = r[:,k+1] + β[k]*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 #计算误差二范数

Conjugate gradient.jl

1.75 KB, 下载次数: 1

评分

参与人数 1威望 +10 宅币 +30 贡献 +10 收起 理由
watermelon + 10 + 30 + 10 赞!

查看全部评分

回复

使用道具 举报

4

主题

7

帖子

392

积分

用户组: 中·技术宅

UID
6285
精华
4
威望
43 点
宅币
208 个
贡献
71 次
宅之契约
0 份
在线时间
4 小时
注册时间
2020-10-2
 楼主| 发表于 2020-12-26 19:34:51 | 显示全部楼层
watermelon 发表于 2020-12-26 17:10
Python写的共轭梯度法来啦!
[mw_shl_code=python,true]# -*- coding: utf-8 -*-
"""

numpy确实方便!

29

主题

331

帖子

1999

积分

用户组: 上·技术宅

UID
3808
精华
11
威望
105 点
宅币
1238 个
贡献
165 次
宅之契约
0 份
在线时间
378 小时
注册时间
2018-5-6
发表于 2020-12-26 17:10:45 | 显示全部楼层
Python写的共轭梯度法来啦!
[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[0, 0] = -1
    b[cnt-1, 0] = -1
    
    for i in range(cnt):
        if i == 0:
            matrix[i, i] = -2
            matrix[i, i+1] = 1
        elif i == cnt-1:
            matrix[i, i] = -2
            matrix[i, i-1] = 1
        else:
            matrix[i, i] = -2
            matrix[i, i-1] = 1
            matrix[i, i+1] = 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)
Passion Coding!

29

主题

331

帖子

1999

积分

用户组: 上·技术宅

UID
3808
精华
11
威望
105 点
宅币
1238 个
贡献
165 次
宅之契约
0 份
在线时间
378 小时
注册时间
2018-5-6
发表于 2020-12-10 22:25:51 | 显示全部楼层
再来一个小弟写的HouseHolder变换的QR分解:
QQ图片20201210222527.jpg
Passion Coding!

1088

主题

2606

帖子

7万

积分

用户组: 管理员

一只技术宅

UID
1
精华
236
威望
474 点
宅币
21344 个
贡献
45935 次
宅之契约
0 份
在线时间
2054 小时
注册时间
2014-1-26
发表于 2020-12-11 01:55:58 | 显示全部楼层
这些数学公式我看了头大。。。

本版积分规则

QQ|申请友链||Archiver|手机版|小黑屋|技术宅的结界 ( 滇ICP备16008837号 )|网站地图

GMT+8, 2021-10-16 16:33 , Processed in 0.040050 second(s), 35 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表