수치 전자기학

**수치 전자기학**은 전자기학의 핵심 방정식인 맥스웰 방정식을 수치적으로 푸는 것을 의미한다. 문제가 복잡하여 접근하기 어려운 경우에 이 방법을 사용하여 근사적인 해를 구할 수 있다. 수치 전자기학에서 사용하는 대부분의 방법들은 문제를 행렬 방정식으로 나타내어 수치적으로 풀어낸다.
# 종류
맥스웰 방정식을 수치적으로 푸는 방법은 방정식을 표현하는 방법에 따라 크게 세 가지 종류로 나뉜다. 종류에는 각각 방정식을 미분형으로 나타내는 방법과 적분형으로 나타내는 방법, 그리고 변분형으로 나타내는 방법이 있다. 미분형 방정식을 사용하는 방법에는 유한 차분법이 있고, 적분형 방정식을 사용하는 방법에는 모멘트법이 있으며, 변분형 방정식을 사용하는 방법에는 유한 요소법이 있다. 유한 차분법과 모멘트법, 그리고 유한 요소법 모두 문제를 행렬 방정식으로 표현한 뒤 방정식을 풀어서 해를 얻기 때문에 행렬 방정식을 효율적으로 푸는 알고리즘 또한 중요하다. 행렬 방정식을 푸는 방법에는 가우스 소거법, 가우스-자이델 방법(Gauss-Seidel), 축차가속완화법(SOR, Successive Over-Relaxation), 경사 하강법(Gradient Descent), 켤레기울기법(Conjugate Gradient) 등이 있다.
# 유한 차분법 (Finite Difference Method)
유한 차분법이란, 미분 방정식을 유한한 차분의 선형 방정식으로 근사하여 해를 구하는 방법을 의미한다. 유한 차분법을 이용하면 미분형으로 나타낸 푸아송 방정식을 수치적으로 풀 수 있다. 푸아송 방정식은 맥스웰 방정식을 다르게 표현한 미분 방정식으로, 전기장 대신 전기 퍼텐셜과 관련된 미분 방정식이다.
## 일계 미분 3점 공식
f(x)를 x_0주위에서 테일러 전개하면 다음과 같이 된다. f(x_0+h)=f(x_0)+hf'(x_0)+\frac{h^2}{2!}f(x_0)+\frac{h^3}{3!}f^{(3)}(x_0)+\cdots f(x_0-h)=f(x_0)-hf'(x_0)+\frac{h^2}{2!}f(x_0)-\frac{h^3}{3!}f^{(3)}(x_0)+\cdots 이 두 식을 서로 빼 주고 정리하면, 일계 미분의 3점 공식을 얻을 수 있다. f'(x_0)=\frac{1}{2h}[f(x_0+h)-f(x_0-h)]+O(h^2)
## 이계 미분 3점 공식
일계 미분 3점 공식을 유도할 때 사용하였던 f(x)의 테일러 전개를 다시 사용하여 이계 미분의 3점 공식을 유도할 수 있다. 두 식을 서로 더해 주고 정리하면, 이계 미분의 3점 공식을 얻을 수 있다. f''(x_0)=\frac{1}{h^2}[f(x_0+h)-2f(x_0)+f(x_0-h)]+O(h^2) 이 식을 통해서 함수의 이계 미분값을 근사할 수 있고, 따라서 2차원 푸아송 방정식을 3점 공식을 이용하여 근사하면 다음과 같아진다. \nabla^2\phi(x,y)|_{(x_0,y_0)}=\phi_{xx}(x_0,y_0)+\phi_{yy}(x_0,y_0)=\frac{1}{h^2}[\phi(x_0+h,y_0)+\phi(x_0,y_0+h)-4\phi(x_0,y_0)+\phi(x_0-h,y_0)+\phi(x_0,y_0-h)]+O(h^2)
## 격자 생성
격자점과 각 점의 좌표를 보여주고 있다. 문제에서 주어진 영역을 우측의 그림과 같이 정사각형(혹은 직사각형)의 격자판으로 나누어준다. 그러면 각 점에서 다음을 만족하게 된다(유전체가 없는 공간을 가정). 여기서 h는 격자점 사이의 거리를 의미한다. \frac{1}{h^2}[\phi(x_{i+1},y_j)+\phi(x_{i-1},y_j)-4\phi(x_i,y_j)+\phi(x_i,y_{j+1})+\phi(x_i,y_{j-1})]\cong\nabla^2\phi(x,y)|_{(x_i,y_j)}=-\frac{\rho(x_i,y_j)}{\epsilon_0} 위 방정식에서 좌항의 \phi(x_i,y_j)는 미지수이고, 우항의 \rho(x_i,y_j)는 이미 알고 있는 값이다. 따라서 주어진 문제는 격자점 개수만큼의 방정식과 미지수 \phi(x_i,y_j)를 가진 연립일차방정식과 같아진다. 연립일차방정식은 행렬 방정식으로 나타낼 수 있기 때문에, 결과적으로 주어진 문제를 푸는 것은 근사적으로 Ax=b 꼴의 행렬 방정식을 푸는 것과 같다. 이때 A는 정사각행렬로, \phi(x_i,y_j)의 계수가 들어가고, x에는 \phi(x_i,y_j)가 들어가며, b에는 \rho(x_i,y_j)가 들어가게 된다. 위의 연립방정식으로부터 행렬을 생성하기 위해서는 각각의 격자점 (x_i,y_j)에 번호를 매겨야 한다. 예를 들어, (x_0,y_0)은 0, (x_0,y_1)은 1, (x_0,y_2)는 2, ... 과 같은 방식으로 번호를 매길 수 있다. 번호를 매기고 나면, 각각의 격자점이 (x_i,y_j)=P_n과 같이 하나의 첨자로 나타내어지게 된다. 따라서 위의 방정식을 다음과 같이 간략하게 나타낼 수 있게 된다. 이때 a, b, c, d는 앞서 매긴 번호와 격자점 사이의 관계에 따라 i가 정해지면 자동으로 결정된다. \phi(P_a)+\phi(P_b)-4\phi(P_i)+\phi(P_c)+\phi(P_d)=-\frac{h^2}{\epsilon_0}\rho(P_i) 이렇게 격자점들을 정렬시키고 나면, x와 b 벡터, 그리고 행렬 A를 다음과 같이 생성할 수 있다. x_i=\phi(P_i), b_i=-\frac{h^2}{\epsilon_0}\rho(P_i), A_{i,j} = \begin{cases} -4, & \text{if }j=i \\\ 1, & \text{if }j=a,b,c,d \\\ 0, & \text{otherwise} \end{cases} 단, 격자점이 주어진 영역의 경계에 있는 경우 방정식이 위에 기술한 것과 조금 달라져 행렬을 수정해야 한다. 이 주제는 경계 조건 항목에서 좀 더 자세하게 다룬다.
## 경계 조건

### 디리클레 경계
디리클레 경계조건은 경계에서 퍼텐셜 값이 상수로 고정되는 조건이다. 따라서 연립방정식으로부터 퍼텐셜을 계산할 필요가 없는 위치이므로, 행렬에서 제외된다. 디리클레 경계와 인접한 격자점의 경우, 디리클레 경계의 퍼텐셜이 미지수가 아니므로, 방정식이 조금 수정되어야 한다. P_n과 인접하여 디리클레 경계 위의 점 P_a가 있을 때 방정식은 다음과 같이 수정된다. 이때 \phi_b(P_a)는 경계에서의 퍼텐셜로, 상수이다. \phi(P_b)-4\phi(P_n)+\phi(P_c)+\phi(P_d)=-\frac{h^2}{\epsilon_0}\rho(P_n) - \phi_b(P_a) P_a가 아닌 다른 점들이 P_n과 인접하여도 같은 논리로 방정식을 수정해주면 된다. 벡터 b와 행렬 A는 디리클레 경계조건이 있을 때 다음과 같이 수정된다. b_i=\begin{cases} -\frac{h^2}{\epsilon_0}\rho(P_i), & \text{if }P_i \text{ is not adjacent to the dirichlet boundaries} \\\ -\frac{h^2}{\epsilon_0}\rho(P_i) - \sum\phi_b(P_j), & \text{otherwise} \end{cases}, A_{i,j} = \begin{cases} -4, & \text{if }j=i \\\ 1, & \text{if }j=k \text{ and } P_k \text{ is not adjacent to the dirichlet boundaries, }k=a,b,c,d \\\ 0, & \text{otherwise} \end{cases} 노이만 경계조건에 따른 추가 격자점과 각 점의 좌표를 보여주고 있다.
### 노이만 경계
노이만 경계조건은 경계에서 퍼텐셜의 일계 미분 값이 상수로 고정되는 조건이다. 이 조건을 적용시키기 위해서는 경계 외부에 인접한 격자점을 추가하여, 경계 격자점에서의 일계 미분값이 일정하도록 만들어야 한다. 우측 그림과 같은 상황에서, 경계에서의 미분값은 일계 미분의 3점 공식에 의해 다음과 같이 주어진다. \phi_y(x_i,y_j) = \frac{1}{2h}[\phi(x_i,y_{j+1})-\phi(x_i,y_{j-1})] 노이만 경계조건에서는 일계 미분 값이 상수로 정해져있으므로, 위의 관계로부터 외부 추가 격자점의 퍼텐셜인 \phi(x_i,y_{j-1})를 다음과 같이 계산할 수 있다. \phi(x_i,y_{j-1})=\phi(x_i,y_{j+1})-2h\phi_y(x_i,y_i) 이렇게 구한 외부 추가 격자점의 퍼텐셜로 부터, 노이만 경계에서의 퍼텐셜을 일반 격자점과 같은 방법으로 계산 가능하다. 노이만 경계에서의 퍼텐셜 방정식은 다음과 같다. \phi(x_{i+1},y_j)+\phi(x_{i-1},y_j)-4\phi(x_i,y_j)+\phi(x_i,y_{j+1})+\phi(x_i,y_{j-1})=-\frac{h^2}{\epsilon_0}\rho(x_i,y_j) 이 방정식을 행렬에 포함시켜 주면 해가 유일하게 존재하는 행렬 방정식을 얻을 수 있다.
## 선형 유전체 내부에서의 퍼텐셜 해
유전체 경계면에서의 격자점과 각 점의 좌표를 보여주고 있다. 선형 유전체란 분극이 전기장에 비례하는 유전체를 의미한다. 이러한 유전체의 경우 유전체가 없는 진공인 상황에서 사용했던 방법을 그대로 사용하여 푸아송 방정식을 풀 수 있다. 앞서 기제된 방정식들의 유전율 항에 진공의 유전율 \epsilon_0대신 선형 유전체 내부에서의 유전율 \epsilon을 대입해주면 유전체 내부에서의 퍼텐셜 해를 구할 수 있다. 선형 유전체의 또 다른 특징은, 유전체 경계면에서 다음과 같은 경계조건을 만족한다는 것이다. \oint_C \epsilon \frac{\partial \phi}{\partial n}dl=-\int_S\rho_fda 여기서 C는 (x_i,y_j) 주위를 둘러싸는 닫힌 곡선을 의미하고, S는 C를 경계로 가지는 곡면을 의미한다. 또한 \frac{\partial \phi}{\partial n}은 곡선에 수직한 방향으로의 퍼텐셜 미분을 의미한다. C를 우측 그림과 같이 각 격자점의 중심을 지나는 정사각형으로 잡고, 격자 사이의 간격이 충분히 작다고 가정하면, 일계 미분의 3점 공식에 의해 위의 식이 다음과 같이 근사된다. 이때 \epsilon_{avg}는 \frac{\epsilon_1+\epsilon_2}{2}를 의미한다. \epsilon_{avg}[\phi(x_{i+1},y_j)+\phi(x_{i-1},y_j)-4\phi(x_i,y_j)]+\epsilon_1\phi(x_i,y_{j+1})+\epsilon_2\phi(x_i,y_{j-1})\cong-\rho_f(x_i,y_j)h^2 이 방정식을 이용하여 선형 유전체 경계면 위에 있는 격자점을 행렬에 포함시켜주면 퍼텐셜을 계산할 수 있다.
# 모멘트법 (Method of Moments)
모멘트법은 푸아송 방정식의 적분형 해를 근사적으로 계산하는 방법이다.
# 유한 요소법 (Finite Elements Method)
유한 요소법은 푸아송 방정식의 변분형 해를 근사적으로 계산하는 방법이다. 이 방법은 공간을 여러 조각으로 잘게 나누어 각 조각에서의 퍼텐셜을 선형 근사하는 방식을 사용한다. 유한 차분법과는 달리 격자점을 자유롭게 선택할 수 있기 때문에 경계면이 곡면일 경우 유한 요소법을 사용하는 것이 유한 차분법을 사용하는 것 보다 효율적이다.
## 변분법과 푸아송 방정식
범함수 I가 다음과 같이 정의 될 때, I를 최소화하는 함수 V(x,y)는 푸아송 방정식의 해가 된다. I(\phi)=\int\bigg[\frac{1}{2}\epsilon(\nabla V)^2-\rho V \bigg]da 이때 피적분 함수의 첫 항은 전기장의 에너지를 의미한다. 따라서, 전하 분포가 없을 때의 푸아송 방정식(라플라스 방정식)의 해는 계의 에너지를 최소화하는 퍼텐셜임을 알 수 있다.
## 퍼텐셜의 선형 근사
유한 요소법의 격자와 각 격자점에서의 퍼텐셜을 나타내고 있다. 유한 요소법에서는 공간을 작은 요소들로 쪼갠 후 각 요소에서의 퍼텐셜을 선형 근사한다. 우측 그림과 같이 삼각형 요소로 공간을 분할하는 경우 (e)번째 요소에서의 선형 퍼텐셜은 다음과 같이 주어진다. V^{(e)}(x,y)=\begin{bmatrix} 1 & x & y \end{bmatrix}\begin{bmatrix} a^{(e)} \\\ b^{(e)} \\\ c^{(e)} \end{bmatrix} 선형 근사된 퍼텐셜이 세 꼭짓점에서 원래의 퍼텐셜과 일치하기 위해서는 a^{(e)}, b^{(e)}, c^{(e)}가 다음 조건을 만족해야 한다. R^{(e)}는 행렬을 축약하기 위해 사용하는 기호이다. \begin{bmatrix} V^{(e)}_1 \\\ V^{(e)}_2 \\\ V^{(e)}_3 \end{bmatrix} = \begin{bmatrix} 1 & x_1^{(e)} & y_1^{(e)} \\\ 1 & x_2^{(e)} & y_2^{(e)} \\\ 1 & x_3^{(e)} & y_3^{(e)} \end{bmatrix} \begin{bmatrix} a^{(e)} \\\ b^{(e)} \\\ c^{(e)} \end{bmatrix} =R^{(e)}\begin{bmatrix} a^{(e)} \\\ b^{(e)} \\\ c^{(e)} \end{bmatrix} R^{(e)} 행렬을 이용하면 요소에서의 퍼텐셜을 다음과 같이 조금 다르게 표현할 수 있다. V^{(e)}(x,y)=\begin{bmatrix} 1 & x & y \end{bmatrix} R^{(e)-1}R^{(e)}\begin{bmatrix} a^{(e)} \\\ b^{(e)} \\\ c^{(e)} \end{bmatrix}=\begin{bmatrix} \alpha_1^{(e)} & \alpha_2^{(e)} & \alpha_3^{(e)} \end{bmatrix} \begin{bmatrix} V^{(e)}_1 \\\ V^{(e)}_2 \\\ V^{(e)}_3 \end{bmatrix} 여기서 \alpha_i^{(e)}(x,y)는 형상함수(shape function)로, 다음 관계를 만족한다. \alpha_i^{(e)}(x_j^{(e)},y_j^{(e)}) = \delta_{ij}=\begin{cases} 1, & \text{if } i = j \\\ 0, & \text{if }i \neq j \end{cases}
## 범함수의 계산
각 요소의 범함수 I^{(e)} 를 형상함수로 나타내면 다음과 같다. I^{(e)}=\sum_i\sum_j\bigg[\frac{1}{2}\epsilon V^{(e)}_iV^{(e)}_j\int\nabla\alpha_i^{(e)}\cdot\nabla\alpha_j^{(e)}da-V^{(e)}_i\rho^{(e)}\int\alpha_i^{(e)}\alpha_j^{(e)}da\bigg] 이때 C^{(e)}_{ij}=\int\nabla\alpha_i^{(e)}\cdot\nabla\alpha_j^{(e)}da 로, T^{(e)}_{ij}=\int\alpha_i^{(e)}\alpha_j^{(e)}da 로 정의하고, \Rho^{(e)} = \rho^{(e)}\begin{bmatrix} 1 \\\ 1 \\\ 1 \end{bmatrix} 라고 정의하면, 요소의 범함수 I^{(e)} 는 다음과 같이 간단해진다. I^{(e)}=\frac{1}{2}\epsilon V^{(e)T}C^{(e)}V^{(e)}-V^{(e)T}T^{(e)}\Rho^{(e)} 여기서 윗 첨자 T 는 전치행렬을 의미한다. 전체 범함수는 각 요소의 피적분 함수의 합이다. I=\sum_e I^{(e)} 전체 범함수를 행렬 형태로 나타내기 위해서는 먼저 각 요소의 관점에서 본 꼭짓점 번호와 전체의 관점에서 본 격자점 번호를 대응시켜야 한다. 이 과정은 유한 차분법에서의 격자 좌표에 번호를 매기는 과정과 유사하다. 꼭짓점 번호를 격자점과 모두 대응시키고 나면, 각 요소 관점에서 만들어졌던 행렬(예를 들어 C^{(e)} )을 각 성분 번호에 대응되는 행렬(C^{(e)}_G )로 바꾸어주어야 한다. 전체 행렬은 C=\sum_eC^{(e)}_G 와 같이 변환된 모든 행렬을 더해주어 만들 수 있다. 이와 같이 모든 행렬을 전체 행렬로 변환시켜 주면, 전체 범함수를 다음과 같이 표현할 수 있다. I=\frac{1}{2}\epsilon V^TCV - V^TT\Rho
## 범함수의 최소화
범함수를 최소화하기 위해서는 다음 조건을 만족해야 한다. \frac{\partial I}{\partial V_i}=0 , V_i 는 i 번째 격자점의 퍼텐셜 이를 모든 i 에 대하여 묶어서 서술하면 아래와 같아진다. \begin{bmatrix} \frac{\partial }{\partial V_1} \\\ \frac{\partial }{\partial V_2} \\\ \vdots \end{bmatrix}I=\epsilon CV - T\Rho=0 이 방정식을 통해서 V 벡터를 구할 수 있다.
## 경계 조건
디리클레 경계 조건을 사용하는 경우, 경계 위에서의 격자점에서는 이미 퍼텐셜이 상수로 고정되어있다. 따라서 경계에서의 격자점들은 행렬 방정식에서 따로 분리해주어야 한다. C , V , T , \Rho 행렬의 원소를 재배열하여 경계에서의 격자점들을 따로 분리해줄 수 있다. 예를 들어서 V 행렬의 경우 경계 격자점들의 퍼텐셜을 나타내는 V_p 와 아직 알지 못하여 구해야 하는 퍼텐셜인 V_f 로 나누어서 V=\begin{bmatrix} V_f \\\ V_p \end{bmatrix} 와 같이 나타낼 수 있다. 마찬가지로, C 행렬도 C=\begin{bmatrix} C_{ff} & C_{fp} \\\ C_{pf} & C_{pp} \end{bmatrix} 와 같이 블록 행렬로 나타낼 수 있다. 이렇게 블록 행렬로 나타내고 나면, 구해야 하는 퍼텐셜 V_f 가 다음과 같이 주어지게 된다. C_{ff}V_f=-C_{fp}V_p+\frac{1}{\epsilon}[T_{ff}\Rho_f+T_{fp}\Rho_p] 이 행렬 방정식은 Ax=b 꼴의 행렬 방정식으로, C_{ff} 가 정사각행렬이므로, 특수한 경우를 제외하고는 풀 수 있다.
## 푸아송 방정식의 해
푸아송 방정식의 수치적 해는 행렬 방정식을 통해 구한 전체 V 벡터를 다시 각 요소에서의 V^{(e)} 로 변환해주어 구할 수 있다. 해는 다음과 같다. V^{(e)}(x,y)=\begin{bmatrix} \alpha_1^{(e)} & \alpha_2^{(e)} & \alpha_3^{(e)} \end{bmatrix} \begin{bmatrix} V^{(e)}_1 \\\ V^{(e)}_2 \\\ V^{(e)}_3 \end{bmatrix}