의사스펙트럼법

편집 역사 토론
수치해석 계산물리 마지막 수정: 2026-08-12 04:49:37

1. 개요[편집]

의사스펙트럼법
Pseudospectral method
별칭유사스펙트럴법 · 스펙트럴 배점법(collocation)
미지수계수 ûk 가 아니라 격자점 값 u(xj)
미분미분행렬 D를 곱한다 (푸리에면 FFT로 O(N log N))
고질병에일리어싱 → 3/2 규칙 (= 2/3 절단)
딴 동네 응용최적제어의 LGL/LGR 배점 + 공변량 사상 정리
이름 충돌의사스펙트럼(pseudospectrum)과 무관

스펙트럴 방법은 계수로 생각하고, 의사스펙트럼법은 격자점 값으로 생각한다. 수학은 같은데 코드가 다르고, 코드가 다르니 버그도 다르다.

의사스펙트럼법(pseudospectral method)은 스펙트럴 방법배점(collocation) 형태로 구현한 것이다. 해를 전역 기저의 계수 u^k\hat{u}_k 로 들고 다니는 대신 배점 노드에서의 값 uj=u(xj)u_j = u(x_j) 로 들고 다니고, 미분은 그 값 벡터에 미분행렬 DD 를 곱하는 행위로, 방정식은 각 노드에서 잔차가 정확히 0이 되라는 조건으로 표현한다.

u(xj)k=0NDjkuk,Lu(xj)=f(xj)  (j=0,,N)u'(x_j) \approx \sum_{k=0}^{N} D_{jk}\,u_k, \qquad \mathcal{L}u(x_j) = f(x_j)\ \ (j = 0,\dots,N)

같은 근사공간을 쓰므로 정확도는 갈레르킨 판본과 본질적으로 같은 지수 수렴이다. 차이는 계산 경로에 있고, 그 차이가 실전에서 전부다. 비선형항 uxuu\,\partial_x u 를 계수 공간에서 다루면 합성곱이 되어 O(N2)O(N^{2}) 이 들지만, 값 공간에서는 성분별 곱셈 한 줄이다. 그래서 나비에-스토크스 방정식이 관여하는 순간 의사스펙트럼법이 사실상 유일한 선택지가 된다.

이름부터 정리하고 가자. 이 문서의 주제는 의사스펙트럼(pseudospectrum, 비정규 행렬의 ε\varepsilon-섭동 고유값 집합 Λε(A)\Lambda_\varepsilon(A))과 아무 관계가 없다. 영어로도 pseudospectral method와 pseudospectrum으로 서로 다른 단어이고, 한국어 번역이 겹쳐 버린 것뿐이다.1

2. 미분행렬 — 이 방법의 전부[편집]

주기 문제라면 푸리에 미분이다. u^kiku^k\hat{u}_k \to ik\hat{u}_k고속 푸리에 변환으로 왕복하면 O(NlogN)O(N\log N) 이고, 대응하는 DD 는 순환 토플리츠 행렬이라 명시적으로 만들 필요조차 없다.

비주기 문제의 표준은 체비셰프 미분행렬이다. 노드는 가우스-로바토(Chebyshev–Gauss–Lobatto) 점

xj=cosjπN,j=0,1,,Nx_j = \cos\frac{j\pi}{N}, \qquad j = 0,1,\dots,N

이고 — 즉 x0=1x_0 = 1 에서 xN=1x_N = -1 까지 내려오며 양 끝에 몰려 있다 — 이 점들에서 보간 다항식을 만들어 미분한 결과가 다음 (N+1)×(N+1)(N+1)\times(N+1) 밀집 행렬이다.

D00=2N2+16,DNN=2N2+16D_{00} = \frac{2N^{2}+1}{6}, \qquad D_{NN} = -\frac{2N^{2}+1}{6} Djj=xj2(1xj2) (1jN1),Dij=cicj(1)i+jxixj (ij)D_{jj} = \frac{-x_j}{2(1-x_j^{2})}\ (1\le j\le N-1), \qquad D_{ij} = \frac{c_i}{c_j}\frac{(-1)^{i+j}}{x_i - x_j}\ (i\ne j)

여기서 c0=cN=2c_0 = c_N = 2, 나머지는 1이다. 노드가 경계에 몰려 있는 것은 보간과 근사의 룽게 현상을 피하기 위한 것이며, 이 분포가 곧 체비셰프 다항식의 극점 분포다.

구현에서 반드시 알아야 할 것 두 가지.

  • 음합 트릭(negative sum trick). DD 는 상수함수를 0으로 보내야 하므로 각 행의 합이 정확히 0이어야 한다. 그런데 위 공식을 그대로 쓰면 NN 이 커질 때 대각 성분이 O(N2)O(N^{2}) 로 크고 비대각 성분들의 소거가 심해서 반올림 오차가 눈에 띄게 쌓인다. 대각을 공식 대신 Dii=jiDijD_{ii} = -\sum_{j\ne i}D_{ij}계산해서 넣으면 상수함수가 기계 정밀도로 소거되고 정확도가 눈에 띄게 좋아진다. 교과서 공식을 그대로 옮겨 적은 코드와 이 한 줄을 넣은 코드의 차이가 큰 NN 에서 유효숫자 몇 자리로 나타난다.
  • DD 는 꽉 차 있고, 거듭제곱은 더 나쁘다. 2계 미분은 D2D^{2} 로 얻는데, DD 를 두 번 곱하면 오차도 함께 증폭된다. 정확도가 중요하면 2계 미분행렬을 별도 공식으로 직접 생성하는 편이 낫다.

경계조건 처리는 놀랄 만큼 단순하다. 디리클레 조건 u(±1)=0u(\pm 1)=0 이면 첫 행·열과 마지막 행·열을 그냥 잘라낸다. 노드값이 곧 미지수이므로 경계값을 고정하는 것이 행을 지우는 것과 같기 때문이다. 계수 공간에서 경계조건을 다루는 타우(tau) 방법과 비교하면 압도적으로 간편하고, 이것이 배점 형태가 널리 쓰이는 실용적 이유 중 하나다.

3. 격자가 만드는 강성[편집]

체비셰프 격자의 최소 간격은 경계 근처에서 Δxminπ2/(2N2)\Delta x_{\min} \approx \pi^{2}/(2N^{2}) 로, 중심부보다 NN 배 촘촘하다. 이 사실이 시간 적분에 곧장 청구서로 온다.

격자1계 미분행렬 고유값 크기2계 미분행렬명시적 시간 적분 제약
푸리에 (균등)O(N)O(N)O(N2)O(N^{2})확산이면 Δt=O(N2)\Delta t = O(N^{-2})
체비셰프 (로바토)O(N2)O(N^{2})O(N4)O(N^{4})확산이면 Δt=O(N4)\Delta t = O(N^{-4})

N=64N=64 짜리 체비셰프 격자에 열방정식을 명시적으로 풀면 Δt107\Delta t \sim 10^{-7} 을 요구받는다. 격자점은 64개인데. 즉 의사스펙트럼 공간 이산화는 태생적으로 강성 방정식을 만들어 낸다. 그래서 실무에서는 확산항을 암시적 룽게-쿠타법이나 후진 차분으로, 대류·비선형항을 명시적으로 다루는 IMEX 조합이 표준이 된다. DD 가 밀집 행렬이라 암시적 풀이 비용이 O(N3)O(N^{3}) 인 것은 감수한다 — 어차피 NN 이 작으니까. “격자를 적게 쓴다”는 스펙트럴의 장점이 “행렬이 꽉 찬다”는 단점을 상쇄하는 구조가 여기서 다시 확인된다.

4. 에일리어싱과 3/2 규칙[편집]

의사스펙트럼법의 대표적 함정. 파수 NN 까지 담는 격자에서 두 신호를 곱하면 최대 2N2N 까지 파수가 생기는데, 격자가 표현 못 하는 고파수 성분이 저파수로 접혀 들어와(alias) 결과를 오염시킨다.

정확한 처방은 격자를 늘리는 것이다. 유지하려는 모드가 NN 개일 때, 이차 비선형항을 계산할 격자점 수 MM

M  32NM \ \ge\ \tfrac{3}{2}N

이면 접힘이 유효 대역 안으로 들어오지 못한다. 절차는 이렇다 — 계수를 3N/23N/2 개로 영 채움(zero padding) 하고, 역변환해 값 공간에서 곱하고, 다시 변환한 뒤 상위 성분을 버린다. 이것이 3/2 규칙이고, 같은 일을 반대편에서 서술하면 ”MM 점 격자에서 상위 1/3 파수를 잘라 낸다”는 2/3 규칙이 된다. 둘은 같은 처방을 패딩 쪽에서 보느냐 절단 쪽에서 보느냐의 차이일 뿐이다. 삼차 비선형항이면 M2NM \ge 2N 이 필요하다. 오르사그(Steven Orszag)가 1971년에 정리했다.

대안으로 위상 이동 제거법(phase-shift dealiasing)이 있다. 격자를 반 칸 어긋나게 이동시켜 한 번 더 계산하고 두 결과를 평균하면 에일리어싱 항의 부호가 반대라 상쇄된다. 격자를 안 늘려도 되지만 변환을 두 배로 해야 하므로 실제 이득은 문제에 따라 갈리고, 이동량을 매 스텝 무작위로 잡으면 에일리어싱 오차가 편향 없이 쌓이게 만들 수 있다.

한 가지 흔한 오해를 짚자. 에일리어싱 제거는 대개 정확도가 아니라 안정성의 문제다. 해가 충분히 매끄러우면 에일리어싱 오차의 점근 차수는 절단 오차와 같아서, 안 지워도 수렴 차수는 그대로다. 문제는 지우지 않았을 때 고파수 끝에 에너지가 쌓여 비선형 불안정으로 폭발하는 경우가 있다는 것이다. 즉 “3/2 규칙을 켰더니 정확해졌다”보다 “3/2 규칙을 껐더니 며칠 뒤에 터졌다” 쪽이 현장의 실제 경험담이다.2

5. 정확도의 천장 — N을 계속 키울 수 없는 이유[편집]

지수 수렴이라니 NN 을 키우면 계속 좋아질 것 같지만, 체비셰프 배점에는 분명한 천장이 있다. 2계 미분행렬 D2D^{2} 의 성분이 O(N4)O(N^{4}) 이므로, 이걸로 세운 선형계의 조건수O(N4)O(N^{4}) 급이다. 반올림 오차는 대략 εmachN4\varepsilon_{\mathrm{mach}}\cdot N^{4} 로 커지고, 배정밀도에서 N200N \approx 200 이면 이미 1016×1.6×10910710^{-16}\times 1.6\times10^{9} \sim 10^{-7} 이다.

그래서 전형적인 오차 곡선은 이렇게 생겼다. NN 을 늘리면 처음에는 지수적으로 떨어지다가, 절단 오차가 반올림 바닥에 닿는 순간 최소점을 찍고 다시 올라간다. 최적 NN 이 존재하고, 그 지점의 정확도가 그 문제에서 얻을 수 있는 전부다. 이것은 알고리즘 결함이 아니라 밀집·고차 연산자에 붙는 원가이며, 유한차분법이 같은 정확도에 훨씬 큰 NN 을 요구하는 대신 조건수는 O(N2)O(N^{2}) 로 얌전한 것과 정확히 반대편의 거래다.

단일 영역에서 NN 을 무한정 키우는 대신 영역을 쪼개고 각 영역의 NN 을 적당히 유지하는 것이 현실적인 처방이고, 그것이 곧 스펙트럴 요소법이다. “고차가 좋다면 무한정 고차로 가자”가 안 되는 이유가 여기 있다.

6. 갈레르킨과의 관계, 그리고 요소로 쪼개기[편집]

갤러킨 방법은 잔차가 모든 기저함수와 직교하도록 요구하고, 배점법은 잔차가 특정 점들에서 0이 되도록 요구한다. 후자는 “델타함수를 시험함수로 쓴 가중잔차법”으로 볼 수 있고, 노드를 구적점에 놓으면 두 정식화가 사실상 같은 이산 시스템을 준다(질량행렬이 대각이 되는 이른바 질량 집중 효과). 르장드르 다항식 기반 스펙트럴 요소법이 LGL 노드를 쓰는 이유가 정확히 이것이다 — 배점의 편의와 갈레르킨의 안정성을 동시에 가져간다.

전역 기저의 한계(복잡 형상 불가, 밀집 행렬, 전역 통신)는 스펙트럴 방법 문서에 정리되어 있고, 그 해법인 스펙트럴 요소법·다영역 분할도 마찬가지다. 불연속이 있으면 깁스 현상이 나오는 것도 여기서 그대로 상속된다.

7. 최적제어의 의사스펙트럼법 — 완전히 다른 무대, 같은 도구[편집]

1990년대 중반부터 최적 제어 문제를 푸는 직접법으로 같은 도구가 쓰이기 시작했다. 발상은 미분방정식 풀이와 똑같다. 상태 x(t)x(t) 와 제어 u(t)u(t) 를 배점 노드에서의 값으로 두고, 동역학 구속조건을 미분행렬로 쓰고, 목적함수는 구적 가중치로 근사한다.

kDikxk=T2f(xi,ui),JT2iwig(xi,ui)\sum_{k} D_{ik}\,x_k = \tfrac{T}{2}\,f(x_i, u_i), \qquad J \approx \tfrac{T}{2}\sum_i w_i\,g(x_i,u_i)

결과는 유한차원 비선형계획(NLP)이고, 이를 SQP나 내점법으로 푼다. 노드 선택에 따라 세 갈래다.

  • 르장드르(LGL) 의사스펙트럼법 — 엘나가르 등(1995). 양 끝점을 포함하는 로바토 노드. 가장 먼저 정립되었고 구현이 직관적이다.
  • 가우스(LG) 의사스펙트럼법 — 끝점을 배점하지 않고 별도 변수로 둔다.
  • Radau(LGR) 의사스펙트럼법 — 한쪽 끝만 포함. 초기값 문제와 궁합이 좋고 공변량 추정이 깨끗해서 최근 구현(GPOPS-II 등)의 주류다.

여기서 이 분야의 핵심 정리가 등장한다. 공변량 사상 정리(covector mapping theorem, 로스–파루). 문제는 이것이다 — 폰트랴긴 최대 원리를 먼저 적용해 얻은 최적성 필요조건(공변량, 즉 수반변수에 대한 미분방정식)을 이산화한 결과와, 먼저 이산화해서 얻은 NLP의 KKT 조건이 일반적으로 일치하지 않는다. “이산화 후 최적화”와 “최적화 후 이산화”가 교환되지 않는다는 것이다.

공변량 사상 정리는 그 간극을 메우는 명시적 사상을 준다. 표준 LGL 정식화에서는 NLP 승수 λ~i\tilde{\lambda}_i 를 구적 가중치로 나눈 λi=λ~i/wi\lambda_i = \tilde{\lambda}_i / w_i 가 이산 공변량에 대응하며, 여기에 끝점에서의 폐합 조건(closure condition)이 추가로 필요하다. 이 사상 덕분에 NLP 풀이 결과에서 곧바로 공변량 추정치를 뽑아 폰트랴긴 조건을 얼마나 만족하는지 검증할 수 있다. 직접법을 쓰면서 간접법의 최적성 증거를 챙기는 셈이고, 이는 실무에서 “이 궤적이 정말 최적인가”를 따질 때 결정적이다.3

이 계보의 유명한 실적이 국제우주정거장의 무추진제 회전 기동(Zero-Propellant Maneuver)이다. 2006~2007년에 CMG(제어 모멘트 자이로)만으로 정거장을 대각도 회전시키는 궤적을 의사스펙트럼 최적제어로 설계해 실제 비행에 적용했고, 추진제를 쓰지 않고 자세 변경을 끝냈다. 스펙트럴 방법이 유체 격자를 벗어나 궤도상에서 값을 한 드문 사례다.4

물론 만능은 아니다. 최적해의 제어 입력이 뱅뱅(bang-bang) 형태로 불연속이면 다항식 근사가 깁스 현상을 일으켜 노드를 늘려도 수렴이 대수적으로 느려진다. 처방은 스위칭 시각을 미지수로 두고 구간을 쪼개는 다구간(multi-phase) 정식화이며, 이 지점에서 의사스펙트럼법은 다시 궤적 최적화 일반론의 문제와 만난다.

8. 관련 문서[편집]

9. Footnotes[편집]

  1. 더 얄궂은 것은 두 주제의 대표적 교과서를 같은 사람이 썼다는 사실이다. 로이드 트레페텐의 Spectral Methods in MATLAB(2000)과 Spectra and Pseudospectra(2005, 엠브리 공저). 그러니 “트레페텐 책 보세요”는 이 동네에서 답이 될 수 없는 조언이다.

  2. 그래서 대형 DNS 코드의 옵션 목록에는 dealiasing 스위치가 거의 항상 있고, 논문에는 “fully dealiased using the 3/2 rule”이 한 줄 박혀 있다. 비용이 대략 (3/2)33.4(3/2)^{3} \approx 3.4 배(3차원)라 아깝지만, 며칠 돌린 계산이 NaN으로 끝나는 것보다는 싸다. 수렴은 신에게 맡기더라도 에일리어싱은 본인이 처리하는 게 맞다.

  3. 직접법 진영이 오랫동안 들었던 비판이 “그래서 그게 진짜 최적이냐, 그냥 NLP 극값 아니냐”였다. 공변량 사상 정리는 그 비판에 대해 “승수를 이렇게 변환하면 폰트랴긴 조건과 대조할 수 있다”고 답한 것이다. 그리고 실제로 대조해 보면 안 맞는 경우가 꽤 나온다 — 그게 이 정리의 쓸모다.

  4. 우주정거장에서 추진제는 화물선으로 올려야 하는 소모품이고, 대각도 자세 변경은 원래 스러스터를 써야 하는 기동이다. 이걸 “중력경사·공력 토크에서 각운동량을 벌어” CMG만으로 해결한다는 발상 자체가 최적제어의 승리에 가깝다. 궤적은 지상에서 미리 최적화해 올려보냈고, 온보드에서 실시간으로 푼 것이 아니다 — 의사스펙트럼 NLP는 아직 그렇게 싸지 않다.