$\bf S$가 positive definite 행렬이고, $\bf C$는 대칭행렬일 때 아래의 일반화된 eigenvalue 문제를 푸는 방법을 알아보자. 
$$\bf  S u = \lambda C u$$

타원을 피팅하는 문제에서 이런 형식의 고유값 문제에 부딛히게 된다. 이 경우 $\bf S$는 scattering matrix이고, $\bf C$는 타원피팅에 걸리는 제한조건때문에 나온다. $\bf S$가 positive definite이므로  $\bf S$의 eigenvalues $\{\sigma_i > 0 \}$는 모두 0보다 크고, eigenvector를 이용하면 

$$ \bf S = R \Lambda R^T, ~~~~\Lambda=\text{diag}(\sigma_1, ...,\sigma_n)$$

처럼 분해할 수 있다. $\bf R$은 $\bf S$의 eigenvector를 열로 가지는 행렬로 orthogonal 행렬이다: $\bf R^{-1}=R^T$.

 

이제 $\bf S$의 제곱근 행렬을 $\bf Q, ~Q^2 =S$라면 

$$ \bf Q= R \Lambda^{1/2} R^T,~~~~\Lambda^{1/2} = \text{diag}( \sqrt{\sigma_1},...,\sqrt{\sigma_n})$$

임을 쉽게 확인할 수 있다. $\bf Q$을 이용하면 구하려는 고유값 문제는 

$$ \bf QQ u = \lambda C ~\to ~  Qu = \lambda Q^{-1} C Q Q^{-1}~~\to ~~ v = \lambda Q^{-1} C Q^{-1} v$$

이므로 $\bf Q^{-1} C Q^{-1}$의 고유값 문제 $(1/\lambda, \bf Qu)$로 단순화됨을 알 수 있다. $\bf Q$의 역행렬이 $\bf Q^{-1} = R \Lambda^{-1/2}R^T$임을 쉽게 체크할 수 있으므로 직접적으로 역행렬을 계산할 필요가 없어진다.

$$\bf \Lambda^{-1/2} = \text{diag}( 1/\sqrt{\sigma_1},...,1/\sqrt{\sigma_n})$$

 
 
 
728x90

'Mathematics' 카테고리의 다른 글

수치적으로 보다 정밀한 이차방정식의 해  (0) 2024.02.23
열방정식의 Green function  (0) 2024.02.12
지구의 나이는?  (0) 2024.02.11
n 차원 공의 부피  (2) 2024.02.07
타원의 둘레길이  (1) 2024.02.06
Posted by helloktk
,

영상처리에서 $2\times2$ 행렬의 고유값과 고유벡터를 수치적으로 구해야 하는 상황이 많이 생긴다. 이 경우 2차인 특성방정식을 풀어서 고유값을 구하게 되는데, 두 고유값의 차이가 매우 큰 경우 단순한 근의 공식을 적용한 결과는 수치적으로 부정확한 값을 줄 수 있다. 이를 피하는 방법을 간단히 살펴보자.

 

이차방정식 

$$x^2 - b x +c =0$$

의 두 해는 

\begin{align}x_1 = \frac{b + d}{2},\quad x_2 = \frac{b - d}{2},\qquad d= \sqrt{b^2 - 4c}\end{align}

로 표현된다.  고유값의 특성방정식으로 보면 계수 $b$는 두 고유값의 합, 판별식 $d$는 두 고유값의 차이다. 그런데 $b ^2 \gg c $ 판별식 $d$는 $b$와 매우 근접한 수치이므로 거의 비슷한 두 수의 뺄셈으로 표현된 $x_2$의 수치계산 결과는 매우 부정확해질 수 있다. 예로 

$$ x^2 - 10^8 x +1=0$$

의 해를 프로그램으로 구하면 

	double b = 1.0e8, c = 1.;
	double d = sqrt(b * b - 4. * c);
	double x1 = (b + d)/2., x2 = (b - d)/2.;
	TRACE("x1=%e, x2=%e", x1, x2);

프로그램을 돌리면 얻게 되는 결과는

$$ \tt  x1=1.000000e+008, \qquad x2=7.450581e-009$$

그런데 테일러 전개를 사용하면 $x_2$는 

$$x_2 = \frac{b}{2}\left( 1 - \sqrt{1-\frac{4c}{b^2}} \right)  \approx \frac{b}{2}\times \frac{2c}{b^2} = \frac{c}{b}=10^{-8}$$ 

이어야 한다. 따라서 $x_2$를 바로 근의 공식으로 구할 경우 25%정도의 오차가 생긴다. 이런 경우에는 근의 공식에 의한 직접적인 계산보다는 이차방정식의 특성을 이용하는 것이 수치적으로 더 정교한 답을 얻을 수 있다. 두 근의 곱이 $x_1 x_2 =c$이므로 $x_2$는 

$$x_2 = c / x_1\approx 10^{-8}$$

을 통해서 구하는 것이 더 수치적으로 정교하다.  일반적으로 두 근의 합  $b$가 양수인 경우는  $x_2$를, 음수인 경우는 $x_1$을 두 근의 곱을 이용해서 계산하는 것이 더 정교한 수치값을 얻을 수 있다.

// Find two roots of a quadratic eq: x * x - b * x + c = 0;
void Quadratic2(double b, double c, double& x1, double& x2) {
    double d = sqrt(b * b - 4. * c);
    x1 = (b + d)/ 2., x2 = (b - d)/ 2.;
    TRACE("x1=%e, x2=%e\r\n", x1, x2);
    if (b > 0) {
        if (x1 != 0) x2 = c / x1;
        TRACE("Improved: x1=%e, x2=%e\r\n", x1, x2);
    } else {
        if (x2 != 0) x1 = c / x2;
        TRACE("Improved: x1=%e, x2=%e\r\n", x1, x2);
    }
}

이 경우 결과는 훨씬 더 정교한 값을 보여준다.

$$\tt Improved: x1=1.000000e+008, \quad x2=1.000000e-008$$

 
 
 
 
 
 
 
 
 
 
 
 
 
 
728x90

'Mathematics' 카테고리의 다른 글

Generalized eigenvalues problem  (0) 2024.03.02
열방정식의 Green function  (0) 2024.02.12
지구의 나이는?  (0) 2024.02.11
n 차원 공의 부피  (2) 2024.02.07
타원의 둘레길이  (1) 2024.02.06
Posted by helloktk
,

1차원 열방정식은 

$$ u_t (x,t)= \kappa u_{x x}(x,t), \qquad -\infty < x < \infty,~t>0 $$

으로 초기조건이 주어지거나 경계조건이 주어야 한다. $u(x,t)$가 해라면 $u(\lambda  x , \lambda^2 t)$ 또한 해가 됨음 쉽게 확인해 볼 수 있다. 즉, 열방정식은 scale 변경에 불변이 특성을 가진다. 이 특성을 이용해서 열방정식의 특별한 해를 구해보자. scale 불변성을 보면 해의 $x,t$에 대한 의존성은 $x/\sqrt{t}\to \lambda x/\sqrt{\lambda^2 t} = x/\sqrt{t}$의 조합형태로 나타나야 될 것으로 보인다.

$$ u(x,t) = h(x/\sqrt{t})$$

그런데 에너지 보존에 의해서 계의 전체 열에너지는 시간에 무관하게 일정한 값을 가져야 한다. 열에너지가 온도에 비례하므로 총 열에너지는 다음 적분값에 비례한다.

$$ \int_{-\infty}^\infty u(x,t) dx = \text{const}$$

따라서 $u = h(x/\sqrt{t})$와 같은 형태의 해는 에너지 보존을 만족시킬 수 없다.

$$ \int _{-\infty}^\infty h(x/\sqrt{t}) dx = \sqrt{t} \int_{-\infty}^\infty h(y) dy \propto \sqrt{t}$$

그렇지만 $u(x,t) dx $가 스케일 불변꼴이면 가능하므로 

$$u(x,t) = \frac{1}{\sqrt{t}} h( x / \sqrt{t})$$

와 같은 형태로 시도하자. 그러면 열방정식의 좌우변의 항이

$$ u_{xx} = t^{-3/2} h'' (x/\sqrt{t}), \qquad u_t = -\frac{1}{2} t^{-3/2}\Big[ h(x/\sqrt{t}) + h'(x/\sqrt{t}) x/\sqrt{t} \Big] $$

이므로 이를 열방정식에 대입하면 ($y= x/\sqrt{t}$)

\begin{align} \kappa h''(y) + \frac{1}{2} h(y) + \frac{1}{2} y h'(y)=0 \\ \to ~ h''(y) + \frac{1}{2\kappa} (y h(y) )' = 0 \\ \to~ h'(y) + \frac{y}{2\kappa } h(y) = \text{const} \to 0 \end{align}

여기서 $\text{const}=0$ 일 때 해만 고려한다. 이 경우

$$ h(y) = A e^{- y^2/4\kappa} \qquad \to \qquad u(x,t) = \frac{A}{\sqrt{t}} e^{- x^2 / 4\kappa t } $$

인 해를 얻는다. $\int _{-\infty}^\infty u(x,t) dx = 1$로 정규화시키면 적분상수 $A$는

$$ A= \frac{1}{\sqrt{4\pi \kappa } }$$

임을 알 수 있다. 이렇게 정규화된 해를 열방정식의 기본해, Green 함수 또는 heat kernel이라고 부른다.

$$K(x,t) = \frac{1}{\sqrt{4\pi\kappa t} }e^{- x^2/4\kappa t}$$

$K(x,t)$는 $t\to 0^+$일 때 $\delta(x)$와 같은 특성을 가지는데, 물리적으로는 $t=0$ 시점에 원점에 $\delta-$함수로 표현되는 열원을 놓았을 때 열의 확산을 기술하는 해이다. 열방정식이 선형이므로  $x-$축을 따라 분포된 각각의 $\{ \delta(x-x_i)| i=1,2,\dots\}$ 열원의 선형결합도 해가 된다. 

$$ u(x,t) = \sum_i K(x-x_i, t) f_i$$

또한 연속적으로 분포하는 열원의 경우에는 합을 적분형태로 변경하면 되므로 $t=0$일 때 온도분포가 $f(x)$로 주어진 경우 다음 표현이

$$ u(x,t) = \int_{-\infty}^\infty K(x-y, t)f(y)dy$$

열방정식을 만족시키고, 초기조건도 만족시킴을 확인할 수 있다.

$$ \lim_{t\to 0} \int_{-\infty}^\infty K(x-y, t) f(y)dy = \int_{-\infty}^\infty \delta(x-y) f(y) dy = f(y)$$

 

해의 몇가지 재미있는 특징을 살펴보자. $t=0$일 때 열원이 유한한 영역에서만 0이 아니더라도 $K(x,t>0)$가 0이 아니기 때문에 매우 짧은 시간에 모든 영역으로 열이 퍼져나간다. 즉 열이 퍼지는 속도가 무한대에 해당한다. 열방정식은 비상대론적인 방정식이다. 두 번째로 $t>0$일 때 $K(x-y,t)$가 $x=y$인 지점에서 벗어나는 경우에 매우 빨리 0으로 변하므로 해는 무한번 미분가능하다. 즉, 처음 주어진 온도 분포가 불연속적이더라도 바로 매끄러운 온도분포로 변하게 만드는 특성이 있다. 이 열방정식을 잡음이 들어있는 신호나 영상에 적용하면 smoothing 효과를 줄 수 있다. 그리고

\begin{align}u(\lambda x, \lambda^2 t) &= \int K(\lambda x - y, \lambda ^2 t) f(y)dy \\ &= \int K(\lambda x - \lambda y' , \lambda^2 t )  {\lambda } f(\lambda y') dy' \\ &= \int K(x-y',t)  {\lambda }f(\lambda y') dy' \end{align}

이므로 스케일링된 해는 초기 온도분포가 $\lambda f(\lambda y)$인 경우와 같음을 알 수 있다. 초기에 국소적으로 가열된 경우라면 $\lambda f(\lambda x)$의 분포는 거의 $\delta$ 함수에 접근할 것이므로 충분히 시간이 지나고 먼 거리에서는 온도분포는 초기조건에 무관하게 됨을 알 수 있고, 이는 smoothing 효과로 이해할 수 있다.

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
728x90

'Mathematics' 카테고리의 다른 글

Generalized eigenvalues problem  (0) 2024.03.02
수치적으로 보다 정밀한 이차방정식의 해  (0) 2024.02.23
지구의 나이는?  (0) 2024.02.11
n 차원 공의 부피  (2) 2024.02.07
타원의 둘레길이  (1) 2024.02.06
Posted by helloktk
,

지구의 나이는?

Mathematics 2024. 2. 11. 22:16

초기의 지구는 매우 뜨거워 모든 암석들이 녹아있는 상태였을 것이다. 시간이 흐르면서 지구는 표면에서 차가운 우주(거의 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 남작에 의해 시도되었다.

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
728x90
Posted by helloktk
,

반지름이 $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
,