CG(conjugate gradient method)+Python+Matlab

没时间敲公式、截个图        

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)