수치 연속법

1. 개요[편집]

수치 연속법 Numerical continuation
대상	$F(u,\mu)=0$, $F:\mathbb{R}^n\times\mathbb{R}\to\mathbb{R}^n$
핵심 기법	의사호장 연속법 (Keller, 1977)
구조	접선 예측자 + 뉴턴 보정자
넘는 것	접힘점(fold) · 분기점 · 가지 전환
도구	AUTO-07p · MatCont · COCO · LOCA

해가 하나뿐인 곳까지는 뉴턴이면 충분하다. 문제는 해가 접혀서 되돌아오는 지점이다.

수치 연속법(numerical continuation)은 파라미터가 들어간 비선형 방정식 $F(u,\mu)=0$ 의 해집합을, 파라미터 $\mu$ 를 조금씩 바꿔가며 곡선(해 가지, solution branch)째로 추적하는 기법이다. 한 점의 해를 구하는 것이 뉴턴-랩슨법의 일이라면, 연속법은 그 해가 파라미터에 따라 어떻게 이동하고, 어디서 사라지고, 어디서 두 갈래로 갈라지는지를 그린다. 동역학계의 정상해 다이어그램, 분기 이론의 실전 계산 도구, 그리고 구조해석의 후좌굴 경로 추적이 전부 같은 알고리즘 위에 서 있다.

암묵 함수 정리가 보장하는 것은 딱 여기까지다. $F_u = \partial F/\partial u$ 가 정칙인 정칙점(regular point) 근방에서는 $u = u(\mu)$ 가 유일하게 잘 정의된다. 연속법의 모든 재미는 이 조건이 깨지는 곳에서 시작한다.

2. 자연 파라미터 연속법이 반드시 실패하는 곳[편집]

가장 소박한 방법은 자연 파라미터 연속법(natural parameter continuation)이다. $\mu_{k+1} = \mu_k + \Delta\mu$ 로 파라미터를 한 칸 올리고, 직전 해 $u_k$ 를 초기 추정으로 삼아 $F(\cdot,\mu_{k+1})=0$ 을 뉴턴으로 푼다. 뉴턴 스텝은

$$F_u(u,\mu_{k+1})\,\delta u = -F(u,\mu_{k+1})$$

이고, 코드 열 줄이면 된다. 잘 굴러가다가, 접힘점(fold point, 극한점·한계점이라고도 한다)에서 예외 없이 죽는다.

접힘점은 곡선이 $\mu$ 축 방향으로 되돌아가는 지점이다. 여기서 $F_u$ 는 계수(rank)가 $n-1$ 로 떨어져 특이해지고, 동시에 $F_\mu$ 가 $F_u$ 의 치역 밖에 놓인다. 즉 뉴턴이 풀어야 할 선형계 자체가 풀리지 않는다. 게다가 접힘점 너머의 $\mu$ 에는 그 근방에 해가 아예 존재하지 않는다. 발산은 알고리즘의 결함이 아니라, 있지도 않은 것을 찾으라고 시킨 결과다.¹

여기서 “그럼 $u$ 의 한 성분을 파라미터로 삼자”는 처방이 나오고(변위 제어), 실제로 접힘점 하나는 넘긴다. 하지만 곡선이 그 성분 방향으로도 되돌아 접히는 순간 똑같이 무너진다. 어느 한 좌표축을 특권화하는 방법은 전부 같은 방식으로 배신한다.

3. 의사호장 연속법[편집]

Keller(1977)의 해법은 관점을 바꾼다. $\mu$ 를 독립변수 자리에서 끌어내려 $u$ 와 대등한 미지수로 승격시키고, 대신 곡선을 따라 잰 거리(호장, arclength)를 새 독립변수로 삼는다. 미지수가 $n+1$ 개가 됐으니 방정식도 하나 더 필요하다. 그 자리를 채우는 것이 의사호장 구속이다.

$$N(u,\mu) = \dot{u}_0^{\top}(u - u_0) + \dot{\mu}_0(\mu - \mu_0) - \Delta s = 0$$

여기서 $(\dot u_0, \dot\mu_0)$ 는 직전 해점에서의 단위 접선벡터, $\Delta s$ 는 보폭이다. 기하적으로는 “접선 방향으로 $\Delta s$ 만큼 간 지점에 세운 초평면 위에서 해를 찾아라”는 뜻이다. 호장법의 원형이 원(球)을 쓰는 것과 달리 구속이 선형이라 야코비 구조가 깔끔하다.

핵심은 확장 야코비다.

$$J_{\text{ext}} = \begin{bmatrix} F_u & F_\mu \\ \dot{u}_0^{\top} & \dot{\mu}_0 \end{bmatrix}$$

단순 접힘점에서 $F_u$ 가 특이해도 $J_{\text{ext}}$ 는 정칙이다. 증명은 짧다. 커널 원소 $(w,\nu)$ 가 있다면 $F_u w + \nu F_\mu = 0$ 인데, 여기에 $F_u$ 의 좌영벡터 $\psi$ 를 곱하면 $\nu\,\psi^{\top}F_\mu = 0$ 이고 접힘점 조건에서 $\psi^{\top}F_\mu \neq 0$ 이므로 $\nu = 0$. 그러면 $w$ 는 $F_u$ 의 영벡터 $\varphi$ 의 상수배인데, 접힘점에서 접선은 정확히 $(\varphi, 0)$ 방향이므로 마지막 행이 $\|\varphi\|^2 c = 0$ 을 강제해 $c=0$. 곡선은 특이점을 정칙 문제로 느끼며 그냥 돌아 나간다.

4. 예측자-보정자와 보폭 조정[편집]

실제 진행은 예측-보정 루프다.

1. 접선 계산. $F_u\dot u + F_\mu\dot\mu = 0$ 을 $J_{\text{ext}}$ 로 풀고 정규화한다. 부호는 직전 접선과의 내적이 양수가 되도록 고른다. 이 한 줄이 없으면 알고리즘은 접힘점에서 왔던 길을 되짚어 돌아가며 같은 구간을 무한 왕복한다.
2. 예측. $(u^{(0)},\mu^{(0)}) = (u_0,\mu_0) + \Delta s\,(\dot u_0,\dot\mu_0)$.
3. 보정. 확장계 $[F;\,N]=0$ 에 뉴턴을 돌린다. 보통 3~5회면 수렴한다.
4. 보폭 적응. 수렴 반복 횟수가 목표보다 적으면 $\Delta s$ 를 키우고, 많으면 줄이며, 실패하면 반으로 접고 재시도한다.

접선 대신 직전 두 해점을 잇는 할선(secant)을 예측자로 쓰는 변형도 흔하다. 접선을 얻으려면 선형계를 한 번 더 풀어야 하는데 할선은 뺄셈 한 번이면 되기 때문이다. 정확도는 접선 쪽이 높지만 보정자가 어차피 뉴턴이라 실전 차이는 크지 않고, 대규모 문제일수록 할선이 이득이다. 반대로 보폭이 지나치게 크면 예측점이 인접한 다른 가지의 수렴 영역으로 넘어가 버려, 뉴턴이 멀쩡히 수렴했는데 정작 다른 곡선 위에 서 있는 사고가 난다. 곡선 위를 걷는 알고리즘이 조용히 옆길로 새는 이 실패 모드가 연속법 디버깅에서 가장 잡기 어려운 부류다.

$n$ 이 큰 문제(이산화된 PDE)에서는 $J_{\text{ext}}$ 를 통째로 다루지 않고, 경계 테두리 제거(bordering) 기법으로 기존 $F_u$ 희소 solver를 두 번 호출해 처리한다. 다만 접힘점 근처에서 $F_u$ 가 특이에 가까워 조건수가 폭발하므로, 대규모 문제에서는 테두리를 유지한 채 크리로프 부분공간법으로 확장계를 직접 푸는 쪽이 안전하다.

5. 분기점 검출과 가지 전환[편집]

곡선을 걷는 동안 각 스텝마다 테스트 함수를 계산해 부호 변화를 감시한다. 부호가 바뀌면 그 구간에 이분법/할선법으로 특이점을 정밀 위치시킨다.

• 접힘점 — $\det F_u$ 의 부호 변화(동등하게 접선의 $\dot\mu$ 성분이 0을 통과). $J_{\text{ext}}$ 는 정칙을 유지한다. → 안장-마디 분기.
• 분기점(branch point) — 이번엔 $\det J_{\text{ext}}$ 자체가 부호를 바꾼다. 접힘점과 분기점을 가르는 것이 정확히 이 차이다.
• 호프 분기 — 켤레복소 고유값 쌍 $\lambda = \pm i\omega$ 가 허수축을 지나는 사건. $\lambda_i + \lambda_j = 0$ 일 때만 특이해지는 이중교대곱(bialternate product) $2F_u \odot I$ 의 행렬식을 테스트 함수로 쓴다. 실수 쌍 $\pm\lambda$ 에서도 걸리므로 걸러내야 한다.

분기점에서는 $J_{\text{ext}}$ 의 커널이 2차원이 된다. 이 2차원 안에서 2차형식(대수적 분기 방정식)을 풀면 교차하는 두 가지의 접선 방향이 모두 나오고, 지금 걷던 방향이 아닌 쪽을 예측자로 삼아 뉴턴을 돌리는 것이 가지 전환(branch switching)이다. 대칭성이 있는 문제에서 대칭 파괴 가지를 찾는 표준 절차이며, 좌굴 문제에서 완전 형상의 자명한 해에서 좌굴 모드로 갈아타는 순간이 바로 이것이다.

주기해도 같은 틀에서 다룬다. 극한 순환은 미지 주기 $T$ 와 위상 고정 조건을 추가한 경계값 문제로 놓고 직교 배열법(orthogonal collocation)으로 이산화하면, 그대로 연속법의 $F(u,\mu)=0$ 형태가 된다. 여기서 더 나아가 분기 조건 자체를 방정식에 넣고 파라미터를 둘로 늘리면 여차원 1 분기 곡선을 2-파라미터 평면에서 추적할 수 있고, 그 곡선들이 만나는 곳에서 첨점·보그다노프-타켄스 같은 여차원 2 점이 나온다. 호모클리닉 궤도 추적도 같은 확장의 연장선이다.

6. 다른 얼굴들 — 호모토피, 호장법, 국소 이론[편집]

6.1. 호모토피 — 초기 추정 제조기[편집]

분기 다이어그램에 아무 관심이 없어도 연속법을 쓸 이유가 있다. 뉴턴에게 줄 초기 추정을 만들어내는 용도다. 풀고 싶은 어려운 문제 $F(u)=0$ 과 해를 이미 아는 쉬운 문제 $G(u)=0$ 를 잇는 호모토피

$$H(u,t) = t\,F(u) + (1-t)\,G(u) = 0$$

를 세우고 $t$ 를 $0$ 에서 $1$ 로 연속법으로 밀면, 각 스텝의 해가 다음 스텝의 초기 추정이 되어 뉴턴이 수렴 반경 안에 계속 머문다. CFD에서 레이놀즈수를 낮은 값부터 천천히 올리는 램핑, 비선형 재료 해석의 하중 증분, 반도체 소자 시뮬레이션의 바이어스 스텝핑이 전부 이름만 다른 같은 기법이다. 이때도 경로가 접히면 자연 파라미터 방식은 멈추므로 의사호장으로 갈아타는 것이 정석이다. 한편 $F$ 가 다항식계라면 베주 수만큼의 시작해에서 출발해 모든 복소 해를 빠짐없이 추적하는 다항식 호모토피 연속법이라는 별도의 분야가 있고, Bertini·HomotopyContinuation.jl 같은 전용 도구가 나와 있다.

6.2. 구조해석의 호장법[편집]

비선형 구조해석에서 쓰는 호장법(Riks·Crisfield)은 이름과 계보가 다를 뿐 사실상 같은 아이디어다. 잔차 $\mathbf{f}_{\text{int}}(\mathbf{u}) - \lambda\mathbf{q} = 0$ 에서 하중계수 $\lambda$ 가 여기의 $\mu$ 이고, 접선강성행렬이 $F_u$ 다. 얕은 아치의 스냅스루는 접힘점, 완전 기둥의 좌굴은 분기점, 좌굴 모드로 갈아타는 것은 가지 전환이다. 구속을 초평면으로 잡으면 Keller, 구(球)로 잡으면 Crisfield일 뿐이다.

6.3. 국소 이론과의 분업, 그리고 도구[편집]

역할 분담도 분명하다. 연속법은 전역적·수치적이다. 가지가 실제로 어디까지 뻗고 어디서 꺾이는지를 알려주지만, 그 분기가 어떤 종류인지는 스스로 말해주지 않는다. 반대로 중심 다양체 정리와 정규형 이론은 국소적·해석적이다. 분기점 근방에서 계를 최소 차원으로 줄이고 표준형 계수(예: 호프의 첫 랴푸노프 계수 부호 → 초임계/아임계 판정)를 준다. 실제 도구들은 둘을 붙여 쓴다 — MatCont는 연속법으로 분기점을 잡은 뒤 그 자리에서 정규형 계수를 수치적으로 계산해 뱉는다.

도구 지형은 대략 이렇다. AUTO-07p(Doedel)는 ODE/경계값 문제 분기 해석의 원조이자 여전히 사실상의 표준이고, MatCont는 MATLAB 환경에서 GUI와 정규형 계산까지 묶은 물건이며, COCO는 연속 문제를 조립식으로 정의하는 프레임워크, LOCA(Trilinos)는 이산화된 대규모 PDE를 병렬로 다루는 쪽이다.² 자코비안 행렬을 해석적으로 줄 수 있느냐가 실사용 품질을 크게 좌우한다.³

7. 관련 문서[편집]

• 분기 이론 · 호프 분기 · 안장-마디 분기 · 극한 순환
• 중심 다양체 정리 · 정규형 이론 · 동역학계
• 뉴턴-랩슨법 · 자코비안 행렬 · 조건수 · 크리로프 부분공간법
• 호장법 · 비선형 구조해석 · 좌굴

8. Footnotes[편집]

1. 그래서 상용 솔버가 “하중 증분에서 수렴 실패”를 뱉었을 때 반사적으로 증분을 잘게 쪼개는 습관은 절반만 맞다. 진짜 접힘점이라면 $\Delta\mu$ 를 $10^{-9}$ 로 줄여도 넘어가지 못한다. 넘어야 할 것은 보폭이 아니라 파라미터화 방식이다. ↩

2. AUTO는 1981년 Doedel의 학위논문에서 출발해 지금까지 살아 있는데, 코드 상당 부분이 여전히 Fortran이고 입력이 c. 로 시작하는 상수 파일이다. 처음 만지면 1980년대로 시간여행한 기분이 들지만, 배열법 기반 주기해 추적의 견고함만큼은 아직도 벤치마크 취급을 받는다. ↩

3. 유한차분으로 야코비를 때우면 스텝당 $n$ 번의 잔차 평가가 추가로 들 뿐 아니라, 접힘점 근처에서 미분 근사 오차가 특이성 판정을 흐려 테스트 함수의 부호 변화를 통째로 놓칠 수 있다. 이 바닥에서 자동 미분이 유난히 대접받는 이유다. ↩