유한 차분법 이란, 미분 방정식을 유한한 차분의 선형 방정식으로 근사하여 해를 구하는 방법을 의미한다. 유한 차분법 을 이용하면 미분형으로 나타낸 푸아송 방정식 을 수치적으로 풀 수 있다. 푸아송 방정식 은 맥스웰 방정식 을 다르게 표현한 미분 방정식으로, 전기장 대신 전기 퍼텐셜 과 관련된 미분 방정식이다.
f
(
x
)
{\displaystyle f(x)}
x
0
{\displaystyle x_{0}}
테일러 전개 하면 다음과 같이 된다.
f
(
x
0
+
h
)
=
f
(
x
0
)
+
h
f
′
(
x
0
)
+
h
2
2
!
f
(
x
0
)
+
h
3
3
!
f
(
3
)
(
x
0
)
+
⋯
{\displaystyle 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
)
−
h
f
′
(
x
0
)
+
h
2
2
!
f
(
x
0
)
−
h
3
3
!
f
(
3
)
(
x
0
)
+
⋯
{\displaystyle 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
)
=
1
2
h
[
f
(
x
0
+
h
)
−
f
(
x
0
−
h
)
]
+
O
(
h
2
)
{\displaystyle f'(x_{0})={\frac {1}{2h}}[f(x_{0}+h)-f(x_{0}-h)]+O(h^{2})}
일계 미분 3점 공식을 유도할 때 사용하였던
f
(
x
)
{\displaystyle f(x)}
테일러 전개 를 다시 사용하여 이계 미분의 3점 공식을 유도할 수 있다. 두 식을 서로 더해 주고 정리하면, 이계 미분의 3점 공식을 얻을 수 있다.
f
″
(
x
0
)
=
1
h
2
[
f
(
x
0
+
h
)
−
2
f
(
x
0
)
+
f
(
x
0
−
h
)
]
+
O
(
h
2
)
{\displaystyle f''(x_{0})={\frac {1}{h^{2}}}[f(x_{0}+h)-2f(x_{0})+f(x_{0}-h)]+O(h^{2})}
이 식을 통해서 함수의 이계 미분값을 근사할 수 있고, 따라서 2차원 푸아송 방정식 을 3점 공식을 이용하여 근사하면 다음과 같아진다.
∇
2
ϕ
(
x
,
y
)
|
(
x
0
,
y
0
)
=
ϕ
x
x
(
x
0
,
y
0
)
+
ϕ
y
y
(
x
0
,
y
0
)
=
1
h
2
[
ϕ
(
x
0
+
h
,
y
0
)
+
ϕ
(
x
0
,
y
0
+
h
)
−
4
ϕ
(
x
0
,
y
0
)
+
ϕ
(
x
0
−
h
,
y
0
)
+
ϕ
(
x
0
,
y
0
−
h
)
]
+
O
(
h
2
)
{\displaystyle \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
{\displaystyle h}
1
h
2
[
ϕ
(
x
i
+
1
,
y
j
)
+
ϕ
(
x
i
−
1
,
y
j
)
−
4
ϕ
(
x
i
,
y
j
)
+
ϕ
(
x
i
,
y
j
+
1
)
+
ϕ
(
x
i
,
y
j
−
1
)
]
≅
∇
2
ϕ
(
x
,
y
)
|
(
x
i
,
y
j
)
=
−
ρ
(
x
i
,
y
j
)
ϵ
0
{\displaystyle {\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}}}}
위 방정식에서 좌항의
ϕ
(
x
i
,
y
j
)
{\displaystyle \phi (x_{i},y_{j})}
ρ
(
x
i
,
y
j
)
{\displaystyle \rho (x_{i},y_{j})}
ϕ
(
x
i
,
y
j
)
{\displaystyle \phi (x_{i},y_{j})}
연립일차방정식 과 같아진다. 연립일차방정식은 행렬 방정식으로 나타낼 수 있기 때문에, 결과적으로 주어진 문제를 푸는 것은 근사적으로
A
x
=
b
{\displaystyle Ax=b}
A
{\displaystyle A}
정사각행렬 로,
ϕ
(
x
i
,
y
j
)
{\displaystyle \phi (x_{i},y_{j})}
x
{\displaystyle x}
ϕ
(
x
i
,
y
j
)
{\displaystyle \phi (x_{i},y_{j})}
b
{\displaystyle b}
ρ
(
x
i
,
y
j
)
{\displaystyle \rho (x_{i},y_{j})}
위의 연립방정식으로부터 행렬을 생성하기 위해서는 각각의 격자점
(
x
i
,
y
j
)
{\displaystyle (x_{i},y_{j})}
(
x
0
,
y
0
)
{\displaystyle (x_{0},y_{0})}
(
x
0
,
y
1
)
{\displaystyle (x_{0},y_{1})}
(
x
0
,
y
2
)
{\displaystyle (x_{0},y_{2})}
(
x
i
,
y
j
)
=
P
n
{\displaystyle (x_{i},y_{j})=P_{n}}
a
{\displaystyle a}
b
{\displaystyle b}
c
{\displaystyle c}
d
{\displaystyle d}
i
{\displaystyle i}
ϕ
(
P
a
)
+
ϕ
(
P
b
)
−
4
ϕ
(
P
i
)
+
ϕ
(
P
c
)
+
ϕ
(
P
d
)
=
−
h
2
ϵ
0
ρ
(
P
i
)
{\displaystyle \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
{\displaystyle x}
b
{\displaystyle b}
A
{\displaystyle A}
x
i
=
ϕ
(
P
i
)
{\displaystyle x_{i}=\phi (P_{i})}
b
i
=
−
h
2
ϵ
0
ρ
(
P
i
)
{\displaystyle b_{i}=-{\frac {h^{2}}{\epsilon _{0}}}\rho (P_{i})}
A
i
,
j
=
{
−
4
,
if
j
=
i
1
,
if
j
=
a
,
b
,
c
,
d
0
,
otherwise
{\displaystyle 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
{\displaystyle P_{n}}
P
a
{\displaystyle P_{a}}
ϕ
b
(
P
a
)
{\displaystyle \phi _{b}(P_{a})}
ϕ
(
P
b
)
−
4
ϕ
(
P
n
)
+
ϕ
(
P
c
)
+
ϕ
(
P
d
)
=
−
h
2
ϵ
0
ρ
(
P
n
)
−
ϕ
b
(
P
a
)
{\displaystyle \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
{\displaystyle P_{a}}
P
n
{\displaystyle P_{n}}
벡터
b
{\displaystyle b}
A
{\displaystyle A}
b
i
=
{
−
h
2
ϵ
0
ρ
(
P
i
)
,
if
P
i
is not adjacent to the dirichlet boundaries
−
h
2
ϵ
0
ρ
(
P
i
)
−
∑
ϕ
b
(
P
j
)
,
otherwise
{\displaystyle 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
=
{
−
4
,
if
j
=
i
1
,
if
j
=
k
and
P
k
is not adjacent to the dirichlet boundaries,
k
=
a
,
b
,
c
,
d
0
,
otherwise
{\displaystyle 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점 공식에 의해 다음과 같이 주어진다.
ϕ
y
(
x
i
,
y
j
)
=
1
2
h
[
ϕ
(
x
i
,
y
j
+
1
)
−
ϕ
(
x
i
,
y
j
−
1
)
]
{\displaystyle \phi _{y}(x_{i},y_{j})={\frac {1}{2h}}[\phi (x_{i},y_{j+1})-\phi (x_{i},y_{j-1})]}
노이만 경계조건에서는 일계 미분 값이 상수로 정해져있으므로, 위의 관계로부터 외부 추가 격자점의 퍼텐셜인
ϕ
(
x
i
,
y
j
−
1
)
{\displaystyle \phi (x_{i},y_{j-1})}
ϕ
(
x
i
,
y
j
−
1
)
=
ϕ
(
x
i
,
y
j
+
1
)
−
2
h
ϕ
y
(
x
i
,
y
i
)
{\displaystyle \phi (x_{i},y_{j-1})=\phi (x_{i},y_{j+1})-2h\phi _{y}(x_{i},y_{i})}
이렇게 구한 외부 추가 격자점의 퍼텐셜로 부터, 노이만 경계에서의 퍼텐셜을 일반 격자점과 같은 방법으로 계산 가능하다. 노이만 경계에서의 퍼텐셜 방정식은 다음과 같다.
ϕ
(
x
i
+
1
,
y
j
)
+
ϕ
(
x
i
−
1
,
y
j
)
−
4
ϕ
(
x
i
,
y
j
)
+
ϕ
(
x
i
,
y
j
+
1
)
+
ϕ
(
x
i
,
y
j
−
1
)
=
−
h
2
ϵ
0
ρ
(
x
i
,
y
j
)
{\displaystyle \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})}
이 방정식을 행렬에 포함시켜 주면 해가 유일하게 존재하는 행렬 방정식을 얻을 수 있다.
유전체 경계면에서의 격자점과 각 점의 좌표를 보여주고 있다. 선형 유전체 란 분극이 전기장에 비례하는 유전체 를 의미한다. 이러한 유전체 의 경우 유전체 가 없는 진공인 상황에서 사용했던 방법을 그대로 사용하여 푸아송 방정식 을 풀 수 있다. 앞서 기제된 방정식들의 유전율 항에 진공의 유전율
ϵ
0
{\displaystyle \epsilon _{0}}
유전율
ϵ
{\displaystyle \epsilon }
유전체 내부에서의 퍼텐셜 해를 구할 수 있다.
선형 유전체의 또 다른 특징은, 유전체 경계면에서 다음과 같은 경계조건을 만족한다는 것이다.
∮
C
ϵ
∂
ϕ
∂
n
d
l
=
−
∫
S
ρ
f
d
a
{\displaystyle \oint _{C}\epsilon {\frac {\partial \phi }{\partial n}}dl=-\int _{S}\rho _{f}da}
여기서
C
{\displaystyle C}
(
x
i
,
y
j
)
{\displaystyle (x_{i},y_{j})}
S
{\displaystyle S}
C
{\displaystyle C}
∂
ϕ
∂
n
{\displaystyle {\frac {\partial \phi }{\partial n}}}
C
{\displaystyle C}
ϵ
a
v
g
{\displaystyle \epsilon _{avg}}
ϵ
1
+
ϵ
2
2
{\displaystyle {\frac {\epsilon _{1}+\epsilon _{2}}{2}}}
ϵ
a
v
g
[
ϕ
(
x
i
+
1
,
y
j
)
+
ϕ
(
x
i
−
1
,
y
j
)
−
4
ϕ
(
x
i
,
y
j
)
]
+
ϵ
1
ϕ
(
x
i
,
y
j
+
1
)
+
ϵ
2
ϕ
(
x
i
,
y
j
−
1
)
≅
−
ρ
f
(
x
i
,
y
j
)
h
2
{\displaystyle \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}}
이 방정식을 이용하여 선형 유전체 경계면 위에 있는 격자점을 행렬에 포함시켜주면 퍼텐셜을 계산할 수 있다.
유한 요소법은 푸아송 방정식 의 변분 형 해를 근사적으로 계산하는 방법이다. 이 방법은 공간을 여러 조각으로 잘게 나누어 각 조각에서의 퍼텐셜을 선형 근사하는 방식을 사용한다. 유한 차분법 과는 달리 격자점을 자유롭게 선택할 수 있기 때문에 경계면이 곡면일 경우 유한 요소법을 사용하는 것이 유한 차분법을 사용하는 것 보다 효율적이다.
범함수
I
{\displaystyle I}
I
{\displaystyle I}
V
(
x
,
y
)
{\displaystyle V(x,y)}
I
(
ϕ
)
=
∫
[
1
2
ϵ
(
∇
V
)
2
−
ρ
V
]
d
a
{\displaystyle I(\phi )=\int {\bigg [}{\frac {1}{2}}\epsilon (\nabla V)^{2}-\rho V{\bigg ]}da}
이때 피적분 함수의 첫 항은 전기장의 에너지 를 의미한다. 따라서, 전하 분포가 없을 때의 푸아송 방정식 (라플라스 방정식 )의 해는 계의 에너지를 최소화하는 퍼텐셜임을 알 수 있다.
유한 요소법의 격자와 각 격자점에서의 퍼텐셜을 나타내고 있다. 유한 요소법에서는 공간을 작은 요소들로 쪼갠 후 각 요소에서의 퍼텐셜을 선형 근사한다. 우측 그림과 같이 삼각형 요소로 공간을 분할하는 경우 (e)번째 요소에서의 선형 퍼텐셜은 다음과 같이 주어진다.
V
(
e
)
(
x
,
y
)
=
[
1
x
y
]
[
a
(
e
)
b
(
e
)
c
(
e
)
]
{\displaystyle V^{(e)}(x,y)={\begin{bmatrix}1&x&y\end{bmatrix}}{\begin{bmatrix}a^{(e)}\\b^{(e)}\\c^{(e)}\end{bmatrix}}}
선형 근사된 퍼텐셜이 세 꼭짓점에서 원래의 퍼텐셜과 일치하기 위해서는
a
(
e
)
{\displaystyle a^{(e)}}
b
(
e
)
{\displaystyle b^{(e)}}
c
(
e
)
{\displaystyle c^{(e)}}
R
(
e
)
{\displaystyle R^{(e)}}
[
V
1
(
e
)
V
2
(
e
)
V
3
(
e
)
]
=
[
1
x
1
(
e
)
y
1
(
e
)
1
x
2
(
e
)
y
2
(
e
)
1
x
3
(
e
)
y
3
(
e
)
]
[
a
(
e
)
b
(
e
)
c
(
e
)
]
=
R
(
e
)
[
a
(
e
)
b
(
e
)
c
(
e
)
]
{\displaystyle {\begin{bmatrix}V_{1}^{(e)}\\V_{2}^{(e)}\\V_{3}^{(e)}\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
)
{\displaystyle R^{(e)}}
V
(
e
)
(
x
,
y
)
=
[
1
x
y
]
R
(
e
)
−
1
R
(
e
)
[
a
(
e
)
b
(
e
)
c
(
e
)
]
=
[
α
1
(
e
)
α
2
(
e
)
α
3
(
e
)
]
[
V
1
(
e
)
V
2
(
e
)
V
3
(
e
)
]
{\displaystyle 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_{1}^{(e)}\\V_{2}^{(e)}\\V_{3}^{(e)}\end{bmatrix}}}
여기서
α
i
(
e
)
(
x
,
y
)
{\displaystyle \alpha _{i}^{(e)}(x,y)}
α
i
(
e
)
(
x
j
(
e
)
,
y
j
(
e
)
)
=
δ
i
j
=
{
1
,
if
i
=
j
0
,
if
i
≠
j
{\displaystyle \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
)
{\displaystyle I^{(e)}}
I
(
e
)
=
∑
i
∑
j
[
1
2
ϵ
V
i
(
e
)
V
j
(
e
)
∫
∇
α
i
(
e
)
⋅
∇
α
j
(
e
)
d
a
−
V
i
(
e
)
ρ
(
e
)
∫
α
i
(
e
)
α
j
(
e
)
d
a
]
{\displaystyle I^{(e)}=\sum _{i}\sum _{j}{\bigg [}{\frac {1}{2}}\epsilon V_{i}^{(e)}V_{j}^{(e)}\int \nabla \alpha _{i}^{(e)}\cdot \nabla \alpha _{j}^{(e)}da-V_{i}^{(e)}\rho ^{(e)}\int \alpha _{i}^{(e)}\alpha _{j}^{(e)}da{\bigg ]}}
이때
C
i
j
(
e
)
=
∫
∇
α
i
(
e
)
⋅
∇
α
j
(
e
)
d
a
{\displaystyle C_{ij}^{(e)}=\int \nabla \alpha _{i}^{(e)}\cdot \nabla \alpha _{j}^{(e)}da}
T
i
j
(
e
)
=
∫
α
i
(
e
)
α
j
(
e
)
d
a
{\displaystyle T_{ij}^{(e)}=\int \alpha _{i}^{(e)}\alpha _{j}^{(e)}da}
P
(
e
)
=
ρ
(
e
)
[
1
1
1
]
{\displaystyle \mathrm {P} ^{(e)}=\rho ^{(e)}{\begin{bmatrix}1\\1\\1\end{bmatrix}}}
I
(
e
)
{\displaystyle I^{(e)}}
I
(
e
)
=
1
2
ϵ
V
(
e
)
T
C
(
e
)
V
(
e
)
−
V
(
e
)
T
T
(
e
)
P
(
e
)
{\displaystyle I^{(e)}={\frac {1}{2}}\epsilon V^{(e)T}C^{(e)}V^{(e)}-V^{(e)T}T^{(e)}\mathrm {P} ^{(e)}}
여기서 윗 첨자
T
{\displaystyle T}
전치행렬 을 의미한다.
전체 범함수는 각 요소의 피적분 함수의 합이다.
I
=
∑
e
I
(
e
)
{\displaystyle I=\sum _{e}I^{(e)}}
전체 범함수를 행렬 형태로 나타내기 위해서는 먼저 각 요소의 관점에서 본 꼭짓점 번호와 전체의 관점에서 본 격자점 번호를 대응시켜야 한다. 이 과정은 유한 차분법 에서의 격자 좌표에 번호를 매기는 과정과 유사하다. 꼭짓점 번호를 격자점과 모두 대응시키고 나면, 각 요소 관점에서 만들어졌던 행렬(예를 들어
C
(
e
)
{\displaystyle C^{(e)}}
C
G
(
e
)
{\displaystyle C_{G}^{(e)}}
C
=
∑
e
C
G
(
e
)
{\displaystyle C=\sum _{e}C_{G}^{(e)}}
I
=
1
2
ϵ
V
T
C
V
−
V
T
T
P
{\displaystyle I={\frac {1}{2}}\epsilon V^{T}CV-V^{T}T\mathrm {P} }
범함수를 최소화하기 위해서는 다음 조건을 만족해야 한다.
∂
I
∂
V
i
=
0
{\displaystyle {\frac {\partial I}{\partial V_{i}}}=0}
V
i
{\displaystyle V_{i}}
i
{\displaystyle i}
이를 모든
i
{\displaystyle i}
[
∂
∂
V
1
∂
∂
V
2
⋮
]
I
=
ϵ
C
V
−
T
P
=
0
{\displaystyle {\begin{bmatrix}{\frac {\partial }{\partial V_{1}}}\\{\frac {\partial }{\partial V_{2}}}\\\vdots \end{bmatrix}}I=\epsilon CV-T\mathrm {P} =0}
이 방정식을 통해서
V
{\displaystyle V}
디리클레 경계 조건을 사용하는 경우, 경계 위에서의 격자점에서는 이미 퍼텐셜이 상수로 고정되어있다. 따라서 경계에서의 격자점들은 행렬 방정식에서 따로 분리해주어야 한다.
C
{\displaystyle C}
V
{\displaystyle V}
T
{\displaystyle T}
P
{\displaystyle \mathrm {P} }
V
{\displaystyle V}
V
p
{\displaystyle V_{p}}
V
f
{\displaystyle V_{f}}
V
=
[
V
f
V
p
]
{\displaystyle V={\begin{bmatrix}V_{f}\\V_{p}\end{bmatrix}}}
마찬가지로,
C
{\displaystyle C}
C
=
[
C
f
f
C
f
p
C
p
f
C
p
p
]
{\displaystyle C={\begin{bmatrix}C_{ff}&C_{fp}\\C_{pf}&C_{pp}\end{bmatrix}}}
블록 행렬 로 나타낼 수 있다. 이렇게 블록 행렬 로 나타내고 나면, 구해야 하는 퍼텐셜
V
f
{\displaystyle V_{f}}
C
f
f
V
f
=
−
C
f
p
V
p
+
1
ϵ
[
T
f
f
P
f
+
T
f
p
P
p
]
{\displaystyle C_{ff}V_{f}=-C_{fp}V_{p}+{\frac {1}{\epsilon }}[T_{ff}\mathrm {P} _{f}+T_{fp}\mathrm {P} _{p}]}
이 행렬 방정식은
A
x
=
b
{\displaystyle Ax=b}
C
f
f
{\displaystyle C_{ff}}
정사각행렬 이므로, 특수한 경우를 제외하고는 풀 수 있다.
푸아송 방정식의 수치적 해는 행렬 방정식을 통해 구한 전체
V
{\displaystyle V}
V
(
e
)
{\displaystyle V^{(e)}}
V
(
e
)
(
x
,
y
)
=
[
α
1
(
e
)
α
2
(
e
)
α
3
(
e
)
]
[
V
1
(
e
)
V
2
(
e
)
V
3
(
e
)
]
{\displaystyle V^{(e)}(x,y)={\begin{bmatrix}\alpha _{1}^{(e)}&\alpha _{2}^{(e)}&\alpha _{3}^{(e)}\end{bmatrix}}{\begin{bmatrix}V_{1}^{(e)}\\V_{2}^{(e)}\\V_{3}^{(e)}\end{bmatrix}}}
![{\displaystyle f'(x_{0})={\frac {1}{2h}}[f(x_{0}+h)-f(x_{0}-h)]+O(h^{2})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/57e9f60299d08c9be4a9e166050cdbf11d60bcd4)
![{\displaystyle f''(x_{0})={\frac {1}{h^{2}}}[f(x_{0}+h)-2f(x_{0})+f(x_{0}-h)]+O(h^{2})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/dc5cbce832f8514e9ceb20aad5e2d3e0139f52da)
![{\displaystyle \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})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b0238e6b5736029af58948369aa023dffdc50cdb)
![{\displaystyle {\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}}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c096f79693a01da002fd690360fef7105df4819b)
![{\displaystyle \phi _{y}(x_{i},y_{j})={\frac {1}{2h}}[\phi (x_{i},y_{j+1})-\phi (x_{i},y_{j-1})]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3233d8ea848f4360fb71be686cb9263c769f6264)
![{\displaystyle \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}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/92654386d09dc168d9b6972dbfc1fc2b1daf0142)
![{\displaystyle I(\phi )=\int {\bigg [}{\frac {1}{2}}\epsilon (\nabla V)^{2}-\rho V{\bigg ]}da}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8f19992eec33b00a8631c1cf055348af47db09f1)
![{\displaystyle I^{(e)}=\sum _{i}\sum _{j}{\bigg [}{\frac {1}{2}}\epsilon V_{i}^{(e)}V_{j}^{(e)}\int \nabla \alpha _{i}^{(e)}\cdot \nabla \alpha _{j}^{(e)}da-V_{i}^{(e)}\rho ^{(e)}\int \alpha _{i}^{(e)}\alpha _{j}^{(e)}da{\bigg ]}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/50e7ef27243f461993ad77f8f0116e5b08b089cb)
![{\displaystyle C_{ff}V_{f}=-C_{fp}V_{p}+{\frac {1}{\epsilon }}[T_{ff}\mathrm {P} _{f}+T_{fp}\mathrm {P} _{p}]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/dd144e0dacc4ec6e3fdff84e3e84838406164684)
