CG(conjugate gradient method)+Python+Matlab

原创已于 2023-04-18 21:32:18 修改 · 1k 阅读

3 ·

本内容遵循CC 4.0 BY-SA版权协议

收录于

Python语言专栏

Matlab语言专栏

于 2021-12-28 09:29:31 首次发布

该博客介绍了使用Matlab和Python实现共轭梯度法来解决线性方程组的过程。通过代码展示了如何构建特定的矩阵A，并给出迭代求解的详细步骤。在不同维度下，如N=100, 200, 400时，比较了迭代次数。最后，验证了所求解与库函数计算结果的一致性。

没时间敲公式、截个图

![在这里插入图片描述](https://img-blog.csdnimg.cn/76ec933a55d5477695bec23feb47eb4f.jpeg) ![在这里插入图片描述](https://img-blog.csdnimg.cn/f4af91e798694ffe9bd71f66dd2bbd05.jpeg) ![在这里插入图片描述](https://img-blog.csdnimg.cn/a73b077672804202aa6892bea0812db5.jpeg) ![在这里插入图片描述](https://img-blog.csdnimg.cn/9c3e25715bcc48b1845b3d88a40c93ff.jpeg) ![在这里插入图片描述](https://img-blog.csdnimg.cn/d3120e8eb5b643ea991606fe8c9264dd.jpeg) ![在这里插入图片描述](https://img-blog.csdnimg.cn/97faeb36d58d40ca8b131cf3d0bea985.jpeg) ![在这里插入图片描述](https://img-blog.csdnimg.cn/e7a1baa8f4604a4484d85c01138323b9.jpeg) ![在这里插入图片描述](https://img-blog.csdnimg.cn/6cdf092d1bb4487390481c9d7ab84ab2.jpeg)

Matlab程序：

python
% ** 文件名:Conjugate Gradient.m
% 
% ** 日 期:2016.12.14
% 
% ** 描 述:共轭梯度法解线性方程组（Conjugate Gradient method）
% 
% ** 函数：113页计算实习3.2  
% 
% ** 参考教材：《数值分析》李乃成，梅立泉，科学出版社
%
clear;
clc;
%%
%获得需要的N阶矩阵A
N=100; %解向量的维数
%N=200;
%N=400;
a0=eye(N);%a0为n阶的单位阵
a1=eye(N-1);%a1为n-1阶的单位阵
a11=diag(a1);%a11是a1的单位阵的元素
a12=diag(a11,1);%a12是一个把a11的值放到其对角阵上一层的 矩阵
A=a12'+a12+(-2)*a0;%A是对角线及其对角线上下（第n-1,n+1条对角线）有元素的矩阵
%%
%获得矩阵b
b=zeros(N,1);
b(1)=-1;
b(N)=-1;
b;
%%  
fprintf('库函数计算结果：');
x=inv(A)*b      %库函数计算结果
x=zeros(N,1);%迭代近似向量
eps=0.0000001;%精度
r=b-A*x;
d=r;
for k=0:N-1
    fprintf('第%d次迭代：',k+1);
    a=(norm(r)^2)/(d'*A*d)
    x=x+a*d;
    rr=b-A*x;    %rr=r(k+1)
    if (norm(rr)<=eps)||(k==N-1)
        break;
    end
    B=(norm(rr)^2)/(norm(r)^2);
    d=rr+B*d;
    r=rr;
end
x
    
%由上述结果可知方程组由共轭梯度法求解所得的解与计算机直接计算所得的解相同，
%故所得的的解释可靠的。
%当n=100,n=200,n=400时分别需要迭代的次数是50次，100次，200次

Python程序：

# -*- coding: utf-8 -*-
"""
  ** 文件名:Conjugate_Gradient_Method.py
 
  ** 日 期:Created on Tue Dec 28 09:56:40 2021

  ** 描 述:共轭梯度法解线性方程组（Conjugate Gradient method）

  ** 函数：113页计算实习3.2  
 
  ** 参考教材：《数值分析》李乃成，梅立泉，科学出版社

  @author:Mr.Yang
"""
import numpy as np
#获得需要的N阶矩阵A
N=5#解向量的维数
#N=200;
#N=400;
a0=np.eye(N)#a0为n阶的单位阵
a1=np.eye(N-1)#a1为n-1阶的单位阵
a11=np.diag(a1)#a11是a1的单位阵的元素
a12=np.diag(a11,1)#a12是一个把a11的值放到其对角阵上一层的 矩阵
A=a12.T+a12+(-2)*a0#A是对角线及其对角线上下（第n-1,n+1条对角线）有元素的矩阵
#%%
#获得矩阵b
b=np.zeros(N)
b[0]=-1
b[N-1]=-1
 
print('Numpy中计算结果：')
x=np.linalg.inv(A)@b#库函数计算结果
print(x)
x=np.zeros(N)#迭代近似向量
eps=0.0000001#精度
r=b-A@x#
d=r
for k in range(0,N):
    print('第%d次迭代：',(k+1))
    a=(np.linalg.norm(r)**2)/(d.T@A@d)
    x=x+a*d
    rr=b-A@x #rr=r(k+1) 残差
    if(np.linalg.norm(rr)<=eps) or (k==N-1):
        break
    B=(np.linalg.norm(rr)**2)/(np.linalg.norm(r)**2)
    d=rr+B*d
    r=rr
print(x)