수치해석
재료공학도를 위한 수치해석
- 1. 학습 목표
- 2. 가우스 소거법
- 3. 손으로 계산하는 예제
- 4. 부분 피벗팅
- 5. Python 구현
- 6. Gauss–Jordan 소거법과의 차이
- 7. 계산량과 주의 사항
- 8. 정리
- 9. 연습 문제
1. 학습 목표
이 자료에서 $[\boldsymbol A]$, $[\boldsymbol x]$는 행렬·벡터의 성분 배열을 뜻하며, 배열의 행렬곱은 $[\boldsymbol A][\boldsymbol x]$처럼 붙여 쓴다. 물리적 벡터·텐서의 단일수축은 $\cdot$로 표시한다.
이번 강의가 끝나면 다음을 할 수 있어야 한다.
- 연립방정식을 확대행렬(augmented matrix)로 나타낼 수 있다.
- 기본 행 연산을 이용하여 확대행렬을 상삼각(upper triangular) 형태로 바꿀 수 있다.
- 소거 계수를 계산하고 전진 소거(forward elimination)를 수행할 수 있다.
- 후진 대입(back substitution)으로 미지수를 구할 수 있다.
- 부분 피벗팅(partial pivoting)이 필요한 이유를 설명할 수 있다.
- 가우스 소거법(Gaussian elimination)을 Python으로 구현하고 잔차로 결과를 확인할 수 있다.
2. 가우스 소거법
가우스 소거법(Gaussian elimination)은 선형 연립방정식
\[[\boldsymbol A][\boldsymbol x]=[\boldsymbol b]\]를 기본 행 연산으로 변형하여 푸는 직접해법이다. 계산 과정은 크게 두 단계로 나뉜다.
- 전진 소거(forward elimination): 계수행렬을 상삼각행렬로 만든다.
- 후진 대입(back substitution): 마지막 식부터 거꾸로 미지수를 구한다.
2.1. 확대행렬
다음 연립방정식을 생각하자.
\[\begin{aligned} a_{11}x_1+a_{12}x_2+a_{13}x_3&=b_1,\\ a_{21}x_1+a_{22}x_2+a_{23}x_3&=b_2,\\ a_{31}x_1+a_{32}x_2+a_{33}x_3&=b_3. \end{aligned}\]계수행렬과 우변 벡터를 결합한 행렬을 확대행렬(augmented matrix)이라고 한다.
\[[\boldsymbol A\mid\boldsymbol b] = \left[ \begin{array}{ccc|c} a_{11}&a_{12}&a_{13}&b_1\\ a_{21}&a_{22}&a_{23}&b_2\\ a_{31}&a_{32}&a_{33}&b_3 \end{array} \right]\]**세로선은 계수행렬과 우변 벡터를 구분하기 위한 표시이며 행렬의 원소가 아니다.
2.2. 기본 행 연산
연립방정식의 해를 바꾸지 않는 기본 행 연산은 세 가지이다.
-
두 행(row)을 서로 교환한다.
\[R_i\leftrightarrow R_j\] -
한 행에 0이 아닌 상수를 곱한다.
\[R_i\leftarrow cR_i, \qquad c\ne0\] -
한 행에 다른 행의 배수를 더한다.
\[R_i\leftarrow R_i+cR_j\]
전진 소거에서는 주로 세 번째 연산을 다음 형태로 사용한다.
\[R_i\leftarrow R_i-m_{ik}R_k, \qquad m_{ik}=\frac{a_{ik}}{a_{kk}}\]여기서 $a_{kk}$는 현재 계산의 기준이 되는 피벗(pivot)이고, $m_{ik}$는 소거 계수(multiplier)이다.
2.3. 전체 알고리듬
$n\times n$ 연립방정식에 대한 과정은 다음과 같다.
- 확대행렬 $[\boldsymbol A\mid\boldsymbol b]$를 만든다.
- 첫 번째 열(column)부터 피벗을 선택한다.
- 필요하면 현재 행과 아래 행을 교환한다.
- 피벗 아래의 원소를 0으로 만든다.
- 다음 열로 이동하여 같은 작업을 반복한다.
- 상삼각 연립방정식이 만들어지면 후진 대입한다.
전진 소거가 끝난 계수행렬은 다음 형태이다.
\[[\boldsymbol U]= \begin{bmatrix} u_{11}&u_{12}&u_{13}\\ 0&u_{22}&u_{23}\\ 0&0&u_{33} \end{bmatrix}\]3. 손으로 계산하는 예제
다음 연립방정식을 풀어보자.
\[\begin{aligned} x+y+z&=6,\\ 2x+3y+z&=11,\\ -x+2y+3z&=12. \end{aligned}\]이 연립방정식의 해는 아래에 이어 계산하여 확인한다.
3.1. 확대행렬 만들기
\[\left[ \begin{array}{ccc|c} 1&1&1&6\\ 2&3&1&11\\ -1&2&3&12 \end{array} \right]\]3.2. 첫 번째 열 소거
첫 번째 피벗은 $a_{11}=1$이다. 두 번째 행에서 첫 번째 원소를 없애기 위한 소거 계수는
\[m_{21}=\frac{a_{21}}{a_{11}}=\frac{2}{1}=\textcolor{red}{2}\]이다. 따라서
\[R_2\leftarrow R_2-\textcolor{red}{2}R_1\]을 수행한다. 세 번째 행의 소거 계수는
\[m_{31}=\frac{a_{31}}{a_{11}}=\frac{-1}{1}=\textcolor{blue}{-1}\]이므로
\[R_3\leftarrow R_3-\textcolor{blue}{(-1)}R_1=R_3+R_1\]을 수행한다. 결과는
\[\left[ \begin{array}{ccc|c} 1&1&1&6\\ 0&1&-1&-1\\ 0&3&4&18 \end{array} \right]\]이다.
3.3. 두 번째 열 소거
두 번째 피벗은 $a_{22}=1$이고 소거 계수는
\[m_{32}=\frac{a_{32}}{a_{22}}=\frac{3}{1}=\textcolor{green}{3}\]이다. 따라서
\[R_3\leftarrow R_3-\textcolor{green}{3}R_2\]를 수행하면
\[\left[ \begin{array}{ccc|c} 1&1&1&6\\ 0&1&-1&-1\\ 0&0&7&21 \end{array} \right]\]이 된다.
3.4. 후진 대입 (back substitution)
마지막 행부터 미지수를 구한다.
\[7z=21 \quad\Longrightarrow\quad z=3\]두 번째 행에 $z=3$을 대입하면
\[y-z=-1 \quad\Longrightarrow\quad y=2\]첫 번째 행에 $y=2$, $z=3$을 대입하면
\[x+y+z=6 \quad\Longrightarrow\quad x=1\]이다. 따라서 최종 해는
\[\boxed{(x,y,z)=(1,2,3)}\]이다.
일반적인 상삼각 연립방정식의 후진 대입은 다음 식으로 나타낼 수 있다.
\[x_i= \frac{ b_i-\displaystyle\sum_{j=i+1}^{n}u_{ij}x_j }{u_{ii}}, \qquad i=n,n-1,\ldots,1\]4. 부분 피벗팅
피벗이 0이면 해당 피벗으로 나눌 수 없다. 다음 연립방정식을 살펴보자.
\[\begin{aligned} 0x+2y&=4,\\ x+y&=3. \end{aligned}\]첫 번째 피벗은 0이지만 두 행을 교환하면 계산할 수 있다.
\[\left[ \begin{array}{cc|c} 0&2&4\\ 1&1&3 \end{array} \right] \quad\xrightarrow{R_1\leftrightarrow R_2}\quad \left[ \begin{array}{cc|c} 1&1&3\\ 0&2&4 \end{array} \right]\]또한 피벗의 절댓값이 매우 작으면 나눗셈 과정에서 반올림 오차가 크게 증폭될 수 있다. 부분 피벗팅(partial pivoting)은 현재 열에서 피벗 후보와 그 아래 원소 중 절댓값이 가장 큰 원소를 찾아 현재 피벗 행과 교환한다.
열 $k$에서 피벗 행은 다음과 같이 선택한다.
\[p=k+\mathop{\mathrm{argmax}}_{i=k,\ldots,n-1}|a_{ik}|\]부분 피벗팅은 모든 행을 단순히 크기순으로 정렬하는 작업이 아니다. 소거가 한 단계 끝날 때마다 아직 처리하지 않은 부분행렬에서 새 피벗을 선택한다.
5. Python 구현
다음 구현은 입력 배열을 복사한 뒤 부분 피벗팅, 전진 소거와 후진 대입을 수행한다.
import numpy as np
def gaussian_elimination(A, b, tol=1e-12):
A = np.asarray(A, dtype=float).copy()
b = np.asarray(b, dtype=float).copy()
if A.ndim != 2 or A.shape[0] != A.shape[1]:
raise ValueError("A는 정사각행렬이어야 합니다.")
if b.ndim != 1 or len(b) != A.shape[0]:
raise ValueError("b의 길이는 A의 행 개수와 같아야 합니다.")
n = len(b)
# 전진 소거
for k in range(n-1):
pivot_row = k+np.argmax(np.abs(A[k:, k]))
if abs(A[pivot_row, k]) < tol:
raise np.linalg.LinAlgError("유일한 해가 없습니다.")
if pivot_row != k:
# partial pivoting
A[[k, pivot_row], :] = A[[pivot_row, k], :]
b[[k, pivot_row]] = b[[pivot_row, k]]
for i in range(k+1, n):
multiplier = A[i, k]/A[k, k]
A[i, k:] = A[i, k:]-multiplier*A[k, k:]
b[i] = b[i]-multiplier*b[k]
if abs(A[-1, -1]) < tol:
raise np.linalg.LinAlgError("유일한 해가 없습니다.")
# 후진 대입
x = np.zeros(n)
for i in range(n-1, -1, -1):
known_terms = A[i, i+1:]@x[i+1:]
x[i] = (b[i]-known_terms)/A[i, i]
return x
5.1. 계산 예제
원래 강의에서 사용한 연립방정식을 풀어보자.
\[\begin{aligned} 2x+3y+z&=9,\\ x-y+z&=1,\\ 3x+11y+5z&=35. \end{aligned}\]A = np.array([
[2.0, 3.0, 1.0],
[1.0, -1.0, 1.0],
[3.0, 11.0, 5.0],
])
b = np.array([9.0, 1.0, 35.0])
x = gaussian_elimination(A, b)
print("solution:", x)
print("NumPy:", np.linalg.solve(A, b))
계산 결과는
\[(x,y,z) = \left( \frac{1}{3}, \frac{23}{12}, \frac{31}{12} \right) \approx (0.3333,1.9167,2.5833)\]이다.
5.2. 잔차 확인
계산한 해가 원래 연립방정식을 얼마나 잘 만족하는지 잔차 벡터로 확인한다.
\[[\boldsymbol r]=[\boldsymbol A][\boldsymbol x]-[\boldsymbol b]\]residual = A@x-b
print("residual:", residual)
print("residual norm:", np.linalg.norm(residual))
부동소수점 연산 때문에 잔차가 정확히 0이 아니라 매우 작은 값으로 나타날 수 있다.
6. Gauss–Jordan 소거법과의 차이
가우스 소거법은 전진 소거로 상삼각행렬을 만든 뒤 후진 대입한다.
\[\left[ \begin{array}{ccc|c} u_{11}&u_{12}&u_{13}&c_1\\ 0&u_{22}&u_{23}&c_2\\ 0&0&u_{33}&c_3 \end{array} \right]\]Gauss–Jordan 소거법은 피벗 위의 원소까지 모두 0으로 만들고 피벗을 1로 만들어 다음과 같은 기약 행 사다리꼴을 구한다.
\[\left[ \begin{array}{ccc|c} 1&0&0&x_1\\ 0&1&0&x_2\\ 0&0&1&x_3 \end{array} \right]\]둘은 관련된 방법이지만 동일한 알고리듬은 아니다. 연립방정식 하나를 푸는 표준 가우스 소거법에서는 모든 비대각 원소를 0으로 만들 필요가 없다.
7. 계산량과 주의 사항
- 피벗이 0이면 행 교환이 필요하다.
- 매우 작은 피벗은 반올림 오차를 증폭시킬 수 있으므로 부분 피벗팅을 사용한다.
- 행렬이 특이행렬이면 유일한 해가 없다.
- 입력 배열을 직접 수정하지 않으려면 복사본을 만들어 계산한다.
- 실제 프로그램에서는 검증된 np.linalg.solve를 사용하는 것이 좋다.
- 학습 목적으로 가우스 소거법을 직접 구현하면 선형대수 알고리듬의 구조를 이해하는 데 도움이 된다.
8. 정리
- 가우스 소거법은 전진 소거와 후진 대입으로 구성된다.
- 확대행렬의 기본 행 연산은 연립방정식의 해를 바꾸지 않는다.
- 소거 계수는 $m_{ik}=a_{ik}/a_{kk}$로 계산한다.
- 부분 피벗팅은 현재 열에서 절댓값이 가장 큰 피벗을 선택한다.
- 계산 결과는 잔차 $[\boldsymbol A][\boldsymbol x]-[\boldsymbol b]$로 확인한다.
- Gauss–Jordan 소거법은 피벗 위까지 소거한다는 점에서 가우스 소거법과 다르다.
9. 연습 문제
강의에서 다룬 행 연산과 계산 과정을 확인하는 기초 문제이다.
-
다음 빈칸을 채워라.
- 계수행렬과 우변 벡터를 결합한 행렬을 ( )이라고 한다.
- 가우스 소거법은 ( )와 후진 대입의 두 단계로 이루어진다.
- 현재 계산에서 기준으로 사용하는 대각 원소를 ( )이라고 한다.
-
다음 연립방정식의 확대행렬을 작성하라.
\[\begin{aligned} 2x+y&=5,\\ x-y&=-1. \end{aligned}\] -
다음 중 기본 행 연산으로 사용할 수 있는 것을 모두 고르라.
- 두 행을 교환한다.
- 한 행에 0을 곱한다.
- 한 행에 0이 아닌 상수를 곱한다.
- 한 행에 다른 행의 배수를 더한다.
-
첫 번째 피벗이 2이고 그 아래에서 없애려는 원소가 6일 때 소거 계수 $m_{21}$을 구하라.
-
다음 확대행렬에 $R_2\leftarrow R_2-2R_1$을 수행하라.
\[\left[ \begin{array}{cc|c} 1&1&3\\ 2&3&8 \end{array} \right]\] -
다음 상삼각 연립방정식을 후진 대입으로 풀어라.
\[\begin{aligned} x+2y&=5,\\ 3y&=6. \end{aligned}\] -
다음 확대행렬에서 첫 번째 단계의 부분 피벗팅을 수행하면 어느 두 행을 교환해야 하는가?
\[\left[ \begin{array}{cc|c} 0.1&1&2\\ 2&3&4 \end{array} \right]\] -
다음 코드의 빈칸을 채워 해를 계산하라.
A = np.array([[1.0, 1.0], [2.0, 3.0]]) b = np.array([3.0, 8.0]) x = gaussian_elimination(____, ____) print(x) -
계산된 해 $\boldsymbol x=(1,2)$와 다음 행렬 및 벡터에 대해 잔차를 구하라.
\[[\boldsymbol A]= \begin{bmatrix} 1&1\\ 2&3 \end{bmatrix}, \qquad [\boldsymbol b]= \begin{bmatrix} 3\\8 \end{bmatrix}.\] -
가우스 소거법과 Gauss–Jordan 소거법의 차이를 한 문장으로 설명하라.