개시하다

물리학 시뮬레이션에서 미분 방정식을 풀는 것은 매우 중요하고 물리 현상은 대체로 다음과 같은 형식으로 나눌 수 있다.
코드는 이쪽에서입니다.
$Poisson$방정식(이번 처리)

\nabla^2 f = \phi (x)

$Laplace$방정식

\nabla^2 f = 0

확산 방정식

\nabla^2 f = C \frac{\partial f}{\partial t}

파동 방정식

\nabla^2 f = C \frac{\partial^2 f}{\partial t^2}

차분법

미분 방정식을 처리할 때 미분을 진행해야 한다. 그러나, 보통 미분을 실제적으로 진행할 수 없다. 따라서 미분의 근사성을 고려할 때 차분법을 사용한다.
$f(x)가 $n회 연속으로 영역을 정의하는 함수인 경우 $Taylor$의 확장을 다음과 같이 정의할 수 있습니다.

f(x+\Delta x) = \sum_{n=0}^{\infty} \frac{f^{(n)}(x)}{n!}\Delta x^n ・・・(A)\\
f(x-\Delta x) = \sum_{n=0}^{\infty} \frac{f^{(n)}(x)}{n!}(-1)^n \Delta x^n ・・・(B)

(A) 변형 후 전진차.

f(x+\Delta x) - f(x) = \Delta xf'(x) + o(\Delta x^2) \\
f'(x) = \frac{f(x+\Delta x) - f(x)}{\Delta x} + o(\Delta x) ・・・(A)'\\

(B) 변형 후 후퇴차.

f(x-\Delta x) - f(x) = -\Delta xf'(x) + o(\Delta x) \\
f'(x) = \frac{f(x) - f(x-\Delta x)}{\Delta x} + o(\Delta x) ・・・(B)'\\

(A) 카페(B)가 변형되면 중심차가 생긴다.

f(x+\Delta x) - f(x-\Delta x) = 2\Delta x f'(x) + o(\Delta x^3) \\
f'(x) = \frac{f(x + \Delta x) - f(x - \Delta x)}{2\Delta x} + o(\Delta x^2)・・・(C)

(A)+(B) 변형 후 2 단계 미분의 중심차

f(x+\Delta x) + f(x-\Delta x) = 2f(x) + \Delta x^2f''(x)+o(\Delta x^3) \\
f''(x) = \frac{f(x+\Delta x) - 2f(x) + f(x-\Delta x)}{\Delta x^2}  + o(\Delta x)・・・(D)

이번 $\Delta×달러는 아주 작습니다. $o(...)$의 부분은 무시할 수 있는 가설 아래 위의 미분 근사 공식을 사용한다.

프로그램에 빠지다

이번에는 간단하게 보기 위해 2차원 $Laplace$방정식을 진행합니다.

\nabla^2 f = \phi (x) \\
\frac{\partial ^2f}{\partial x^2} + \frac{\partial ^2f}{\partial y^2} = \phi (x, y)

공간 분할

우선, 차분법을 사용하기 위해서는 공간을 격자 모양으로 나누어야 합니다. 따라서 $0\leq×\leq(m+1)\Deltax,\0\leqy\leq(n+1)\Deltay달러($m,\n달러는 자연수),\leqj\leqn+1$$(i\Deltax,\j\Deltay)는 $에서만 정의됩니다. 시뮬레이션을 해결하기 위해서는 경계 조건이 필요합니다. 따라서 $i=0, m+1달러 또는 $j=0, n+1달러의 격자에 있는 점의 값은 $f달러로 알고 있습니다.
또한 $i, j$에 따라 격자에 있는 점을 유일하게 구할 수 있기 때문에 $r{i,j}=(i-1)+m(j-1)$. 또한, 간단하도록 $f(i\Deltax,\j\Deltay)=f{r{i,j}달러입니다.

위 그림의 옅은 파란색 부분은 경계 조건이고 격자 위의 점은 알 수 없는 점이다.

공식 근사

(D) 의 근사 공식을 사용하여

\frac{\partial ^2f_{r_{i,j}}}{\partial x^2} \simeq \frac{f_{r_{i+1,j}} - 2f_{r_{i,j}} + f_{r_{i-1,j}}}{\Delta x^2} \\
\frac{\partial ^2f_{r_{i,j}}}{\partial y^2} \simeq \frac{f_{r_{i,j+1}} - 2f_{r_{i,j}} + f_{r_{i,j-1}}}{\Delta y^2} 

이 두 공식을 더하면

\phi _{r_{i,j}} = \frac{\partial ^2f_{r_{i,j}}}{\partial x^2} + \frac{\partial ^2f_{r_{i,j}}}{\partial y^2} 
\simeq \frac{f_{r_{i+1,j}} - 2f_{r_{i,j}} + f_{r_{i-1,j}}}{\Delta x^2} + \frac{f_{r_{i,j+1}} - 2f_{r_{i,j}} + f_{r_{i,j-1}}}{\Delta y^2} \\
= \frac{1}{\Delta x^2}f_{r_{i+1,j}} + \frac{1}{\Delta x^2}f_{r_{i-1,j}} + \Biggl(\frac{-2}{\Delta x^2}+\frac{-2}{\Delta y^2}\Biggl)f_{r_{i,j}} + \frac{1}{\Delta y^2}f_{r_{i,j-1}} + \frac{1}{\Delta y^2}f_{r_{i,j+1}} 

얻다{i,j}=0,1,2,...,mn-1달러는 모든 물건을 요구할 수 있습니다. 또한, 아래는 $c입니다.1=\frac{1}{\Delta x^2},c_2=\frac{1}{\Delta y^2},c_3=\Bigl(\rac{-2}{Deltax^2]+\rac{2-}{Deltay^2}\Bigl)$로 설정합니다.

연립 방정식의 구축

위의 모델에 따라 다음 방정식을 입식으로 풀고 $\boldsymbol {f} (f 0, f 1,..., f {mn-1}) 를 결정하십시오.

\begin{eqnarray}
\left\{ \begin{array}{l}

c_1 f_{r_{2,1}} + c_1 f_{r_{0,1}} +c_3 f_{r_{1,1}} + c_2 f_{r_{1,0}} + c_2 f_{r_{1,2}} = \phi _{r_{1,1}} \\

c_1 f_{r_{3,1}} + c_1 f_{r_{1,1}} +c_3 f_{r_{2,1}} + c_2 f_{r_{2,0}} + c_2 f_{r_{2,2}} = \phi _{r_{2,1}} \\
. \\
. \\
. \\
c_1 f_{r_{m,n}} + c_1 f_{r_{m-2,n}} +c_3 f_{r_{m - 1,n}} + c_2 f_{r_{m-1,n-1}} + c_2 f_{r_{m-1,n+1}} = \phi _{r_{m - 1,n}} \\

c_1 f_{r_{m+1,n}} + c_1 f_{r_{m-1,n}} +c_3 f_{r_{m,n}} + c_2 f_{r_{m,n-1}} + c_2 f_{r_{m,n+1}} = \phi _{r_{m,n}} 

\end{array} \right.
\end{eqnarray}

그러나 $i=0, m+1달러 또는 $j=0, n+1달러라는 점에서 $f는 경계 조건보다 알고 있기 때문에 오른쪽에 있는 이미 알고 있는 쪽으로 옮깁니다.
그리고

\begin{eqnarray}
\left\{ \begin{array}{l}

c_1 f_{r_{2,1}} +c_3 f_{r_{1,1}} + c_2 f_{r_{1,2}} = \phi _{r_{1,1}} - c_1 f_{r_{0,1}} - c_2 f_{r_{1,0}} \\

c_1 f_{r_{3,1}} + c_1 f_{r_{1,1}} +c_3 f_{r_{2,1}} + c_2 f_{r_{2,2}} = \phi _{r_{2,1}} - c_2 f_{r_{2,0}}\\
. \\
. \\
. \\
c_1 f_{r_{m,n}} + c_1 f_{r_{m-2,n}} +c_3 f_{r_{m - 1,n}} + c_2 f_{r_{m-1,n-1}} = \phi _{r_{m - 1,n}} - c_2 f_{r_{m-1,n+1}}\\

c_1 f_{r_{m-1,n}} +c_3 f_{r_{m,n}} + c_2 f_{r_{m,n-1}} = \phi _{r_{m,n}} - c_1 f_{r_{m+1,n}} - c_2 f_{r_{m,n+1}}

\end{array} \right.
\end{eqnarray}

매트릭스 형식으로 변형가능,

L\boldsymbol{f} = \boldsymbol{\psi}

이 밖에 이 행렬과 벡터 요소들은

\begin{eqnarray}
L_{r_{i,j},r_{i',j'}}=\left\{ \begin{array}{ll}
c_1 & (|r_{i,j} - r_{i',j'}| = 1)\\
c_2 & (|r_{i,j} - r_{i',j'}| = m)\\
c_3 & (|r_{i,j} - r_{i',j'}| = 0)\\
0 & (otherwise)\\
\end{array} \right.
\end{eqnarray}

\begin{eqnarray}
\boldsymbol{\psi_{r_{i,j}}}= \phi_{r_{i,j}} - \left\{ \begin{array}{ll}
c_1 f_{r_{0,j}} & (i = 0)\\
c_1 f_{r_{m,j}} & (i = m)\\
c_2 f_{r_{i,0}} & (j = 0)\\
c_2 f_{r_{i,n}} & (j = n)\\
0 & (otherwise)\\
\end{array} \right.
\end{eqnarray}

오른쪽 길이는 조건에 따라 필요한 값을 줄여서 원소의 값을 얻을 수 있습니다. 이것은 $Poisson 방정식의 입식을 완성합니다. 그리고 $L의 역행렬을 찾아서 $\boldsymbol {f}달러를 얻을 수 있습니다.

실제 문제를 해결하다

위의 모델링을 통해 메쉬 크기 $\Deltax,\Deltay$, 메쉬 수 $m, n$, 경계 조건 $f${r_{0,j}},f_{r_{m+1,j}},f_{r_{i,0}},f_{r{i, n+1}달러만 주면 정의된 공간에서 $\boldsymbol {f}달러를 얻을 수 있습니다.
이번에는 여기. 중의 열전도 온도 분포 프로젝트를 모델링하는 데 주력한다.
열원 분포는 $s(x)$이고 열전도율은 $\lambda이며 온도 분포는 $T일 때의 관계입니다.

\nabla^2T = - \frac{s}{\lambda} = 0

이번에는 특정한 열원이 없다. 즉, $s(x)=0$의 열확산을 시뮬레이션한 것이다.

LT = \boldsymbol{\psi}

$L^{-1}달러를 왼쪽에서 양쪽으로 걸면 안정 상태의 온도 분포 $T를 발견할 수 있습니다.

위의 경계 조건과 열원이 없는 상태에서 안정적인 온도 분포를 확정한다.

진실을 추구해야 이런 느낌이 든다.