线性方程组数值迭代法

import numpy as np
from typing import Tuple, Optional

class IterativeLinearSolvers:
    """
    Collection of iterative methods for solving linear systems Ax = b
    """
    
    @staticmethod
    def jacobi(A: np.ndarray, b: np.ndarray, x0: Optional[np.ndarray] = None,
               tol: float = 1e-10, max_iter: int = 1000) -> Tuple[np.ndarray, int, list]:
        """
        Jacobi iteration method
        x^(k+1) = D^(-1)(b - (L+U)x^(k))
        """
        n = len(b)
        x = np.zeros(n) if x0 is None else x0.copy()
        x_new = np.zeros(n)
        residuals = []
        
        D = np.diag(np.diag(A))
        R = A - D
        
        for k in range(max_iter):
            x_new = np.linalg.solve(D, b - R @ x)
            residual = np.linalg.norm(x_new - x)
            residuals.append(residual)
            
            if residual < tol:
                return x_new, k + 1, residuals
            
            x = x_new.copy()
        
        print(f"Warning: Jacobi did not converge in {max_iter} iterations")
        return x, max_iter, residuals
    
    @staticmethod
    def gauss_seidel(A: np.ndarray, b: np.ndarray, x0: Optional[np.ndarray] = None,
                     tol: float = 1e-10, max_iter: int = 1000) -> Tuple[np.ndarray, int, list]:
        """
        Gauss-Seidel iteration method
        x^(k+1) = (D+L)^(-1)(b - Ux^(k))
        """
        n = len(b)
        x = np.zeros(n) if x0 is None else x0.copy()
        x_new = x.copy()
        residuals = []
        
        for k in range(max_iter):
            for i in range(n):
                sum1 = np.dot(A[i, :i], x_new[:i])
                sum2 = np.dot(A[i, i+1:], x[i+1:])
                x_new[i] = (b[i] - sum1 - sum2) / A[i, i]
            
            residual = np.linalg.norm(x_new - x)
            residuals.append(residual)
            
            if residual < tol:
                return x_new, k + 1, residuals
            
            x = x_new.copy()
        
        print(f"Warning: Gauss-Seidel did not converge in {max_iter} iterations")
        return x, max_iter, residuals
    
    @staticmethod
    def sor(A: np.ndarray, b: np.ndarray, omega: float = 1.5, 
            x0: Optional[np.ndarray] = None, tol: float = 1e-10, 
            max_iter: int = 1000) -> Tuple[np.ndarray, int, list]:
        """
        Successive Over-Relaxation (SOR) method
        x^(k+1) = (D + ωL)^(-1)(ωb - [ωU + (ω-1)D]x^(k))
        omega: relaxation parameter (0 < omega < 2)
        """
        n = len(b)
        x = np.zeros(n) if x0 is None else x0.copy()
        x_new = x.copy()
        residuals = []
        
        for k in range(max_iter):
            for i in range(n):
                sum1 = np.dot(A[i, :i], x_new[:i])
                sum2 = np.dot(A[i, i+1:], x[i+1:])
                x_gs = (b[i] - sum1 - sum2) / A[i, i]
                x_new[i] = omega * x_gs + (1 - omega) * x[i]
            
            residual = np.linalg.norm(x_new - x)
            residuals.append(residual)
            
            if residual < tol:
                return x_new, k + 1, residuals
            
            x = x_new.copy()
        
        print(f"Warning: SOR did not converge in {max_iter} iterations")
        return x, max_iter, residuals
    
    @staticmethod
    def conjugate_gradient(A: np.ndarray, b: np.ndarray, x0: Optional[np.ndarray] = None,
                          tol: float = 1e-10, max_iter: int = 1000) -> Tuple[np.ndarray, int, list]:
        """
        Conjugate Gradient method (for symmetric positive definite matrices)
        """
        n = len(b)
        x = np.zeros(n) if x0 is None else x0.copy()
        residuals = []
        
        r = b - A @ x
        p = r.copy()
        rs_old = r @ r
        
        for k in range(max_iter):
            Ap = A @ p
            alpha = rs_old / (p @ Ap)
            x = x + alpha * p
            r = r - alpha * Ap
            rs_new = r @ r
            
            residual = np.sqrt(rs_new)
            residuals.append(residual)
            
            if residual < tol:
                return x, k + 1, residuals
            
            beta = rs_new / rs_old
            p = r + beta * p
            rs_old = rs_new
        
        print(f"Warning: CG did not converge in {max_iter} iterations")
        return x, max_iter, residuals
    
    @staticmethod
    def steepest_descent(A: np.ndarray, b: np.ndarray, x0: Optional[np.ndarray] = None,
                        tol: float = 1e-10, max_iter: int = 1000) -> Tuple[np.ndarray, int, list]:
        """
        Steepest Descent method (for symmetric positive definite matrices)
        """
        n = len(b)
        x = np.zeros(n) if x0 is None else x0.copy()
        residuals = []
        
        for k in range(max_iter):
            r = b - A @ x
            residual = np.linalg.norm(r)
            residuals.append(residual)
            
            if residual < tol:
                return x, k + 1, residuals
            
            Ar = A @ r
            alpha = (r @ r) / (r @ Ar)
            x = x + alpha * r
        
        print(f"Warning: Steepest Descent did not converge in {max_iter} iterations")
        return x, max_iter, residuals
    
    @staticmethod
    def gmres(A: np.ndarray, b: np.ndarray, x0: Optional[np.ndarray] = None,
              tol: float = 1e-10, max_iter: int = 1000, restart: int = 50) -> Tuple[np.ndarray, int, list]:
        """
        Generalized Minimal Residual (GMRES) method with restart
        """
        n = len(b)
        x = np.zeros(n) if x0 is None else x0.copy()
        residuals = []
        total_iter = 0
        
        for _ in range(max_iter // restart):
            r = b - A @ x
            beta = np.linalg.norm(r)
            residuals.append(beta)
            
            if beta < tol:
                return x, total_iter, residuals
            
            V = np.zeros((n, restart + 1))
            H = np.zeros((restart + 1, restart))
            V[:, 0] = r / beta
            
            for j in range(restart):
                total_iter += 1
                w = A @ V[:, j]
                
                for i in range(j + 1):
                    H[i, j] = w @ V[:, i]
                    w = w - H[i, j] * V[:, i]
                
                H[j + 1, j] = np.linalg.norm(w)
                
                if H[j + 1, j] > 1e-12:
                    V[:, j + 1] = w / H[j + 1, j]
                
                e1 = np.zeros(j + 2)
                e1[0] = beta
                y, _, _, _ = np.linalg.lstsq(H[:j+2, :j+1], e1, rcond=None)
                
                residual = np.linalg.norm(e1 - H[:j+2, :j+1] @ y)
                residuals.append(residual)
                
                if residual < tol:
                    x = x + V[:, :j+1] @ y
                    return x, total_iter, residuals
            
            e1 = np.zeros(restart + 1)
            e1[0] = beta
            y, _, _, _ = np.linalg.lstsq(H, e1, rcond=None)
            x = x + V[:, :restart] @ y
        
        print(f"Warning: GMRES did not converge in {max_iter} iterations")
        return x, total_iter, residuals
    
    @staticmethod
    def richardson(A: np.ndarray, b: np.ndarray, alpha: float = 0.1,
                   x0: Optional[np.ndarray] = None, tol: float = 1e-10,
                   max_iter: int = 1000) -> Tuple[np.ndarray, int, list]:
        """
        Richardson iteration method
        x^(k+1) = x^(k) + α(b - Ax^(k))
        alpha: relaxation parameter
        """
        n = len(b)
        x = np.zeros(n) if x0 is None else x0.copy()
        residuals = []
        
        for k in range(max_iter):
            r = b - A @ x
            residual = np.linalg.norm(r)
            residuals.append(residual)
            
            if residual < tol:
                return x, k + 1, residuals
            
            x = x + alpha * r
        
        print(f"Warning: Richardson did not converge in {max_iter} iterations")
        return x, max_iter, residuals
    
    @staticmethod
    def bicgstab(A: np.ndarray, b: np.ndarray, x0: Optional[np.ndarray] = None,
                 tol: float = 1e-10, max_iter: int = 1000) -> Tuple[np.ndarray, int, list]:
        """
        BiConjugate Gradient Stabilized (BiCGSTAB) method
        """
        n = len(b)
        x = np.zeros(n) if x0 is None else x0.copy()
        residuals = []
        
        r = b - A @ x
        r_hat = r.copy()
        rho = alpha = omega = 1.0
        v = p = np.zeros(n)
        
        for k in range(max_iter):
            rho_prev = rho
            rho = r_hat @ r
            
            if abs(rho) < 1e-12:
                print("BiCGSTAB breakdown: rho too small")
                return x, k, residuals
            
            if k == 0:
                p = r.copy()
            else:
                beta = (rho / rho_prev) * (alpha / omega)
                p = r + beta * (p - omega * v)
            
            v = A @ p
            alpha = rho / (r_hat @ v)
            s = r - alpha * v
            
            residual = np.linalg.norm(s)
            residuals.append(residual)
            
            if residual < tol:
                x = x + alpha * p
                return x, k + 1, residuals
            
            t = A @ s
            omega = (t @ s) / (t @ t)
            x = x + alpha * p + omega * s
            r = s - omega * t
        
        print(f"Warning: BiCGSTAB did not converge in {max_iter} iterations")
        return x, max_iter, residuals


# Example usage and testing
if __name__ == "__main__":
    # Create a test system: Ax = b
    n = 100
    np.random.seed(42)
    
    # Create a diagonally dominant matrix (ensures convergence)
    A = np.random.rand(n, n)
    A = A + A.T  # Make symmetric
    A = A + n * np.eye(n)  # Make diagonally dominant
    
    # True solution
    x_true = np.random.rand(n)
    b = A @ x_true
    
    solver = IterativeLinearSolvers()
    
    print("Testing Iterative Methods for Linear Systems")
    print("=" * 60)
    print(f"System size: {n}x{n}\n")
    
    methods = [
        ("Jacobi", lambda: solver.jacobi(A, b)),
        ("Gauss-Seidel", lambda: solver.gauss_seidel(A, b)),
        ("SOR (ω=1.5)", lambda: solver.sor(A, b, omega=1.5)),
        ("Conjugate Gradient", lambda: solver.conjugate_gradient(A, b)),
        ("Steepest Descent", lambda: solver.steepest_descent(A, b)),
        ("GMRES", lambda: solver.gmres(A, b, restart=20)),
        ("Richardson (α=0.01)", lambda: solver.richardson(A, b, alpha=0.01)),
        ("BiCGSTAB", lambda: solver.bicgstab(A, b))
    ]
    
    for name, method in methods:
        x_sol, iters, residuals = method()
        error = np.linalg.norm(x_sol - x_true)
        print(f"{name:25s} | Iterations: {iters:4d} | Error: {error:.2e} | Final Residual: {residuals[-1]:.2e}")
    
    print("\n" + "=" * 60)
    print("Note: For non-symmetric systems, use GMRES or BiCGSTAB")
    print("For symmetric positive definite systems, CG is optimal")