Python-线性最优解迭代方法

定义矩阵乘法

def mult(h, x):
    result = []
    for x in h:
        summ = 0
        for index, y in enumerate(x):
            summ += y * x[index]
        result.append(summ)
    return result

创建希尔伯特矩阵

def create_hobert(n):
    h=[]
    for x in range(1, n + 1):
        row = []
        for y in range(1, n + 1):
            row.append(1 / (x + y - 1))
        h.append(row)
    return h

雅克比迭代法

def jacobi(A,B,sigma,N):
    n = len(A)
    x0 = []
    x = []
    for i in range(0,n):
        x0.append(0)
        x.append(0)
    for k in range(1,N+1):
        R = 0
        for i in range(0,n):
            sum_ax = 0
            for j in range(0,n):
                sum_ax = sum_ax + A[i][j] * x0[j]
            x[i] = x0[i] + ((B[i] - sum_ax)/A[i][i])
            if abs(x[i] - x0[i]) > R:
                R = abs(x[i] - x0[i])
        if R <= sigma:
            print("精确度等于",sigma,"时,jacobi迭代法需要迭代",k,"次收敛")
            return (x,k)
        for i in range(0,n):
            x0[i] = x[i]
    return (x,k)

Hauss-Seidel迭代法

def gauss_seidel(A,B,sigma,N):
    n = len(A)
    x0 = []
    x = []
    for i in range(0,n):
        x0.append(0)
        x.append(0)
    for k in range(1,N+1):
        R = 0
        for i in range(0,n):
            sum_ax = 0
            for j in range(0,n):
                if j >= i:
                    sum_ax = sum_ax + A[i][j] * x0[j]
                else:
                    sum_ax = sum_ax + A[i][j] * x[j]
            x[i] = x0[i] + ((B[i] - sum_ax)/A[i][i])
            if abs(x[i] - x0[i]) > R:
                R = abs(x[i] - x0[i])
        if R <= sigma:
            print("精确度等于",sigma,"时,gauss_seidel迭代法需要迭代",k,"次收敛")
            return (x,k)
        for i in range(0,n):
            x0[i] = x[i]
    return (x,k)

SOR迭代法

def sor(A,B,sigma,N,omega):
    n = len(A)
    x0 = []
    x = []
    for i in range(0,n):
        x0.append(0)
        x.append(0)
    for k in range(1,N+1):
        R = 0
        for i in range(0,n):
            sum_ax = 0
            for j in range(0,n):
                if j >= i:
                    sum_ax = sum_ax + A[i][j] * x0[j]
                else:
                    sum_ax = sum_ax + A[i][j] * x[j]
            x[i] = x0[i] + omega * ((B[i] - sum_ax)/A[i][i])
            if abs(x[i] - x0[i]) > R:
                R = abs(x[i] - x0[i])
        if R <= sigma:
            print("精确度等于",sigma,"时,松弛因子为",omega,"时,sor迭代法需要迭代",k,"次收敛")
            print(x)
            return (x,k)
        for i in range(0,n):
            x0[i] = x[i]
    return (x,k)

测试希尔伯特矩阵H

if __name__ == "__main__":
    n = 10
    H = create_hobert(n)
    x = [1 for x in range(n)]
    b = mult(H,x)
    print('雅克比迭代法:',jacobi(H, b, 0.1, 200))
    print('SOR迭代法:',sor(H, b, 0.001, 20000, 1.5))
    # print('Guass-Seidel迭代法:',gauss_seidel(H,b,0.00001,20))

原创声明,本文系作者授权云+社区发表,未经许可,不得转载。

如有侵权,请联系 yunjia_community@tencent.com 删除。

编辑于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏人工智能

Apache Spark中的决策树

原文地址:https://dzone.com/articles/decision-trees-in-apache-spark

8568
来自专栏ATYUN订阅号

深度学习图像识别项目(中):Keras和卷积神经网络(CNN)

在下篇文章中,我还会演示如何将训练好的Keras模型,通过几行代码将其部署到智能手机上。

1.9K5
来自专栏机器学习人工学weekly

机器学习人工学weekly-2018/5/20

Dynamic Control Flow in Large-Scale Machine Learning

1036
来自专栏贾志刚-OpenCV学堂

OpenCV中图像算术操作与逻辑操作

在图像处理中有两类最重要的基础操作分别是图像点操作与块操作,简单点说图像点操作就是图像每个像素点的相关逻辑与几何运算、块操作最常见就是基于卷积算子的各种操作、实...

43510
来自专栏CreateAMind

代码开放--模仿学习3篇论文及官方效果视频

https://sites.google.com/view/one-shot-imitation

1032
来自专栏IT派

如何用200行Python代码换张脸

在这篇文章中我将介绍如何写一个简短(200行)的 Python 脚本,来自动地将一幅图片的脸替换为另一幅图片的脸。

992
来自专栏数据结构与算法

HDU3853

Akemi Homura is a Mahou Shoujo (Puella Magi/Magical Girl). Homura wants to hel...

2595
来自专栏用户2442861的专栏

基于级联分类器的多目标检测

原文地址:http://blog.csdn.NET/ariesjzj/article/details/8639208

2511
来自专栏AI深度学习求索

CAM实践:基于pytorch的使用方法

注意:如果为了快一点,不使用网络的图片以及文件的话,记得更改图片地址和已下载文件地址哦

3415
来自专栏C语言及其他语言

【每日一题】问题 1266: 马拦过河卒

棋盘上A点有一个过河卒,需要走到目标B点。卒行走的规则:可以向下、或者向右。同时在棋盘上C点有一个对方的马,该马所在的点和所有跳跃一步可达的点称为对方马的控制点...

1392

扫码关注云+社区