Gaussian Elimination for Solving Linear Systems

Task: Implement the Gaussian Elimination Method

Your task is to implement the Gaussian Elimination method, which transforms a system of linear equations into an upper triangular matrix. This method can then be used to solve for the variables using backward substitution.

Write a function gaussian_elimination(A, b) that performs Gaussian Elimination with partial pivoting to solve the system (Ax = b).

The function should return the solution vector (x).

Example:

Input:

A = np.array([[2,8,4], [2,5,1], [4,10,-1]], dtype=float)
b = np.array([2,5,1], dtype=float)

print(gaussian_elimination(A, b))

Output:

[11.0, -4.0, 3.0]

Reasoning:

The Gaussian Elimination method transforms the system of equations into an upper triangular matrix and then uses backward substitution to solve for the variables.

import numpy as np

def gaussian_elimination(A, b):
    """
    Solves the system Ax = b using Gaussian Elimination with partial pivoting.
    
    :param A: Coefficient matrix
    :param b: Right-hand side vector
    :return: Solution vector x
    """
    # Create copies to avoid modifying the original matrices
    A = A.copy()
    b = b.copy()
    
    n = len(b)
    
    # Augment the matrix A with vector b
    augmented = np.column_stack((A, b))
    
    # Forward elimination with partial pivoting
    for i in range(n):
        # Partial pivoting: find the row with the largest absolute value in the current column
        max_row = i
        max_val = abs(augmented[i, i])
        
        for k in range(i + 1, n):
            if abs(augmented[k, i]) > max_val:
                max_val = abs(augmented[k, i])
                max_row = k
        
        # Swap the rows if necessary
        if max_row != i:
            augmented[[i, max_row]] = augmented[[max_row, i]]
        
        # Check if the matrix is singular
        if abs(augmented[i, i]) < 1e-10:
            raise ValueError("Matrix is singular or nearly singular")
        
        # Eliminate entries below the pivot
        for j in range(i + 1, n):
            factor = augmented[j, i] / augmented[i, i]
            augmented[j, i:] -= factor * augmented[i, i:]
    
    # Back substitution to solve for x
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = augmented[i, -1]
        
        for j in range(i + 1, n):
            x[i] -= augmented[i, j] * x[j]
            
        x[i] /= augmented[i, i]
    
    return x.tolist()

Test Results3/3

官方题解
 

import numpy as np

def partial_pivoting(A_aug, row_num, col_num):
    rows, cols = A_aug.shape
    max_row = row_num
    max_val = abs(A_aug[row_num, col_num])
    for i in range(row_num, rows):
        current_val = abs(A_aug[i, col_num])
        if current_val > max_val:
            max_val = current_val
            max_row = i
    if max_row != row_num:
        A_aug[[row_num, max_row]] = A_aug[[max_row, row_num]]
    return A_aug

def gaussian_elimination(A, b):
    rows, cols = A.shape
    A_aug = np.hstack((A, b.reshape(-1, 1)))

    for i in range(rows-1):
        A_aug = partial_pivoting(A_aug, i, i)
        for j in range(i+1, rows):
            A_aug[j, i:] -= (A_aug[j, i] / A_aug[i, i]) * A_aug[i, i:]

    x = np.zeros_like(b, dtype=float)
    for i in range(rows-1, -1, -1):
        x[i] = (A_aug[i, -1] - np.dot(A_aug[i, i+1:cols], x[i+1:])) / A_aug[i, i]
    return x