【数值计算】花式解线性方程组

作者：娟紫恋蓝_610 | 来源：互联网 | 2023-10-12 21:56

Assignment1LU分解codes参考文献Assignment2迭代法

Assignment1 LU分解
- codes
- 参考文献
Assignment2 迭代法
- Problem
- 数据生成
- norm
- slow_Jacobi
- 清真_Jacobi
- Gauss-Seidel
- SOR超松弛法
- 参考文献
Assignment3CGQR
- Problem
- Conjugate Gradient Method 共轭梯度法
- QR Method
- 参考文献

Assignment1: LU分解
老师让只交.py，

于是很多东西直接在注释里写了

codes

# -*-encoding:utf-8-*-
# Numerical Computation & Optimization
# homework1: LU decomposition
# cww97
# 2017/09/27
# from cww97jh@gmail.com
# to num_com_opt@163.com
from numpy import *
n = 100


def lu(A):  #
    """ Doolittle decomposition(course book page45) I just copy the formula from the course book :param A: the Coefficient matrix :return: L, U the lower & upper matrix """
    L = mat(eye(n))  # low matrix
    U = mat(zeros((n, n)))   # upper matrix
    for k in range(n):  # k-th step
        for i in range(k, n):  # calc k-th row of U
            U[k, i] = A[k, i] - sum(L[k, j]*U[j, i] for j in range(k))
        for i in range(k+1, n):  # k-th col of L
            L[i, k] = (A[i, k] - sum(L[i, j]*U[j, k] for j in range(k))) / U[k, k]
    return L, U


def calc(A, B):
    """ calc the equation Ax = b(course book page46) I just copy the formula from the course book :param L: lower triangle matrix :param U: upper triangle matrix :param B: B = A * X :return: the x of A * X = B """
    L, U = lu(A)
    # first, calc L * Y = B, get Y
    Y = mat(zeros((n, 1)))
    for k in range(n):
        Y[k, 0] = B[k, 0] - sum(L[k, j]*Y[j, 0] for j in range(k))
    # then, calc U * X = Y, get X
    X = mat(zeros((n, 1)))
    for k in range(n-1, -1, -1):
        X[k, 0] = (Y[k, 0] - sum(U[k, j]* X[j, 0] for j in range(k+1, n))) / U[k, k]
    return X


""" sample data on course book, for debug A = mat([[2, 1, 5], [4, 1, 12], [-2, -4, 5]]) B = mat([[11], [27], [12]]) """
if __name__ == '__main__':
    # Generate a random matrix M &
    M = mat(random.randint(1, 100, size=(n, n)))
    # A = M + I(Identity matrix)
    A = M + mat(eye(n))
    print('---------matrix A:----------\n', A,)
    # Generate a vector x = (1,2,··· ,100).T
    X = mat([[i+1] for i in range(n)])
    print('---------vector X:----------\n', X.T, '.T')
    # Generate the vector b as b = A * x
    B = A * X
    print('-------B = A * X:-----------\n', B.T, '.T')
    x = calc(A, B)
    print('calc the equation A * x = B, x:\n', x.T, '.T')

""" finished this, when I need to calc some equation, I think I may more likely to use this: x = np.linalg.solve(A, b) However, If I use this, this homework may cannot pass >_

参考文献

python 矩阵操作

Assignment2: 迭代法
TAGS: 数值计算

数值计算与优化指导老师：王祥丰

陈伟文 10152510217

Problem

这里写图片描述

数据生成

def get_data():
    # Generate a random matrix M &
    A = mat(zeros((N, N)))
    for i in range(N):
        A[i, i] = 2
    for i in range(N-1):
        A[i, i + 1] = -1
        A[i + 1, i] = -1
    print('---------matrix A:----------\n', A,)
    # Generate a vector b = (1,1,··· ,1).T
    b = mat(ones((N, 1)))
    print('---------vector X:----------\n', b.T, '.T')
    return A, b

norm

vector norm

这里写图片描述

matrix norm

这里写图片描述

slow_Jacobi

Wiki

核心算法

这里写图片描述

根据最后一个公式，写出如下代码

def norm(x):
    ans = 0
    for t in x:
        ans += square(t[0, 0])
    return sqrt(ans)


def jacobi(A, b):
    n = len(b)
    x0 = mat(zeros((n, 1)))
    x1 = mat(zeros((n, 1)))
    d = 1
    ti = 0
    while d > eps:
        for i in range(n):
            tmp = 0
            for j in range(n):
                if j != i:
                    tmp += A[i, j] * x0[j, 0]
            x1[i, 0] = 1./A[i, i] * (b[i, 0] - tmp)
            #print(x1.T)
        x0 = x1
        d = 1.*norm(A * x0 - b) / norm(b)
        ti += 1
        print('time = %d, d = %lf' % (ti, d))
    return x1

跑了好久好久好久，，，我问主席跑了多久

这里写图片描述为啥他能跑这么快

于是乎我瞄了一眼他的代码，发现，他没像我这样一个值一个值的计算，而是使用numpy自带的矩阵运算
好的，我傻了，又自造车轮，太好了这里写图片描述

于是我决定把上面的函数命名为slow_jacobi，来纪念我这个zz的操作

不过好歹等了十几分钟还是有结果出来的

这里写图片描述

运行速度慢了，不过迭代次数少了，不行，这么慢的代码我不能忍，重写

清真_Jacobi

用这个式子

Xk+1=D−1(b−Rxk)

def jacobi(A, b):  # 这次吸取教训，多用矩阵运算
    # X^(k+1) = D^−1 (b − R * x^k)
    n = len(b)
    x0 = mat(zeros((n, 1)))
    x1 = mat(zeros((n, 1)))
    D = mat(diag(diag(A).tolist()))
    R = A - D
    norm_b = linalg.norm(b)
    d = 1; ti = 0
    while d > eps:
        x1 = D.I * (b - R * x0)
        x0 = x1; ti += 1
        d = linalg.norm(A * x0 - b) / norm_b
        print('time = %d, delta = %.8lf' % (ti, d))
    return x0

迭代次数多了，不过10秒内出解

这里写图片描述

不要自造车轮，不要自造车轮，不要自造车轮

Gauss-Seidel

又是一波紧张刺激的抄公式~~（偷懒直接贴ppt了）~~

这里写图片描述

核心迭代还是抄一下吧

xk+1=BGxk+fG

BG=(D−L)−1，fG=(D−L)−1b

这个迭代可以用x迭代x

def gauss_seidel(A, b):
    # $x ^ {k + 1} = B_G x ^ k + f_G$
    # $B_G = (D - L) ^ {-1}, f_G = (D - L) ^ {-1}b$
    n = len(b)
    D = mat(diag(diag(A).tolist()))
    U = mat(zeros((n, n))) - triu(A, 1)
    L = mat(zeros((n, n))) - tril(A, -1)
    bg = (D - L).I * U
    fg = (D - L).I * b
    norm_b = linalg.norm(b)
    x = mat(zeros((n, 1)))
    d = 1; ti = 0
    while d > eps:
        x = bg * x + fg
        ti += 1
        d = linalg.norm(A * x - b) / norm_b
        print('time = %d, delta = %.8lf' % (ti, d))
    return x

迭代次数少了很多

这里写图片描述

SOR超松弛法

抄ppt

这里写图片描述

xk+1=Bwxk+fw

Bw=(D−wL)−1[(1−w)D+wU]

fw=w(D−wL)−1b

关于 ω 的取值（ppt上写成了 w ），这篇论文，下载之后发现一共只有四页

这里写图片描述

一个跟实验报告一样的论文，还是c语言实现的

互动百科这么说

这里写图片描述

看来这个 ω 取值，在 1 之间取吧 2 ，就是看脸

def sor(A, b, w):
    # $x^{k+1} = B_w x^k + f_w$
    # $B_w = (D - wL)^{-1}[(1-w)D + wU]$
    # $f_w = w(D - wL)^{-1}b$
    n = len(b)
    D = mat(diag(diag(A).tolist()))
    U = mat(zeros((n, n))) - triu(A, 1)
    L = mat(zeros((n, n))) - tril(A, -1)
    bw = (D - w * L).I * ((1-w)*D + w*U)
    fw = w * (D - w*L).I * b

    norm_b = linalg.norm(b)
    x = mat(zeros((n, 1)))
    d = 1; ti = 0
    while d > eps:
        x = bw * x + fw
        ti += 1
        d = linalg.norm(A * x - b) / norm_b
    print('w = %.2lf, time = %d' % (w, ti))
    return x

这个时候我已经不关心是否能出正解了，肯定是正解

为了测试下 ω 的取值对迭代次数的影响

我又写了如下循环

    w = 1.
    while w <= 2:
        w += 0.1
        x = sor(A, b, w)

output:

w = 1.10, time = 11597
w = 1.20, time = 9448
w = 1.30, time = 7630
w = 1.40, time = 6071
w = 1.50, time = 4719
w = 1.60, time = 3536
w = 1.70, time = 2489
w = 1.80, time = 1554
w = 1.90, time = 696

是越大越快吗，我又改循环：

    w = 1.8
    while w <2:
        w += 0.01
        x = sor(A, b, w)

得到如下输出

w = 1.81, time = 1466
w = 1.82, time = 1378
w = 1.83, time = 1292
w = 1.84, time = 1205
w = 1.85, time = 1120
w = 1.86, time = 1035
w = 1.87, time = 950
w = 1.88, time = 866
w = 1.89, time = 781
w = 1.90, time = 696
w = 1.91, time = 609
w = 1.92, time = 520
w = 1.93, time = 423
w = 1.94, time = 293
w = 1.95, time = 303
w = 1.96, time = 404
w = 1.97, time = 505
w = 1.98, time = 808
w = 1.99, time = 1515

我觉得没有必要继续枚举了

我只能说对于本题， ω 取 1.94 的时候速度最快

参考文献

[1] Lecture-3课程ppt

[2] wiki_jacobi

[3] 数值计算——线性方程组的迭代法

[4] 胡枫 ,于福溪《超松弛迭代法中松弛因子 ω的选取方法》青海师范大学学报 2006.01.13 42-46

[5] 互动百科_松弛法

[6] numpy文档_上三角矩阵

[7] numpy文档_norm

Assignment3：CG&QR
陈伟文 10152510217

2017/10/19

Problem

Calculate the equation Ax = b through Conjugate Gradient Method
and QR Method respectively

Conjugate Gradient Method 共轭梯度法

wiki，中文wiki只有寥寥几句，英文的比较详细

贴ppt

这里写图片描述

观察了一下每次循环的逻辑：

先x，用到了上一步的x和p，
r,用到上一步的r和p
p，用到刚刚算出的r和上一步的p

很显然可以发现，x，r，p只需要一个变量就可以完成循环

def conjugate_gradient(A, b):
    n = len(b)
    norm_b = linalg.norm(b)
    # generate x0, r0, p0
    x = mat(zeros((n, 1)))
    r = b - A * x; p = r  # 1
    # here we go
    d = 1; ti = 0
    while d > eps:  #2
        ap = A * p  # 3
        a = (float)(1.*(r.T * r) / (p.T * ap))
        x += a * p  # 4
        r -= a * ap  # 5
        p = r + (float)(1.*(r.T * ap) / (p.T * ap)) * p  #6
        # for counting
        ti += 1
        d = linalg.norm(A * x - b) / norm_b
        print('time = %d, d = %f' % (ti, d))
    return x

我已经不关心方程解出来的正确性了（肯定是对的），看迭代次数吧

这里写图片描述

QR Method

先贴公式

先算Q
这里写图片描述

然后算R

这里写图片描述

（tips：这里的R特别容易算错）

最后： RX=QTb

def QR(A, b):
    n = len(b)
    # calc Q
    R = mat(zeros((n, n)))
    R[0, 0] = linalg.norm(A[:, 0])
    Q = A.copy()
    Q[:, 0] = A[:, 0] / linalg.norm(A[:, 0])
    for j in range(1, n):
        for i in range(j):
            Q[:, j] -= float(A[:, j].T * Q[:, i]) * Q[:, i]
        R[j, j] = linalg.norm(Q[:, j])
        Q[:, j] /= linalg.norm(Q[:, j])
    # calc R
    for i in range(n):
        for j in range(i+1, n):
            R[i, j] = float(A[:, j].T * Q[:, i])
    # calc x from Rx = Q^T * b
    b = Q.T * b
    x = mat(zeros((n, 1)))
    for i in range(n-1, -1, -1):
        x[i] = (b[i] - R[i, i+1:] * x[i+1:]) / R[i, i]
    return x

输出取了整：

calc the equation A * x = B, x:
 [[   50.    99.   147.   194.   240.   285.   329.   372.   414.   455.
    495.   534.   572.   609.   645.   680.   714.   747.   779.   810.
    840.   869.   897.   924.   950.   975.   999.  1022.  1044.  1065.
   1085.  1104.  1122.  1139.  1155.  1170.  1184.  1197.  1209.  1220.
   1230.  1239.  1247.  1254.  1260.  1265.  1269.  1272.  1274.  1275.
   1275.  1274.  1272.  1269.  1265.  1260.  1254.  1247.  1239.  1230.
   1220.  1209.  1197.  1184.  1170.  1155.  1139.  1122.  1104.  1085.
   1065.  1044.  1022.   999.   975.   950.   924.   897.   869.   840.
    810.   779.   747.   714.   680.   645.   609.   572.   534.   495.
    455.   414.   372.   329.   285.   240.   194.   147.    99.    50.]] .T

参考文献

[1] numpyt.dot 计算矩阵内积

[2] dot常见error

[3] numpy线性代数

[4] numpy行列操作

推荐阅读

main
将图像平移到画布中心python_python – Tkinter画布缩放移动/平移

高级缩放示例.就像谷歌地图一样.它仅缩放图块,但不缩放整个图像.因此,缩放的瓷砖占据了恒定的记忆,并且不会为大型缩放图像调整大小的图像.对于简化的缩放示例lookhere.在Win ... [详细]

蜡笔小新 2024-11-20 15:47:04
hash
Go从入门到精通系列视频之go编程语言密码学哈希算法（二）

Go从入门到精通系列视频之go编程语言密码学哈希算法（二） ... [详细]

蜡笔小新 2024-11-21 10:55:36
sum
视觉Transformer综述

本文综述了视觉Transformer在计算机视觉领域的应用，从原始Transformer出发，详细介绍了其在图像分类、目标检测和图像分割等任务中的最新进展。文章不仅涵盖了基础的Transformer架构，还深入探讨了各类增强版Transformer模型的设计思路和技术细节。 ... [详细]

蜡笔小新 2024-11-22 19:53:16
python
Python算法实践：多维缩放技术的应用

本文介绍了多维缩放（MDS）技术，这是一种将高维数据映射到低维空间的方法，通过保持原始数据间的关系，以便于可视化和分析。文章详细描述了MDS的原理和实现过程，并提供了Python代码示例。 ... [详细]

蜡笔小新 2024-11-21 20:04:27
main
UVALive 8201 - BBP 公式计算圆周率

在1995年，Simon Plouffe 发现了一种特殊的求和方法来表示某些常数。两年后，Bailey 和 Borwein 在他们的论文中发表了这一发现，这种方法被命名为 Bailey-Borwein-Plouffe (BBP) 公式。该问题要求计算圆周率 π 的第 n 个十六进制数字。 ... [详细]

蜡笔小新 2024-11-21 18:32:57
main
Ryanair Expands Frankfurt Operations, Challenges Lufthansa's Dominance

Irish budget airline Ryanair announced plans to significantly increase its route network from Frankfurt Airport, marking a direct challenge to Lufthansa, Germany's leading carrier. ... [详细]

蜡笔小新 2024-11-21 13:09:01
python
AI炼金术：KNN分类器的构建与应用

本文介绍了如何使用Python及其相关库（如NumPy、scikit-learn和matplotlib）构建KNN分类器模型。通过详细的数据准备、模型训练及新样本预测的过程，展示KNN算法的实际操作步骤。 ... [详细]

蜡笔小新 2024-11-21 11:40:55
main
c语言二元插值,二维线性插值c语言

c语言二元插值,二维线性插值c语言 ... [详细]

蜡笔小新 2024-11-20 12:20:16
main
线段树详解与实现

本文详细介绍了线段树的基本概念及其在编程竞赛中的应用，并提供了一个具体的线段树实现代码示例。 ... [详细]

蜡笔小新 2024-11-19 21:26:45
chat
诚邀资深后端开发工程师加盟，共谋家乡发展

一家位于长沙的知名网络安全企业，现面向全国诚聘高级后端开发工程师，特别欢迎具有一线城市经验的技术精英回归故乡，共创辉煌。 ... [详细]

蜡笔小新 2024-11-22 19:22:15
python
Ubuntu 14.04 环境下搭建 Caffe（仅限 CPU）

本文详细介绍了如何在 Ubuntu 14.04 系统上搭建仅使用 CPU 的 Caffe 深度学习框架，包括环境准备、依赖安装及编译过程。 ... [详细]

蜡笔小新 2024-11-22 16:43:30
python
为何 TypeScript 如此流行而 Python 类型注解却鲜少使用？

本文探讨了Python类型注解使用率低下的原因，主要归结于历史背景和投资回报率（ROI）的考量。文章不仅分析了类型注解的实际效用，还回顾了Python类型注解的发展历程。 ... [详细]

蜡笔小新 2024-11-22 14:02:28
python
基于OpenCV和Python的边缘检测与四点变换实现

本文介绍了如何利用OpenCV库进行图像的边缘检测，并通过Canny算法提取图像中的边缘。随后，文章详细说明了如何识别图像中的特定形状（如矩形），并应用四点变换技术对目标区域进行透视校正。 ... [详细]

蜡笔小新 2024-11-22 11:31:52
main
解决iOS应用推送通知错误：未找到有效aps-environment权限

在尝试加载支持推送通知的iOS应用程序的Ad Hoc构建时，遇到了‘no valid aps-environment entitlement found for application’的错误提示。本文将探讨此错误的原因及多种可能的解决方案。 ... [详细]

蜡笔小新 2024-11-21 19:26:31
python
Markdown 编辑技巧详解

本文介绍如何使用 Typora 编辑器高效编写 Markdown 文档，包括代码块的插入方法等实用技巧。Typora 官方网站：https://www.typora.io/ 学习资源：https://www.markdown.xyz/ ... [详细]

蜡笔小新 2024-11-20 23:42:54