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차원에서 최대가 된 후 다시 감소한다.

타원을 표현하는 방정식은 

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

처럼 쓸 수 있다. 매개변수로 표현하면

$$ 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]  \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}$$

임을 확인할 수 있다.

지구 표면의 온도는 1년을 단위로 거의 주기적으로 변한다. 그럼 땅속의 온도는 시간과 깊이에 따라 어떻게 변할까? 지표면이 태양으로부터 받은 열은 일부는 반사되고 일부는 땅속으로 전달된다. 땅속에서 온도의 변화는 열방정식에 의해서 표현할 수 있다. 땅속의 온도분포 $u$가 지표면에서 깊이 $x$와 시간 $t$에만 의존한다면 $u(x,t)$가 만족하는 열방정식은(거리척도를 적당하게 잡아 계수를 단순화시킨다)

$$ u_t(x,t) = u_{xx} (x,t)$$

로 주어진다. $u(x,t)$가 $t$에 대해서 주기가 $1$년인 주기함수이므로 Fourier 급수를 이용해서 방정식을 풀도록 하자. 

\begin{gather} u(x,t) = \sum _{n=-\infty}^{\infty} C_n (x) e^{i 2\pi n t} \\ C_n(x)= \int_0^1 u(x,t) e^{-i 2\pi n t} dt \end{gather}

Fourier 계수 $C_n(x)$를 두 번 미분하면

\begin{align} \frac{d^2 C_n(x)}{dx^2} & = \int_0^1 u_{xx}(x,t) e^{-i 2 \pi n t} dt \\ &= \int_0^1  \frac{ \partial u (x,t)}{\partial t} e^{-i 2 \pi n t} dt \\ &= i 2\pi n \int_0^1 u(x,t) e^{-i 2 \pi n t} dt \\ &= i 2 \pi n C_n(x)\end{align}

이므로 $C_n$은 다음의 방정식을 만족해야 한다.

$$ \frac{d^2 C_n (x)}{dx^2}  = i 2\pi n C_n (x) $$

깊이 $x$가 증가할 때 온도가 발산하지 않는 조건을 고려하면 이 방정식의 해는 

$$ C_n (x) = \left\{ \begin{matrix}  A_n e^{- \sqrt{\pi n} (1+i)x}~~~n \ge 0 \\A_n e^{ -\sqrt{\pi |n|}(1-i)x}~~~n<0 \end{matrix} \right.$$

로 쓸 수 있음을 쉽게 알 수 있다. 따라서 열방정식의 해 $u(x,t)$는

$$ u(x,t) = \sum _{n=- \infty}^{\infty} A_n e^{- \sqrt{\pi  |n| }x} e^{ i (2\pi n t - \text{sign}(n) \sqrt{\pi |n| } x) }   $$

처럼 쓰인다. 온도는 깊이에 따라 감쇄를 하여 계절에 따른 온도변화가 점점 작아진다. 그리고 깊이에 따른 위상이 추가되므로 지표면에서의 온도변화와 다른 양상을 가지게 된다. 이를 구체적으로 보기 위해서 계절에 따른 지표면에서 온도 $u(x=0,t)$을 간단히 시간 $t$에 대한 사인함수로 근사하자. 이 경우 평균온도가 0인데, 평균온도가 0이 아니 경우는 여기서 구한 해에  평균온도만큼을 더해주면 된다.

$$u(0, t) = \sin (2\pi t)$$

지표면에서 Fourier 계수는 

\begin{align} C_n(0) =A_n  &= \int_0^1 u(0,t)   e^{-i 2\pi n t} dt \\ &=\left\{  \begin{matrix} \pm \frac{1}{2i} ~~~n=\pm 1 \\ 0~~\text{otherwise}  \end{matrix} \right. \end{align}

따라서 해는 

\begin{align} u(x,t) &= \frac{1}{2i} e^{-\sqrt{\pi}(1+i)x} e^{i2 \pi t} -\frac{1}{2i} e^{-\sqrt{\pi} (1-i)x }e^{-i 2\pi t} \\ &= e^{-\sqrt{\pi} x} \sin (2 \pi t - \sqrt{\pi}x) \end{align}

해를 보면 온도는 깊이에 따라 감쇄를 하여 온도변화가 점점 사라지고, 시간에 대해서는 깊이에 따른 위상변화가 생긴다. 특히 깊이 $x= \sqrt{\pi}$에서는 위상이 $\pi$ 만큼 변해서 지표면에서의 시간에 대한 온도변화와 완전히 반대로 행동한다. 즉, 겨울에는 따뜻하고 여름에는 시원해진다. 물론 이 깊이는 땅의 열확산계수($\kappa$)에 따라 달라진다. 열확산계수를 고려하려면 $x \to x/\sqrt{\kappa}$을 사용하면 된다. 땅의 열확산계수가 $\kappa \sim 0.1\times 10^{-6} \rm{m^2 / s=3.15 m^2/yr}$ 정도이므로 $x=\sqrt{\kappa \pi} \sim 3 \rm m$이다. 즉, 땅 속 깊이 $x\sim 3\rm m$ 정도이면 온도변화는 지표면의 $e^{-\pi}=0.043$배 정도로 줄어들고 겨울이 여름보다 상대적으로 더 따뜻하게 된다.

이는 물리적으로 쉽게 이해를 할 수 있는 현상으로 깊이 들어갈수록 온도차가 작아져서 열전달이 느려지므로 상대적으로 빨리 변하는 표면에서의 온도변화에 맞추지 못하여 변화가 지연되어 나타나는 것으로 볼 수 있다.