목록으로

Chapter 7: 선형시스템

진행률: 7 / 8
예상 시간: 25분난이도: 중급

선형시스템

학습 목표

  • 연립방정식의 다양한 해법 마스터
  • LU 분해와 효율적인 계산
  • 반복법과 수치적 안정성

선형시스템 Ax = b

유일해

det(A) ≠ 0
rank(A) = n

무수히 많은 해

rank(A) = rank([A|b]) < n

해 없음

rank(A) < rank([A|b])

가우스 소거법

1

전진 소거

상삼각 행렬로 변환

2

후진 대입

아래에서 위로 해 계산

계산 복잡도: O(n³)

LU 분해

PA = LU

  • • P: 순열 행렬 (피벗팅)
  • • L: 하삼각 행렬
  • • U: 상삼각 행렬

장점

  • • 여러 b에 대해 효율적
  • • 행렬식 계산 용이
  • • 역행렬 계산 효율적

해법 과정

  1. 1. Ly = Pb (전진 대입)
  2. 2. Ux = y (후진 대입)

반복법

야코비 방법

x_i^(k+1) = (b_i - Σ a_ij x_j^(k)) / a_ii

가우스-자이델

새 값을 즉시 사용
더 빠른 수렴

SOR 방법

과이완 계수 ω 사용
최적 수렴

선형시스템 솔버 구현

import numpy as np
from scipy.linalg import lu, solve_triangular

# 가우스 소거법
def gaussian_elimination(A, b):
    """가우스 소거법으로 Ax = b 해결"""
    n = len(A)
    Ab = np.column_stack([A.copy(), b.copy()])
    
    # 전진 소거
    for i in range(n):
        # 피벗팅
        max_row = i + np.argmax(np.abs(Ab[i:, i]))
        Ab[[i, max_row]] = Ab[[max_row, i]]
        
        # 소거
        for j in range(i+1, n):
            factor = Ab[j, i] / Ab[i, i]
            Ab[j, i:] -= factor * Ab[i, i:]
    
    # 후진 대입
    x = np.zeros(n)
    for i in range(n-1, -1, -1):
        x[i] = (Ab[i, -1] - np.dot(Ab[i, i+1:n], x[i+1:])) / Ab[i, i]
    
    return x

# LU 분해 솔버
def lu_solver(A, b):
    """LU 분해를 이용한 선형시스템 해결"""
    P, L, U = lu(A)
    
    # Ly = Pb
    y = solve_triangular(L, P @ b, lower=True)
    
    # Ux = y
    x = solve_triangular(U, y)
    
    return x

# 야코비 반복법
def jacobi_method(A, b, x0=None, tol=1e-10, max_iter=1000):
    """야코비 반복법"""
    n = len(A)
    x = x0 if x0 is not None else np.zeros(n)
    
    for iteration in range(max_iter):
        x_new = np.zeros(n)
        
        for i in range(n):
            sum_ax = sum(A[i, j] * x[j] for j in range(n) if j != i)
            x_new[i] = (b[i] - sum_ax) / A[i, i]
        
        # 수렴 확인
        if np.linalg.norm(x_new - x) < tol:
            print(f"수렴 (반복: {iteration+1})")
            return x_new
        
        x = x_new
    
    print("최대 반복 횟수 도달")
    return x

# 가우스-자이델 방법
def gauss_seidel(A, b, x0=None, tol=1e-10, max_iter=1000):
    """가우스-자이델 반복법"""
    n = len(A)
    x = x0 if x0 is not None else np.zeros(n)
    
    for iteration in range(max_iter):
        x_old = x.copy()
        
        for i in range(n):
            sum1 = sum(A[i, j] * x[j] for j in range(i))
            sum2 = sum(A[i, j] * x_old[j] for j in range(i+1, n))
            x[i] = (b[i] - sum1 - sum2) / A[i, i]
        
        # 수렴 확인
        if np.linalg.norm(x - x_old) < tol:
            print(f"수렴 (반복: {iteration+1})")
            return x
    
    print("최대 반복 횟수 도달")
    return x

# 조건수 분석
def condition_analysis(A):
    """행렬의 조건수 분석"""
    cond = np.linalg.cond(A)
    
    print(f"조건수: {cond:.2e}")
    
    if cond < 10:
        print("매우 안정적")
    elif cond < 100:
        print("안정적")
    elif cond < 1000:
        print("주의 필요")
    else:
        print("불안정 (ill-conditioned)")
    
    return cond

# 예제
A = np.array([[10, -1, 2], 
              [-1, 11, -1], 
              [2, -1, 10]])
b = np.array([6, 25, -11])

print("가우스 소거법:", gaussian_elimination(A, b))
print("LU 분해:", lu_solver(A, b))
print("야코비 방법:", jacobi_method(A, b))
print("가우스-자이델:", gauss_seidel(A, b))
print("\n조건수 분석:")
condition_analysis(A)