유한차분 시간영역법(FDTD) 유도 과정
1. 서론
이 페이지에서는 시간영역에서 맥스웰 방정식을 해석하기 위해 사용되는 유한차분 시간영역법(FDTD) 기반의 주요 업데이트 방정식에 대한 수학적 유도 과정을 상세히 설명합니다. 이 페이지의 목적은 기본 이론 페이지와는 다릅니다. 여기서는 간략한 개요를 제시하는 것이 아니라, 연속적인 형태의 컬 방정식들이 어떻게 컴퓨터에서 처리 가능한 이산 대수식으로 변환되는지를 명확히 보여주는 데 중점을 둡니다. 또한 이 부분이 바로 코드 작성 시 오류가 발생하기 가장 쉬운 FDTD 방법의 영역이기도 합니다. 부호 오류, 인덱스 배치 실수, 혹은 시간 단계 간 불일치가 발생할 경우 즉시 수치적으로 불안정하거나 물리적으로 부적절한 결과가 나타나게 됩니다.
따라서 우리는 미분 형태의 암페어 법칙과 패러데이 법칙으로부터 출발하여 이를 직교좌표계 성분들로 전개한 뒤 2차원 광학 문제에 적용합니다. 공간 미분항들은 중심 유한차분법을 이용해 치환하고, 시간 미분항들은 표준적인 리프플로그 방식을 사용해 이산화합니다. 이후 도출된 식들을 재배열하여 전기장 및 자기장 성분들에 대한 명시적 업데이트 공식을 얻게 됩니다.
유도 과정을 명확히 하기 위해 여기서는 방정식을 2차원 경우로 단순화하였으며, 이때 \(\partial/\partial y = 0\)와 자기 전도도는 0으로 설정하였다 \(\sigma_m = 0\). 이를 통해 수식의 전개 과정을 보다 쉽게 이해할 수 있을 뿐만 아니라 실제로 차원 축소가 어떻게 이루어지는지도 명확히 보여준다. 다만 OghmaNano의 전체 FDTD 해석기에서는 전자기 방정식이 3차원 공간에서 해석된다. 따라서 여기서 제시된 2차원 유도 과정은 실제 해석기가 2차원에만 국한된다는 의미가 아니라, 해당 방법론을 명확히 설명하기 위한 예시로 활용된 것이다.
이 페이지에서는 전기장 성분들이 정수 단위의 시간 단계 \(n\Delta t\)에서 계산되는 반면, 자기장 성분들은 반정수 단위의 시간 단계 \(\left(n+\tfrac{1}{2}\right)\Delta t\)에서 계산됩니다. 이는 시간영역유한차분법(FDTD)에서 사용되는 전형적인 Yee 격자 기반의 점프프로그 구조입니다. 상대적으로 간략한 개요 페이지와는 달리, 여기서는 인덱싱 방식이 보다 명확하게 설명되어 있어 이 유도 과정을 수학적 설명뿐만 아니라 구현 가이드로도 활용할 수 있습니다.
2. 기초 방정식
우리는 전기전도도 \(\sigma\), 유전율 \(\epsilon\), 자기투자율 \(\mu\)를 갖는 선형 매질에서 전기장과 자기장에 관한 맥스웰의 회전 방정식으로부터 시작하며, 일반성을 위해 자기전도도 \(\sigma_m\)도 함께 고려합니다.
\[ \nabla \times \boldsymbol{H} = \sigma \boldsymbol{E} + \epsilon \frac{\partial \boldsymbol{E}}{\partial t}, \qquad \nabla \times \boldsymbol{E} = - \sigma_m \boldsymbol{H} - \mu \frac{\partial \boldsymbol{H}}{\partial t}. \]
행렬식 형태로 표현한 경우, 컬 연산식은 다음과 같다
\[ \nabla \times \boldsymbol{H} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \dfrac{\partial}{\partial x} & \dfrac{\partial}{\partial y} & \dfrac{\partial}{\partial z} \\ H_x & H_y & H_z \end{vmatrix}, \qquad \nabla \times \boldsymbol{E} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \dfrac{\partial}{\partial x} & \dfrac{\partial}{\partial y} & \dfrac{\partial}{\partial z} \\ E_x & E_y & E_z \end{vmatrix}. \]
이것들은 연속적인 시간 및 공간에서 성립하는 방정식들입니다. 이 페이지의 나머지 부분에서는 이를 시간에 따라 계산을 진행할 수 있는 명시적인 유한차분 형식으로 변환하는 과정을 설명합니다.
3. 상태밀도 방정식의 성분 형태
암페어 법칙을 직교좌표계 성분 형태로 전개하면 다음과 같다
\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = \frac{\partial H_z}{\partial y} - \frac{\partial H_y}{\partial z}, \]
\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = \frac{\partial H_x}{\partial z} - \frac{\partial H_z}{\partial x}, \]
\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y}. \]
패러데이 법칙을 전개하면 다음과 같다
\[ -\sigma_m H_x - \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_z}{\partial y} - \frac{\partial E_y}{\partial z}, \]
\[ -\sigma_m H_y - \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_x}{\partial z} - \frac{\partial E_z}{\partial x}, \]
\[ -\sigma_m H_z - \mu \frac{\partial H_z}{\partial t} = \frac{\partial E_y}{\partial x} - \frac{\partial E_x}{\partial y}. \]
현시점에서는 차원성에 관해 아무런 가정도 이루어지지 않았습니다. 위 식들은 전체 3차원 커얼 방정식의 성분 형태에 불과합니다.
4. 2차원 축소
유도 과정을 용이하게 하기 위해, 이제 \(y\) 방향으로는 어떠한 변화도 없다고 가정하여 문제를 2차원으로 축소시키면 \(\partial / \partial y = 0\)가 된다. 이 단계는 단지 대수적 계산을 간소화하고 업데이트 방정식의 구조를 보다 쉽게 이해할 수 있도록 하기 위함이다. 동일한 절차를 3차원으로도 직접 확장할 수 있으며 실제로 OghmaNano는 맥스웰 방정식의 완전한 3차원 형태를 해석한다.
이러한 가정을 적용하면 전기장 방정식은 다음과 같이 단순화된다.
\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = - \frac{\partial H_y}{\partial z}, \]
\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = - \frac{\partial H_z}{\partial x} + \frac{\partial H_x}{\partial z}, \]
\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x}. \]
자기 전도도를 0으로 설정하면 \(\sigma_m = 0\), 자기장 방정식은 다음과 같이 된다.
\[ \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_y}{\partial z}, \]
\[ \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} - \frac{\partial E_x}{\partial z}, \]
\[ \mu \frac{\partial H_z}{\partial t} = - \frac{\partial E_y}{\partial x}. \]
이들은 이제 이산화될 연속형 방정식들입니다. 이 문제는 본래 3차원 문제이지만 2차원으로 축소된 형태임을 명확히 밝혀두는 것이 중요합니다. 이는 모든 FDTD 시뮬레이션이 본질적으로 2차원적이라는 주장을 하기 위함이 아니라, 한 방향의 변화가 무시되는 일반적인 경우에 적용 가능한 명확한 업데이트 방정식을 도출하기 위함입니다.
5. Yee 격자, 인덱싱 및 표기법
위 방정식들을 수치적으로 해결하기 위해 공간과 시간을 모두 이산화한다. \(x\) 방향과 \(z\) 방향에서의 계산 격자 간격을 각각 \(\Delta x\)와 \(\Delta z\)라 하고, 시간 간격은 \(\Delta t\)라 하자. \(x\) 방향과 \(z\) 방향에 대해서는 정수 인덱스 \(i\)와 \(k\)를 사용한다.
\[ x_i = i \Delta x, \qquad z_k = k \Delta z. \]
시간은 정수 \(n\)로 인덱싱되며, 이에 따라
\[ t^n = n \Delta t. \]
Yee 격자 구조에서는 전기장과 자기장이 공간적으로나 시간적으로 동일한 위치에 저장되지 않고 서로 어긋나게 배치됩니다. 전기장 성분들은 정수 단위의 시간 레벨에서 계산되는 반면 \(n\Delta t\), 자기장 성분들은 반정수 단위의 시간 레벨에서 계산됩니다 \(\left(n+\tfrac{1}{2}\right)\Delta t\). 이러한 특성이 바로 시간영역유한차분법(FDTD)의 점프프로그 구조가 형성되는 원인입니다.
공간적으로도 각 전기장 및 자기장 성분들은 서로 오프셋되어 있어서, 모든 중심 유한차분 기반의 컬 연산이 인접한 격자점들 사이에서 자연스럽게 정의될 수 있다. 실제 코드 구현 시에는 각 전기장 및 자기장 성분마다 고유한 공간 인덱스 쌍이 할당되는데, 예컨대 \(E_x^n(i,k)\), \(H_y^{n+\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2})\) 등과 같다. 수식의 가독성을 높이기 위해 아래 단계들에서는 반복적으로 등장하는 ‘관찰용’ 인덱스들을 항상 명시하지는 않았으나, 그 의미는 언제나 동일하다.
- \(i\)은 \(x\) 방향을 나타냅니다.
- \(k\)은 \(z\) 방향을 인덱싱합니다.
- \(n\)은 정수 단위의 전기장 시간 레벨을 나타내며
- \(n+\tfrac{1}{2}\)은 자기장 값의 반단계 시간 레벨을 지정합니다.
예를 들어, \(\partial H_y / \partial z\)에 대한 중심 유한차분 근사식은 다음과 같습니다.
\[ \frac{\partial H_y}{\partial z} \approx \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}, \]
마찬가지로 \(\partial H_z / \partial x\)에 대한 중심차분 근사식은 다음과 같다
\[ \frac{\partial H_z}{\partial x} \approx \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]
이와 같이 표기법을 명확히 정의한 뒤에는 아래의 유도 과정을 Yee 격자에 대한 정확한 업데이트 방식으로 무오류적으로 이해할 수 있다.
6. 전기장 업데이트 방정식의 유도
전기장은 정수 시간 단계에서 계산되며 \(n\Delta t\), 자기장은 반정수 시간 단계에서 계산됩니다 \(\left(n+\tfrac{1}{2}\right)\Delta t\). 업데이트 방정식을 유도하기 위해 시간 미분항은 중심 유한차분법으로 근사하고, 전도도 항은 이전 및 현재 전기장 값의 평균을 이용해 계산합니다. 이를 통해 전도도에 대한 수치적으로 안정적인 사다리꼴 근사가 가능해집니다.
이러한 평균화 과정은 매우 중요합니다. 만약 전도도 항을 부적절하게 처리할 경우 업데이트 계수들이 기존의 감쇠 형태를 유지하지 못하게 되며, 일반적인 유한차분 시간영역법(FDTD) 방정식들과의 일치성도 상실하게 됩니다.
6.1 \(E_x\)에 대한 업데이트 방정식
시작점으로 삼기
\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = - \frac{\partial H_y}{\partial z}. \]
\(E_x^{n+1}\)과 \(E_x^n\)의 평균값을 이용하여 전도도 항을 계산하고, 중심 유한차분법을 활용하여 시간 미분항을 산출한다.
\[ \sigma \frac{E_x^{n+1}(i,k) + E_x^n(i,k)}{2} + \epsilon \frac{E_x^{n+1}(i,k) - E_x^n(i,k)}{\Delta t} = - \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]
\(E_x^{n+1}\) 항들을 좌변으로 모아보자:
\[ \left( \frac{\sigma}{2} + \frac{\epsilon}{\Delta t} \right) E_x^{n+1}(i,k) = - \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} + \left( \frac{\epsilon}{\Delta t} - \frac{\sigma}{2} \right) E_x^n(i,k). \]
따라서
\[ E_x^{n+1}(i,k) = \frac{ \left( \dfrac{\epsilon}{\Delta t} - \dfrac{\sigma}{2} \right) E_x^n(i,k) - \dfrac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} }{ \dfrac{\epsilon}{\Delta t} + \dfrac{\sigma}{2} }. \]
분자와 분모에 \(2\Delta t\)을 곱하면 익숙한 계수 형태가 도출된다:
\[ E_x^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_x^n(i,k) - \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]
6.2 \(E_y\)에 대한 업데이트 방정식
시작점으로 삼기
\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = - \frac{\partial H_z}{\partial x} + \frac{\partial H_x}{\partial z}. \]
Yee 격자 상에서 두 미분항을 이산화한다:
\[ \sigma \frac{E_y^{n+1}(i,k) + E_y^n(i,k)}{2} + \epsilon \frac{E_y^{n+1}(i,k) - E_y^n(i,k)}{\Delta t} = - \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x} + \frac{ H_x^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_x^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]
식을 재배열하면 다음과 같다
\[ E_y^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_y^n(i,k) + \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \left[ - \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x} + \frac{ H_x^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_x^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} \right]. \]
2차원 모델에서는 전기장 성분 중 오직 \(E_z\)만이 두 개의 서로 다른 자기장 미분값으로부터 영향을 받으므로 코딩 시 가장 실수하기 쉬운 부분입니다. 특히 \(x\) 및 \(z\) 항에 포함된 부호는 커울 방정식의 전개 결과에 따라 결정되므로 절대로 변경해서는 안 됩니다.
6.3 \(E_z\)에 대한 업데이트 방정식
시작점으로 삼는다
\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x}. \]
이산화 과정을 통해 다음과 같은 결과가 도출된다
\[ \sigma \frac{E_z^{n+1}(i,k) + E_z^n(i,k)}{2} + \epsilon \frac{E_z^{n+1}(i,k) - E_z^n(i,k)}{\Delta t} = \frac{ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_y^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]
따라서
\[ E_z^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_z^n(i,k) + \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \frac{ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_y^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]
이제 전기장에 대한 이산화 업데이트 방정식들이 표준 감쇠 형태로 도출되었습니다. 다음 단계로는 반정수 시간 단계에서 업데이트되는 자기장 방정식들을 유도해야 합니다.
7. 자기장 업데이트 방정식의 유도
자기장은 반정수 시간 단계에서 업데이트된다. \(\sigma_m = 0\)에 따르면 자기장 방정식에는 감쇠 항이 존재하지 않으므로 업데이트 공식은 보다 간단한 전진 형태를 띤다. 이 부분의 유도 과정은 사다리꼴 형태의 전도도 항이 포함되지 않기 때문에 전기장의 경우에 비해 대수적으로 더 간결하다.
7.1 \(H_x\)에 대한 업데이트 방정식
시작점으로 삼기
\[ \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_y}{\partial z}. \]
\(n-\tfrac{1}{2}\)과 \(n+\tfrac{1}{2}\) 사이의 시간 간격 및 \(z\) 방향에 대해 중심차분법을 적용하여 이산화하면 다음과 같은 결과가 도출된다.
\[ \mu \frac{ H_x^{n+\frac{1}{2}}(i,k+\frac{1}{2}) - H_x^{n-\frac{1}{2}}(i,k+\frac{1}{2}) }{\Delta t} = \frac{ E_y^n(i,k+1) - E_y^n(i,k) }{\Delta z}. \]
유도 과정에서 자주 사용되는 간략한 반지수 표기법을 사용하면 동일한 결과를 다음과 같이 표현할 수 있다.
\[ \mu \frac{ H_x^{n+\frac{1}{2}} - H_x^{n-\frac{1}{2}} }{\Delta t} = \frac{ E_y^n\!\left(k+\frac{1}{2}\right) - E_y^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]
따라서
\[ H_x^{n+\frac{1}{2}} = H_x^{n-\frac{1}{2}} + \frac{\Delta t}{\mu} \frac{ E_y^n\!\left(k+\frac{1}{2}\right) - E_y^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]
7.2 \(H_y\)에 대한 업데이트 방정식
시작점으로 삼기
\[ \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} - \frac{\partial E_x}{\partial z}. \]
이산화 과정을 통해 다음과 같은 결과가 도출된다
\[ \mu \frac{ H_y^{n+\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2}) - H_y^{n-\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2}) }{\Delta t} = \frac{ E_z^n(i+1,k+\frac{1}{2}) - E_z^n(i,k+\frac{1}{2}) }{\Delta x} - \frac{ E_x^n(i+\frac{1}{2},k+1) - E_x^n(i+\frac{1}{2},k) }{\Delta z}. \]
다시 말해, 축약된 반지수 표기법에서
\[ \mu \frac{ H_y^{n+\frac{1}{2}} - H_y^{n-\frac{1}{2}} }{\Delta t} = \frac{ E_z^n\!\left(i+\frac{1}{2}\right) - E_z^n\!\left(i-\frac{1}{2}\right) }{\Delta x} - \frac{ E_x^n\!\left(k+\frac{1}{2}\right) - E_x^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]
따라서
\[ H_y^{n+\frac{1}{2}} = H_y^{n-\frac{1}{2}} + \frac{\Delta t}{\mu} \left[ \frac{ E_z^n\!\left(i+\frac{1}{2}\right) - E_z^n\!\left(i-\frac{1}{2}\right) }{\Delta x} - \frac{ E_x^n\!\left(k+\frac{1}{2}\right) - E_x^n\!\left(k-\frac{1}{2}\right) }{\Delta z} \right]. \]
7.3 \(H_z\)에 대한 업데이트 방정식
시작점으로 삼는다
\[ \mu \frac{\partial H_z}{\partial t} = - \frac{\partial E_y}{\partial x}. \]
이산화 과정을 거치면 다음과 같은 결과가 도출된다
\[ \mu \frac{ H_z^{n+\frac{1}{2}}(i+\frac{1}{2},k) - H_z^{n-\frac{1}{2}}(i+\frac{1}{2},k) }{\Delta t} = - \frac{ E_y^n(i+1,k) - E_y^n(i,k) }{\Delta x}. \]
동일하게, 간결한 표기법으로 표현하면
\[ \mu \frac{ H_z^{n+\frac{1}{2}} - H_z^{n-\frac{1}{2}} }{\Delta t} = - \frac{ E_y^n\!\left(i+\frac{1}{2}\right) - E_y^n\!\left(i-\frac{1}{2}\right) }{\Delta x}. \]
따라서
\[ H_z^{n+\frac{1}{2}} = H_z^{n-\frac{1}{2}} - \frac{\Delta t}{\mu} \frac{ E_y^n\!\left(i+\frac{1}{2}\right) - E_y^n\!\left(i-\frac{1}{2}\right) }{\Delta x}. \]
이러한 치환 과정을 마치면 표준적인 리프플로그 방식의 유한차분 시간영역법(FDTD) 업데이트 방정식들이 바로 도출된다.
8. 최종 업데이트 방정식
계산 결과를 종합하면 2차원 유한차분 시간영역법(FDTD) 업데이트 방정식들이 도출된다.
\[ E_x^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_x^n(i,k) - \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}, \]
\[ E_y^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_y^n(i,k) + \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \left[ - \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x} + \frac{ H_x^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_x^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} \right], \]
\[ E_z^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_z^n(i,k) + \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \frac{ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_y^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}, \]
\[ H_x^{n+\frac{1}{2}}(i,k+\tfrac{1}{2}) = H_x^{n-\frac{1}{2}}(i,k+\tfrac{1}{2}) + \frac{\Delta t}{\mu} \frac{ E_y^n(i,k+1) - E_y^n(i,k) }{\Delta z}, \]
\[ H_y^{n+\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2}) = H_y^{n-\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2}) + \frac{\Delta t}{\mu} \left[ \frac{ E_z^n(i+1,k+\tfrac{1}{2}) - E_z^n(i,k+\tfrac{1}{2}) }{\Delta x} - \frac{ E_x^n(i+\tfrac{1}{2},k+1) - E_x^n(i+\tfrac{1}{2},k) }{\Delta z} \right], \]
\[ H_z^{n+\frac{1}{2}}(i+\tfrac{1}{2},k) = H_z^{n-\frac{1}{2}}(i+\tfrac{1}{2},k) - \frac{\Delta t}{\mu} \frac{ E_y^n(i+1,k) - E_y^n(i,k) }{\Delta x}. \]
이와 같은 형태로 표현하면 리프플로그 방식의 작동 원리가 명확해진다. 시간 단계 \(n+1\)에서의 전기장은 시간 단계 \(n+\tfrac{1}{2}\)에서의 자기장으로부터 계산되며, 시간 단계 \(n+\tfrac{1}{2}\)에서의 자기장은 시간 단계 \(n\)에서의 전기장으로부터 계산된다. 이에 따라 전기장과 자기장은 시간적으로 번갈아가며 진화하게 된다.
9. 구현에 관한 참고사항
실제 코드 구현 시에는 이전 시간 단계의 전기장 값 및 컬 연산 결과에 곱해지는 계수들이 사전에 계산되어 저장되는 경우가 많으며, 특히 \(\epsilon\), \(\mu\), \(\sigma\)가 격자 셀마다 달라질 때 그러합니다. 예를 들어 전기장 업데이트 과정은 일반적으로 다음과 같은 형태를 띱니다.
\[ E^{n+1} = C_a E^n + C_b (\nabla \times H), \]
여기서 지역 계수 \(C_a\)과 \(C_b\)는 재료 파라미터와 시간 간격을 나타냅니다. 이를 통해 구현 과정이 간결해지며 매 시간 간격마다 동일한 계수들을 재계산할 필요가 없게 됩니다.
앞서 설명한 유도 과정에서는 시간영역유한차분법(FDTD)의 핵심 업데이트 방식만을 도출하였을 뿐입니다. 실제 시뮬레이션을 수행하려면 광원, 경계 조건, 완벽 정합층(PML)과 같은 흡수층, 그리고 쿠랑 안정성 조건을 만족하는 시간 간격도 반드시 필요합니다. 이러한 추가 요소들은 앞서 도출된 리프플로그 구조 자체를 변경하지는 않지만, 각 격자 지점에서 업데이트 방정식이 적용되는 방식에는 영향을 미칩니다.
마지막으로, 위에서 설명한 표기법은 구현에 활용할 수 있을 정도로 충분히 명확하게 기술되었으나 각 구성요소의 정확한 공간 배치는 코드 내에서 사용되는 격자 규칙에 따라 달라집니다. 가장 중요한 요건은 표기 형태 자체가 아니라 일관성입니다. 즉, 컬 연산자와 공간상의 배치 방식, 그리고 부호들이 전체 시뮬레이션 과정에서 동일한 Yee 격자 구조에 부합해야 합니다.