时间戳

首页 > 解决方案 > Python - 用 LU 分解求解线性方程组的数值问题

python - Python - 用 LU 分解求解线性方程组的数值问题

问题描述

信息

操作系统：Manjaro Linux

开发环境：Spyder4 (Python 3.7)

库：Numpy

问题

你好，

我写了一些函数来解决这三种方法的线性方程组：

  1. LU分解
  2. 雅可比法
  3. 高斯-赛德尔法

该程序运行完美。然而，LU分解的结果让我很困扰。例如，如果我的矩阵 A 和向量 b 是

A = np.array([[10., -1., 2., 0.],
          [-1., 11., -1., 3.],
          [2., -1., 10., -1.],
          [0., 3., -1., 8.]])
b = np.array([6., 25., -11., 15.])

然后我的结果是

  1. 陆：[0.95034091 1.67840909 -0.9125 1.875]
  2. 雅可比：[ 1. 2. -1. 1.]
  3. 高斯-赛德尔：[ 1. 2. -1. 1.]

正如你所看到的，LU 给了我一个稍微不同的结果。是否存在舍入或截断错误的问题？任何建议或帮助将不胜感激。

干杯!

编辑1：

我修复了这个问题，它只是 LU 分解期间代码中的一个错误。感谢大家的反馈。


共享代码如下：

def lu_decomposition(matrix_in, b_in, n):

lower = [[0 for x in range(n)]
            for y in range(n)];
upper = [[0 for x in range(n)]
            for y in range(n)];

# Doolittle's Method
# Decomposing matrix into upper and lower matrices
for i in range(n):
    # Upper triangle
    for k in range(i, n):
        # Sigma from j to i for each row k of U
        sum = 0
        for j in range(i):
            sum += lower[i][j] * upper[j][k]
        # Evaluate U for row k
        upper[i][k] = matrix_in[i][k] - sum
        
    # Lower Triangle
    for k in range(i,n):
        if(i == k): # Entry of a diagonal element
            lower[i][i] = 1 
        else:
            # Sigma from j to i for each column k of L
            sum = 0
            for j in range(i):
                sum += (lower[k][j] * upper[j][i])
            # Evaluate L for column k
            lower[k][i] = int( (matrix_in[k][i] - sum)/ upper[i][i])
            
# Perform forward substitution Ly=b
y = [0 for x in range(n)]
lower_inv = np.linalg.inv(lower)
y = np.dot(lower_inv, b_in)

# Perform back substitution Ux=y
x_sol = [0 for x in range(n)]
upper_inv = np.linalg.inv(upper)
x_sol = np.dot(upper_inv, y)

# printing results
# setw is for displaying nicely 
print("Lower Triangular\t\tUpper Triangular"); 

# Displaying the result : 
for i in range(n): 
      
    # Lower 
    for j in range(n): 
        print(lower[i][j], end = "\t");  
    print("", end = "\t"); 

    # Upper 
    for j in range(n): 
        print(upper[i][j], end = "\t"); 
    print(""); 
            
print("Here's the solution for your matrix: ")
print(x_sol)

标签： pythonpython-3.xnumpynumerical-methods

解决方案

您的代码中有很多事情需要解决：

  • Python 不需要；要结束一行，python 知道行何时结束。Python 只使用 ; 将两行代码放在同一物理行上。像“x=7;print(x)”。
  • 如果您不在列表推导中使用变量，则习惯上使用 _。例如，[0 for _ in range(10)]。当然，python 有更好的方法“[0]*10”。
  • 我注意到你导入了 numpy，但你没有使用 numpy 数组的操作。这将极大地改进您的代码（并使其更快）并使其更易于阅读。例如，您可以一次写入和读取整行或整列。“matrix[0,:]”（第一行），“matrix[:,0]”（第一列）。
  • 不需要在 python 中包含对象的长度，因为 python 会自动存储它的对象的长度，它总是可以使用 len() 内置函数来检索。
  • 对您的代码来说并不是很重要（因为您可能会这样做是为了练习），但您可能已经知道，lu 分解已经存在，例如 scipy.linalg.lu()。事实上，快速检查 linalg.lu(A,True) 会发现 L 和 U 中都有错误。
  • 看到你手动生成 L 和 U 矩阵然后继续在 L 和 U 上使用 np.linalg.inv() 真的很奇怪。如果你愿意使用 np.linalg.inv() 函数，那么你的答案是一行，“np.linalg.inv(A) @ b”。通常，人们更容易找到 L 和 U 来手动求解 X。要找到 L 和 U 然后使用 numpy 的反函数，有点达不到目的。
  • 尽管有时它会有所帮助，但 python 并不要求您在创建对象之前在内存中分配空间。Python 自动管理您的内存创建和删除（无需创建空的零列表）。
  • np.dot() 只是在 python 中访问“@”运算符的手动方式。

一些例子：

import numpy as np
from scipy import linalg

def lu_dec_solver_v1(A,b):
    return np.linalg.inv(A) @ b

def lu_dec_solver_v2(A,b):
    L, U  = linalg.lu(A,True)
    L_inv = np.linalg.inv(L)
    U_inv = np.linalg.inv(U)
    return U_inv @ (L_inv @ b)

def lu_dec_solver_v3(A,b):
    U = A.copy()
    L = np.identity(len(A))
    for n in range(0,len(A)-1):
        for m in range(n+1,len(A)):
            L[m,n]  = U[m,n]/U[n,n]
            U[m,:] += -L[m,n]*U[n,:]
    L_inv = np.linalg.inv(L)
    U_inv = np.linalg.inv(U)
    return U_inv @ (L_inv @ b)

推荐阅读