1. 선형 시스템과 수치해석
기계공학도로 대부분 배워온 구조와 유체, 전자기장 등이 있을 것이다. 대부분의 공학적 분야는 이런 선형 시스템을 분석하는 것인데, 여기서 선형 시스템은 입력에 비례하는 응답이 선형적인 것으로 간주되는 시스템이다. 만약 트러스 구조 해석이나 전기 회로(이거는 선대로 풀어본 적은 없지만)라면 시스템은 불연속적이다. 그때 그때 상태가 기록되면 각각이 방정식이 된다. 불연속(discrete)한 시스템은 선형 대수 방정식으로 바로 넘어가고, 연속(linear)한 시스템은 미분 방정식으로 설명된다. 그러나 수치 해석은 이산(discrete) 변수로만 처리할 수 있어 미방을 대수 연립방정식으로 근사한다. 여기에 쓰이는게 유한 요소, 경계 요소 방법 등등이다. 근사하는 방법은 다를 수 있어도 결국 문제는 Ax = b를 푸는 것이다.
A는 시스템의 특성, x는 시스템의 응담, b는 입력이다.
결국 문제는 풀리지만 우리는 공학도니까, 최소한의 계산 비용으로 방정식을 푸는 알고리즘을 써먹는게 중요하다.
이제부터 알아보자.
2. 해를 찾는 방법 - 직접방법
선형 대수 방정식을 푸는 법은 두 가지가 있다. 직접 Direct method, 반복 Iterative method 방법이다.
직접 방법은 있는 방정식을 더 쉽게 풀 수 있도록 바꾼다. 알고리즘 수업을 들었거나 수치해석을 배웠다면 가우스 소거를 배웠을 것이다. 이게 대표적인 직접 방법으로, 세 가지 기본 연산을 사용한다.
- 두 방정식 교환하기 → 행렬식 \( |A| \) 부호 변경함
- 방정식에 0 아닌 상수 곱하기 → 같은 상수로 행렬식 값 곱함
- 방정식에 0 아닌 상수 곱하고 다른 방정식에서 빼기 → 행렬식 값 변경 X
이 기본 연산들은 해에 변경을 주지 않으면서 더 간단한 계수 행렬 A를 얻을 수 있도록 해준다.
(반복 Iterative method 는 해 x의 추측으로 시작해 특정한 수렴 조건까지 계속 해를 수정해나간다. 계수 행렬이 매우 크고 sparse할 경우 더 좋다고 한다)
유명한 Direct method는 다음과 같다.
| 방법 | 초기 형태 | 최종 형태 |
| 가우스 소거 | \( Ax = b\) | \( Ux = c\) |
| LU 분해 | \( Ax = b\) | \( LUx = b \) |
| 가우스-조르단 소거 | \( Ax = b\) | \( Ix = c \) |
U는 상삼각행렬 upper triangular matrix , L은 하삼각행렬 lower triangular matrix , I는 항등행렬 identity matrix이다.
상삼각행렬로 방정식이 정리된 경우(가우스 소거) 후진대입으로 해를 찾을 수 있다.
하삼각행렬로 방정식이 정리된 경우 전진대입, 그리고 항등행렬로 정리된 경우 그 자체가 해이다.
2.1 가우스 소거 Gauss Elimination
Eigen::VectorXd gaussElimin(const Eigen::MatrixXd &a_, const Eigen::VectorXd &b_){
int n = b_.rows();
Eigen::MatrixXd a = a_;
Eigen::VectorXd b = b_;
// Elimination
for(size_t k=0; k<n-1; k++){ // k번째 행은 이미 상삼각행렬로 변환됨 -> k번째 행이 pivot 행이다. 처음(0) ~ 마지막에서 두 번째(n-2) 행
for(size_t i=k+1; i<n; i++) { // 변환될 행으로 pivot 행의 바로 다음(k+1) ~ 마지막(n-1) 행
if (a(i, k) != 0) {
float lambda = a(i, k) / a(k, k);
a(i, Eigen::seq(k+1, Eigen::last)) = a(i, Eigen::seq(k+1, Eigen::last)) - lambda * a(k, Eigen::seq(k+1, Eigen::last));
b(i) = b(i) - lambda * b(k);
}
}
}
// Back substitution
for(int i= n-1; i > -1; i--){ // 후진대입이니까 마지막(n-1) ~ 처음(0) 행까지
Eigen::VectorXd substituteA, substituteB;
substituteA = a(i, Eigen::seq(i+1, Eigen::last));
substituteB = b(Eigen::seq(i+1, Eigen::last));
b(i) = ( b(i) - substituteA.dot(substituteB) ) / a(i,i);
}
return b;
}넘파이의 슬라이싱 기능이 진짜 좋다는걸 느끼게 됐다 :)
해보니까 소거 부분에서 열 요소를 k~n-1으로 하나 k+1 ~ n-1으로 하나 결과는 당연히 같지만 손으로 푸는거랑 조금 간격이 있어 당황스러웠다. 알고리즘적으로는 람다 상수배를 한걸 빼줬으니 당연히 0열이 0이 되겠거니 하고 아예 0열을 안보는데, 손으로 푼거랑은 조금 달라서!
2.2 LU 분해 factorization
앞선 가우스 소거는 계수행렬 A를 기본연산으로 편리하게 만든 후 사용했다.
LU분해 decomposition, 혹은 LU 인수분해 factorization 는 계수행렬 A를 상삼각행렬 U와 하삼각행렬 L 로 분해한다. 이 과정에서 LU는 고유하지 않다. 즉, L과 U의 조합은 무한하다.
| 이름 | 제약 조건 |
| Doolittle 분해 | \( L_{ii} =1,\; i=1, \; 2\; ,\cdots, n \) |
| Crout 분해 | \( U_{ii} =1,\; i=1, \; 2\; ,\cdots, n \) |
| Choleski 분해 | \( \bf L = U^T \) |
\( Ax = b\)를 \( LUx=b \)로 분해했을 때, 풀이법은 다음과 같다.
- \(Ux = y\)로 치고 \(Ly = b\)를 전진대입으로 y에 대해 푼다.
- \(Ux=y\)가 있을 때 y를 알고 있으므로 후진대입으로 x에 대해 푼다.
2.2.1 Doolittle 분해
Eigen::MatrixXd LUdecomp_doolittle(const Eigen::MatrixXd &a_){
Eigen::MatrixXd a = a_;
int n = a.rows();
for(size_t k=0; k<n-1; k++){
for(size_t i=k+1; i<n; i++){
if( a(i,k) != 0.0){
float lambda = a(i,k) / a(k,k);
a(i, Eigen::seq(k+1, Eigen::last)) -= lambda * a(k, Eigen::seq(k+1, Eigen::last));
a(i,k) = lambda;
}
}
}
return a;
}- 분해할 때 U는 가우스 소거할 때의 상삼각행렬과 동일하다.
- \( A=LU\)에 따라 L의 대각선이 아닌 요소는 가우스 소거 중에 사용되는 피벗 계산에 사용되는 lambda다.
이렇게 분해하고 나면
$$ A = [L/U] = \begin{bmatrix}U_{11} & U_{12} & U_{13}\\
L_{21} & U_{22} & U_{23} \\
L_{31} & L_{32} & U_{33} \end{bmatrix} $$
이 되어 L과 U의 혼합이 된다.
이제 L과 U를 알았으니
- \(Ly = b\)에서 전진대입해 y를 풀고
- \(Ux=y\)에서 y를 알고 있으므로 후진대입으로 x에 대해 푼다.
Eigen::VectorXd LUsolve(const Eigen::MatrixXd &a_, const Eigen::VectorXd &b_){
Eigen::MatrixXd a = a_;
Eigen::VectorXd b = b_;
int n = b_.rows();
// forward substitute
for(size_t k=1; k<n; k++) b(k) -= a(k, Eigen::seq(0,k-1)).dot(b(Eigen::seq(0,k-1)));
b(n-1) /= a(n-1, n-1);
// backward substitute
for(int k=n-2; k>-1; k--) b(k) = (b(k) - a(k, Eigen::seq(k+1, n-1)).dot(b(Eigen::seq(k+1, n-1))) ) / a(k,k);
return b;
}2.2.2 Choleski 분해
Choleski 분해 \( A=LL^T \)는 행렬 A가 대칭행렬 Symmetric matrix 일때만 성립한다.
또, 분해 과정에서 A 원소들의 제곱근을 취하는데(try문 내부) 음수의 제곱근을 피하기 위해서 A 행렬은 양의 정부호 행렬 positive definite matrix 이어야 한다.
공돌이의 수학정리노트 블로그에서 아주아주아주 잘 설명되어 있다.
https://angeloyeo.github.io/2021/12/20/positive_definite.html
양의 정부호 행렬 (positive definite matrix) - 공돌이의 수학정리노트 (Angelo's Math Notes)
angeloyeo.github.io
하여튼 cholesky 분해는 모든 요소에 접근하는 doolittle 분해와 다르게 \( A^T = A,\; A=LL^T \)의 대칭성을 써먹기 때문에 알고리즘의 시간 비용이 적은 편이다. 써먹기 위해서는 계수 행렬 A가 대칭이어야 한다는 것 뿐!
Eigen::MatrixXd choleskiDecomp(const Eigen::MatrixXd &a_){
int n = a_.rows();
Eigen::MatrixXd a = a_;
for(size_t k=0; k<n; k++){
try
{
a(k,k) = std::sqrt(a(k,k) - a(k,Eigen::seq(0,k-1)).dot(a(k,Eigen::seq(0,k-1))) );
if(a(k,k) < 0) throw std::invalid_argument( "received negative value" );
}
catch ( const std::invalid_argument& e )
{
std::cout << "Matrix is not positive definite" << std::endl;
}
for(size_t i=k+1; i<n; i++) a(i,k) = (a(i,k) - a(i,Eigen::seq(0,k-1)).dot(a(k,Eigen::seq(0,k-1)))) / a(k,k);
for(size_t k=1; k<n; k++){
for(size_t i=0; i<k; i++) a(i,k) = 0.0;
}
}
return a;
}
해를 구하는 choleskiSolve는 L을 받아서 LU분해 풀이법에 사용한다.
- \(Ly = b\)에서 전진대입해 y를 풀고
- \(L^Tx=y\)에서 y는 1단계에서 기록된 b에 반영되어있으니 후진대입으로 x에 대해 풀면, 마지막 b가 연립방정식의 해 x가 된다.
Eigen::VectorXd choleskiSolve(const Eigen::MatrixXd &L, const Eigen::VectorXd &b_){
int n = L.rows();
Eigen::VectorXd b = b_;
// SOL for Ly=b
for(size_t k=0; k<n; k++) b(k) = (b(k) - L(k,Eigen::seq(0,k-1)).dot(b(Eigen::seq(0,k-1))) ) / L(k,k);
// SOL for L_transpose x = y
for(int k=n-1; k>-1; k--) b(k) = (b(k) - L(Eigen::seq(k+1,n-1), k).dot(b(Eigen::seq(k+1,n-1))) ) / L(k,k);
return b;
}2.3 돌려보기
헤더 파일은 위 함수들을 넣은 것들이다.
지금은 비공개로 채워나가는데 나중에 완성되어가면 공개로 돌릴 것 같다.
졸업준비랑 기말고사가 다가오다보니 좀 늦었다.
'Study > 공학수치해석' 카테고리의 다른 글
| 1.1 연립 선형 대수 방정식 Intro (0) | 2023.11.08 |
|---|---|
| 시작 글 (0) | 2023.11.06 |