# 수치적 불안정성의 개념

수치해석에서 말하는 불안정성은 계산 과정에서 발생하는 오차가 계속해서 증폭되어 결과적으로 의미 있는 해석이 불가능해지는 현상을 지칭한다. 실제 연산은 이산화된 컴퓨터 아키텍처, 유효 자릿수가 제한된 부동소수점 표현, 내부적으로 발생하는 반올림(rounding)과 잘림(chopping) 등으로 인해 완벽하게 이루어지지 않는다. 이때 약간의 오차가 누적되거나 특정 조건에서 급격히 커져서 최종 계산 결과가 왜곡되는 문제가 생길 수 있다. 이러한 문제를 방지하거나 예측하기 위해서는 먼저 수치적 불안정성의 근원, 형태, 그리고 그 해석 방법을 이해해야 한다.

#### 연산 오차와 정밀도

컴퓨터에서 사용하는 부동소수점 표현은 실제 수를 유한한 비트로 근사 표현하기 때문에 언제나 오차를 수반한다. 예를 들어 실수 $x$를 표현할 때 컴퓨터에서 사용할 수 있는 인접 부동소수점 수를 $\mathrm{fl}(x)$라고 하면 $x$와 $\mathrm{fl}(x)$의 오차는 근본적으로 불가피하다. 실제로 연산을 수행할 때도 두 수의 덧셈, 뺄셈, 곱셈, 나눗셈 등이 이루어질 때마다 다음과 같은 반올림 현상이 발생한다.

$$
\mathrm{fl}(a ; \circ ; b) = (a ; \circ ; b)\bigl(1 + \epsilon\bigr)
$$

위 식에서 $\circ$는 사칙연산을 나타내고, $\epsilon$은 매우 작지만 완전히 0이 아닌 값을 의미한다. 이러한 과정이 누적되는 것이 궁극적으로 수치적 불안정성을 일으키는 토대가 된다. 정밀도가 큰(예: 배정밀도, 쿼드러플 정밀도 등) 연산 환경일수록 이러한 오차가 상대적으로 작지만, 알고리즘적 특성과 문제의 조건수에 따라 오차가 악화되는 상황은 여전히 발생할 수 있다.

#### 조건수(Condition Number)와 민감도(Sensitivity)

수치적 불안정성을 이해하려면 먼저 문제 자체가 얼마나 민감하게 반응하는지를 알아야 한다. 민감도는 주어진 문제의 입력이 미세하게 변할 때 해가 얼마나 크게 달라지는지를 측정하는 개념이다. 예를 들어 선형시스템 $\mathbf{A}\mathbf{x} = \mathbf{b}$에서 행렬 $\mathbf{A}$가 임의의 소규모 변화에 대해 해 $\mathbf{x}$가 크게 요동한다면, 문제 자체가 이미 민감하다고 해석할 수 있다. 이처럼 문제의 민감도는 조건수(condition number)로 측정되며, 보통 노름 $|\cdot|$에 대하여 다음과 같이 정의된다.

$$
\kappa(\mathbf{A}) = |\mathbf{A}||\mathbf{A}^{-1}|
$$

$\kappa(\mathbf{A})$가 클수록 $\mathbf{x}$의 오차가 입력 $\mathbf{b}$ 또는 행렬 $\mathbf{A}$의 작은 변화에 의해 더 크게 증폭될 가능성이 있다. 문제 자체가 민감한 경우를 일컫어 ill-conditioned라 하고, 그렇지 않은 경우를 well-conditioned라고 한다. 문제 자체가 ill-conditioned라면 아무리 안정적인 알고리즘으로 계산해도 일정 수준 이상의 오차 증가를 피하기 어렵다.

#### 전방오차(Forward Error)와 후방오차(Backward Error)

해석적으로 문제를 다룰 때는 실제 해(정확해)와 수치해(근사해)의 차이를 정량화해야 한다. 전방오차는 최종적으로 구해진 수치해 $\mathbf{\tilde{x}}$가 참해 $\mathbf{x}$와 얼마나 다른지를 직접 측정하는 방식이다.

$$
\text{forward error} = |\mathbf{x} - \mathbf{\tilde{x}}|
$$

후방오차는 수치해를 얻은 결과를 다른 시각으로 해석한다. 즉, 실제로는 다음을 만족하는 근사 문제

$$
\mathbf{\tilde{x}} = \mathbf{b} + \Delta \mathbf{b}
$$

에 대한 정확해가 $\mathbf{\tilde{x}}$였다고 본다. 이때 $\Delta \mathbf{A}$와 $\Delta \mathbf{b}$가 매우 작다면, 수치해가 원래 문제와 거의 동일한 문제에 대한 '정확한 해'로서 해석될 수 있다. 이를 두고 그 알고리즘이 후방안정(backward stability)을 가진다고 말한다. 어떤 알고리즘의 후방안정성 여부는 문제의 조건수와 함께 수치적 안정성을 판단하는 핵심 기준이 된다.

#### 알고리즘의 안정성(Stability)

알고리즘 자체가 불안정하면, 아무리 컴퓨터의 정밀도가 커도 연산 과정에서 발생하는 오차가 지속적으로 증폭된다. 안정적 알고리즘(stable algorithm)은 오차가 문제가 되는 순간에 내부적으로 자동 억제되어 결과적 오차가 커지지 않도록 설계된 알고리즘이다. 반면 불안정(unstable)한 알고리즘은 작은 입력 오차나 연산 오차가 큰 결과 편차로 이어진다. 특히 ill-conditioned 문제에서 불안정 알고리즘을 사용하면, 결과가 전혀 유의미하지 않을 정도로 오차가 커질 수 있다. 이 때문에 실제로는 조건수가 낮은 문제라도, 혹은 그 반대 경우라도 적절한 알고리즘 선택이 필수적이다.

알고리즘의 안정성을 판단하기 위해서는 대부분의 경우 전방오차와 후방오차 모두를 고려해야 한다. 예컨대 Gauss 소거법에서 피벗팅(pivoting) 전략을 어떻게 쓰느냐에 따라 매우 다른 정도의 안정성을 보인다. 부분 피벗팅(Partial Pivoting)을 적용하면 오차 증폭을 효과적으로 억제할 수 있지만, 피벗팅을 전혀 하지 않는 알고리즘은 극단적인 경우 수치적 불안정성을 쉽게 초래할 수 있다. 또한 계수 행렬의 구조, 대칭성 여부, 스펙트럼 분포 등에 따라 필요한 안정화 기법이 달라지기도 한다.

#### 절대오차와 상대오차

수치해석에서 해의 정확도를 평가할 때 오차를 나타내는 지표로 절대오차와 상대오차가 있다. 문제의 입력 혹은 출력값을 $x$라 하고, 근사값을 $\tilde{x}$라 할 때 절대오차는 다음과 같이 정의된다.

$$
\text{Absolute Error} = |,x - \tilde{x},|
$$

상대오차는 절대오차를 참값의 크기로 나눈 비율로 정의한다.

$$
\text{Relative Error} = \frac{|,x - \tilde{x},|}{|,x,|}
$$

상대오차는 참값 $x$의 크기에 비례하여 오차 정도를 판단하고자 할 때 사용하며, 특히 $x$가 매우 클 때 혹은 매우 작을 때 오차가 얼마나 중요한지를 상대적으로 살펴보는 데 유용하다. 반면 $x$가 0에 매우 가깝거나 실제로 0일 경우에는 상대오차의 정의가 곤란해지므로, 이때는 절대오차가 더 적합한 지표가 될 수 있다.

이러한 절대오차와 상대오차는 문제의 민감도, 알고리즘의 안정성, 그리고 부동소수점에서의 실제 구현을 종합적으로 고려해야 한다. 문제에 따라 어떤 오차 척도가 더 의미 있는지는 다를 수 있지만, 일반적으로 수치해석에서는 상대오차가 좀 더 널리 사용된다. 이는 부동소수점 시스템에서 큰 값에 대해서는 작은 절대오차도 중요한 결과 왜곡을 일으킬 수 있고, 작은 값에 대해서는 비교적 큰 절대오차라도 전체 문제에서의 영향이 적을 수 있기 때문이다.

#### 큰수와 작은수 연산에서의 취약성

컴퓨터 부동소수점 연산에서 유효숫자 자릿수는 한정적이므로, 매우 큰 수와 매우 작은 수가 함께 연산될 때 극단적인 손실(loss of significance)이 발생할 가능성이 커진다. 예를 들어 큰 실수 $M$에 비해 매우 작은 실수 $\delta$를 더하거나 빼는 연산은 다음과 같은 문제를 야기한다.

$$
\mathrm{fl}(M \pm \delta) = \begin{cases} M, & \text{if } \delta \text{ is too small}, \ M \pm \delta, & \text{otherwise}. \end{cases}
$$

만일 $\delta$가 부동소수점 정밀도로는 표현하기 어려울 만큼 작은 수라면, 실제로는 $M \pm \delta \approx M$가 되어 $\delta$ 자체가 반영되지 않는 현상이 생긴다. 이처럼 실제 수학적 연산과 부동소수점 연산의 결과가 달라지는 것을 catastrophic cancellation(치명적 소거)이라고 부르기도 한다. 특히 뺄셈에서 대략 같은 크기를 가진 두 수가 서로 소거될 때, 남는 유효숫자가 현저히 줄어드는 문제가 대표적이다.

예를 들어 $x \gg y > 0$인 상황에서 $x - (x - y)$와 같이, 실제로는 $y$가 되어야 하는 연산이 부동소수점 환경에서는 $y$보다 훨씬 작은 값으로 계산되거나 0이 되어버릴 수도 있다. 이 현상은 알고리즘의 안정성과도 깊이 관련된다. 어떤 알고리즘 내에서 두 수의 뺄셈이 빈번히 일어날 때, 특히 차이가 작은 수끼리 뺄셈이 자주 반복되면 오차가 크게 누적될 위험이 있다. 이러한 위험은 문제의 스케일 조정, 적절한 재배열(summation ordering), 혹은 특정 수학적 변형을 통해 완화할 수 있다.

#### 다항식 평가와 불안정성

다항식의 값 $p(x)$을 직접 계산하기 위해 계수와 $x$의 여러 거듭제곱을 순차적으로 곱해가면서 더하는 방식이 반드시 안정적이지는 않다. 예를 들어 다항식

$$
p(x) = a\_n x^n + a\_{n-1} x^{n-1} + \cdots + a\_1 x + a\_0
$$

를 한 번에 계산하려면, 큰 $n$에서 시작해 $x^n, x^{n-1}, \ldots$ 등을 구하여 차례로 곱셈과 덧셈을 수행한다. 이 과정에서 $|x| \gg 1$이거나 $|x| \ll 1$인 경우에 지수 항이 극단적인 크기 차이를 보이게 되면, 덧셈이나 뺄셈 과정에서 소수점 이하 자리들이 소거될 위험이 있다. Horner’s Method와 같이 알고리즘을 재구성하면 동일한 결과를 좀 더 안정적으로 계산할 수 있는데, 이는 곱셈과 덧셈의 순서를 바꾸어 오차를 억제하기 때문이다.

Horner의 방법은 다음과 같은 형태로 다항식을 재배열하여 계산한다.

$$
p(x) = a\_n x^n + a\_{n-1} x^{n-1} + \cdots + a\_1 x + a\_0 = a\_n \bigl(x(\cdots (x(a\_2 x + a\_1) + a\_0)\cdots)\bigr)
$$

실제로는 괄호와 순서가 좀 더 체계적으로 배치되어, 큰 지수 항들을 바로 계산하지 않고 여러 번 나누어 곱셈과 덧셈을 수행하여 불필요한 반올림 오차가 누적되는 것을 억제한다. 이는 일종의 후방안정성을 목표로 한 재배열 방법이라 볼 수 있으며, 수치적인 안정성 증대에 매우 효과적이다.

#### 예시: 취소오차(cancellation error)의 간단한 시연

부동소수점 연산 환경에서, 예를 들어 $a \approx b$인 두 실수에 대하여 $a - b$를 구하면 유효 자릿수가 많이 소실될 수 있음을 보인다. 다음은 간단한 예시 코드를 Python으로 작성한 것이며, 매우 근접한 두 값을 연산했을 때 발생하는 차이를 확인할 수 있다.

```python
import math

a = 1.000000000000001
b = 1.000000000000000
c = a - b

print("a - b =", c)
print("정확 계산(이론상 a - b) =", (1.000000000000001 - 1.000000000000000))
```

대부분의 이중정밀도(double precision) 연산 환경에서 $c$ 값이 $1 \times 10^{-15}$ 정도로 계산되기는 하지만, 실제로는 더 정확하게 표현하고 싶은 $\delta$가 잘려나갈 수 있다. $a$와 $b$가 더 큰 수라면(예: 1e10 정도의 스케일) 이 차이는 더욱 심각해질 수 있다. 이것이 곧 취소오차의 전형적인 예이며, 문제의 조건수와 알고리즘의 특성을 고려하여 적절한 대응책을 마련해야 한다.

#### 가우스 소거와 피벗팅의 중요성

선형시스템 $\mathbf{A}\mathbf{x}=\mathbf{b}$를 풀기 위한 기본 알고리즘인 가우스 소거법(Gaussian Elimination)은 매우 널리 사용되고 있지만, 피벗팅(pivoting) 전략이 적용되지 않거나 부적절하게 적용되면 수치적 불안정성이 쉽게 발생한다. 가우스 소거의 단계별 연산 과정을 살펴보면, 특정 단계에서 분할(divisor)이 되는 피벗(pivot) 원소가 매우 작으면 작은 수로 나누는 일이 빈번해져 연산 오차가 크게 증폭될 가능성이 커진다. 또한, 연산 도중 일부 행(또는 열)을 교환하여 피벗 위치를 적절히 선택하는 과정을 누락하면, 작은 피벗 원소가 그대로 유지되어 오차가 누적될 확률이 훨씬 높아진다.

부분 피벗팅(Partial Pivoting)은 각 단계에서 현재 열을 기준으로 가장 큰 절댓값을 갖는 행을 찾아 피벗 행과 교환하는 방식이다. 완전 피벗팅(Complete Pivoting)은 행뿐만 아니라 열까지 교환하여 피벗 원소를 최대로 만드는 전략이다. 완전 피벗팅이 이론적으로 더 안전하지만, 연산량이 증가하고 실제 구현에서 행 교환에 비해 열 교환은 기억장치 접근을 복잡하게 할 수 있기 때문에, 실무에서는 부분 피벗팅이 주로 사용된다. 부분 피벗팅만으로도 대부분의 경우 후방안정성을 확보하기에 충분하고, ill-conditioned 사례에서도 큰 오차 폭증을 어느 정도 억제할 수 있다.

#### 스케일링(Scaling)과 행, 열 정규화

가우스 소거법을 적용하기 전에 계수 행렬 $\mathbf{A}$의 각 행 또는 열을 적절히 정규화(normalization)하는 방식으로, 알고리즘 수행 과정에서 상대적으로 더 균일한 크기의 연산을 유지할 수 있다. 예를 들어 행별로 가장 큰 계수의 절댓값이 1이 되도록 만들어두면, 뺄셈 과정에서 유의미한 자릿수가 과도하게 손실되는 일을 줄일 수 있다. 실제 수학적 해석에서는 스케일링이 문제의 조건수를 근본적으로 바꾸지는 않지만, 부동소수점 연산 관점에서는 알고리즘적 불안정성을 어느 정도 완화하는 실용적 기법이 된다.

스케일링과 피벗팅은 상호 보완적인 역할을 하며, 때로는 스케일 팩터(scale factor)를 추적 관리하면서 연산 과정 전반에 걸쳐 동적으로 갱신하기도 한다. 이 과정에서 특히 주의해야 할 점은, 지나친 스케일링 또는 불필요한 행과 열 교환이 새로운 round-off 오차를 유발할 수도 있다는 것이다. 따라서 실제 구현 시에는 문제 크기, 계수 행렬의 구조, 기억장치 접근 방식 등을 포괄적으로 고려하여 가장 효율적인 수준의 스케일링과 피벗팅을 설계하게 된다.

#### 잔차(Residual)와 잔차 기반 보정

선형시스템을 풀고 난 후 해 $\mathbf{\tilde{x}}$의 유효성을 확인하기 위해 잔차 $\mathbf{r} = \mathbf{b} - \mathbf{A}\mathbf{\tilde{x}}$를 살펴본다. 이 잔차가 이상적으로는 0에 가까워야 하지만, 실제 수치 계산에서는 부동소수점 표현 오차로 인해 약간의 편차가 발생한다. 만약 잔차가 예측보다 크게 나타난다면(예: $|\mathbf{r}|\gg \epsilon|\mathbf{b}|$, 어떤 작동 임계치 $\epsilon$ 설정), 알고리즘적 불안정성 또는 문제 자체가 상당히 ill-conditioned인지 진단할 필요가 있다.

잔차 기반 보정 기법(iterative refinement)은 이미 구한 근사해 $\mathbf{\tilde{x}}$를 이용하여, 다시 한 번 $\mathbf{r}$를 계산한 후 이를 해석적으로 보완하여 새로운 해를 업데이트하는 방식이다. 간단한 형태의 보정 알고리즘은 대략 다음과 같은 수식을 따른다.

$$
\begin{align}
\mathbf{r} &= \mathbf{b} - \mathbf{A}\mathbf{\tilde{x}}
\\
\Delta \mathbf{x} &= \mathbf{A}^{-1} \mathbf{r} \quad (\text{또는} \ \mathbf{A}\Delta \mathbf{x} = \mathbf{r} \ \text{를 다시 풀어} \ \Delta \mathbf{x}\ \text{얻음})
\\
\mathbf{\tilde{x}} &\leftarrow \mathbf{\tilde{x}} + \Delta \mathbf{x}
\end{align}
$$

연산 과정마다 $\mathbf{r}$를 다시 계산하여, 충분히 작아질 때까지(또는 최대 반복 횟수에 도달할 때까지) 이 과정을 반복한다. 다만 $\mathbf{A}^{-1}$을 매번 직접 구하는 대신, 원래 가우스 소거 결과의 LU 분해 등을 활용하여 $\mathbf{A}\Delta \mathbf{x}=\mathbf{r}$를 푼다. 이 방법은 해가 이미 어느 정도 근사 상태에 있을 때 효과적이고, 이상적으로는 잔차가 점차 줄어들어 더욱 정확한 근사해를 얻게 된다. 그러나 문제가 지나치게 ill-conditioned이거나 첫 단계 해가 잘못 구해졌으면, 잔차 기반 보정의 효과가 미미할 수 있다.

#### 고유값 문제에서의 수치적 불안정성

고유값 문제(eigenvalue problem)는 선형대수 전반에 걸쳐 다양한 응용을 갖는다. $\mathbf{A}$가 $n\times n$ 행렬일 때, 고유값 $\lambda$와 고유벡터 $\mathbf{v}$는 다음을 만족한다.

$$
\mathbf{A}\mathbf{v} = \lambda \mathbf{v}
$$

고유값 문제를 직접 해석적으로 푸는 것이 아닌, QR 알고리즘, 멱반복법(power method) 등 수치 알고리즘을 통해 구할 때도 안정성 문제가 대두된다. 예를 들어 $\mathbf{A}$가 대칭 실대각화가 가능한 행렬인 경우, QR 알고리즘의 내재적 안정성(inherent numerical stability) 덕분에 비교적 안정적으로 고유값과 고유벡터를 얻을 수 있다. 반면 비대칭 행렬이나 조건수가 큰 사례에서는, 매 반복 단계에서 발생하는 반올림 오차가 고유값 분리(separation)에 따라 크게 증폭될 수 있다.

멱반복법의 경우 특정 지배적인 고유값(dominant eigenvalue)과 그에 대응하는 고유벡터를 근사하는 간단한 알고리즘이지만, 관심 있는 고유값이 지배적이지 않은 상황에서는 오차 증폭 문제, 근접 고유값 간의 구분 문제 등의 불안정성이 존재한다. 또한 실수 행렬의 복소 고유값 문제에서는 실수 연산을 복소수 연산으로 확장하는 과정에서 연산 비용뿐만 아니라 수치적인 위험도도 증가한다. 이처럼 고유값 문제도 기본적인 선형 연산의 반복으로 구성되므로, 이미 다루었던 수치불안정 요인들이 동일하게 적용됨을 염두에 두어야 한다.

#### 다항식 근 찾기와 수치적 불안정

다항식의 근을 수치적으로 구하는 문제에서도 계산 불안정성이 자주 부각된다. 예를 들어 고차 다항식의 해를 직접 해석적(analytic) 공식으로 푸는 것은 매우 복잡하거나 불가능한 경우가 많다. 따라서 일반적으로 뉴턴-랩슨(Newton-Raphson)이나 이분법(bisection) 등의 수치 기법을 사용한다. 그러나 특정 기법에서 초기값 설정이 잘못되거나, 근이 다중근(double root)인 경우, 혹은 서로 근접한 근들이 존재하면 알고리즘이 불안정하게 수렴하거나 무한 루프에 빠질 수도 있다.

특히 뉴턴-랩슨 방법은 다음과 같은 반복식을 가진다.

$$
x\_{k+1} = x\_k - \frac{f(x\_k)}{f'(x\_k)}
$$

이때 $f'(x\_k)$가 0에 근접하면 분모가 아주 작은 값이 되어 큰 오차가 발생하거나 발산할 수 있다. 또한 두 근이 서로 매우 가까울 때, 부동소수점 정밀도가 충분히 크지 않으면 두 근을 구분하지 못하고 하나의 근으로 합쳐져 버리는 문제가 생길 수 있다. 불안정성은 주로 미분계수 $f'(x)$가 작아지는 영역, 혹은 오차가 근 해석 상으로 더 크게 증폭되는 지점에서 두드러진다.

#### 간단한 근 찾기 예시

다음은 Python 코드로 뉴턴-랩슨 방법을 적용해보는 예시를 보인다. 여기에서는 $f(x)=x^3 - 2x - 5$와 같은 삼차 다항식을 대상으로 한다.

```python
def f(x):
    return x**3 - 2*x - 5

def df(x):
    return 3*x**2 - 2

def newton_raphson(f, df, x0, tol=1e-12, max_iter=1000):
    x = x0
    for i in range(max_iter):
        fx = f(x)
        dfx = df(x)
        if abs(dfx) < 1e-20:  # 분모가 너무 작아지면 중단
            print("경고: 미분계수가 거의 0에 근접")
            break
        x_next = x - fx/dfx
        if abs(x_next - x) < tol:
            break
        x = x_next
    return x

init_guess = 2.0
root_approx = newton_raphson(f, df, init_guess)
print("구해진 근의 값:", root_approx)
print("함수값(검증):", f(root_approx))
```

위 코드는 단순한 형태의 뉴턴-랩슨 반복을 수행하며, 미분계수가 지나치게 작아질 경우 경고 메시지를 띄우도록 했다. 다항식에 따라 적절한 초기값을 주면 안정적으로 수렴하지만, 만약 초기값 근방에서 $f'(x)$가 0에 가까운 지점이라면 반복 과정에서 큰 오차가 발생하고, 그 오차가 계속 누적되어 엉뚱한 방향으로 튈 가능성이 있다.

#### 비선형 방정식과 고차원 문제에서의 안정성

단일 변수 다항식이 아니라 여러 변수를 포함하는 비선형 방정식을 풀거나, 최적화 문제를 해결하는 과정에서도 수치적 불안정성은 동일한 원인으로 발생한다. 뉴턴-랩슨 방법을 다변수로 확장한 뉴턴-KKT(N-KKT) 기법 또는 준뉴턴(quasi-Newton) 기법 등에서도, 야코비 행렬 또는 헤시안 행렬이 치환 가능한 형태로 잘 조건화되어 있지 않다면, 작은 반올림 오차가 여러 단계를 거치며 과도하게 증폭될 수 있다.

예컨대 다변수 뉴턴 방법의 반복식을 살펴보면,

$$
\mathbf{x}\_{k+1} = \mathbf{x}\_k - \mathbf{J}^{-1}(\mathbf{x}\_k) , \mathbf{f}(\mathbf{x}\_k)
$$

여기에서 $\mathbf{J}(\mathbf{x}\_k)$는 $\mathbf{f}$의 야코비 행렬이다. $\mathbf{J}(\mathbf{x}\_k)$가 심각하게 ill-conditioned면, $\mathbf{J}^{-1}(\mathbf{x}\_k)$를 구하는 과정(또는 선형시스템 $\mathbf{J}(\mathbf{x}\_k)\mathbf{p}=-\mathbf{f}(\mathbf{x}\_k)$를 푸는 과정)에서 불안정성이 크게 나타난다. 이러한 문제를 방지하기 위해서는 선형 해법을 사용할 때처럼 피벗팅, 적절한 전처리(preconditioning), 혹은 새로운 접근 방식을 모색할 필요가 있다.

#### 미분 근사와 수치적 불안정

수치해석에서 미분을 근사하는 대표적인 방식은 전진 차분(forward difference), 후진 차분(backward difference), 중앙 차분(central difference)이다. 예를 들어 스칼라 함수 $f(x)$에 대해 전진 차분을 사용한 1차 미분 근사는

$$
f'(x) \approx \frac{f(x + h) - f(x)}{h}
$$

여기서 $h$는 아주 작은 양의 수로, 이 이론적 식만 봤을 때는 $h$가 작을수록 근사 정확도가 좋아야 할 것 같지만, 실제로는 $h$가 너무 작아지면 $f(x + h)$와 $f(x)$가 부동소수점 오차에 의해 소거되어, $f(x + h) - f(x)$에서 유효 숫자가 사라지는 문제가 발생할 수 있다. 다시 말해 근사 공식이 오차를 더 크게 유발할 수 있다. 따라서 수치적으로 최적화된 $h$를 찾는 것이 중요하다. 이 문제는 고차 차분 혹은 고차원 도함수 근사에서도 동일하게 일어난다.

미분이 들어가는 문제, 예를 들어 미분방정식이나 최적화 문제가 자주 등장하는 실제 응용에서는 이처럼 근사 과정에서 발생하는 수치적 불안정성에 대한 주의를 기울여야 한다. 일부 해석적 기법(예: 기호적 미분)을 활용하거나, 자동미분(automatic differentiation) 라이브러리를 사용하여 오차 누적을 줄이는 방식을 채택할 수도 있다.

#### 간단한 전진 차분 미분 예시

아래의 Python 예시는 $f(x)=\sin(x)$의 도함수를 전진 차분으로 구해보고, $h$가 지나치게 작거나 큰 경우의 문제를 보여준다.

```python
import numpy as np

def f(x):
    return np.sin(x)

def forward_diff(f, x, h):
    return (f(x + h) - f(x)) / h

x_val = 1.0
for h in [1e-1, 1e-3, 1e-5, 1e-7, 1e-9, 1e-11]:
    approx = forward_diff(f, x_val, h)
    true_val = np.cos(x_val)
    print("h =", h, "근사미분 =", approx, "오차 =", abs(approx - true_val))
```

출력 결과를 살펴보면, $h$를 너무 크게 잡으면 근사식 자체의 오차가 커지고, $h$를 너무 작게 잡으면 부동소수점 반올림 오류로 인한 취소가 커진다. 실제로는 그 중간 어딘가에서 오차가 최소가 되는 $h$가 존재하며, 이것은 수학적으로 $\sqrt{\epsilon\_{\text{machine}}}$ 정도의 크기와 관련이 있다고 알려져 있다(단, 문제의 스케일과 $f'(x)$의 변동성도 함께 고려해야 한다).

#### 강직성(Stiffness)이 있는 미분방정식의 불안정성

초기값 문제 형태의 상미분방정식(ODE)을 수치적으로 풀 때, 특히 계(system)가 강직성(stiff)을 가지면 수치해의 불안정성이 더욱 심각하게 나타날 수 있다. 강직성은 물리적으로 빠른 시간 스케일과 느린 시간 스케일이 공존하거나, 계수 행렬(선형화했을 때)이 고도로 조건이 나쁜 경우 등에서 나타나며, 수치 알고리즘 선택에 매우 민감하다. 예를 들어 초기값 문제

$$
\mathbf{y}'(t) = \mathbf{f}\bigl(t, \mathbf{y}(t)\bigr), \quad \mathbf{y}(t\_0) = \mathbf{y}\_0
$$

를 전진 오일러(Forward Euler)와 같은 명시적(explicit) 방법으로 풀면, 시간구간에서 매우 작은 스텝 크기($\Delta t$)를 강제로 사용하지 않으면 해가 발산하거나 진동하면서 물리적으로 의미 없는 결과를 얻을 수 있다. 이는 문제 자체가 강직성을 가져서, 고유치의 실부(Real Part)가 매우 음수인 방향으로 급격히 감소하는 모드(Mode)가 존재하기 때문이다. 명시적 방법은 이러한 급격한 스케일 차이를 안정적으로 처리하기가 어렵다.

#### 명시적 방법의 안정영역(Region of Absolute Stability)

한 단계 방법으로 ODE를 풀이할 때, 방법 자체가 가지는 안정영역(absolute stability region)이 문제가 된다. 단순화된 스칼라 방정식

$$
y' = \lambda y
$$

에 대해 수치 방법이 적용될 때, $\lambda$의 실제 부분이 음수이면 연속해(정확해)는 $t$가 증가함에 따라 0으로 수렴한다. 하지만 명시적 방법은 고유치가 이론적으로 안정 영역 바깥에 위치하면, 반올림 오차가 증폭되어 결국 발산할 수 있다. 전진 오일러의 경우

$$
y\_{n+1} = y\_n + \Delta t , \lambda , y\_n = (1 + \Delta t ,\lambda) , y\_n
$$

이므로, $|1 + \Delta t ,\lambda|<1$이 되어야 수치해가 안정적으로 0에 수렴한다. 따라서 $\lambda<0$인 상황에서 $\Delta t$가 충분히 작아야 $|1+\Delta t,\lambda|<1$을 만족할 수 있다. 즉, $\Delta t$에 대한 제약이 커질수록 계산 비용이 증가하고, 문제 크기가 커지면 이러한 제약은 수치적 불안정성을 막기 위해 더욱 엄격해진다.

#### 암시적 방법의 안정성과 계산 비용

암시적(implicit) 방법은 한 단계마다 비교적 더 많은 연산을 필요로 하지만, 일반적으로 명시적 방법에 비해 훨씬 넓은 안정영역을 가진다. 예를 들어 후진 오일러(Backward Euler) 방법은 다음과 같은 형태를 가지며

$$
y\_{n+1} = y\_n + \Delta t ,\lambda,y\_{n+1},
$$

결국 이를 $y\_{n+1}$에 대해 풀면

$$
y\_{n+1} = \frac{y\_n}{1 - \Delta t,\lambda}.
$$

이때 안정조건은 $|1-\Delta t,\lambda|>0$ 정도에 불과하여(즉, 극단적으로 $\Delta t,\lambda$가 실수 음수라면 큰 값을 가지는 한이 있어도 발산을 초래하지 않는다) 명시적 방법보다 상대적으로 훨씬 안정적이다. 그러나 실제 다변수 계에서 암시적 방법은 매 스텝마다 선형시스템 $\bigl(\mathbf{I} - \Delta t,\mathbf{J}\bigr)\Delta \mathbf{y} = \cdots$ 형태를 풀어야 하므로, $\mathbf{J}$가 큰 차원이라면 계산 비용이 상당하다. 강직성 ODE를 다루는 경우에는 보통 암시적 혹은 반-암시적(semi-implicit) 방법을 사용하며, 이 과정에서의 수치 연산 불안정성을 줄이기 위해서는 적절한 피벗팅, 사전조건화(Preconditioning) 등을 적용해야 한다.

#### 예시: 간단한 강직성 시스템

간단한 선형 강직성 계를 한 번 살펴보자.

$$
\begin{cases} y\_1' = -1000,y\_1 \ y\_2' = -y\_2 \end{cases}
$$

초기 조건을 $y\_1(0) = 1,\ y\_2(0)=1$로 두면, $y\_1(t)$는 매우 빠르게 0으로 수렴하고, $y\_2(t)$는 상대적으로 완만하게 감쇠한다. 만약 전진 오일러 방법으로 동일한 $\Delta t$를 적용하면, $y\_1$의 급격한 감쇠 모드에 의해 $\Delta t$를 극도로 작게 선택해야 안정적으로 수렴할 수 있다. $\Delta t$가 조금만 커도 $y\_1$이 진동하거나 발산하는 결과가 나오고, 이는 실제 물리해와 전혀 다른 거동이다.

간단한 Python 코드를 통해 전진 오일러와 후진 오일러의 차이를 시연할 수 있다.

```python
import numpy as np
import matplotlib.pyplot as plt

def stiff_ode_explicit(dt, t_end):
    t_vals = []
    y_vals = []
    y = np.array([1.0, 1.0])
    lam1, lam2 = -1000.0, -1.0
    t = 0.0
    while t < t_end:
        t_vals.append(t)
        y_vals.append(y.copy())
        # 전진 오일러
        y = y + dt * np.array([lam1*y[0], lam2*y[1]])
        t += dt
    return np.array(t_vals), np.array(y_vals)

def stiff_ode_implicit(dt, t_end):
    t_vals = []
    y_vals = []
    y = np.array([1.0, 1.0])
    lam1, lam2 = -1000.0, -1.0
    t = 0.0
    while t < t_end:
        t_vals.append(t)
        y_vals.append(y.copy())
        # 후진 오일러
        # y_{n+1} = y_n + dt * lam * y_{n+1}  -->  y_{n+1} - dt*lam*y_{n+1} = y_n
        # (1 - dt*lam)*y_{n+1} = y_n  -->  y_{n+1} = y_n / (1 - dt*lam)
        y[0] = y[0] / (1.0 - dt*lam1)
        y[1] = y[1] / (1.0 - dt*lam2)
        t += dt
    return np.array(t_vals), np.array(y_vals)

if __name__ == "__main__":
    dt = 0.001
    t_end = 0.1
    t_exp, y_exp = stiff_ode_explicit(dt, t_end)
    t_imp, y_imp = stiff_ode_implicit(dt, t_end)

    plt.plot(t_exp, y_exp[:,0], label="Explicit y1")
    plt.plot(t_exp, y_exp[:,1], label="Explicit y2")
    plt.plot(t_imp, y_imp[:,0], '--', label="Implicit y1")
    plt.plot(t_imp, y_imp[:,1], '--', label="Implicit y2")
    plt.legend()
    plt.show()
```

이 예시를 실행해보면, 전진 오일러의 경우 $\Delta t$가 조금만 커져도 $y\_1$에서 수치 폭발이 발생할 위험이 훨씬 크다. 후진 오일러는 비교적 큰 $\Delta t$에서도 감쇠 거동이 유지되어 안정적으로 0으로 수렴한다. 이는 강직성 문제에서 명시적 방법을 무턱대고 적용하면 큰 수치적 불안정성을 겪게 됨을 보여준다.

#### 고차 ODE 방법과 안정성

룽게-쿠타(Runge-Kutta) 계열의 고차 정확도 방법들도, 계 자체가 강직성을 가지면 명시적 방식을 적용하기 어렵다. 다단계(Multistep) 방법이나 가변 스텝(variable step) 전략을 도입해 해결하기도 하지만, 결국 안정성이 확보되지 않으면 스텝 크기를 과도하게 줄이는 결과로 이어진다. 암시적 룽게-쿠타, 암시적 다단계 방법(BDF: Backward Differentiation Formula) 등은 안정성이 높아서 강직성 문제에 적합하지만, 각 스텝에서 비선형 시스템을 풀거나 큰 행렬 연산을 동반하기 때문에 계산 비용이 만만치 않다. 이 과정에서의 수치적 불안정성은 대부분 앞서 논의한 선형시스템이나 비선형 방정식 풀이 과정에서 기인하므로, 철저한 전처리나 피벗팅, 혹은 문제 스케일링 작업이 함께 이뤄져야 한다.

#### 부분미분방정식(PDE)에서의 수치적 불안정

편미분방정식(PDE)을 수치적으로 풀 때도 불안정성이 부각된다. 대표적으로 열방정식(Heat Equation) $u\_t = \alpha u\_{xx}$을 명시적 유한차분(Explicit Finite Difference)으로 이산화하면, 시간 스텝 크기와 공간 격자 간격이 특정한 안정조건(Courant-Friedrichs-Lewy, CFL 조건)을 만족해야 한다. 이를 무시하고 스텝을 크게 잡으면 수치해가 발산하거나 가짜 진동이 발생한다. 반면 암시적 유한차분(Implicit Finite Difference) 또는 Crank-Nicolson 방법을 쓰면 안정 범위가 넓어지지만, 매 단계 선형시스템을 풀어야 하므로 계산량이 증가한다.

PDE에서의 불안정성은 미분방정식 자체의 물리적 성질(확산, 전파, 반응-확산 등)과 격자 스케일, 시간 스케일의 부적절한 결합 때문에 발생한다. 동시에, 큰 차원의 선형시스템 또는 비선형 방정식을 반복해서 풀어야 하므로, 알고리즘적 불안정성이 더욱 쉽게 누적된다. 이를 피하기 위해서는 주어진 PDE의 물리적·수학적 해석을 고려하여 적절한 이산화 방법과 수치 방법(명시적 vs 암시적 vs 혼합형), 스케일링, 필수적인 전처리 등을 종합적으로 적용해야 한다.

\--- 대신, 실무 사례에서의 수치적 불안정 대처

수치적 불안정성은 다양한 분야에서 실제 문제를 풀 때 본질적으로 나타난다. 예컨대 전산 유체역학(CFD), 전기회로 해석, 구조해석(FEM), 영상처리, 신호처리 등에서 대규모의 선형/비선형 연산이 반복적으로 수행된다. 이때 문제 자체가 큰 조건수를 가지고 있으면, 아무리 정밀도가 높은 계산 환경을 사용해도 오차가 급격히 확대되기 쉽다. 반면 문제는 비교적 잘 조건화되어 있어도, 알고리즘 선택이 부적절하거나 부동소수점 제한을 고려하지 않으면 불안정한 결과를 얻을 수 있다.

수치적 불안정성을 완화하기 위해 실제로는 다음과 같은 처리를 자주 고려한다.

#### 문제 스케일링(Scaling)과 전처리(Preconditioning)

고려 대상인 방정식이나 함수가 극단적으로 큰 계수나 아주 작은 계수를 포함한다면, 알맞게 스케일링을 실시하여 알고리즘이 다루는 수의 범위를 균등하게 만든다. 선형시스템 풀이에서는 전처리(Preconditioning) 기법을 통해 $\mathbf{A}\mathbf{x} = \mathbf{b}$를 $\mathbf{M}^{-1}\mathbf{A}\mathbf{x} = \mathbf{M}^{-1}\mathbf{b}$의 형태로 변환함으로써, $\mathbf{M}^{-1}\mathbf{A}$의 조건수가 크게 낮아지도록 설계한다. 이때 $\mathbf{M}$은 적절히 선택된 전처리 행렬로, $\mathbf{A}$와 비슷한 구조를 가지면서도 역행렬을 구하기 쉽거나, 최소한 해당 시스템을 빠르게 풀기 위한 도구를 제공한다.

#### 알고리즘의 선택과 단계별 오류 확인

선형대수 문제에서는 LU 분해, QR 분해, SVD 등 다양한 분해법이 존재하고, 각각의 안정성과 계산 비용이 다르다. Gauss 소거법과 LU 분해가 후방안정적인 것으로 잘 알려져 있지만, ill-conditioned 문제에서는 QR나 SVD가 더 좋은 안정성을 보이기도 한다. 비선형 문제에서는 뉴턴 방법 외에도 다양한 로버스트(robust) 기법을 고려하고, 중간 단계마다 잔차나 에너지 함수 등을 관측하여 수렴 거동을 주시한다. 이처럼 각 단계별로 “오차가 갑자기 커지지 않는지”, “잔차가 의도와 달리 증가하고 있지 않은지” 등을 실시간 모니터링해 문제를 조기 발견하는 방식이 효과적이다.

#### 혼합정밀도(Mixed Precision) 기법

최근 컴퓨팅 환경에서는 단정밀도(float32), 반정밀도(float16) 등을 활용하여 연산 속도를 높이되, 최종적으로 필요한 단계에서만 배정밀도(float64) 연산을 사용해 정확도를 보장하는 혼합정밀도 방법이 널리 사용된다. 이는 GPU나 특수 하드웨어에서 매우 유용하지만, 그만큼 수치적 불안정성을 방지하기 위한 후방안정 알고리즘이나 적절한 반복보정(iterative refinement) 기법이 전제되어야 한다. 잘못된 곳에서 저정밀도를 사용하면 오차가 걷잡을 수 없이 증폭될 수 있기 때문이다.

#### 병렬 및 분산 계산에서의 주의점

대규모 문제를 슈퍼컴퓨팅 환경에서 병렬화, 분산화하여 푸는 일이 일반화됨에 따라, 연산 순서가 매 실행마다 달라지는 문제가 발생할 수 있다. 특히 덧셈과 뺄셈의 순서가 달라지면 부동소수점 합산 결과가 달라져서, 최종 오차가 다르게 나타난다. 병렬 계산에서는 그 순서를 완전히 고정하기 어렵기 때문에, 알고리즘 자체가 후방안정적이어야 하고, 부분 합산 과정을 철저히 관리할 필요가 있다. 예를 들어 트리(tree) 구조로 합산을 수행하거나 Kahan Summation 알고리즘 등을 적용하여, 누적 오차가 커지지 않도록 주의한다.

#### 잘못된 해석을 피하기 위한 추가 검증

수치 결과만 의존해 엔지니어링적 결론을 도출할 때, 결과가 갑자기 폭증하거나 물리적으로 의미 없는 값을 보일 경우, 단순히 “모델이 불안정하다”라고 판단하기 전에 먼저 수치 알고리즘의 불안정성과 문제 조건수를 의심해볼 필요가 있다. 필요하다면 계수 행렬이나 중간 계산을 재점검하고, 독립적인 풀(풀이) 방식으로 검증 또는 근사적 유효범위를 확인하는 절차를 거치는 것이 안전하다.