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

 
 
 
 
 
 
 
 
 
 
 
728x90

'Mathematics' 카테고리의 다른 글

열방정식의 Green function  (0) 2024.02.12
지구의 나이는?  (0) 2024.02.11
타원의 둘레길이  (1) 2024.02.06
왜 땅속은 여름에 시원하고 겨울에는 따뜻한가?  (1) 2024.02.03
Rejection Sampling  (0) 2023.02.20
Posted by helloktk
,

타원을 표현하는 방정식은 

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

임을 확인할 수 있다.

 
 
 
 
 
 
 
 
 
 
728x90
Posted by helloktk
,

지구 표면의 온도는 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$배 정도로 줄어들고 겨울이 여름보다 상대적으로 더 따뜻하게 된다.

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

 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
728x90

'Mathematics' 카테고리의 다른 글

n 차원 공의 부피  (2) 2024.02.07
타원의 둘레길이  (1) 2024.02.06
Rejection Sampling  (0) 2023.02.20
Fourier Series를 이용한 Isoperimetric Inequality 증명  (0) 2023.02.01
Brachistochrone inside the Earth  (0) 2023.01.25
Posted by helloktk
,

Rejection Sampling

Mathematics 2023. 2. 20. 22:08

절대온도 $T$일 때 질량 $m$인 이상기체의 속력은 Maxwell-Boltzmann 분포를 따른다.

$$P(v)dv = 4\pi \left(\frac{m}{2\pi kT}\right)^{3/2}v^2 e^{-mv^2/2kT}dv$$

속력 $v$는 0부터 $\infty$까지의 값을 갖는다. 그럼 이 속력분포를 가지는 $N$개의 속력 샘플을 만들고 싶은 경우 어떻게 할까? 우선 무한대 범위까지 모두 처리할 수 없으므로 적당한 속력 범위 내의 샘플만 만들자. 계산하기 쉽도록 하기 위해서 $x^2 = mv^2 /kT$인 차원없는 변수 $x$를 사용한다. 확률이 보존되어야 하므로 $P(v)dv= f(x)dx$에서 $x$에 대한 확률밀도함수는

$$ f(x) = \sqrt{\frac{2}{\pi}} x^2 e^{-x^2/2}$$

로 주어진다. 가능한 샘플 $x$는 다음의 단계로 얻어진다.

1. 주어진 $x$ 범위 $[x_\text{min}, x_\text{max}]$ 사이에서 랜덤하게 한 값 $x$를 취하고,

2. $[0,1]$ 사이에서 랜덤한 값 $p$을 취한 후 

3. $ p> f(x)$이면 $x$는 버리고 다시 과정 1로 돌아간다. 

 

다음 예는 $x$를 $[0,4]$사이에서 500 구간으로 나누어서 만든 50000개의 샘플의 분포와 Maxwell-Boltzmann 분포를 비교한 것이다.

그림의 붉은점은 accepted된 값이고 푸른점은 rejected된 값이다.

    const CSize sz = CSize(500, 500);
    raster.SetDimensions(sz, 24);
    srand(unsigned(time(0)));
    int hist[500] = {0};
    int required = 50000;
    const double xmin  = 0, xmax = 4;
    while (required) {
        double x = double(rand()) / RAND_MAX; //[0,1]
        double v = (xmax - xmin) * x + xmin;  //[xmin,xmax]
        double p = double(rand()) / RAND_MAX; //[0,1]
        double maxwell = MaxwellDist(v);
        CPoint q = CPoint(int(x * (sz.cx - 1)), int(p * (sz.cy - 1)));
        if (p <= maxwell) {
            raster.SetPixel(q.x, sz.cy - 1 - q.y, 0xFF0000);
            hist[q.x]++;
            required--;
        }
        else
            raster.SetPixel(q.x, sz.cy - 1 - q.y, 0xFF);	
    }
    int s = 0;
    for (int i = 0; i < sz.cx; i++) 
        s += hist[i];
    FILE *fp = fopen("maxwell.txt", "w");
    for (int i = 0; i < sz.cx; i++)
        fprintf(fp, "%f, %f\n", (xmax - xmin) * double(i) /(sz.cx - 1) + xmin, \
                                double(hist[i]) / s);
    fclose(fp);
728x90
Posted by helloktk
,

면적이 $A$인 단순 평면 도형이 있을 때 그 둘레의 길이 $L$과의 사이에는 다음 부등식이 성립한다:

$$ L^2 \ge 4\pi  A$$

정 $n$ 각형에 대해서 이 부등식을 체크해보면

$$ \frac{L^2}{4\pi A} = \frac{\tan(\pi/n)}{\pi/n}=\left\{\begin{array}{cc} 1.654 & \text{정삼각형}\\ 1.273&\text{정사각형} \\ 1.156&\text{정오각형}\\ 1.034&\text{정십각형}\end{array} \right.$$

증명: 평면에서 둘레의 길이가 $L$인 한 폐곡선을 곡선을 따라 움직인 길이를 매개변수로 사용해서 표현하자: $(x(s), y(s)),~s\in [0,L]$. 이 폐곡선을 복소평면으로 옮기면 $z(s)= x(s) + i y(s)$로 표시할 수 있다. $z(s)$는 $s$에 대해서 주기가 $L$ 복소함수로 볼 수 있다. Fourier series를 사용하기 편리하게 하기 위해서 매개변수를 바꿔 $z(t) = x(Lt/2\pi) + i y(Lt/2\pi)$로 표현하자. 그러면 $t$에 대한 주기가 $2\pi$가 되고 다음과 같이 Fourier 전개를 쓸 수 있다. (영상처리에서는 물체의 윤곽선에 대한 Fourier  계수를 이용해서 물체의 모양을 분류하는데 이용이 된다)

$$ z(t) = \sum_{k=-\infty}^\infty a_k e^{i k t}$$

먼저 $L^2$항을 뽑아내기 위해 다음의 적분을 고려하자. $s=Lt/2\pi$로 놓으면 $z(s)$은 곡선의 길이를 매개변수로 사용한 것이어서  $|dz/ds|=1$이므로

$$  \int_0^{2 \pi} \left|  \frac{dz}{dt} \right|^2 dt = \int_0^{2\pi} \left|\frac{dz}{ds} \right| ^2 \left( \frac{ds}{dt}\right)^2 dt = \frac{L^2}{4\pi^2} \int_0^{2 \pi} dt = \frac{L^2}{2\pi}$$ 

좌변에 Fourier 전개를 대입하면 (note: $ \int_0^{2\pi} e^{i (k-j)t} dt = 2\pi \delta _{kj}$)

\begin{align} \int_0^{2\pi}\left|  \frac{dz}{dt} \right|^2 dt= 2 \pi \sum_{k=-\infty}^\infty k^2 |a_k|^2 \end{align}

을 얻는다. 따라서

$$ \frac {L^2 }{4\pi^2 } = \sum_{k=-\infty}^\infty k^2 |a_k|^2 $$

다음으로 도형의 면적으로 표현하자.

\begin{align} A &= \frac {1}{2} \int_0^ {2\pi } \left(x \frac{dy}{dt}- y \frac{dx}{dt}  \right)dt \\ &= \frac{1}{2} \text {Re} \int_0^{2\pi} i z \overline{\frac{dz}{dt}}  dt= \pi \sum _{k=-\infty}^\infty k |a_k|^2 \end{align}

그런데 $$ \sum_k k^2 |a_k| ^2 \ge   \sum _k k |a_k|^2$$

이므로

$$ \frac{L^2}{4\pi} \ge A$$

을 얻는다. 그리고 등호 성립하려면 Fourier 계수가 $k=0, 1$이외에는 모두 0이어야 된다. 즉

$$ z(t) = a_0 + a_1 e^{i t}$$

인 경우인데, 이는 복소평면에서  중심이 $a_0$, 반지름이 $|a_1|$인 원을 나타낸다. 변분법을 이용해서도 원이 주어진 둘레길이를 가지는 평면도형 중에서 가장 넓은 면적의 도형임을 쉽게 보일 수 있다. 

https://kipl.tistory.com/187

 

주어진 길이의 폐곡선으로 가둘 수 있는 최대 면적의 도형은?

평면에서 주어진 길이의 폐곡선으로 가둘 수 있는 최대의 면적은 얼마이고 그 모양은 어떤 것일까? 이 문제는 변분법을 이용하면 쉽게 해결할 수 있다. 평면에서 폐곡선은 매개변수 $t$의 함수로

kipl.tistory.com

 

728x90
Posted by helloktk
,