This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.

<span class="iconify mr-1 svg-loading" data-icon="mdi:arrow-left"></span>======전자기학에서====== 경계 $\Gamma$로 둘러싸인 공간 $\Omega$에서 푸아송 방정식 $$\nabla^2 u = -f$$ 이 정의되어 있다. 디리클레 경계조건 $\left. u \right|_\Gamma = u_0$ 혹은 노이만 경계조건 $\left. \partial u / \partial n \right|_\Gamma = \nabla u \cdot \hat{n} = 0$이 만족된다고 가정한다. 이러한 $u$를 찾는 것은 다음 적분을 최소화하는 해 $u = \arg\min_{w} J\left[w\right]$를 찾는 것과 동치이다. $$J\left[w\right] \equiv \frac{1}{2} \int_\Omega \left| \nabla w \right|^2 - \int_\Omega f w ~dA.$$ 우리의 추정해인 $w$가 그러한 답 $u$에서 벗어난 정도가 $v$라고 하자. 경계조건을 만족해야 하므로 디리클레라면 $\left. v \right|_\Gamma = 0$일 것이다. 적분을 해보면 \begin{eqnarray*} J\left[ u+v \right] &=& \frac{1}{2} \int_\Omega \left| \nabla u \right|^2 dA - \int_\Omega fu ~dA\\ &&+\int_\Omega \nabla u \cdot \nabla v ~dA - \int_\Omega fv ~dA\\ && + \frac{1}{2} \int_\Omega \left| \nabla v \right|^2 dA \end{eqnarray*} 여기에서 세 번째와 네 번째 항을 따로 적어보자. 발산정리를 활용하면 다음처럼 고쳐 쓸 수 있다: \begin{eqnarray*} \delta J &=& \int_\Omega \nabla u \cdot \nabla v ~dA - \int_\Omega fv ~dA\\ &=& \int_\Omega \nabla \cdot \left( v\nabla u \right) dA - \int_\Omega v \nabla^2 u ~dA - \int_\Omega fv ~dA\\ &=& \int_\Gamma v \nabla u \cdot \hat{n} dl + \int_\Omega v \left(\nabla^2 u - f\right) ~dA. \end{eqnarray*} 두 항 모두 0임을 알 수 있다. 따라서 다음 결론을 얻는다: \begin{eqnarray*} J\left[ u+v \right] &=& \frac{1}{2} \int_\Omega \left| \nabla u \right|^2 dA - \int_\Omega fu ~dA + \frac{1}{2} \int_\Omega \left| \nabla v \right|^2 dA\\ &\ge& \frac{1}{2} \int_\Omega \left| \nabla u \right|^2 dA - \int_\Omega fu ~dA\\ &=& J\left[ u \right]. \end{eqnarray*} =====예제===== 푸아송 방정식 $\nabla^2 V = -\rho_0 = const.$가 2차원 공간 $\Omega \equiv [-1,1] \times [-1,1]$에서 정의되어 있다. 경계조건은 $V(\pm 1, y) = V(x, \pm 1)=0$이다. 시험 함수 $\Phi$를 아래의 기저로 나타내자: $$u_{mn} = \left(1-x^2 \right) \left(1-y^2 \right) \left( x^{2m} y^{2n} + x^{2n} y^{2m} \right).$$ 여기에서 $m,n = 0,1,2,\ldots$이다. 1차 근사로서 $m=n=0$인 하나의 기저만 사용해보자. 그러면 시험 해는 $\Phi = a_1 u_1$이며, 변분법을 활용하기 위해 이를 적분에 집어 넣으면 \begin{eqnarray*} J &=& \frac{1}{2} \int_\Omega \left| \nabla \Phi \right|^2 dA - \int_\Omega \rho_0 \Phi ~dA\\ &=& \frac{a_1^2}{2} \int_{-1}^1 \int_{-1}^1 \left[ 4x^2 \left(1-y^2\right)^2 + \left(1-x^2\right) 4y^2 \right] dx dy - \rho_0 a_1 \int_{-1}^1 \int_{-1}^1 \left(1-x^2\right) \left(1-y^2\right) dx dy\\ &=& \frac{128}{45} a_1^2 - \frac{16}{9} \rho_0 a_1. \end{eqnarray*} $a_1$에 대한 이 이차함수의 최소를 구해보면 $a_1 = \frac{5}{16} \rho_0$임을 얻는다. $m=n=1$인 기저까지 활용해서 $\Phi = a_1 \left(1-x^2 \right) \left(1-y^2 \right) + a_2 \left(1-x^2 \right) \left(1-y^2 \right) \left( x^2 y^2 + x^2 y^2 \right)$로도 연습하고 해의 모양이 얼마나 바뀌었는지 확인해볼 수 있다. 답은 $a_1 = \frac{1295}{4432} \rho_0$와 $a_2 = \frac{525}{8864} \rho_0$이다. ======양자역학에서====== =====변분을 통해 에너지 준위 찾기===== 1차원의 시간 무관 슈뢰딩거 방정식 $H \psi_n = E_n \psi_n$을 고려하자 ($n=0,1,2,\ldots$). 임의의 파동함수 $\psi$에 대해 에너지의 기대값은 $$\left< H \right> = \frac{\int \psi^\ast H \psi ~dx}{\int \psi^\ast \psi ~dx}$$ 인데 이 값은 언제나 바닥 상태 에너지 $E_0$보다 크거나 같으며, 바닥 상태가 유일하면 $\psi = \psi_0$일 때에 성립한다. 그 이유는 아래와 같다: 먼저 $\psi$를 에너지 고유벡터들로 분해하여 쓰자: $\psi = \sum_n c_n \psi_n$. 그러면 \begin{eqnarray} \left< H \right> &=& \frac{\int \psi^\ast H \psi ~dx}{\int \psi^\ast \psi ~dx} = \frac{\sum_n |c_n|^2 E_n}{\sum_n |c_n|^2}\\ &=& E_0 + \frac{\sum_n |c_n|^2 (E_n - E_0)}{\sum_n |c_n|^2} \ge E_0. \end{eqnarray} 그러니까 개념적으로는 시험 삼아 파동함수를 적어보고, 그 모양을 잡아주는 맺음변수들 $\alpha_1, \alpha_2, \ldots, \alpha_r$을 이리저리 바꿔가면서 에너지의 기대값을 최소화하면 될 것이다. 수식으로 적어보면 $\psi = \psi(x;\alpha_1, \alpha_2, \ldots, \alpha_r)$에 대해서 $$E = \frac{\int \psi^\ast H \psi ~dx}{\int \psi^\ast \psi ~dx}$$ 가 정의되고 이를 편미분해서 $$\frac{\partial E}{\partial \alpha_1} = \frac{\partial E}{\partial \alpha_2} = \ldots = \frac{\partial E}{\partial \alpha_r} = 0$$ 으로 하면 될 것이다. =====라그랑주 곱수의 사용===== 실제 계산에서는 라그랑주 곱수를 사용해서 $E$의 분모에 표현된 규격화 조건을 유지하면서 분자를 최소화하는 쪽이 더 간편하다. 변분의 표현으로는 $$\delta \left( \int \psi^\ast H \psi ~dx - \epsilon \int \psi^\ast \psi ~dx \right) = 0,$$ $$\int \delta\psi^\ast (H-\epsilon) \psi ~dx + \int \psi^\ast (H-\epsilon) \delta \psi ~dx = 0.$$ 여기에서 $H \psi = \epsilon \psi$가 나오므로 라그랑주 곱수의 물리적 의미는 에너지 고유값이다. =====기저함수로 전개===== 라그랑주 곱수를 사용해 범함수 $$G[\psi] = \int \psi^\ast H \psi ~dx - \epsilon \int \psi^\ast \psi ~dx$$ 를 적고 직교정규화된 $N$개의 기저함수로 파동함수를 표현하자. 즉 기저함수 $b_1, b_2, \ldots, b_N$은 $$\int b_i^\ast b_j ~dx = \delta_{ij}$$ 을 만족하며 $\psi \approx \sum_i c_i b_i$로 표현한다. 그러면 범함수는 사실상 $c_1, c_2, \ldots, c_N$의 함수여서 $$G[c_1, \ldots, c_N] = \sum_{ij} c_i^\ast c_j H_{ij} - \epsilon \sum_{ij} c_i^\ast c_j \delta_{ij}$$ 이며 여기서 $$H_{ij} \equiv \int b_i^\ast H b_j ~dx$$ 로 정의된다. $H$가 에르미트 성질을 가지므로 $H_{ij} = H_{ji}^\ast$를 만족해야 한다. 이제 $\partial G / \partial c_i = 0$을 풀어야한다. 위에 적힌 $G[c_1, \ldots, c_N]$을 미분해보면 이는 $$\sum_j (H_{ij} - \epsilon \delta_{ij}) c_j = 0,$$ 행렬식으로 표현하면 $$(H - \epsilon I) \vec{c} = 0$$ 임을 의미한다. 즉 변분의 해 $c_1, c_2, \ldots, c_N$을 찾는 일은 해밀토니안 행렬의 고유값 문제 $H \vec{c} = \epsilon \vec{c}$로 귀착된다. =====예제===== 입자가 $0<x<L$의 상자 안에 갇혀있고 상자에는 $V(x) = V_0 \frac{x}{L} \left( \frac{x}{L}-1 \right)$의 퍼텐셜이 걸려있다. 해밀토니안은 다음처럼 주어지고 $$H = -\frac{-\hbar^2}{2m} \frac{d^2}{dx^2} + V_0 \frac{x}{L} \left( \frac{x}{L}-1 \right).$$ 경계조건은 $\psi(0) = \psi(L) = 0$이다. 기저함수를 평면파 $$b_n = \sqrt{\frac{2}{L}} \sin \frac{n\pi x}{L}$$ 로 고르면, 해밀토니안 행렬은 다음처럼 원소를 가진다. $$H_{nk} = \left\{ \begin{array}{ll} \frac{\hbar^2 \pi^2 n^2}{2mL^2}-V_0 & (n=k)\\ \frac{4V_0 nk [1+(-1)^{k+n}]}{\pi^2 (n-k)^2 (n+k)^2} & (n \neq k) \end{array} \right.$$ 이제 코드를 적어보자. <code:python> from __future__ import print_function,division from math import pi from numpy import array, empty from numpy.linalg import eigvalsh # symmetric or hermitian hbar = 1.0546e-34 # Planck constant / (2*pi) m, eV = 9.1094e-31, 1.6022e-19 # electron mass and 1eV = 1.6022e-19J L = 5.2918e-11 # box size N = 10 # number of basis functions V0 = 10**4*eV H = empty([N,N], float) for n in range(N): for k in range(N): n1,k1 = n+1,k+1 # the physical indices start from 1 if n==k: # diagonal elements H[n,k] = hbar**2*pi**2*n1**2/(2*m*L**2) \ - V0*(pi**2*n1**2+3)/(6*pi**2*n1**2) else: # off-diagonal elements H[n,k] = 4*V0*n1*k1*(1+(-1)**(k1+n1))/(pi**2*(n1-k1)**2*(n1+k1)**2) e = eigvalsh(H) print(e/eV) # output in units of eV </code> 위 코드의 변분법으로는 $N=10$개의 기저함수만을 사용해도 $E_0 \approx -2128.8767 eV$를 얻는다. 같은 상수값들을 가지고 [[:전산물리학:사격법]]을 사용할 경우 $E_0 \approx -2128.879 eV$를 얻는다. 분명 변분법의 결과가 더 크기는 하지만 그 차이는 매우 작다. ======함께 보기====== * [[수학:범함수]] * [[전산물리학:모멘트법]] * [[전산물리학:유한요소법]] ======참고문헌====== * http://www.fisica.uniud.it/~giannozz/Corsi/MQ/LectureNotes/mq.pdf * https://ocw.mit.edu/courses/16-920j-numerical-methods-for-partial-differential-equations-sma-5212-spring-2003/b08aec3e6f8cd87891d84d05b6d7a5eb_lecs13_14_notes.pdf * Karl E. Lonngren, Sava V. Savov, and Randy J. Jost, //Fundamentals of Electromagnetics with MATLAB//, (SciTech Publishing, 2007). * Matthew N.O. Sadiku, //Computational Electromagnetics with MATLAB//, (CRC Press, Boca Raton, FL, 2019).