Geometry & Recognition

초기의 지구는 매우 뜨거워 모든 암석들이 녹아있는 상태였을 것이다. 시간이 흐르면서 지구는 표면에서 차가운 우주(거의 0K이다)로 열이 배출되어 점차 식게 된다. 처음에서 표면과 내부의 온도가 같았지만 열의 배출과정에서 표면의 온도는 내부의 온도보다 낮게 되어 지구 내부에서 온도의 gradient가 만들어진다. 실험적으로 지표면 근처에서 이 기울기를 측정하면 지구 나이를 추정할 수 있다.

지구 내부에서 온도의 분포는 열방정식에 의해서 구할 수 있다. 우선 지구가 매우 크므로 평평한 판(그리고 매우 두꺼운)으로 근사를 하자. 그러면 열이 식는 과정은 1차원 열방정식으로 설명할 수 있다. 지구 내부 깊이 $x>0$인 지점에서 온도 분포 $u(x, t)$는 다음의 열방정식에 의해서 결정된다.

\begin{align}  u_t(x, t) = \kappa u_{xx} (x,t), \quad u(x, t=0)=f(x), \qquad x>0\end{align}

여기서 $f(x)$는 초기 지구 내부 온도로 일정한 값으로 생각할 수 있다. 지표면 밖(우주)은 무한히 큰 차가운 열원으로 생각할 수 있으므로 온도가 $0\rm K$으로 고정되어 있다고 볼 수 있다. 즉, $u(x,t)$은 모든 시간에서 다음의 경계조건을 만족해야 한다.

$$ u(x=0, t)=0$$

이 $x\ge 0$에서 열방정식은 전체 실수구간에서의 열방정식을 푸는 문제로 변환하면 보다 쉽게 해결할 수 있다. 단, $x=0$에서 경곗값을 유지해야 하므로 $t=0$일 때 초기값은 원 문제의 기함수로 확장을 쓰면 된다. 물리적으로 $t=0$ 일 때 $x=0$을 기준으로 $T$와 $-T$(물론 음의 절대온도는 없다!)를 가지는 영역을 열접촉을 시키고 난 후의 온도 분포를 찾는 문제로 환원시킨 것과 같다. 따라서 다음의 열방정식의 해 $v(x,t)$를 우선 고려하자.

\begin{align} v_t (x,t) = \kappa v_{xx}(x,t),     \quad v(x,0)= {f}_o(x)=\left\{ \begin{matrix} f(x) & x>0 \\ -f(-x) & x <0\end{matrix}\right. \end{align}

이 방정식의 해는 다음처럼 써진다.(Fourier 변환과 convolution을 사용하면 쉽게 구할 수 있다). $\hat{v}(s,t)$와 $\hat{f_o}(s)$를 각각 $v(x,t)$와 ${f_o}(x)$의 공간에 대한 Fourier 변환이라면 열방정식은 

$$\hat{v}_t (s,t)= -\kappa s^2 \hat{v}(s,t),\quad \hat{v}(s, 0) = \hat{f_o}(s)$$

$t$에 대해 적분하면 다음처럼 두 함수의 곱으로 표현된다.

\begin{align} \hat{v}(s,t) =  e^{- \kappa s^2 t /2} \hat{f_o}(s) = \hat{K}(s,t) \hat{f_o}(s)\end{align}

따라서 $v(x,t)$는 heat kernel $K(x,t)$와 ${f}_o(s)$의 convolution으로 주어지는데, $e^{-\kappa s^2 t}$의 역변환이 $K(x,t)=\frac{1}{\sqrt{4\pi\kappa t}} e^{- x^2/4\kappa t}$이므로 

\begin{align} v(x,t) &= (K * {f}_o)(x,t)\\ &=\frac{1}{\sqrt{ 4\pi\kappa t}} \int_{-\infty}^\infty e^{- (x-y)^2/4\kappa t } {f}_o(y) dy \\ &= \frac{1}{\sqrt{4\pi\kappa t}} \int_0^\infty \Big( e^{-(x-y)^2/4\kappa t} - e^{-(x+y)^2 / 4\kappa t} \Big) f(y)dy \end{align}

항상 $v(0,t) = 0$을 만족함을 쉽게 확인할 수 있다. 따라서 $x\ge 0$에서 $u(x,t) = v(x,t)$이다. 초기온도 분포가 일정한 ($f(x)=T$) 경우는 적분을 정리하면 error function으로 써짐을 알 수 있다.

$$u(x,t) = \frac{T}{\sqrt{\pi \kappa t}} \int_0^x e^{-s^2 /4\kappa t } ds = T\times  \text{Erf}(x/2\sqrt{\kappa t}) \quad (x\ge0)$$

현시점($t=t_0$)에서 $x=0$에서 온도 기울기 $\beta$를 구하면 

$$\beta= u_x (0, t_0) = \frac{T}{\sqrt{\pi \kappa t_0}}  \quad \rightarrow \quad t_0 = \frac{T^2}{\pi \kappa \beta^2} $$

지구 초기 온도로 $T = 3900 \rm K$,  지각의 평균 열확산도로 $\kappa = 1.2 \times 10^{-2} \rm cm^2/s$, 그리고 현시점의 지표면 근처에서 실험적으로 측정된 온도 기울기가 $\beta = 3.6\times 10^{-4} \rm K/cm$ 정도이므로

$$ t _0 =3.1131 \times 10^{15} \rm s =  9.7 \times 10^7 yr$$

이 추정에 의하면 지구나이가 1억 년이 안된다. 방사선 동위원소법에 의해서 지구 나이가 45억 년쯤 되는 것으로 나오는데 이와 큰 차이가 있다. 구형 지구를 고려하면 3차원 열방정식 문제가 되는데 이 방법을 쓰더라도 평평한 지구 모형을 쓴 결과와 큰 차이가 나지 않는다. 이는 지구가 생성된 이후 추가적인 열의 공급 없이 식는 과정만 고려한 모델인데 실제로는 지구내부에 있는 방사성 원소가 붕괴할 때 발생하는 열이 지속적으로 공급되고 있으므로 이에 대한 고려를 해야 좀 더 정확한 나이를 추정할 수 있을 것이다. 열방정식을 이용한 지구 나이 추정은 영국의 수리물리학자이자 공학자인 Kelvin 남작에 의해 시도되었다.

반지름이 $r$인 $n$차원 공의 부피 $V_n(r)$은 $r^n$에 비례하므로 

$$ V_n(r) = K_n r^n$$

처럼 쓸 수 있다. 계수 $K_n$는 반지름이 1인 공의 부피를 나타낸다. $K_n$을 구하기 위해서 다음의 $n$ 차원 gaussian 적분을 고려하자.

$$ I_n = \int_{-\infty}^\infty dx_1 \int_{-\infty}^\infty dx_2 \cdots \int_{-\infty}^\infty dx_n e^{-x_1^2 - x_2^2 - \cdots - x_n^2}$$

이 적분은 

$$I_n= \Big( \int_{-\infty}^\infty  e^{-x^2} dx\Big)^n = (\sqrt{\pi})^{n}$$

임을 쉽게 알 수 있다.  $n$ 공간에서 보면 적분인자가 중심에서 $(x_1, x_2,\cdots, x_n)$까지의 거리의 제곱 $x_1^2 + x_2^2 +\cdots + x_n^2$에만 의존하는 구대칭을 가진다. 따라서 $I_n$은 공간을 미소구각으로 나눈 후 $r$로 적분만 하면 된다. 이 경우 적분요소는 

$$ dx_1 dx_2 \cdots dx_n =( \text{구각면적}) \times dr= \frac{dV_n}{dr} dr = nK_n r^{n-1} dr$$ 

$$ I_n = \int_0^\infty e^{-r^2 }  nK _n r^{n-1} dr = \frac{n K_n }{2} \int_0^\infty e^{-t} t^{n/2-1} dt =K_n \frac{n}{2} \left(n/2-1\right)! = { K_ n } \times  (n/2)!$$

따라서 

$$ K_n = \frac{\pi^{n/2}}{(n/2)!}=\frac{\pi^{n/2}}{\Gamma(n/2+1)}$$

$(-1/2)!=\sqrt{\pi}$, $(1/2)!= \sqrt{\pi}/2$임을 고려하면

\begin{align} K_0 &= 1 \\ K_1 &= 2 \\  K_2 &= \pi=3.1415 \\ K_3 &= \frac{4\pi }{3}= 4.18879 \\ K_4 &= \frac{\pi^2}{2} = 4.9348\\ K_5 &= \frac{8\pi^2}{15} =  5.26379 \\ K_6 &= \frac{\pi^3}{6}  = 5.16771 \\ K_7 &= \frac{16\pi^3}{105}=4.72477 \\ K_ 8 &= \frac{\pi^4}{24}=4.05871 \\ K_{2n}& = \frac{\pi^n }{ n!},\qquad K_{2n+1} = \frac{2(2\pi)^n}{(2n+1)!!} \end{align} 

이 결과를 보면 차원이 증가하면 단위공의 부피값이 증가하다가 5차원에서 최대가 된 후 다시 감소한다.

$n$ 차원 단위공은 한변의 길이가 $2$인 $n$ 차원 cube에 담기게 되는데 큐브에서 공이 차지하는 부피는 

$$ \frac{K_n}{2^n}$$이다. 차원이 커지면 이 값은 빠르게 0으로 수렴한다. Stirling 근사를 쓰면

$$ n! \sim \sqrt{2\pi n} \left( \frac{n}{e}\right)^n , \qquad n \gg 1$$

이므로 

$$\frac{K_n }{ 2^n} \sim \frac{1}{\sqrt{2\pi n}} \left(  \frac{2n}{\pi e}\right)^{-n}$$

타원은 두 지점(초점)에서 거리의 합이 일정한 점들의 자취로 표준형의 타원 방정식은

$$  \frac{x^2}{a^2} + \frac{y^2}{b^2}=1$$

처럼 쓸 수 있다. 이 곡선이 감싸는 면적은 단순하게 $\text{area}= \pi ab$로 주어지지만 둘레길이를 정확히 구하기 위해서는 적분이 필요하다. 타원방정식을 매개변수 $\theta$를 써서 표현하면

$$ x= a \cos \theta, \quad y=b \sin \theta, \quad 0\le \theta \le 2\pi$$

이므로 둘레길이는

$$L = \int _0^{2\pi} \sqrt{ \dot{x}^2 + \dot{y}^2 } d \theta  = \int_0^{2\pi} \sqrt{ a^2 \sin ^2  \theta + b^2 \cos ^2 \theta}d \theta$$

로 써진다.

$b\ge a$인 경우만 고려하여 적분인자를

$$ \sqrt{b^2 - (b^2 - a^2)\sin^2 \theta} d \theta = b \sqrt{1-e^2 \sin^2  \theta}d \theta$$

처럼 쓰자. 여기서 $e= \sqrt{b^2-a^2}/{b}$는 이심률(eccentricity)로 타원이 눌려진 정도를 나타낸다. 타원은 장축 반지름($b$)과 모양을 결정하는 이심률 ($e$)을 써서 표현할 수 있다. 대칭성에 의해서 둘레길이는 1 사분면에 걸친 길이의 4배를 하면 되므로

$$ L= 4 b   \int _0^{\pi/2} \sqrt{1-e^2 \sin ^2  \theta }d \theta $$

적분 부분은 장축 반지름이 1이고 이심률이 $e$ 타원 둘레길이의 $1/4$를 나타낸다. 그런데 이 적분은 닫힌 꼴을 구할 수 없고 무한수열의 합으로 표현해야 한다(참고: complete elliptic integral of the second kind). $x= e^2 \sin ^2 \theta$로 놓고 적분인자를 이항전개하면

\begin{align} (1-x)^{1/2} &= \sum_{n=0}^\infty \left( \begin{matrix}1/2\\n\end{matrix}\right) (-x)^n \\ &= 1-\frac{1}{2}x-\frac{1}{8}x^2 - \frac{1}{16}x^3-\frac{5}{128}x^4-\frac{7}{256}x^5-\cdots \\ &= 1-\frac{1}{2}x -\sum_{n=2}^\infty \frac{(2n-3)!!}{2^n n!}x^n \end{align}

그리고

$$ \int_0^{\pi/2} \sin^{2n} \theta  d \theta = \frac{\pi}{2}\frac{1}{2^{2n}} \left(\begin{matrix}  2n \\ n \end{matrix}\right)$$

임을 이용하면 다음의 결과를 얻는다:

\begin{align} \int_0^{\pi/2} \sqrt{1-e^2 \sin^2  \theta }  d \theta  &= \frac{\pi}{2} \left( 1- \frac{1}{2\cdot 2} e^2 -\sum_{n=2}^\infty \frac{(2n-1)!!(2n-3)!!}{ ((2n)!!)^2}  e^{2n}   \right)\\  &= \frac{\pi}{2}\left[1-  \Big(\frac{1}{2}\Big)^2 e^2 -\Big(\frac{1\cdot 3}{2\cdot 4}\Big)^2 \frac{e^4}{3} -\Big(\frac{1\cdot 3\cdot 5}{2\cdot 4 \cdot 6}\Big)^2 \frac{e^6}{5}-\cdots\right]    \end{align}

물론 이심률이 $e=0$인 원의 경우에는 $L=4b\times \frac{\pi}{2}=2\pi b$임을 확인할 수 있고, 이심률이 $e=1$인 경우(완전히 눌린 타원으로 끝이 연결된 두 선분)는  적분값이 1 이므로 둘레길이가 $L=4b$임도 확인할 수 있다.

참고:  $\sin \theta = (e^{i\theta}-e^{-i \theta})/2i$임과 이항정리를 쓰면

$$\sin^{2n} \theta = \frac{1}{(2i)^{2n}} \sum _{k=0}^{2n} \left( \begin{matrix}2n\\ k\end{matrix}\right)e^{-i 2(n-k) \theta }(-1)^k $$

이므로 적분을 하면

$$ \int_0^{\pi/2} \sin^{2n}\theta d\theta = \frac{(-1)^n}{2^{2n}} \left[ \sum _{k=0, k\ne n}^{2n} \left( \begin{matrix} 2n \\ k \end{matrix}\right)   \frac{(-1)^n - (-1)^k }{2i(n-k)} +\left(\begin{matrix}2n \\n\end{matrix}\right) \frac{\pi}{2}(-1)^n\right]$$

인데, $\left( \begin{matrix}2n\\ k\end{matrix}\right)=\left( \begin{matrix}2n\\ 2n-k\end{matrix}\right)$임을 고려하면 $\sum$ 내부의 항들은 서로 상쇄되어 위의 결과를 얻는다.  

그렇지만 위의 무한수열은 이심률이 작은 경우를 제외하고는 수렴이 빠르게 이루어지지 않아 실질적인 계산에 쓰기에는 문제가 있다. 좀 더 빠르게 수렴하는 타원 둘레길이 표현을 찾아보자. 이심률을 쓰는 식은 장축과 단축에 대해서 비대칭적이다. 좀 더 대칭적인 표현을 구하기 위해서 다음과 같이 두축의 차이와 두축의 합의 비를 이용하자. 

$$a = \frac{a+b}{2} (1- \sqrt{h}), ~b = \frac{a+b}{2} (1 + \sqrt{h})$$

로 정의하면

$$ h = \Big( \frac{a-b }{ a+b } \Big)^2$$

임을 알 수 있고, $a,b$에 대해서 대칭적이다.  그리고,

\begin{align} a^2  \sin ^2 \theta + b^2 \cos ^2 \theta &= \frac{a^2+b^2}{2} - \frac{a^2 - b^2 }{2} \cos(2\theta) \\ &= \frac{(a+b)^2 }{4} (1 + h +2\sqrt{h} \cos 2   \theta) \\ &=\frac{(a+b)^2 }{4}   (1+z)(1+\bar{z}) \end{align}

여기서 $z= \sqrt{h} e^{i 2\theta}$, $\bar{z}= \sqrt{h}e^{ - i 2\theta}$로 정의했다. 이항정리를 이용해서 적분인자를 전개하면 

\begin{align} L &= 2(a+b) \int_0^{\pi/2} \sqrt{ (1+z)(1+\bar{z})} d \theta    \\ &= 2(a+b) \sum_{m=0}^\infty \sum _{n=0}^\infty \left(\begin{matrix}1/2 \\ m \end{matrix} \right) \left(\begin{matrix}1/2\\n\end{matrix} \right)  h^{(m+n)/2} \int_0 ^{\pi/2}  e^{i 2\theta (m-n)}d \theta \end{align}

을 얻는다. 그리고

$$ \int _0^{\pi/2} e^{i 2\theta (m-n)} d \theta = \frac{\pi}{2} \delta_{mn}$$

이므로 다음과 같은 타원 둘레길이에 대한 식을 얻을 수 있다.

\begin{align} L& = \pi (a+b) \sum_{n=0}^\infty  \left( \begin{matrix} 1/2 \\ n\end{matrix} \right)^2 h^n \\ &= \pi (a+b) \left[ 1 + \frac{1}{4}h+ \frac{1}{64}h^2 + \frac{1}{256} h^3 + \dots \right] , \qquad h = \left( \frac{a-b}{a+b}\right)^2 \end{align}

여기서 $\left(\begin{matrix}1/2\\n\end{matrix}\right) = (1/2)(1/2-1)\cdots(1/2-n+1)/n!$로

$$\left(  \begin{matrix} 1/2 \\n \end{matrix}\right) = \left(\begin{matrix} 2n\\n\end{matrix} \right) \frac{1}{(1-2n)(-4n)^n} $$

이다. 완전히 납작한 타원($h=1$ or $e=1$)인 경우 $L=4b$이므로 

$$ \sum_{n=0}^\infty  \left( \begin{matrix} 1/2 \\ n\end{matrix} \right)^2=\frac{4}{\pi}$$

임을 확인할 수 있다.