共轭梯度法python实现

python代码

        先上代码,如果代码有不懂的地方,可以继续看后面的数学理论依据,算法参考西安交大数值分析教材(李乃成,梅立泉著)

import numpy as np
import matplotlib.pyplot as plt

class conjugate_gradient_method():
    """
    共轭梯度法:A是n阶对称正定矩阵
    输入:
    A:系数矩阵
    b:b就是b
    n:矩阵阶数
    x0:初始向量
    precision:精度
    输出:
    求解结果
    """
    def __init__(self,A,b,x0,n,precision):
        self.A=np.array(A)
        self.n=n
        self.precision=precision
        self.b=b
        self.x0=np.array(x0)
        self.x=[]
        self.x.append(list(x0))
        self.iter=0

    #计算向量范数
    def Norm2(self,p): #p为输入向量
        sum_of_p=sum([i**2 for i in p])
        return sum_of_p**0.5
    
    #共轭梯度法求解
    def calculate(self):
        r0=self.b-np.dot(self.A,self.x0)
        d0=r0
        r_pre=r0
        d_pre=d0
        for k in range(self.n):
            temp=np.dot(d_pre.T,self.A)
            alpha=np.dot(r_pre.T,r_pre)/np.dot(temp,d_pre)
            x=np.array(self.x[-1])+np.dot(alpha,d_pre)
            r=self.b-np.dot(self.A,x)
            self.x.append(list(x))
            if self.Norm2(r)<=self.precision:
                self.iter=k+1
                break
            else:
                beta=(self.Norm2(r)**2)/(self.Norm2(r_pre)**2)
                d=r+beta*d_pre
            r_pre = r
            d_pre = d
        print("计算结果为:",self.x[-1])

    #各个分量的差的绝对值求和
    def c_error(self,a,b):
        result=[]
        for j,k in zip(a,b):
            temp=np.array(j)-np.array(k)
            temp_new=sum([abs(i) for i in temp])
            result.append(np.log10(temp_new)) #计算log10
        return(np.array(result))

    def picture_show(self):
        horizontal=np.linspace(0,self.iter-1,len(self.x)-1)
        #print(horizontal)
        vertical_1=self.x[0:len(self.x)-1]
        vertical_2=[self.x[-1]]*(len(self.x)-1)
        vertical=self.c_error(vertical_1,vertical_2)
        #print(vertical)
        plt.plot(horizontal,vertical)
        plt.title("error")
        plt.xlabel("iter")
        plt.ylabel("log10-error")
        plt.show()

#############测试################
if __name__ == "__main__":
    n=200
    #构造向量b
    b=[0]*n
    b[-1]=-1
    b[0]=-1
    #构造矩阵A
    A=np.zeros((n,n))
    for i in range(n):
        for j in range(n):
            if i==0:
                A[i][0]=-2
                A[i][1]=1
            elif i==n-1:
                A[i][n-1]=-2
                A[i][n-2]=1
            elif i==j:
                A[i][j]=-2
                A[i][j-1]=1
                A[i][j+1]=1
            
    print(A)
    print(b)
    x0=[0]*n #初始值向量
    precision=0.001
    my_cgm=conjugate_gradient_method(A,b,x0,n,precision) #初始化
    my_cgm.calculate() #使用共轭梯度法求解
    print("迭代次数:", my_cgm.iter)
    my_cgm.picture_show() #显示收敛速度

共轭梯度法原理

共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次函数极小化的问题,从任意给定的初始点出发,沿一组关于矩阵A的共轭方向进行线性搜索,在无舍入误差的假定下,最多迭代n次(n为矩阵的阶数),就可以求得二次函数的极小点,也就求得了线性方程组Ax=b的解。

共轭梯度法在形式上具有迭代法的特征,给定初始向量x(0),由迭代格式

 产生的迭代序列x(1),x(2),x(3),… 在无舍入误差的假定下,最多经过n次迭代就可求得f(x)的最小点。

 

算法伪代码

 

  • 6
    点赞
  • 35
    收藏
    觉得还不错? 一键收藏
  • 4
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值