线性代数-Python-04:线性系统+高斯消元的实现

这篇具有很好参考价值的文章主要介绍了线性代数-Python-04:线性系统+高斯消元的实现。希望对大家有所帮助。如果存在错误或未考虑完全的地方,请大家不吝赐教,您也可以点击"举报违法"按钮提交疑问。

1 线性系统

线性系统 csdn,线性代数python,python,线性代数,开发语言
线性系统 csdn,线性代数python,python,线性代数,开发语言
线性系统 csdn,线性代数python,python,线性代数,开发语言

2 高斯-jordon消元法的实现

2.1 Matrix

from .Vector import Vector

class Matrix:

    def __init__(self, list2d):
        self._values = [row[:] for row in list2d]

    @classmethod
    def zero(cls, r, c):
        """返回一个r行c列的零矩阵"""
        return cls([[0] * c for _ in range(r)])

    @classmethod
    def identity(cls, n):
        """返回一个n行n列的单位矩阵"""
        m = [[0]*n for _ in range(n)]
        for i in range(n):
            m[i][i] = 1;
        return cls(m)

    def T(self):
        """返回矩阵的转置矩阵"""
        return Matrix([[e for e in self.col_vector(i)]
                       for i in range(self.col_num())])

    def __add__(self, another):
        """返回两个矩阵的加法结果"""
        assert self.shape() == another.shape(), \
            "Error in adding. Shape of matrix must be same."
        return Matrix([[a + b for a, b in zip(self.row_vector(i), another.row_vector(i))]
                       for i in range(self.row_num())])

    def __sub__(self, another):
        """返回两个矩阵的减法结果"""
        assert self.shape() == another.shape(), \
            "Error in subtracting. Shape of matrix must be same."
        return Matrix([[a - b for a, b in zip(self.row_vector(i), another.row_vector(i))]
                       for i in range(self.row_num())])

    def dot(self, another):
        """返回矩阵乘法的结果"""
        if isinstance(another, Vector):
            # 矩阵和向量的乘法
            assert self.col_num() == len(another), \
                "Error in Matrix-Vector Multiplication."
            return Vector([self.row_vector(i).dot(another) for i in range(self.row_num())])

        if isinstance(another, Matrix):
            # 矩阵和矩阵的乘法
            assert self.col_num() == another.row_num(), \
                "Error in Matrix-Matrix Multiplication."
            return Matrix([[self.row_vector(i).dot(another.col_vector(j)) for j in range(another.col_num())]
                           for i in range(self.row_num())])

    def __mul__(self, k):
        """返回矩阵的数量乘结果: self * k"""
        return Matrix([[e * k for e in self.row_vector(i)]
                       for i in range(self.row_num())])

    def __rmul__(self, k):
        """返回矩阵的数量乘结果: k * self"""
        return self * k

    def __truediv__(self, k):
        """返回数量除法的结果矩阵:self / k"""
        return (1 / k) * self

    def __pos__(self):
        """返回矩阵取正的结果"""
        return 1 * self

    def __neg__(self):
        """返回矩阵取负的结果"""
        return -1 * self

    def row_vector(self, index):
        """返回矩阵的第index个行向量"""
        return Vector(self._values[index])

    def col_vector(self, index):
        """返回矩阵的第index个列向量"""
        return Vector([row[index] for row in self._values])

    def __getitem__(self, pos):
        """返回矩阵pos位置的元素"""
        r, c = pos
        return self._values[r][c]

    def size(self):
        """返回矩阵的元素个数"""
        r, c = self.shape()
        return r * c

    def row_num(self):
        """返回矩阵的行数"""
        return self.shape()[0]

    __len__ = row_num

    def col_num(self):
        """返回矩阵的列数"""
        return self.shape()[1]

    def shape(self):
        """返回矩阵的形状: (行数, 列数)"""
        return len(self._values), len(self._values[0])

    def __repr__(self):
        return "Matrix({})".format(self._values)

    __str__ = __repr__

2.2 Vector

import math
from ._globals import EPSILON
class Vector:

    def __init__(self, lst):
        """
        __init__ 代表类的构造函数
        双下划线开头的变量 例如_values,代表类的私有成员
        lst是个引用,list(lst)将值复制一遍,防止用户修改值
        """
        self._values = list(lst)

    def underlying_list(self):
        """返回向量的底层列表"""
        return self._values[:]

    def dot(self, another):
        """向量点乘,返回结果标量"""
        assert len(self) == len(another), \
            "Error in dot product. Length of vectors must be same."
        return sum(a * b for a, b in zip(self, another))

    def norm(self):
        """返回向量的模"""
        return math.sqrt(sum(e**2 for e in self))

    def normalize(self):
        """
        归一化,规范化
        返回向量的单位向量
        此处设计到了除法: def __truediv__(self, k):
        """
        if self.norm() < EPSILON:
            raise ZeroDivisionError("Normalize error! norm is zero.")
        return Vector(self._values) / self.norm()
        # return 1 / self.norm() * Vector(self._values)
        # return Vector([e / self.norm() for e in self])

    def __truediv__(self, k):
        """返回数量除法的结果向量:self / k"""
        return (1 / k) * self

    @classmethod
    def zero(cls, dim):
        """返回一个dim维的零向量
        @classmethod 修饰符对应的函数不需要实例化,不需要 self 参数,但第一个参数需要是表示自身类的cls参数,可以来调用类的属性,类的方法,实例化对象等。
        """
        return cls([0] * dim)

    def __add__(self, another):
        """向量加法,返回结果向量"""
        assert len(self) == len(another), \
            "Error in adding. Length of vectors must be same."
        # return Vector([a + b for a, b in zip(self._values, another._values)])
        return Vector([a + b for a, b in zip(self, another)])

    def __sub__(self, another):
        """向量减法,返回结果向量"""
        assert len(self) == len(another), \
            "Error in subtracting. Length of vectors must be same."
        return Vector([a - b for a, b in zip(self, another)])

    def __mul__(self, k):
        """返回数量乘法的结果向量:self * k"""
        return Vector([k * e for e in self])

    def __rmul__(self, k):
        """
        返回数量乘法的结果向量:k * self
        self本身就是一个列表
        """
        return self * k

    def __pos__(self):
        """返回向量取正的结果向量"""
        return 1 * self

    def __neg__(self):
        """返回向量取负的结果向量"""
        return -1 * self

    def __iter__(self):
        """返回向量的迭代器"""
        return self._values.__iter__()

    def __getitem__(self, index):
        """取向量的第index个元素"""
        return self._values[index]

    def __len__(self):
        """返回向量长度(有多少个元素)"""
        return len(self._values)

    def __repr__(self):
        """打印显示:Vector([5, 2])"""
        return "Vector({})".format(self._values)

    def __str__(self):
        """打印显示:(5, 2)"""
        return "({})".format(", ".join(str(e) for e in self._values))

2.3 线性系统

from .Matrix import Matrix
from .Vector import Vector


class LinearSystem:

    def __init__(self, A, b):

        assert A.row_num() == len(b), "row number of A must be equal to the length of b"
        self._m = A.row_num()
        self._n = A.col_num()
        assert self._m == self._n  # TODO: no this restriction

        self.Ab = [Vector(A.row_vector(i).underlying_list() + [b[i]])
                   for i in range(self._m)]

    def _max_row(self, index_i, index_j, n):

        best, ret = abs(self.Ab[index_i][index_j]), index_i
        for i in range(index_i + 1, n):
            if abs(self.Ab[i][index_j]) > best:
                best, ret = abs(self.Ab[i][index_j]), i
        return ret

    def _forward(self):

        n = self._m
        for i in range(n):
            # Ab[i][i]为主元
            max_row = self._max_row(i, i, n)
            self.Ab[i], self.Ab[max_row] = self.Ab[max_row], self.Ab[i]

            # 将主元归为一
            self.Ab[i] = self.Ab[i] / self.Ab[i][i]  # TODO: self.Ab[i][i] == 0?
            for j in range(i + 1, n):
                self.Ab[j] = self.Ab[j] - self.Ab[j][i] * self.Ab[i]

    def _backward(self):

        n = self._m
        for i in range(n - 1, -1, -1):
            # Ab[i][i]为主元
            for j in range(i - 1, -1, -1):
                self.Ab[j] = self.Ab[j] - self.Ab[j][i] * self.Ab[i]

    def gauss_jordan_elimination(self):

        self._forward()
        self._backward()

    def fancy_print(self):

        for i in range(self._m):
            print(" ".join(str(self.Ab[i][j]) for j in range(self._n)), end=" ")
            print("|", self.Ab[i][-1])

3 行最简形式

线性系统 csdn,线性代数python,python,线性代数,开发语言

4 线性方程组的结构

线性系统 csdn,线性代数python,python,线性代数,开发语言
线性系统 csdn,线性代数python,python,线性代数,开发语言

5 线性方程组-通用高斯消元的实现

5.1 global

# 包中的变量,但是对包外不可见,因此使用“_”开头
EPSILON = 1e-8


def is_zero(x):
    return abs(x) < EPSILON


def is_equal(a, b):
    return abs(a - b) < EPSILON

5.2 Vector-引入is_zero

import math
from ._globals import is_zero
class Vector:

    def __init__(self, lst):
        """
        __init__ 代表类的构造函数
        双下划线开头的变量 例如_values,代表类的私有成员
        lst是个引用,list(lst)将值复制一遍,防止用户修改值
        """
        self._values = list(lst)

    def underlying_list(self):
        """返回向量的底层列表"""
        return self._values[:]

    def dot(self, another):
        """向量点乘,返回结果标量"""
        assert len(self) == len(another), \
            "Error in dot product. Length of vectors must be same."
        return sum(a * b for a, b in zip(self, another))

    def norm(self):
        """返回向量的模"""
        return math.sqrt(sum(e**2 for e in self))

    def normalize(self):
        """
        归一化,规范化
        返回向量的单位向量
        此处设计到了除法: def __truediv__(self, k):
        """
        if is_zero(self.norm()):
            raise ZeroDivisionError("Normalize error! norm is zero.")
        return Vector(self._values) / self.norm()
        # return 1 / self.norm() * Vector(self._values)
        # return Vector([e / self.norm() for e in self])

    def __truediv__(self, k):
        """返回数量除法的结果向量:self / k"""
        return (1 / k) * self

    @classmethod
    def zero(cls, dim):
        """返回一个dim维的零向量
        @classmethod 修饰符对应的函数不需要实例化,不需要 self 参数,但第一个参数需要是表示自身类的cls参数,可以来调用类的属性,类的方法,实例化对象等。
        """
        return cls([0] * dim)

    def __add__(self, another):
        """向量加法,返回结果向量"""
        assert len(self) == len(another), \
            "Error in adding. Length of vectors must be same."
        # return Vector([a + b for a, b in zip(self._values, another._values)])
        return Vector([a + b for a, b in zip(self, another)])

    def __sub__(self, another):
        """向量减法,返回结果向量"""
        assert len(self) == len(another), \
            "Error in subtracting. Length of vectors must be same."
        return Vector([a - b for a, b in zip(self, another)])

    def __mul__(self, k):
        """返回数量乘法的结果向量:self * k"""
        return Vector([k * e for e in self])

    def __rmul__(self, k):
        """
        返回数量乘法的结果向量:k * self
        self本身就是一个列表
        """
        return self * k

    def __pos__(self):
        """返回向量取正的结果向量"""
        return 1 * self

    def __neg__(self):
        """返回向量取负的结果向量"""
        return -1 * self

    def __iter__(self):
        """返回向量的迭代器"""
        return self._values.__iter__()

    def __getitem__(self, index):
        """取向量的第index个元素"""
        return self._values[index]

    def __len__(self):
        """返回向量长度(有多少个元素)"""
        return len(self._values)

    def __repr__(self):
        """打印显示:Vector([5, 2])"""
        return "Vector({})".format(self._values)

    def __str__(self):
        """打印显示:(5, 2)"""
        return "({})".format(", ".join(str(e) for e in self._values))

5.3 LinearSystem

from .Matrix import Matrix
from .Vector import Vector
from ._globals import is_zero


class LinearSystem:

    def __init__(self, A, b):

        assert A.row_num() == len(b), "row number of A must be equal to the length of b"
        self._m = A.row_num()
        self._n = A.col_num()
        # assert self._m == self._n  # TODO: no this restriction

        self.Ab = [Vector(A.row_vector(i).underlying_list() + [b[i]])
                   for i in range(self._m)]
        self.pivots = []

    def _max_row(self, index_i, index_j, n):

        best, ret = abs(self.Ab[index_i][index_j]), index_i
        for i in range(index_i + 1, n):
            if abs(self.Ab[i][index_j]) > best:
                best, ret = abs(self.Ab[i][index_j]), i
        return ret

    def _forward(self):

        i, k = 0, 0
        while i < self._m and k < self._n:
            # 看Ab[i][k]位置是否可以是主元
            max_row = self._max_row(i, k, self._m)
            self.Ab[i], self.Ab[max_row] = self.Ab[max_row], self.Ab[i]

            if is_zero(self.Ab[i][k]):
                k += 1
            else:
                # 将主元归为一
                self.Ab[i] = self.Ab[i] / self.Ab[i][k]
                for j in range(i + 1, self._m):
                    self.Ab[j] = self.Ab[j] - self.Ab[j][k] * self.Ab[i]
                self.pivots.append(k)
                i += 1

    def _backward(self):

        n = len(self.pivots)
        for i in range(n - 1, -1, -1):
            k = self.pivots[i]
            # Ab[i][k]为主元
            for j in range(i - 1, -1, -1):
                self.Ab[j] = self.Ab[j] - self.Ab[j][k] * self.Ab[i]

    def gauss_jordan_elimination(self):
        """如果有解,返回True;如果没有解,返回False"""
        self._forward()
        self._backward()

        for i in range(len(self.pivots), self._m):
            if not is_zero(self.Ab[i][-1]):
                return False
        return True

    def fancy_print(self):

        for i in range(self._m):
            print(" ".join(str(self.Ab[i][j]) for j in range(self._n)), end=" ")
            print("|", self.Ab[i][-1])

5.4 main

from playLA.Matrix import Matrix
from playLA.Vector import Vector
from playLA.LinearSystem import LinearSystem


if __name__ == "__main__":

    A = Matrix([[1, 2, 4], [3, 7, 2], [2, 3, 3]])
    b = Vector([7, -11, 1])
    ls = LinearSystem(A, b)
    ls.gauss_jordan_elimination()
    ls.fancy_print()
    print()
    # [-1, -2, 3]

    A2 = Matrix([[1, -3, 5], [2, -1, -3], [3, 1, 4]])
    b2 = Vector([-9, 19, -13])
    ls2 = LinearSystem(A2, b2)
    ls2.gauss_jordan_elimination()
    ls2.fancy_print()
    print()
    # [2, -3, -4]

    A3 = Matrix([[1, 2, -2], [2, -3, 1], [3, -1, 3]])
    b3 = Vector([6, -10, -16])
    ls3 = LinearSystem(A3, b3)
    ls3.gauss_jordan_elimination()
    ls3.fancy_print()
    print()
    # [-2, 1, -3]

    A4 = Matrix([[3, 1, -2], [5, -3, 10], [7, 4, 16]])
    b4 = Vector([4, 32, 13])
    ls4 = LinearSystem(A4, b4)
    ls4.gauss_jordan_elimination()
    ls4.fancy_print()
    print()
    # [3, -4, 0.5]

    A5 = Matrix([[6, -3, 2], [5, 1, 12], [8, 5, 1]])
    b5 = Vector([31, 36, 11])
    ls5 = LinearSystem(A5, b5)
    ls5.gauss_jordan_elimination()
    ls5.fancy_print()
    print()
    # [3, -3, 2]

    A6 = Matrix([[1, 1, 1], [1, -1, -1], [2, 1, 5]])
    b6 = Vector([3, -1, 8])
    ls6 = LinearSystem(A6, b6)
    ls6.gauss_jordan_elimination()
    ls6.fancy_print()
    print()
    # [1, 1, 1]

    A7 = Matrix([[1, -1, 2, 0, 3],
                 [-1, 1, 0, 2, -5],
                 [1, -1, 4, 2, 4],
                 [-2, 2, -5, -1, -3]])
    b7 = Vector([1, 5, 13, -1])
    ls7 = LinearSystem(A7, b7)
    ls7.gauss_jordan_elimination()
    ls7.fancy_print()
    print()

    A8 = Matrix([[2, 2],
                 [2, 1],
                 [1, 2]])
    b8 = Vector([3, 2.5, 7])
    ls8 = LinearSystem(A8, b8)
    if not ls8.gauss_jordan_elimination():
        print("No Solution!")
    ls8.fancy_print()
    print()

    A9 = Matrix([[2, 0, 1],
                 [-1, -1, -2],
                 [-3, 0, 1]])
    b9 = Vector([1, 0, 0])
    ls9 = LinearSystem(A9, b9)
    if not ls9.gauss_jordan_elimination():
        print("No Solution!")
    ls9.fancy_print()
    print()

线性系统 csdn,线性代数python,python,线性代数,开发语言文章来源地址https://www.toymoban.com/news/detail-769031.html

到了这里,关于线性代数-Python-04:线性系统+高斯消元的实现的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

  • MATLAB数值分析学习笔记:线性代数方程组的求解和高斯-赛德尔方法

    迭代法是前面介绍的消元法的有效替代,线性代数方程组常用的迭代法有 高斯-赛德尔方法 和 雅克比迭代法, 下面会讲到二者的不同之处,大家会发现两者的实现原理其实类似,只是方法不同,本篇只重点介绍高斯-赛德尔方法。 看了我之前的笔记的同学应该已经对迭代法不

    2024年02月05日
    浏览(60)
  • 04 MIT线性代数-矩阵的LU分解 Factorization into A=LU

    目的: 从矩阵的角度理解高斯消元法, 完成 LU 分解得到 A = LU U 为上三角阵(Upper triangular matrix),  L 为下三角阵(Lower triangular matrix), 通过分解得到对角阵 D (diagonal matrix) 设定一组消元矩阵,其中 E31 为单位阵 I ,其它两个消元矩阵如下: row3 -5 newrow2 = row3 -5( row2 -2 row1 )= row3 -

    2024年02月07日
    浏览(41)
  • 线性代数1:线性方程和系统

    Digital Collection (staedelmuseum.de) 图片来自施泰德博物馆         通过这些文章,我希望巩固我对这些基本概念的理解,同时如果可能的话,通过我希望成为一种基于直觉的数学学习方法为其他人提供额外的清晰度。如果有任何错误或机会需要我进一步阐述,请分享,我可以进

    2024年02月06日
    浏览(45)
  • 基于python的线性代数运算

    前言:这是学校多元统计分析课程布置的实验(包括基于python的线性代数运算、线性回归分析实验、聚类分析、因子分析和主成分分析),这里分享出来,注解标注的比较全,供大家参考。 使用Python语言开发完成以下运算。 1、已知有两个矩阵A和B,如下所示: ①求A+B、A-B;

    2023年04月08日
    浏览(68)
  • 线性代数Python计算:矩阵对角化

    线性变换 T T T 的矩阵 A ∈ P n × n boldsymbol{A}in P^{ntimes n} A ∈ P n × n 的对角化,即寻求对角阵 Λ boldsymbol{Lambda} Λ ,使得 A boldsymbol{A} A ~ Λ boldsymbol{Lambda} Λ ,需分几步走: (1)解方程 det ⁡ ( λ I − A ) = 0 det(lambdaboldsymbol{I}-boldsymbol{A})=0 det ( λ I − A ) = 0 ,得根 λ 1 , λ

    2024年02月08日
    浏览(47)
  • Python处理矩阵和线性代数问题

    如未作说明,下列方法均调用自 linalg 矩阵分解 cholesky , qr ,奇异值分解 svd 求特征值 eigvals ,共轭对称阵特征值 eigvalsh(a[, UPLO]) 特征值和 特征向量 eig ,共轭对称的特征值和向量 eigh(a[, UPLO]) 特征数字 范数 norm ,迹 trace 条件数 cond ,行列式 det ,符号 slogdet 通过SVD方法求秩

    2024年02月05日
    浏览(36)
  • 线性代数Python计算:线性方程组的最小二乘解

    给定ℝ上无解线性方程组 A x = b boldsymbol{Ax}=boldsymbol{b} Ax = b ,构造 A T A boldsymbol{A}^text{T}boldsymbol{A} A T A 及 A T b boldsymbol{A}^text{T}boldsymbol{b} A T b ,然后调用博文《线性方程组的通解》定义的mySolve函数,解方程组 A T A x = A T b boldsymbol{A}^text{T}boldsymbol{Ax}=boldsymbol{A}^text{T

    2023年04月08日
    浏览(59)
  • 【Python · PyTorch】线性代数 & 微积分

    本文采用Python及PyTorch版本如下: Python:3.9.0 PyTorch:2.0.1+cpu 本文为博主自用知识点提纲,无过于具体介绍,详细内容请参考其他文章。 线性代数是数学的一个分支,它的研究对象是向量、向量空间(线性空间)、线性变换及有限维的线性方程组。线性代数已被广泛地应用于

    2024年02月08日
    浏览(49)
  • python-np.linalg-线性代数

    np.linalg 是NumPy库中用于线性代数运算的子模块。 1. 矩阵和向量的乘法: np.dot() 2. 矩阵的逆: np.linalg.inv(A) 矩阵必须是方阵且可逆,否则会抛出LinAlgError异常。 3. 矩阵的转置: np.transpose(A) 4. 矩阵的行列式: np.linalg.det(A) 5. 矩阵的特征值和特征向量: np.linalg.eig() linalg模块中,

    2024年04月23日
    浏览(40)
  • Python在高等数学和线性代数中的应用

    Python数学实验与建模学习 目录 1. SymPy工具库 1.1 符号运算基础 1.2 用SymPy做符号函数画图  2. 高等数学的符号解 2.1 极限 2.2 导数  2.3 级数求和  2.4 泰勒展开  2.5 不定积分和定积分  2.6 代数方程  2.7 微分方程  3. 高等数学问题的数值解 3.1 一重积分 3.1.1 梯形计算 3.1.2 辛普森

    2024年01月25日
    浏览(53)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包