码迷,mamicode.com
首页 > 编程语言 > 详细

多线性方程组迭代算法——Gauss-Seidel迭代算法的Python实现

时间:2019-06-29 01:13:50      阅读:129      评论:0      收藏:0      [点我收藏+]

标签:youdao   tar   eal   art   index   ros   mda   ash   div   

多线性方程组(张量)迭代算法的原理请看这里:原理部分请留言,不方便公开分享

Jacobi迭代算法里有详细注释:多线性方程组迭代算法——Jacobi迭代算法的Python实现

import numpy as np
import time

1.1 Gauss-Seidel迭代算法

def GaussSeidel_tensor_V2(A,b,Delta,m,n,M):
    start=time.perf_counter()
    find=0
    X=np.ones(n)
    d=np.ones(n)
    m1=m-1
    m2=2-m
    for i in range(M):
        print(X,X)
        x=np.copy(X)
        #迭代更新
        for j in range(n):
            a=np.copy(A)
            for k in range(m-2):
                a=np.dot(a,x)
            for k in range(n):
                d[k]=a[k,k]
                a[k,k]=m2*a[k,k]            
            x[j]=(b[j]-np.dot(a[j],x))/(m1*d[j])
        #判断是否满足精度要求
        if np.max(np.fabs(X-x))<Delta:
            find=1
            break  
        X=np.copy(x)
    end=time.perf_counter()
    print(时间:,end-start)
    print(迭代,i)
    return X,find,i,end-start

1.2张量A的生成函数和向量b的生成函数:

def Creat_A(m,n):#生成张量A
    size=np.full(m, n)
    X=np.ones(n)
    while 1:
        #随机生成给定形状的张量A
        A=np.random.randint(-49,50,size=size)
        #判断Dx**(m-2)是否非奇异,如果是,则满足要求,跳出循环
        D=np.copy(A)
        for i1 in range(n):
            for i2 in range(n):
                if i1!=i2:
                    D[i1,i2]=0
        for i in range(m-2):
                D=np.dot(D,X)
        det=np.linalg.det(D)
        if det!=0:
            break
    #将A的对角面张量扩大十倍,使对角面占优
    for i1 in range(n):
        for i2 in range(n):
            if i1==i2:
                A[i1,i2]=A[i1,i2]*10
    print(A:)
    print(A)
    return A

#由A和给定的X根据Ax**(m-1)=b生成向量b
def Creat_b(A,X,m):
    a=np.copy(A)
    for i in range(m-1):
        a=np.dot(a,X)
    print(b:)
    print(a)
    return a

1.3 对称张量S的生成函数:

def Creat_S(m,n):#生成对称张量B
    size=np.full(m, n)
    S=np.zeros(size)
    print(S,S)
    for i in range(4):
        #生成n为向量a
        a=np.random.random(n)*np.random.randint(-5,6)
        b=np.copy(a)
        #对a进行m-1次外积,得到秩1对称张量b
        for j in range(m-1):
            b=outer(b,a)
        #将不同的b叠加得到低秩对称张量S
        S=S+b
    print(S:)
    print(S)
    return S
def outer(a,b):
    c=[]
    for i in b:
        c.append(i*a)
    return np.array(c)
    return a

1.4 实验二

def test_2():
    Delta=0.01#精度
    m=3#A的阶数
    n=3#A的维数
    M=200#最大迭代步数
    X_real=np.array( [2,3,4])
    A=Creat_A(m,n) 
    b=Creat_b(A,X_real,m)
    GaussSeidel_tensor_V2(A,b,Delta,m,n)

 

多线性方程组迭代算法——Gauss-Seidel迭代算法的Python实现

标签:youdao   tar   eal   art   index   ros   mda   ash   div   

原文地址:https://www.cnblogs.com/Fengqiao/p/Gauss-Seidel-tensor.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!