슈뢰딩거 방정식을 수치 해법으로 풀기 (1): 슈팅 방법으로 알아보는 고유값 문제, 무한 포텐셜 우물 문제 풀기

지난 세 번의 포스팅을 통해 미분 방정식을 수치 해법으로 푸는 과정을 알아봤습니다. 수치 해법으로 미분 방정식을 푸는것에 대해서 잘 모르신다면, 아래 포스팅을 먼저 보고 나서 이 포스팅을 읽는 것을 추천합니다. 

https://studyingrabbit.tistory.com/53

수치 미분 방정식 풀이법과 이를 이용한 고전 역학 문제 해결 : 오일러 방법(Euler method)을 이용하

https://studyingrabbit.tistory.com/54

수치 미분 방정식 풀이법을 이용하여 최종값 조건이 주어진 문제 풀기 : 슈팅 방법

https://studyingrabbit.tistory.com/55

룽게-쿠타 방법(Runge-Kutta method) 을 활용하여 보다 빠르고 정확하게 미분 방정식 문제를 수치 해법

가장 기본적인 수치 해법인 오일러 방법과 룽게-쿠타 방법을 소개하였고, 이를 이용하여 고전 역학의 미분 방정식 문제를 풀어 보았습니다. 경계 조건이 초기값으로 주어진 경우에 대해서는 별다른 추가 과정이 없이 오일러 방법이나 룽게-쿠타 방법을 통해서 임의의 시간 $t$에서의 해를 쉽게 구할 수 있는 반면, 경계 조건이 최종값으로 주어진 경우에는 슈팅 방법 (shooting method)와 같이, 주어지지 않은 초기값을 조금씩 바꿔 가면서 원하는 최종값을 갖는 해를 구할 수 있다는 것을 실제 예를 통해서 알아보았습니다. 

 

이번 포스팅에서 풀어 볼 미분 방정식의 양자 역학의 "뉴턴 방정식"에 해당하는 슈뢰딩거 방정식 입니다. 좀 더 정확히는 시간에 무관한 슈뢰딩거 방정식을 풀어 볼 텐데, 질량이 $m$인 입자가 1차원 공간에서 $V(x)$라는 포텐셜장에 있는 경우, 이 입자의 파동 방정식을 기술하는 슈뢰딩거 방정식은 

$$-\frac{1}{2} \psi''(x) + V(x) \psi(x) = E \psi(x)$$

이 됩니다. $x$에 대한 미분을 $'$을 통하여 표시하였습니다. 또한 편의를 위해 $m = 1, \hbar = 1$로 놓았습니다.

 

이번 포스팅에서는 입자가 포텐셜에 의해 구속된 상태 (bound state)만을 고려할 것인데, 이 경우에는 $V(x)$가 $x=0$ 근처에서 최소값을 갖고, $x$가 원점에서 멀어질 수록 크기가 증가하게 됩니다. 또한 $x \rightarrow \pm \infty$ 일때, $V(x) = \infty$ 가 됩니다. 이 경우에는 파동함수 $\psi(x)$ 는 주로 $x=0$ 주변에서만 $0$이 아닌 값을 갖게 되고, $x$값이 원점에서 멀어질 수록 파동함수는 $0$에 가까워 집니다. 따라서 계산의 편의상, $|x| \gt L$ 영역에서는 파동 함수의 값이 $0$이라고 가정하도록 하겠습니다. (이는 $V(x)$는 $|x| \gt L$ 영역에서는 무한대 값을 갖는다고 가정하는 것과 같습니다) 즉 일반적인 구속 상태의 포텐셜에 의해서 $\psi(L) = \psi(-L) = 0$이라는 경계 조건이 생긴다고 할 수 있습니다. 실제로는  $\psi(\infty) = \psi(-\infty) = 0$ 가 맞는데, 수치 계산에서는 무한대 $\infty$를 다룰 수 없기 때문에, 편의상 적당히 큰 숫자인 $L$를 무한대 대신으로 사용하는 것 입니다. $L$이 크면 클수록 해석적인 해와 수치 해법으로 구한 해가 유사해 질 것 입니다. 

 

하나 더 고려해야 할 것은, 슈뢰딩거 방정식은 $\psi(x)$에 대한 선형 방정식으로 $\psi(x)$가 위 방정식의 해가 된다면, 이 함수를 실수배 한 $c\psi(x)$ 역시 방정식의 해가 됩니다. 따라서 양자 역학적 응용을 위해서는 $\psi(x)$를 정규화(Normalization) 시키는 것이 필요하지만, 정규화의 과정은 다소 당연하게 할 수 있기도 하고, 이는 우리의 주된 관심사가 아니기 때문에 정규화에 대한 고려는 하지 않도록 하겠습니다. 

 

그리하여, 우리가 풀어야할 문제의 최종 형태는 아래와 같이 요약 됩니다. 

$$\psi''(x) = 2 \Big(V(x) - E \Big) \psi(x)$$

$$x \in [x_a = -L, x_b = L], \psi(x_a=-L) = \psi(x_b=L) = 0$$

 

위 경계 조건은 이전 포스팅에서 정의했던 "최종값 문제" 입니다. 2계 미분 방정식이기 때문에 경계 조건은 2개가 있어야하는데 $x=x_a$에서의 함수값과 미분값이 모두 주어졌다면 룽게-쿠타 방법을 통해 한 번에 풀면 되지만, $x=x_a, x=x_b$에서의 함수값만 주어졌기 때문에, 슈팅 방법을 통해서 경계 조건을 만족 시키는 $x=x_a$에서의 미분값 $\psi'(x_a)$를 찾는 문제로 생각할 수 있습니다. 

 

그러나, 위와 같은 해석은 옳지 않습니다! 슈뢰딩거 방정식이 뉴턴 방정식처럼, 경계 조건이 주어진 경우 단 하나의 해 $\psi(x)$를 구하는 방정식이라면 위와 같은 해석이 옳지만, 슈뢰딩거 방정식은 단순한 미분 방정식이 아니라 "경계값" 문제 입니다. 위 슈뢰딩거 방정식에서는 문제를 풀기 전에 $E$의 구체적인 값이 숫자로 주어지는 것이 아니라, 경계 조건을 만족하는 $\psi(x)$에 대응되는 $E$값을 찾아 주어야 합니다. 문제의 고유값 중 하나가 먼저 제시된 경우라면, 바로 위 문단에서 쓴 것과 같이 슈팅 방법을 이용하면 되지만, 보통의 슈뢰딩거 방정식에서는 $E$의 값이 미리 주어지지 않습니다.

 

또한, 더 앞에서 이미 언급한 바와 같이 위 방정식은 $\psi(x)$에 대한 선형 방정식이기 때문에 $\psi(x)$가 해가 된다면 임의의 실수배 $c\psi(x)$ 역시 해가 됩니다. 따라서 일반적인 슈팅 방법에서 가정하는 $x=x_a$에서의 미분값 $\psi'(x_a)$는 의미가 없습니다. $x=x_a$에서의 미분값을 $1$이라고 가정했을 때 $\psi(x)_1$가 해로 얻어졌다면, $x=x_a$에서의 미분값을 $2$로 가정하는 경우, 해를 찾을 수 있고 이 때의 해는 $2 \cdot \psi(x)_1$이 됩니다. $E$가 주어진 문제의 고유값이라면 초기값의 미분값을 어떤 값으로 가정하는가와는 상관 없이 해를 구할 수 있고, $E$가 고유값이 아니라면 초기값의 미분값과는 상관 없이 해를 구할 수 없습니다. 

 

이와 같은 슈뢰딩거 방정식의 성질에 따라서, 슈뢰딩거 방정식을 풀 때에 사용되는 슈팅 방법의 변수는 초기 미분값 $\psi(x_a)$가 아니라 에너지값 $E$가 되어야 합니다. 즉, $E$값을 조금씩 변화시켜 가면서 얻어지는 수치해가 경계조건을 만족하는지를 확인하는 것 입니다. 앞선 문단에서 언급한 바와 같이 슈뢰딩거 방정식의 선형성에 따라서 $\psi(x_a)$는 어떤 값으로 가정하든 상관이 없습니다. 물론 $\psi(x_a)$값이 $0$이 되어서는 안되는 경우가 있습니다. 

 

(또 다시) 그리하여, 우리가 풀어야할 슈팅 방법의 최종 과정은 아래와 같습니다. 이전 포스팅의 고전 역학 문제를 풀 때의 과정과 비교해 보면 좋습니다. 

 

(1) $\psi(x_a) = 0, \psi'(x_a) = 1$을 초기 조건으로 하고, 특정한 $E$에 대한 슈뢰딩거 방정식을 푼다. 룽게-쿠타 방정식을 이용한다.

(2) 위에서 얻은 해의 $\psi(x_b)$에서의 값이 $0$인지 확인한다. 

(3) $\psi(x_b) = 0$이라면 (1)에서 가정한 $E$가 실제로 시스템의 고유값이 되고, (1)에서 얻은 $\psi(x)$가 해가 된다. 그렇지 않다면, (1)의 과정으로 돌아가서 새로운 $E$에 대해서 과정을 반복한다

 

입니다. 다시 설명하자면, 일반적인 미분 방정식에서는 고정 된 $\psi(x_a), E$ 값에 대해서 $\psi'(x_a)$를 바꾸어 가면서 $x=x_b$에서의 경계조건을 만족하는지에 대해서 다루었다면, 슈뢰딩거 방정식에서는 고정된 $\psi(x_a), \psi'(x_a)$ 값에 대해서 $E$를 바꾸어 가면서 $x=x_b$에서의 경계조건을 만족하는지를 탐색하게 됩니다. 슈뢰딩거 방정식이 뉴턴 방정식과 같이 초기값이 정해지면 모든 결과가 정해지는 방정식이 아니라, 문제에 꼭 맞는 고유값을 구해야하는 고유값 문제이기 때문입니다. 이 설명이 이해가 안된다면, 미분 방정식 혹은 슈팅 방법을 통해서 문제를 푸는 것에 대해서 이해가 부족한 것이라기 보다는 슈뢰딩거 방정식 그 자체에 대한 물리적 혹은 수학적인 이해가 부족한 것이라고 할 수 있습니다. 이 경우에는 이 포스팅을 계속해서 읽는것 보다는 슈뢰딩거 방정식 그 자체에 대해서 공부를 하는 편이 더 좋습니다. 

 

무한 포텐셜 우물 문제 풀기

 

서론이 조금은 길었습니다. 그럼 이제 실제로 문제를 풀어 보도록 하겠습니다. 역시나 처음에 나오는 문제는 무한 포텐셜 우물에 갖혀 있는 입자에 대한 슈뢰딩거 방정식 풀이 입니다. 포텐셜 우물의 너비를 $2\pi$라고 한다면, 위치 에너지는 

$$\begin{align}
    V(x) = \Bigg\{
        \begin{array}{ll}
        0, & 0  \leq x \leq 2\pi,\\
        \infty, & \text{otherwise}
        \end{array}
\end{align}$$

입니다. 위 문단에서 우리가 가정한 일반적인 조건에서는 $x \in [-L, L]$로 $x$의 공간을 정의한다고 했지만, 이 문제의 경우에는 그냥 $x \in [0, 2\pi]$를 사용하는 것이 편리합니다. 그래서 이번 문제는 간단하게, 

$$\psi''(x) = -2E \psi(x)$$

$$x \in [0, 2\pi], \psi(0) = \psi(2\pi) = 0, \psi'(0) = 1$$

가 됩니다. 경계 조건의 마지막 $\psi'(0) = 1$ 는 정확히는 경계조건은 아니고, 정규화 되지 않은 파동함수 중에서 $\psi'(0) = 1$인 파동함수를 찾겠다는 뜻 입니다. 위 문제를 서로 다른 $E$값에 따라서 룽게-쿠타 방법을 통해서 $\psi(x)$를 계산해 보도록 하겠습니다. 

 

우선 테스트로 $E=0.1$에 대한 $\psi(x)$를 구하면 위와 같습니다. 편의상 이 파동함수를 $\psi_{E=0.1}(x)$라고 표시하겠습니다. 그래프에서 볼 수 있는 것과 같이 $\psi_{E=0.1}(2\pi)$값은 $0$보다 한참이나 큽니다. 따라서 $E = 0.1$은 고유값이 아닙니다. 당연히 $\psi_{E=0.1}(x)$ 역시 슈뢰딩거 방정식의 해가 아닙니다. 

$E = 0.15$에 대한 $\psi_{E=0.15}(x)$를 붉은색으로 그렸습니다. 이번에는 $\psi_{E=0.15}(2\pi) \lt 0$입니다. 이를 통해서 $0.1 \lt E \lt 0.15$ 인 $E$에서 고유값이 있다는 것을 알 수 있습니다. Trial and Error를 통해서 $E$ 값을 구해보면, 

위와 같이 $E = 0.125 = \frac{1}{8}$에서 정확히 $\psi_{E=0.125}(2\pi) = 0$가 됨을 얻을 수 있습니다. 즉 이 문제의 고유값은 $E_0 = 0.125$, 이에 대응되는 고유 함수는 $\psi_0(x) = \psi_{E=0.125}(x)$를 구했습니다. 여기서 $E$와 $\psi(x)$의 아래첨자 $_0$은 그라운드 상태를 나타냅니다. 

 

무한 포텐셜 우물의 구속된 상태는 무한개 존재합니다. 따라서 첫 번째 여기 상태를 구해보겠습니다. 같은 방식으로 $E$를 키워가면서 $\psi_{E}(2\pi) = 0$가 되는 $E$를 찾으면 됩니다. 

$E = 0.4, 0.5, 0.6$ 에 대한 $\psi_{E=0.4}(x), \psi_{E=0.5}(x), \psi_{E=0.6}(x)$을 차레로 구한 것 입니다. $E = 0.5$에서 정확히 슈팅에 성공(즉, $\psi_{E=0.5}(2\pi) = 0$) 하였습니다. 따라서 첫 번째 여기 상태의 에너지는 $E_1 = 0.5$가 됩니다. 

 

이와 같은 방식으로 계속 해서 $E$를 증가시키면서 슈팅에 성공하는 $E$를 찾으면 됩니다. 매번 위와 같은 그래프를 그리는 것은 다소 귀찮으니, $E$ vs $\psi_{E}(2\pi)$ 를 그리고 이 그래프와 $x$과의 교점, 즉 $\psi_{E}(2\pi) = 0$를 만족하는 $E$를 찾으면 됩니다. 

$0 \lt E \lt 5$ 구간을 에너지 간격 $0.01$마다 $\psi_E(2\pi)$를 계산한 그래프 입니다. 이 그래프와 $x$축 ($E$축)과의 교점에 해당하는 $E$값이 고유값에 해당하고, 그 때의 파동 함수가 고유 함수에 해당합니다. 위 그래프의 $x$축 근처를 조금 확대해서 보면, 

와 같고, 여기서 붉은색으로 동그라미 친 값이 $\psi_E(2\pi) = 0$ 을 만족시키는 $E$값이 됩니다. 그 값은

$$E = 0.12, 0.5, 1.13, 2.0, 3.13, 4.51$$

입니다. 

 

이 문제는 당연히 해석적인 풀이가 가능한데, 양자 역학 시간에 배웠듯 고유값은 $\frac{n^2}{8}$로 주어집니다. 따라서 정확한 고유값은 

$$E = 0.125, 0.5, 1.125, 2.0, 3.125, 4.5$$

입니다. 이 값과 위에서 수치적구한 고유값이 거이 일치 함을 확인할 수 있습니다. 위 계산 과정에서 에너지의 간격이 $0.01$이었기 때문에, 차이가 생기게 된 것 입니다. 

 

이분법(bisection method)를 이용하여 해를 정확하게 구하기

 

여러번 언급한 바와 같이, 임의의 초기값 문제에 대한 미분 방정식을 룽게-쿠타 방법과 같이 수치적인 해법으로 풀 수 있는 방법을 마련해 놓았다면, 최종적으로 풀어야 하는 문제는 $\psi_E(2\pi) = 0$ 를 만족하는 $E$를 찾는 방정식 문제로 바뀌게 됩니다. 위에서는 $E$ vs $\psi_E(2\pi)$ 의 개형을 보여주기 위해서 $E$를 작은값 부터 큰값으로 특정한 간격을 두고 모든 $\psi_E(2\pi)$를 구했습니다. 이 방법은 직관적이기는 하지만, 부정확하고 비효율적인 방법입니다. 앞에서 본 것과 같이 $E$이 간격이 크다면, 해가 부정확해 지고, $E$의 간격이 좁다면 해가 정확해 지지만 계산량이 커진다는 단점이 있습니다. 

 

따라서, 수치 방정식을 푸는 일반적인 방법을 도입하여 해 $E$를 찾는 것이 효율적입니다. 수치 방정식 해법 중에서 직관적이면서도 간단한 방법으로는 이분법(bisection method)가 있습니다. 

이 그림은 위키백과의 사진을 그대로 따왔습니다

이분법은 위 그림 하나로 설명이 될 것 같습니다. 연속함수 $F(x)$의 $F(x) = 0$의 해를 찾고 있습니다. $F(x)$는 해 주변에서 위와 같은 개형을 갖게 될 것 입니다. 만일 $F(a_1) \gt 0, F(b_1) \lt 0$ 이라면, $F(x)$의 연속성에 의해서 반드시 구간 $[a_1, b_1]$에 적어도 하나 이상의 해가 존재하게 됩니다. 이 구간을 점점 좁혀가면서 해를 정확하게 찾는것이 이분벙 방법의 핵심인데요, $b_2 = \frac{a_1 + b_1}{2}$에 대해서 함수값 $F(b_2)$를 구합니다. 만일 이 값이 음수라면, 해는 $[a_1, b_2]$ 구간에 있는 것이고, 이 값이 양수라면, 해는 $[b_2, b_1$ 구간에 있는 것이 됩니다. 이렇게 한 번 과정을 거치게 되면 해가 존재하는 구간의 크기는 $1/2$로 줄어들게 됩니다. 이 과정을 $n$ 번 반복한다면, 해가 존재해야 하는 구간의 크기는 $\frac{1}{2^n}$ 으로 줄어들게 되고, $n$이 큰 수가 되면 해가 존재하는 구간이 거의 0에 가깝기 때문에 그 구간의 평균값이 해의 값이라고 근사할 수 있습니다. 이 때 근사의 에러는 $\frac{1}{2^{n+1}}$ 크기 정도가 되니, $n$을 10으로 한다면 대략 $0.001$의 오차로 해를 구할 수 있습니다. 위에서 모든 $E$에 대해서 계산한 것에 비해서 정확도가 훨씬 증가하고 계산량을 훨씬 감소합니다. 

 

그렇다면, 이분법을 이용하여 위 $E$ vs $\psi_E(2\pi)$ 그래프를 통해서 얻은 6개의 해를 정확도 높게 구해 보도록 하겠습니다. 첫번째 해의 값은 0.12로 얻었는데, 정확한 값은 $[0.1, 0.15]$ 범위에는 존재할 것 입니다. 이 값을 이분법의 초기 영역으로 하여 정확한 $E$값을 구하면, $E = 0.12502$ 를 얻게 됩니다. 참값 $0.125$와 비교하면 $0.00002$의 오차가 발생합니다. 이 오차는 룽게-쿠타 방법을 푸는데서 생긴 오차입니다. 같은 방법으로 여기 상태의 에너지를 차례로 구할 수 있습니다. 그 결과를 소숫점 넷째짜리 까지 구하면,

$$E = 0.1250, 0.5001, 1.1252, 2.0004, 3.1256, 4.5009$$

입니다. 해석적인 방법으로 얻은 참값과 비교하면 에러의 크기는 모두 $0.0001$ 수준입니다. 즉, 수치적으로 얻은 해는 소숫점 아래 세 번째 자리 까지 정확합니다. 

 

고유 함수 $\psi_E(x)$ 도 구해보자

 

위에서 구한 고유값에 대한 고유 함수를 구해서 그리면, 결과는 아래 그래프와 같습니다.

위 그래프를 보면 고유 함수 마다 크기(최대값)이 각각 다른데, 고유 함수를 정규화 하지 않았기 때문입니다. $x=0$에서의 기울기를 $1$로 두고 계산을 했기 때문에 $x=0$에서의 개형이 거의 동일 합니다. 고유값 $E$에 의해서 함수의 곡률이 정해지기 때문에 $E$값이 큰 함수일수록 그래프가 급격하게 바뀌게 되고, 잘 아는대로 고유 값이 커질수록 노드의 갯수에 하나씩 커지게 됩니다. 이전에 포스팅한 바 있는 적분 방법론을 통해서 위 고유 함수를 정규화 할 수 있습니다. 이것은 독자분께 숙제로 남기겠습니다. 

 

여기서 끝나긴 아쉽다.

 

프로그램 코드를 만들어 놨으니, 임의의 $V(x)$에 대한 슈뢰딩거 방정식의 고유값과 고유함수를 구할 수 있습니다. 다음 포스팅에서는 조금 더 복잡한 형태의 포텐셜에 대해서 문제를 풀어 보도록 하겠습니다.