fork download 
  1. import numpy as np
  2.  
  3. # =========================
  4. # Input Handling for Matrix A Block
  5. # =========================
  6. # This block defines a function to handle user input for the number of variables and the coefficient matrix A.
  7. # It prompts the user for n, reads the matrix rows, validates the input length, and returns A as a NumPy array.
  8. # Exceptions are raised for invalid inputs to be handled in the calling code.
  9. def take_A_input():
  10.     try:
  11.         n = int(input("dimension of matrix: "))
  12.     except ValueError:
  13.         raise ValueError("Invalid input: Number of variables must be an integer.")
  14.  
  15.     print("\nEnter the coefficient matrix (A):")
  16.     A = []
  17.     for i in range(n):
  18.         try:
  19.             row = list(map(float, input(f"Row {i+1}: ").split()))
  20.             if len(row) != n:
  21.                 raise ValueError(f"Row {i+1} must have exactly {n} elements.")
  22.             A.append(row)
  23.         except ValueError as e:
  24.             raise ValueError(f"Invalid input in row {i+1}: {e}")
  25.  
  26.     return np.array(A, dtype=float), n  # Return A and n for later use in dimensions.
  27.  
  28. # =========================
  29. # Input Handling for Vector B Block
  30. # =========================
  31. # This block defines a function to handle user input for the constants vector B, given the dimension n.
  32. # It prompts for each element of B, validates them as floats, and returns B as a NumPy array.
  33. # Exceptions are raised for invalid inputs.
  34. def take_B_input(n):
  35.     print("\nEnter the constants vector (B):")
  36.     B = []
  37.     for i in range(n):
  38.         try:
  39.             val = float(input(f"B[{i+1}]: "))
  40.             B.append(val)
  41.         except ValueError:
  42.             raise ValueError(f"Invalid input for B[{i+1}]: Must be a float.")
  43.  
  44.     return np.array(B, dtype=float)
  45.  
  46. # =========================
  47. # LU Decomposition Block
  48. # =========================
  49. # This block performs LU decomposition on matrix A without pivoting.
  50. # It decomposes A into a lower triangular matrix L (with 1s on diagonal) and an upper triangular matrix U such that A = L * U.
  51. # Raises an error if a zero pivot is encountered during the process.
  52. def lu_decomposition(A):
  53.     n = len(A)
  54.     L = np.zeros((n, n))  # Initialize L as zero matrix; diagonals will be set to 1.
  55.     U = A.copy().astype(float)  # Copy A to U to modify it in place without affecting original A.
  56.  
  57.     try:
  58.         for i in range(n):
  59.             L[i, i] = 1  # Set diagonal of L to 1, as per standard LU decomposition.
  60.             for j in range(i+1, n):
  61.                 if U[i, i] == 0:
  62.                     raise ZeroDivisionError("Zero pivot encountered during LU decomposition")  # Check for zero pivot to avoid division by zero.
  63.                 factor = U[j, i] / U[i, i]  # Compute the multiplier factor for elimination.
  64.                 L[j, i] = factor  # Store the factor in L below the diagonal.
  65.                 U[j, i:] -= factor * U[i, i:]  # Eliminate entries below pivot in U and update the row.
  66.     except Exception as e:
  67.         raise RuntimeError(f"Error during LU decomposition: {e}")
  68.  
  69.     return L, U
  70.  
  71. # =========================
  72. # LU Solve Block
  73. # =========================
  74. # This block solves the linear system Ax = B using LU decomposition.
  75. # First, it performs forward substitution to solve Ly = B for y (assuming L diagonal is 1).
  76. # Then, backward substitution to solve Ux = y for x.
  77. # Returns the solution vector x.
  78. def lu_solve(A, B):
  79.     try:
  80.         L, U = lu_decomposition(A)  # Decompose A into L and U.
  81.     except Exception as e:
  82.         raise RuntimeError(f"Error in decomposition during solve: {e}")
  83.  
  84.     n = len(B)
  85.  
  86.     # Forward substitution: Solve Ly = B
  87.     y = np.zeros_like(B, dtype=float)  # Initialize y vector.
  88.     try:
  89.         for i in range(n):
  90.             y[i] = B[i] - np.dot(L[i, :i], y[:i])  # Compute y[i] using previously solved y values.
  91.     except Exception as e:
  92.         raise RuntimeError(f"Error during forward substitution: {e}")
  93.  
  94.     # Backward substitution: Solve Ux = y
  95.     x = np.zeros_like(y, dtype=float)  # Initialize x vector.
  96.     try:
  97.         for i in range(n-1, -1, -1):
  98.             x[i] = (y[i] - np.dot(U[i, i+1:], x[i+1:])) / U[i, i]  # Compute x[i] using previously solved x values; divide by pivot.
  99.     except ZeroDivisionError:
  100.         raise ZeroDivisionError("Division by zero during backward substitution")
  101.     except Exception as e:
  102.         raise RuntimeError(f"Error during backward substitution: {e}")
  103.  
  104.     return x
  105.  
  106. # =========================
  107. # Matrix Inverse Block
  108. # =========================
  109. # This block computes the inverse of matrix A using LU decomposition.
  110. # It solves Ax = e_i for each standard basis vector e_i to find each column of the inverse.
  111. # Returns the inverse matrix.
  112. def inverse(A):
  113.     try:
  114.         n = len(A)
  115.         inv_A = np.zeros_like(A, dtype=float)  # Initialize inverse matrix with zeros.
  116.         for i in range(n):
  117.             e = np.zeros(n)  # Create standard basis vector e_i with 1 at position i.
  118.             e[i] = 1
  119.             inv_A[:, i] = lu_solve(A, e)  # Solve for the i-th column of the inverse.
  120.     except Exception as e:
  121.         raise RuntimeError(f"Error during inverse calculation: {e}")
  122.     return inv_A
  123.  
  124. # =========================
  125. # Determinant Calculation Block
  126. # =========================
  127. # This block computes the determinant of matrix A using LU decomposition.
  128. # The determinant is the product of the diagonal elements of U (since det(A) = det(L) * det(U) and det(L) = 1).
  129. def determinant(A):
  130.     try:
  131.         _, U = lu_decomposition(A)  # Only need U for determinant.
  132.         det = np.prod(np.diag(U))  # Product of diagonal elements of U gives the determinant.
  133.     except Exception as e:
  134.         raise RuntimeError(f"Error during determinant calculation: {e}")
  135.     return det
  136.  
  137. # =========================
  138. # Main Execution Block
  139. # =========================
  140. # This is the entry point of the program.
  141. # It first handles input for A, performs LU decomposition, computes determinant and inverse, and prints them.
  142. # Then, it prompts the user if they want to solve for roots; if yes, takes B input and solves Ax = B.
  143. # All operations are wrapped in try-except to catch and display errors gracefully.
  144. def main():
  145.     print("=== LU Decomposition Solver ===\n")
  146.     try:
  147.         A, n = take_A_input()
  148.     except Exception as e:
  149.         print("Input Error for A:", e)
  150.         return
  151.  
  152.     print("\n--- LU Decomposition ---")
  153.     try:
  154.         L, U = lu_decomposition(A)
  155.         print("Lower Triangular Matrix L:\n", L)
  156.         print("Upper Triangular Matrix U:\n", U)
  157.     except Exception as e:
  158.         print("Error:", e)
  159.  
  160.     print("\n--- Determinant of A ---")
  161.     try:
  162.         detA = determinant(A)
  163.         print(detA)
  164.     except Exception as e:
  165.         print("Error:", e)
  166.  
  167.     print("\n--- Inverse of A ---")
  168.     try:
  169.         invA = inverse(A)
  170.         print(invA)
  171.     except Exception as e:
  172.         print("Error:", e)
  173.  
  174.     # Prompt user for optional solving of roots
  175.     solve_roots = input("\nDo you want to know the roots (yes/no): ").strip().lower()
  176.     if solve_roots == 'yes':
  177.         try:
  178.             B = take_B_input(n)
  179.         except Exception as e:
  180.             print("Input Error for B:", e)
  181.             return
  182.  
  183.         print("\n--- Solve using LU Decomposition ---")
  184.         try:
  185.             x_lu = lu_solve(A, B)
  186.             print("Solution:", x_lu)
  187.         except Exception as e:
  188.             print("Error:", e)
  189.  
  190. if __name__ == "__main__":
  191.     main()
Success #stdin #stdout 0.8s 41396KB
comments (?)
stdin 
Standard input is empty
stdout 
=== LU Decomposition Solver ===

dimension of matrix: Input Error for A: EOF when reading a line
https://ideone.com/wz8oYO
language:
Python 3 (python  3.12)
created:
16 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!