수치 전자기학 문서 원본 보기

{{위키데이터 속성 추적}}
'''수치 전자기학'''은 전자기학의 핵심 방정식인 [[맥스웰 방정식]]을 수치적으로 푸는 것을 의미한다. 문제가 복잡하여 접근하기 어려운 경우에 이 방법을 사용하여 근사적인 해를 구할 수 있다. 수치 전자기학에서 사용하는 대부분의 방법들은 문제를 행렬 방정식으로 나타내어 수치적으로 풀어낸다.

== 종류 ==
[[맥스웰 방정식]]을 수치적으로 푸는 방법은 방정식을 표현하는 방법에 따라 크게 세 가지 종류로 나뉜다. 종류에는 각각 방정식을 미분형으로 나타내는 방법과 적분형으로 나타내는 방법, 그리고 변분형으로 나타내는 방법이 있다. 미분형 방정식을 사용하는 방법에는 [[유한차분법|유한 차분법]]이 있고, 적분형 방정식을 사용하는 방법에는 모멘트법이 있으며, 변분형 방정식을 사용하는 방법에는 [[유한요소법|유한 요소법]]이 있다.

[[유한차분법|유한 차분법]]과 모멘트법, 그리고 [[유한요소법|유한 요소법]] 모두 문제를 행렬 방정식으로 표현한 뒤 방정식을 풀어서 해를 얻기 때문에 행렬 방정식을 효율적으로 푸는 알고리즘 또한 중요하다. 행렬 방정식을 푸는 방법에는 [[가우스 소거법]], [[가우스-자이델 방법]](Gauss-Seidel), [[축차가속완화법]](SOR, Successive Over-Relaxation), [[경사 하강법]](Gradient Descent), [[켤레기울기법]](Conjugate Gradient) 등이 있다.

== 유한 차분법 (Finite Difference Method) ==
[[유한차분법|유한 차분법]]이란, 미분 방정식을 유한한 차분의 선형 방정식으로 근사하여 해를 구하는 방법을 의미한다. [[유한차분법|유한 차분법]]을 이용하면 미분형으로 나타낸 [[푸아송 방정식]]을 수치적으로 풀 수 있다. [[푸아송 방정식]]은 [[맥스웰 방정식]]을 다르게 표현한 미분 방정식으로, [[전기장]] 대신 [[전기 퍼텐셜]]과 관련된 미분 방정식이다.

=== 일계 미분 3점 공식 ===
<math>f(x)</math>를 <math>x_0</math>주위에서 [[테일러 전개]]하면 다음과 같이 된다.

<math>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</math>

<math>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</math>

이 두 식을 서로 빼 주고 정리하면, 일계 미분의 3점 공식을 얻을 수 있다.

<math>f'(x_0)=\frac{1}{2h}[f(x_0+h)-f(x_0-h)]+O(h^2)</math>

=== 이계 미분 3점 공식 ===
일계 미분 3점 공식을 유도할 때 사용하였던 <math>f(x)</math>의 [[테일러 전개]]를 다시 사용하여 이계 미분의 3점 공식을 유도할 수 있다. 두 식을 서로 더해 주고 정리하면, 이계 미분의 3점 공식을 얻을 수 있다.

<math>f''(x_0)=\frac{1}{h^2}[f(x_0+h)-2f(x_0)+f(x_0-h)]+O(h^2)</math>

이 식을 통해서 함수의 이계 미분값을 근사할 수 있고, 따라서 2차원 [[푸아송 방정식]]을 3점 공식을 이용하여 근사하면 다음과 같아진다.

<math>\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)</math>

=== 격자 생성 ===
[[파일:Grid points.jpg|섬네일|격자점과 각 점의 좌표를 보여주고 있다.]]
문제에서 주어진 영역을 우측의 그림과 같이 정사각형(혹은 [[직사각형]])의 격자판으로 나누어준다. 그러면 각 점에서 다음을 만족하게 된다(유전체가 없는 공간을 가정). 여기서 <math>h</math>는 격자점 사이의 거리를 의미한다.

<math>\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}</math>

위 방정식에서 좌항의 <math>\phi(x_i,y_j)</math>는 미지수이고, 우항의 <math>\rho(x_i,y_j)</math>는 이미 알고 있는 값이다. 따라서 주어진 문제는 격자점 개수만큼의 방정식과 미지수 <math>\phi(x_i,y_j)</math>를 가진 [[연립 일차 방정식|연립일차방정식]]과 같아진다. 연립일차방정식은 행렬 방정식으로 나타낼 수 있기 때문에, 결과적으로 주어진 문제를 푸는 것은 근사적으로 <math>Ax=b</math> 꼴의 행렬 방정식을 푸는 것과 같다. 이때 <math>A</math>는 [[정사각행렬]]로, <math>\phi(x_i,y_j)</math>의 계수가 들어가고, <math>x</math>에는 <math>\phi(x_i,y_j)</math>가 들어가며, <math>b</math>에는 <math>\rho(x_i,y_j)</math>가 들어가게 된다.

위의 연립방정식으로부터 행렬을 생성하기 위해서는 각각의 격자점 <math>(x_i,y_j)</math>에 번호를 매겨야 한다. 예를 들어, <math>(x_0,y_0)</math>은 0, <math>(x_0,y_1)</math>은 1, <math>(x_0,y_2)</math>는 2, ... 과 같은 방식으로 번호를 매길 수 있다. 번호를 매기고 나면, 각각의 격자점이 <math>(x_i,y_j)=P_n</math>과 같이 하나의 첨자로 나타내어지게 된다. 따라서 위의 방정식을 다음과 같이 간략하게 나타낼 수 있게 된다. 이때 <math>a</math>, <math>b</math>, <math>c</math>, <math>d</math>는 앞서 매긴 번호와 격자점 사이의 관계에 따라 <math>i</math>가 정해지면 자동으로 결정된다.

<math>\phi(P_a)+\phi(P_b)-4\phi(P_i)+\phi(P_c)+\phi(P_d)=-\frac{h^2}{\epsilon_0}\rho(P_i)</math>

이렇게 격자점들을 정렬시키고 나면, <math>x</math>와 <math>b</math> 벡터, 그리고 행렬 <math>A</math>를 다음과 같이 생성할 수 있다.

<math>x_i=\phi(P_i)</math>, <math>b_i=-\frac{h^2}{\epsilon_0}\rho(P_i)</math>, <math>A_{i,j} = \begin{cases} -4, & \text{if }j=i \\ 1, & \text{if }j=a,b,c,d \\ 0, & \text{otherwise} \end{cases}</math>

단, 격자점이 주어진 영역의 경계에 있는 경우 방정식이 위에 기술한 것과 조금 달라져 행렬을 수정해야 한다. 이 주제는 경계 조건 항목에서 좀 더 자세하게 다룬다.

=== 경계 조건 ===

==== 디리클레 경계 ====
[[디리클레 경계 조건|디리클레 경계조건]]은 경계에서 퍼텐셜 값이 상수로 고정되는 조건이다. 따라서 연립방정식으로부터 퍼텐셜을 계산할 필요가 없는 위치이므로, 행렬에서 제외된다. 디리클레 경계와 인접한 격자점의 경우, 디리클레 경계의 퍼텐셜이 미지수가 아니므로, 방정식이 조금 수정되어야 한다. <math>P_n</math>과 인접하여 디리클레 경계 위의 점 <math>P_a</math>가 있을 때 방정식은 다음과 같이 수정된다. 이때 <math>\phi_b(P_a)</math>는 경계에서의 퍼텐셜로, 상수이다.

<math>\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)</math>

<math>P_a</math>가 아닌 다른 점들이 <math>P_n</math>과 인접하여도 같은 논리로 방정식을 수정해주면 된다.

벡터 <math>b</math>와 행렬 <math>A</math>는 디리클레 경계조건이 있을 때 다음과 같이 수정된다.

<math>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}</math>, <math>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}</math>
[[파일:노이만 경계조건과 추가 격자점.jpg|섬네일|노이만 경계조건에 따른 추가 격자점과 각 점의 좌표를 보여주고 있다.]]

==== 노이만 경계 ====
[[노이만 경계 조건|노이만 경계조건]]은 경계에서 퍼텐셜의 일계 미분 값이 상수로 고정되는 조건이다. 이 조건을 적용시키기 위해서는 경계 외부에 인접한 격자점을 추가하여, 경계 격자점에서의 일계 미분값이 일정하도록 만들어야 한다. 우측 그림과 같은 상황에서, 경계에서의 미분값은 일계 미분의 3점 공식에 의해 다음과 같이 주어진다.

<math>\phi_y(x_i,y_j) = \frac{1}{2h}[\phi(x_i,y_{j+1})-\phi(x_i,y_{j-1})]</math>

노이만 경계조건에서는 일계 미분 값이 상수로 정해져있으므로, 위의 관계로부터 외부 추가 격자점의 퍼텐셜인 <math>\phi(x_i,y_{j-1})</math>를 다음과 같이 계산할 수 있다.

<math>\phi(x_i,y_{j-1})=\phi(x_i,y_{j+1})-2h\phi_y(x_i,y_i)</math>

이렇게 구한 외부 추가 격자점의 퍼텐셜로 부터, 노이만 경계에서의 퍼텐셜을 일반 격자점과 같은 방법으로 계산 가능하다. 노이만 경계에서의 퍼텐셜 방정식은 다음과 같다.

<math>\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)</math>

이 방정식을 행렬에 포함시켜 주면 해가 유일하게 존재하는 행렬 방정식을 얻을 수 있다.

=== 선형 유전체 내부에서의 퍼텐셜 해 ===
[[파일:유전체 경계면.jpg|섬네일|유전체 경계면에서의 격자점과 각 점의 좌표를 보여주고 있다.]]
[[선형 유전체]]란 분극이 전기장에 비례하는 [[유전체 (물리학)|유전체]]를 의미한다. 이러한 [[유전체 (물리학)|유전체]]의 경우 [[유전체 (물리학)|유전체]]가 없는 진공인 상황에서 사용했던 방법을 그대로 사용하여 [[푸아송 방정식]]을 풀 수 있다. 앞서 기제된 방정식들의 [[유전율]] 항에 진공의 [[유전율]] <math>\epsilon_0</math>대신 선형 유전체 내부에서의 [[유전율]] <math>\epsilon</math>을 대입해주면 [[유전체 (물리학)|유전체]] 내부에서의 퍼텐셜 해를 구할 수 있다.

선형 유전체의 또 다른 특징은, 유전체 경계면에서 다음과 같은 경계조건을 만족한다는 것이다.

<math>\oint_C \epsilon \frac{\partial \phi}{\partial n}dl=-\int_S\rho_fda</math>

여기서 <math>C</math>는 <math>(x_i,y_j)</math> 주위를 둘러싸는 닫힌 곡선을 의미하고, <math>S</math>는 <math>C</math>를 경계로 가지는 곡면을 의미한다. 또한 <math>\frac{\partial \phi}{\partial n}</math>은 곡선에 수직한 방향으로의 퍼텐셜 미분을 의미한다. <math>C</math>를 우측 그림과 같이 각 격자점의 중심을 지나는 정사각형으로 잡고, 격자 사이의 간격이 충분히 작다고 가정하면, 일계 미분의 3점 공식에 의해 위의 식이 다음과 같이 근사된다. 이때 <math>\epsilon_{avg}</math>는 <math>\frac{\epsilon_1+\epsilon_2}{2}</math>를 의미한다.

<math>\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</math>

이 방정식을 이용하여 선형 유전체 경계면 위에 있는 격자점을 행렬에 포함시켜주면 퍼텐셜을 계산할 수 있다.

== 모멘트법 (Method of Moments) ==
모멘트법은 [[푸아송 방정식]]의 적분형 해를 근사적으로 계산하는 방법이다.

== 유한 요소법 (Finite Elements Method) ==
유한 요소법은 [[푸아송 방정식]]의 [[변분법|변분]]형 해를 근사적으로 계산하는 방법이다. 이 방법은 공간을 여러 조각으로 잘게 나누어 각 조각에서의 퍼텐셜을 선형 근사하는 방식을 사용한다. [[유한차분법|유한 차분법]]과는 달리 격자점을 자유롭게 선택할 수 있기 때문에 경계면이 곡면일 경우 유한 요소법을 사용하는 것이 유한 차분법을 사용하는 것 보다 효율적이다.

=== 변분법과 푸아송 방정식 ===
[[범함수]] <math>I</math>가 다음과 같이 정의 될 때, <math>I</math>를 최소화하는 함수 <math>V(x,y)</math>는 푸아송 방정식의 해가 된다.

<math>I(\phi)=\int\bigg[\frac{1}{2}\epsilon(\nabla V)^2-\rho V \bigg]da</math>

이때 피적분 함수의 첫 항은 [[전위 에너지|전기장의 에너지]]를 의미한다. 따라서, 전하 분포가 없을 때의 [[푸아송 방정식]]([[라플라스 방정식]])의 해는 계의 에너지를 최소화하는 퍼텐셜임을 알 수 있다.

=== 퍼텐셜의 선형 근사 ===
[[파일:유한요소법 격자.jpg|섬네일|유한 요소법의 격자와 각 격자점에서의 퍼텐셜을 나타내고 있다.]]
유한 요소법에서는 공간을 작은 요소들로 쪼갠 후 각 요소에서의 퍼텐셜을 선형 근사한다. 우측 그림과 같이 삼각형 요소로 공간을 분할하는 경우 (e)번째 요소에서의 선형 퍼텐셜은 다음과 같이 주어진다.

<math>V^{(e)}(x,y)=\begin{bmatrix} 1 & x & y \end{bmatrix}\begin{bmatrix} a^{(e)} \\ b^{(e)} \\ c^{(e)} \end{bmatrix}</math>

선형 근사된 퍼텐셜이 세 꼭짓점에서 원래의 퍼텐셜과 일치하기 위해서는 <math>a^{(e)}</math>, <math>b^{(e)}</math>, <math>c^{(e)}</math>가 다음 조건을 만족해야 한다. <math>R^{(e)}</math>는 행렬을 축약하기 위해 사용하는 기호이다.

<math>\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}</math>

<math>R^{(e)}</math> 행렬을 이용하면 요소에서의 퍼텐셜을 다음과 같이 조금 다르게 표현할 수 있다.

<math>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}</math>

여기서 <math>\alpha_i^{(e)}(x,y)</math>는 형상함수(shape function)로, 다음 관계를 만족한다.

<math>\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}</math>

=== 범함수의 계산 ===
각 요소의 범함수 <math>I^{(e)} </math>를 형상함수로 나타내면 다음과 같다.

<math>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] </math>

이때 <math>C^{(e)}_{ij}=\int\nabla\alpha_i^{(e)}\cdot\nabla\alpha_j^{(e)}da </math>로, <math>T^{(e)}_{ij}=\int\alpha_i^{(e)}\alpha_j^{(e)}da </math>로 정의하고, <math>\Rho^{(e)} = \rho^{(e)}\begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} </math>라고 정의하면, 요소의 범함수 <math>I^{(e)} </math>는 다음과 같이 간단해진다.

<math>I^{(e)}=\frac{1}{2}\epsilon V^{(e)T}C^{(e)}V^{(e)}-V^{(e)T}T^{(e)}\Rho^{(e)} </math>

여기서 윗 첨자 <math>T </math>는 [[전치행렬]]을 의미한다.

전체 범함수는 각 요소의 피적분 함수의 합이다.

<math>I=\sum_e I^{(e)} </math>

전체 범함수를 행렬 형태로 나타내기 위해서는 먼저 각 요소의 관점에서 본 꼭짓점 번호와 전체의 관점에서 본 격자점 번호를 대응시켜야 한다. 이 과정은 [[유한차분법|유한 차분법]]에서의 격자 좌표에 번호를 매기는 과정과 유사하다. 꼭짓점 번호를 격자점과 모두 대응시키고 나면, 각 요소 관점에서 만들어졌던 행렬(예를 들어 <math>C^{(e)} </math>)을 각 성분 번호에 대응되는 행렬(<math>C^{(e)}_G </math>)로 바꾸어주어야 한다. 전체 행렬은 <math>C=\sum_eC^{(e)}_G </math>와 같이 변환된 모든 행렬을 더해주어 만들 수 있다. 이와 같이 모든 행렬을 전체 행렬로 변환시켜 주면, 전체 범함수를 다음과 같이 표현할 수 있다.

<math>I=\frac{1}{2}\epsilon V^TCV - V^TT\Rho </math>

=== 범함수의 최소화 ===
범함수를 최소화하기 위해서는 다음 조건을 만족해야 한다.

<math>\frac{\partial I}{\partial V_i}=0 </math>, <math>V_i </math>는 <math>i </math>번째 격자점의 퍼텐셜

이를 모든 <math>i </math>에 대하여 묶어서 서술하면 아래와 같아진다.

<math>\begin{bmatrix} \frac{\partial }{\partial V_1} \\ \frac{\partial }{\partial V_2} \\ \vdots \end{bmatrix}I=\epsilon CV - T\Rho=0 </math>

이 방정식을 통해서 <math>V </math> 벡터를 구할 수 있다.

=== 경계 조건 ===
디리클레 경계 조건을 사용하는 경우, 경계 위에서의 격자점에서는 이미 퍼텐셜이 상수로 고정되어있다. 따라서 경계에서의 격자점들은 행렬 방정식에서 따로 분리해주어야 한다.

<math>C </math>, <math>V </math>, <math>T </math>, <math>\Rho </math> 행렬의 원소를 재배열하여 경계에서의 격자점들을 따로 분리해줄 수 있다. 예를 들어서 <math>V </math> 행렬의 경우 경계 격자점들의 퍼텐셜을 나타내는 <math>V_p </math>와 아직 알지 못하여 구해야 하는 퍼텐셜인 <math>V_f </math>로 나누어서 <math>V=\begin{bmatrix} V_f \\ V_p \end{bmatrix} </math>와 같이 나타낼 수 있다.

마찬가지로, <math>C </math> 행렬도 <math>C=\begin{bmatrix} C_{ff} & C_{fp} \\ C_{pf} & C_{pp} \end{bmatrix} </math>와 같이 [[블록 행렬]]로 나타낼 수 있다. 이렇게 [[블록 행렬]]로 나타내고 나면, 구해야 하는 퍼텐셜 <math>V_f </math>가 다음과 같이 주어지게 된다.

<math>C_{ff}V_f=-C_{fp}V_p+\frac{1}{\epsilon}[T_{ff}\Rho_f+T_{fp}\Rho_p] </math>

이 행렬 방정식은 <math>Ax=b</math> 꼴의 행렬 방정식으로, <math>C_{ff} </math>가 [[정사각행렬]]이므로, 특수한 경우를 제외하고는 풀 수 있다.

=== 푸아송 방정식의 해 ===
푸아송 방정식의 수치적 해는 행렬 방정식을 통해 구한 전체 <math>V </math> 벡터를 다시 각 요소에서의 <math>V^{(e)} </math>로 변환해주어 구할 수 있다. 해는 다음과 같다.

<math>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}</math>

== 같이 보기 ==
* [[계산물리학]]
* [[물리광학]]
== 참고 문헌 ==
* {{서적 인용|url=https://books.google.co.kr/books?hl=ko&lr=&id=2tOOCgAAQBAJ&oi=fnd&pg=PP1&dq=sadiku+numerical+electromagnetics&ots=KLHtLbd4Nx&sig=k36eClSk7Fi2vA8_x7ObXCAEVz4#v=onepage&q=sadiku%20numerical%20electromagnetics&f=false|제목=Numerical Techniques in Electromagnetics with MATLAB |판=3 |성=Sadiku|이름=Matthew N. O.|날짜=2011-05-31|출판사=CRC Press|언어=en|isbn=9781439883044}}

[[분류:수치해석학]]
[[분류:전자기학]]
[[분류:전자공학]]
[[분류:편미분 방정식]]