[작성자:] snowall

이번엔 왜 6번이냐고 물으신다면 몇번까지 했는지 기억이 나지 않기 때문…이라고.

# Elementary Numerical analysis 6

# based on Python

# (C) 2013. Keehwan Nam, Dept. of physics, KAIST.

# snowall@gmail.com / snowall@kaist.ac.kr

# Relaxation method

# http://snowall.tistory.com/2561 이 글의 개선된 버전이 되겠다.

# 라플라스 방정식이나 푸아송 방정식을 풀기 위해서 사용하는 가장 대표적인 방법은 뭐니뭐니해도 가우스 법칙이다. 가우스 법칙은 표면에서 필드를 적분하면 내부에 존재하는 전하량에 비례한다는 사실을 말해준다. 여기서 완화법(relaxation method)이 왜 작동하는지 설명할 수 있다.

# 라플라스 방정식의 해가 조화 함수(Harmonic function)이라고도 불리운다는 사실을 기억하자. 조화함수의 가장 중요한 특성 중 하나는 평균값 정리인데, 이 정리는 조화 함수가 잘 정의되는 임의의 점에서 성립한다. 주어진 임의의 점에서 거리가 일정한 점의 집합, 즉 어떤 ‘구면'(spherical surface)을 정의할 수 있다. 이 구면에서(3차원이나 2차원 표면이 아니라도) 조화함수의 함수값의 평균값은 구면의 중심점에서의 함수값과 같다.

# 그런데, 우리가 해결하고자 하는 문제의 공간은 격자로 잘게 쪼개져 있다. 따라서, ‘구면’은 바로 옆에 있는 점들로 정의된다. 자, 이제 2차원 구면을 생각해 보자.

# 함수 f(x, y)가 잘 정의되고, x0, y0에서 h만큼 떨어져 있는 ‘구면’은 다음과 같이 정의된다.

# 구면 = {(x0+h, y0), (x0-h, y0). (x0, y0+h), (x0, y0-h)}

# 2차원이고 네모난 격자이므로 이렇게 점 4개가 ‘구면’을 정의한다. 이 점 위에서의 함수값들의 평균값은 다음과 같이 정의된다.

# f_mean(x0, y0) = (f(x0+h, y0)+f(x0-h, y0)+f(x0, y0+h)+f(x0, y0-h))/4

# 적분은 덧셈으로 바뀌었고, 4개 있으니까 4로 나눈 것이다. 아주 쉽다.

# 조화 함수의 특성으로부터, f_mean(x0, y0) = f(x0, y0) 이어야 한다.

# 하지만, f(x,y)를 알고 있다면 당연히 이게 성립하는걸 확인할 수 있겠지만 모르는 마당에 어떻게 확인할 수 있을까?

# 라플라스 방정식이나 푸아송 방정식을 공부해 본 사람이라면 누구나 알고 있겠지만, 이 방정식들의 해는 만약 존재한다면 유일하게 존재한다. 즉, 수단과 방법을 가리지 않고 그런 해를 찾아냈다면 그 해는 우리가 찾아 헤메이던 바로 그것이다.

# 만약 함수값 f(x0, y0)을 위에서 계산한 f_mean(x0, y0)으로 대신한다면, 그래도 아무것도 모르는 것 보다는 답에 조금 더 가까울 것이다. f_mean(x0, y0)을 그렇게 두고, f_mean(x0+h, y0)을 계산하고, f_mean(x0, y0+h)를 계산하고, 그렇게 계산할 수 있다. 이렇게 h씩 움직이면서 모든 점을 한바퀴 돌고 나면, 아마 원래 알고 있던 함수보다는 조금 더 진짜 함수에 가까워진 모양이 될 것이다. 이 짓을 더이상 함수값이 변하지 않을 때 까지 반복하자.

# 이 방법은 왜 정답에 수렴할까? 왜냐하면, 경계조건이 있기 때문이다. 경계에서 주어진 함수값은 정해져 있고, 따라서 경계 근처에서 평균을 내면 어떻게든 경계의 함수값을 따라갈 수밖에 없다. 경계 근처에서 라플라스, 푸아송 방정식을 만족한다면, 그 옆에서도, 그 옆에서도, 어떤 위치에서도 다 만족할 수 밖에 없다. 충분히 오래 반복한다면 답은 반드시 나올 것이다.

# 이제, 예를 들어서 답을 찾아보자.

import numpy as np

f = []

for i in range(52):

f.append([])

for j in range(52):

f[-1].append(0)

# 0으로 가득 찬 52 x 52 짜리 행렬을 만들었다.

for i in range(52):

f[0,i] = sin(float(i))

f[-1,i] = sin(float(51-i))

f[i,0] = sin(float(51-i))

f[i,-1] = sin(i)

# 네모난 경계에 경계조건을 주었다.

desire = 1.e-10 # 내가 원하는 오차의 한계

error = 1.

while error>desire:

error = 0.

for i in range(1, 51):

for j in range(1, 51):

tmp_error= f[i, j]

f[i, j] = (f[i-1, j]+f[i+1, j]+f[i, j+1]+f[i, j-1])/4. # 문제를 해결하는 핵심 코드.

error+= abs(tmp_error-f[i, j])

result = open(“solution.txt”, “w”)

result.write(f)

result.close()

# 이게 루틴의 끝이다.

# 심지어 문제를 푸는 루틴이 경계조건을 부여하는 루틴보다 간단하다!

# 푸아송 방정식이라면?

while error>desire:

error = 0.

for i in range(1, 51):

for j in range(1, 51):

tmp_error= f[i, j]

f[i, j] = -h*h*rho(x[i], y[j])+(f[i-1, j]+f[i+1, j]+f[i, j+1]+f[i, j-1])/4. # 문제를 해결하는 핵심 코드.

error+= abs(tmp_error-f[i, j])

# -h*h*rho(x[i], y[i])이 추가된 것 외에 똑같다. rho는 푸아송 방정식에서 inhomogenius인 부분을 나타내는 함수이다. 저게 0이면 라플라스 방정식과 똑같다. x[i]와 y[j]는 i, j에서의 실제로 주어진 x, y값이다. 이것도 속도를 빠르게 하고 싶으면 미리 rho(x, y)를 rho[i, j]로 만들어 놓고 시작하면 된다.

# 완화법의 문제는 느리다는 점이다. 이것을 해결하기 위해서, 해가 대충 알려져 있는 경우 f(x, y)를 위에서 쓴 것처럼 0으로 배치하지 말고 대충 알려진 근사함수로 미리 넣고 시작하면 보다 빠르게 수렴할 것이다. 사실 미리 넣어주는 함수를 아무 함수나 집어넣어도 함수가 무한대로 뻗치지 않는 한 무조건 수렴한다. 그러나 빠르게 값을 얻고 싶다면 미리미리 비슷해 보이는 함수에서 시작하는 것이 조금 더 빠른 답을 얻는데 도움이 된다.

# 미리 넣고 시작하는 근사 함수를 얻기 위해서, 격자를 우선 성기게 잘라서 대충 답을 구하고, 이 답을 잘게 쪼갠 격자에 대입해서 시작하면 보다 빠르게 답을 얻을 수 있다. 이런 방법을 다중격자(multigrid) 방법이라고 한다.

# 이 방법을 확장한 Successive Over-relaxation (SOR) method 라는 방법이 있다. 이것은 주변의 평균값을 더할 때 가중치를 주어서 더 빠르게 수렴하도록 하는 방법이다.

# 라플라스 방정식과 푸아송 방정식 등, 타원꼴 편미분방정식은 이 방법으로 언제나 풀 수 있다. 고차원에서도 성립하므로, 2차원보다 높은 차원에 대해서도 손쉽게 확장할 수 있다. 문제는 격자를 쪼개는 것이 사각형이 아니면 매우 곤란해진다는 점이다. 이 문제를 해결하기 위해 유한요소법(Finite element method, FEM)이 도입되었다. 그러나 다음 시간에는 유한차 시간영역 법(Finite difference time domain, FDTD)을 다룰 것이다. FEM은 완화법을 보다 일반적으로 확장한 것이다. 어쨌든, 다음 시간에.

*내일 숙제 내는 사람들 베껴가기 없기.

자발적 대칭성 붕괴(Spontaneous symmetry breaking)는 자연에서 흔히 일어나는 대칭성 붕괴 과정의 하나이다. 이것이 무엇인지 설명하기 위해서는 우선 대칭성을 이해해야 하고, 그 붕괴가 무엇을 의미하는지 알아야 한다. 그리고 그것이 자발적인 것과 그렇지 않은 것이 어떤 차이가 있는지 알아봐야 할 것이다.

대칭성이란 일반적으로 변환에 대하여 변하지 않는 성질을 뜻한다. 물리학적으로는, 라그랑지안이나 해밀토니안 등 주어진 물리계를 설명하는 방정식이 변환을 작용하더라도 그 모양이 변하지 않는 것으로 정의한다. 그리고 물리계를 설명하는 라그랑지안에 대칭성이 있다면 그에 대응하는 보존량이 존재한다는 것이 알려져 있다. 실제로, 공간에 대한 이동이 물리 현상을 바꾸지 않는다는 사실(Space uniformity, translation symmetry)로부터 운동량 보존법칙이 성립하고, 보는 방향에 따라 물리 현상이 달라지지 않다는 사실(Rotation symmetry)으로부터 각운동량 보존법칙을 유도할 수 있다.

대칭성 붕괴는 따라서 다루려는 물리계에 대한 해밀토니안이 변환에 대하여 대칭성을 갖지 못하는 경우에 나타난다. 가령, 수소원자를 설명하는 해밀토니안은 회전 연산자에 대하여 대칭적이다. 따라서 각운동량이 단일한 고유값으로 나타나며, 모든 각운동량에 대해 같은 에너지를 갖는 겹친 상태로 관찰된다. 그러나 수소원자가 강한 자기장 안에 들어가게 되면 해밀토니안이 더이상 회전에 대해 대칭적이지 않게 된다. 왜냐하면, 자기장과 수소원자의 각운동량이 상호작용하게 되는데, 이 상호작용에서 나타나는 에너지가 자기장과 각운동량이 서로 어느 방향을 향하고 있느냐에 따라 달라지기 때문이다. 따라서 이 경우 대칭성이 붕괴되고, 서로 겹쳐져 있던 상태들은 다른 에너지 상태로 분리되어 분광선에 나타나게 된다.

물리계를 설명하는 해밀토니안이 본질적으로 공간에 대해 대칭성을 갖고 있는 경우더라도, 대칭적인 상태가 바닥상태가 아닐 수 있다. 가령, 강자성체의 경우 각각의 스핀 각운동량이 무작위적으로 배열된 상태보다 어떤 방향을 기준으로 동일하게 정렬된 상태가 더 낮은 에너지를 갖는다. 따라서, 이 경우 바닥상태는 대칭성이 깨진 상태이고 이 물리계는 대칭성이 깨진 상태에서 안정화된다. 이러한 현상은 우리 주변에서 흔히 발견할 수 있는데, 특히 가장 극적인 사례는 물질 입자의 질량을 만들어 내는 힉스 매커니즘이다.

우리 우주를 기술하는 입자 물리학의 표준모형에서, 표준 모형을 그대로 적용하면 질량을 나타내는 항은 대칭성을 깨트리기 때문에 자연스럽게 도입되지 않는다. 여기에 힉스 입자와의 상호작용을 도입하면, 에너지가 높은 상태에서는 라그랑지안 자체의 대칭성을 깨트리지 않으면서 질량에 해당하는 항을 자연스럽게 유도할 수 있다. 즉, 힉스 입자의 진공기대값이 0이 아닌 유한한 값으로 수렴하면서 힉스 입자와 상호작용하는 모든 입자는 힉스 입자의 진공기대값에 결합상수를 곱한 크기만큼의 질량을 얻게 되는 것이다.

자발적 대칭성 붕괴의 결과로 골드스톤 입자가 도입된다. 골드스톤 정리에 따르면, 연속 대칭성이 깨지는 경우 깨진 대칭성의 차원에 해당하는 골드스톤 입자가 나타나는데, 이는 대칭성이 깨지면서 해밀토니안이 대칭적이더라도 그 해밀토니안이 기술하는 물리계의 바닥상태가 대칭적이지 않기 때문이다. 대표적인 골드스톤 입자가 양성자와 중성자의 SU(2) 아이소스핀 대칭성의 붕괴에서 나타나는 파이온이다. (원래 골드스톤 입자는 질량이 없지만, 파이온은 질량을 갖고 있는데, 이건 어쨌든 사정이 있다. 이런 경우를 Pseudo-Goldstone boson이라 부른다.)

이렇듯, 자발적 대칭성 붕괴 현상은 우리 우주의 곳곳에서 찾을 수 있으며, 물리를 이해하는데 매우 중요한 틀이다.

2번이 왜 이제 올라오느냐도 역시 의미 없는 질문이다.

—

# Elementary Numerical analysis 2

# based on Python

# (C) 2013. Keehwan Nam, Dept. of physics, KAIST.

# snowall@gmail.com / snowall@kaist.ac.kr

# Differential equation

# 추천 참고자료: http://www.swarthmore.edu/NatSci/echeeve1/Ref/NumericInt/FrameNumInt.html

# 미분방정식은 미분 변수가 1개인 상미분방정식(Ordinary differential equation:ODE)과 2개 이상인 편미분방정식(Partial differential equation:PDE)으로 나눠진다.

# 일단 1개인 미방을 어떻게 푸는지 알아보자. 연립 미방은 아직 다루지 않고, 나중에 다룬다.

# ODE는 선형 미방과 비선형 미방으로 나눠진다. 선형미방은 미방의 해가 두개 이상 있을 때, 두 해의 선형결합도 해가 되는 경우이다.

# 즉, f1과 f2가 미방의 연산자 L에 대해서 L(f1)=L(f2)=0일 때 L(f1+f2)=0이 만족되는 경우이다. 물론 그 역을 선형 미방이라고 하지는 않는다.

# 비선형 미방은 풀기 골치아프므로 일단 선형 미방을 어떻게 푸는지 알아보자.

# 선형 미방은 다시 최대 차수인 도함수에 따라 n차 미방으로 구분된다. 즉, 가령 2차 도함수가 최대차수인 도함수라면 2차 미분방정식이 된다.

# 1차 이상의 미분을 포함하는 고차 미분방정식은 여러개의 중간 변수를 도입하여 1차 연립 미분 방정식으로 간주할 수 있다. 따라서, 1차 미분방정식을 잘 풀 수 있다면 고차 미분방정식은 손쉽게 풀린다. 그러므로 1차 미방에 집중하자.

# 1차 미분방정식을 이산화(discretize)하게 되면 차분방정식(difference dquation)이 된다. 차분방정식은 고등학교 때 배우는 계차수열의 한 형태이다.

# 그럼, 이제 차분방정식을 풀어보자. 가장 간단하게, 오일러 풀이법이 있다.

# 오일러 풀이법에서는 dy/dt + a*y(t) = f(t) 라는 문제를 풀게 된다. a는 값을 알고 있는 상수이고, f(t)는 이미 알고 있는 함수이다.

# 오일러 풀이법은 정말 간단하다. 주어진 시간 t로부터 h만큼의 시간이 지난 후의 값인 y(t+h)를 다음과 같이 근사한다.

# y(t+h)=y(t)+(dy/dt)h

# 아주 간단한 1차 근사가 되겠다. 그럼, 이제 y(t+h)를 알았으니, 같은 방법으로 한번 더 계산하면 y(t+2h)도 알아낼 수 있다. 이 작업을 반복하면, y(t=t0)를 알고 있을 때 우리는 아주 많은 문제를 풀 수 있게 된다.

# y(t+h) = y(t)+(dy/dt)h = y(t) + (- a*y(t) + f(t))*h

# 잘 생각해보면, 이것은 아주 간단한 반복작업이 된다. 이를 위한 파이썬 코드를 쓴다면, 다음과 같을 것이다.

import numpy as np

a=2.

f = lambda x, 3.*np.exp(-4.*x)

y0 = 0. # t=t0일때의 값이라고 하자.

h = 0.1 # time slice가 된다. 이 값은 작을수록 더 정확한 답을 얻게 되지만, 계산이 느려진다.

y = [y0] # 일단 여기서 시작해 보자.

t = 0.

tmax = 10. # 얻고싶은 시간의 마지막 순간. 즉, 이 시간 이후에는 관심이 없고 그 전까지만 관심이 있다는 뜻이다.

while t y.append(y[-1]+(f(t)-a*y[-1])*h)

t+=h

# while 안에 있는 코드가 가장 중요하다. 잘 보면, 별거 없다는 사실을 알 수 있을 것이다.

# 이제 2차 미방을 1차 미방 2개로 쪼개는 기법을 아주 간단히 이해하고 넘어가자.

# 주어진 2차 미방이 d(dy/dt)/dt + a*dy/dt + b*y(t) = f(t) 라고 하면, 일단 dy/dt = z 라고 하자. 그리고 z를 대입하면

# dz/dt + a*z + b*y = f(t)

# dy/dt = z

# 이렇게 2개의 미분방정식으로 쪼개졌다. 한번 풀어보도록 하자.

# 일단 y0와 z0=dy/dt(t=0)는 알려져 있다고 하자.

y = [y0]

z = z0.

# z는 우리가 문제를 풀기를 원하는 값이 아니라 중간 저장 변수이므로 그냥 이렇게 두면 된다. 만약 2차 미분방정식인데 1차 미분도 같이 구해야 하는 경우(뉴턴 방정식에서 속도와 위치를 둘 다 구하는 경우, 회로 방정식에서 전하와 전류를 둘 다 구하는 경우) 등등에서는 z를 리스트로 선언해서 써도 된다.

while t y.append(y[-1]+z*h)

z+=(-a*z-b*y+f(t))*h

t+=h

# 2차 미방이지만 위와 같이 딱 1줄 더 늘어난 아름다운 코드로 만들 수 있다.

# 오일러 풀이법은 가장 간단하지만 h가 상당히 작아야 정확한 답을 얻을 수 있다. 미분이나 적분을 할 때 함수값의 대칭성을 이용해서 같은 구간 h에 대해서도 보다 정확한 계산을 할 수 있었는데, 여기서도 마찬가지 테크닉을 사용할 수 있다.

# 룽게-쿠타(Runge-Kutta) 풀이법은 한발 앞서 간 도함수와 평균을 내서 다시 구한다. 즉, 기울기 계산을 두번한다. 말로는 알아듣기 어려울 것이고, 직접 코드를 보는 것이 보다 이해하기 쉬울 것이다.

# y(t+h) = y(t)+(dy/dt)h = y(t) + (- a*y(t) + f(t))*h

# dy/dt + a*y(t) = f(t)

y = [y0]

while t y1= -a*y[-1]+f(t) # 일단 기울기 계산

y.append(y1*h+y[-1]) # 일단 값 집어넣고

y2 = -a*y[-1]+f(t+h) # 한칸 간 곳에서의 기울기 계산

y[-1] = y[-2]+0.5*(y1+y2)*h # 실제 값은 두 기울기의 평균만큼만 증가하는 것으로.

t+=h

# 계산 시간을 아주 조금이라도 더 절약하고 싶으면 0.5*h를 미리 h2같은걸로 정의해놓고 0.5*h대신 h2를 곱하는 것으로 바꾸면 조금 빨라진다.

# y2를 계산하는 부분은 그 다음줄로 한번에 집어넣을 수도 있지만, 읽기가 힘들어지므로 굳이 그렇게 할 필요는 없다.

# 2차 룽게-쿠타 방법 뿐만 아니라 더 높은 차수의 룽게-쿠타 방법도 가능하다. 위와 같이 기울기 가중치를 어떻게 주느냐를 잘 생각하면 되는데, 아이디어는 비슷하므로 크게 어려울 것 없다.

# 조금 더 일반화하면 predictor-corrector 방법이 된다.

# 미분방정식이 dy/dt+f(y(t), t)=0으로 주어진 경우에 사용하는데, 위의 미분방정식은 f(y, t) = a*y-f(t)로 주어진 하나의 사례에 해당한다.

# Predictor-corrector 방법은 그럼 RK방법을 그대로 갖고 오면 된다.

y = [y0]

while t y1= f(y[-1], t) # 일단 기울기 계산

y.append(y1*h+y[-1]) # 일단 값 집어넣고

y2 = f(y[-1], t+h) # 한칸 간 곳에서의 기울기 계산

y[-1] = y[-2]+0.5*(y1+y2)*h # 실제 값은 두 기울기의 평균만큼만 증가하는 것으로.

t+=h

# 물리학에서 풀어야 하는 수많은 미분방정식들은 대체로 운동방정식이다. 운동방정식은 2차 미분방정식이므로, 2차 도함수를 알고 있다는 사실로부터 뭔가를 직접 해결할 수 있다.

# 그래서 나온 방법이 벌레(Verlet) 방법이다. Verlet은 사람 이름이고, 프랑스 사람이라 t가 묵음이다.

# 벌레 알고리즘에 대해서는 이 글을 참고해 볼 수 있다.

# http://www.physics.udel.edu/~bnikolic/teaching/phys660/numerical_ode/node5.html

# 일단, 테일러 전개로부터 다음과 같은 사실을 알 수 있다.

# y(t-h) = y(t)-v(t)*h+a(t)*h*h/2

# 여기서 v=dy/dt, a=dv/dt 이다.

# -h가 아니라 +h에 대해서 다시 쓰면

# y(t+h) = y(t)+v(t)*h+a(t)*h*h/2

# 이렇게 될 것이다.

# 두개를 더해보자.

# y(t-h)+y(t+h) = 2y(t)+a(t)*h*h

# v가 사라졌다!

# 이제 다음과 같이 y(t+h)를 추정할 수 있다.

# y(t+h) = 2y(t)-y(t-h)+a(t)*h*h

# 물리 문제를 해결할 때, 이 공식은 매우 기가막힌 방법이다. 2차 미분방정식을 풀어야 하는데, 속력을 구하지 않고 힘으로부터 곧바로 위치를 얻어내기 때문이다.

# 또한, 두개를 빼면 속력을 얻을 수 있다.

# v(t) = (y(t+h)-y(t-h))/2h

# 이 공식은 미분에서 공부했던 3점 공식과 같다.

# 이렇게 계산하는 방법을 original Verlet 알고리즘이라고 한다. 그런데, 오리지널 벌레는 문제가 있다. y(t=0)일때를 알고 있어도, y(t+h)를 구하기 위해서는 y(t-h)를 알아야 한다. 즉, 최초 2점을 알아야 한다는 문제가 있다. 그래서 이 방법을 개선하기 위해서 velocity Verlet 알고리즘이 등장한다. 일단 공식부터 보고 시작하자.

# y(t+h) = y(t)+v(t)*h+0.5*a(t)*h*h

# v(t+h) = v(t)+0.5*(a(t+h)+a(t))*h

# 왠지 앞에서 보았던 오일러 풀이법으로 2차 미분방정식을 푸는 것과 비슷해 보이지만, 그보다 더 정확하다. 그리고 이 경우 알지도 못하는 과거인 y(t-h)를 몰라도 문제를 해결할 수 있으므로 고민이 한층 덜어졌다.

# 이 문제에서 a(t) 가 알려져 있다고 하면, 코드를 다음과 같이 짜볼 수 있을 것이다.

y=[y0]

v=v0 # 2차 미방이므로 y0와 v0는 알고 있어야 한다.

while t y.append(y[-1]+v*h+0.5*a(t)*h*h)

v+=0.5*(a(t+h)+a(t))*h

t+=h

# 이외에도 많은 알고리즘이 있지만 대체로 이런 아이디어들의 변형과 개선을 도모한 방법들이다. 실제로 자신이 풀려고 하는 문제에 적합한 알고리즘을 선정하고, 그 문제에 적합하게 알고리즘을 개조/개선해서 적용하는 것이 좋다.

# 수치해석적인 방법의 장점은 선형 미분방정식 뿐만 아니라 비선형 미분방정식의 경우에도 적용할 수 있다는 것이다. 사실, 선형 제차 미분방정식의 해는 매우 쉽게 해석적으로 구할 수 있다. 그러나 비제차 미분방정식이나 비선형 미분방정식의 경우 해가 알려진 것이 매우 적으나, 실용적 관점에서 수치해석을 이용하여 문제를 해결할 수 있다.

# 수치해석으로 비선형 미분방정식을 해결할 때 주의해야 하는 것은 오차의 전파이다. 비선형 효과에 의해 오차가 예상보다 커질 수 있으므로 오차를 줄이기 위해서 기존의 공식을 적절히 변형하는 것이 필요할 수 있다.

http://news.khan.co.kr/kh_news/khan_art_view.html?artid=201309012128025&code=990100

강신주 교수님을 싫어하지는 않지만, 칼럼에 동의할 수 없는 부분이 있어서 적어둔다.

설악산에 케이블카를 설치하는 것이 환경을 파괴한다는 점은 동의한다. 그리고 설악산에 케이블카를 설치하는 것이 얻는 것 보다는 잃는 것이 더 많을 것이라고 생각한다. 또한 나는 설악산에 케이블카를 설치하는 것에도 반대하는 편이다.

하지만 산을 오르는 것에 자신의 걸음으로만 올라가야 반드시 그 기쁨을 느낄 수 있다는 점에 있어서는 동의할 수 없다. 어떤 사람이 기쁨을 느끼는 행위는 다른 사람의 통제나 지시를 받을 수 없다. 그 자신에게는 자신의 다리로 걸어가는 것이 더 희열을 느낄 수 있을지 모르지만, 또 대다수의 사람들이 그렇게 생각할 수도 있겠지만 그럼에도 불구하고 그런 사람이 많다는 점이 누구나 그래야 한다는 결론을 짓는 근거가 되지는 않는다.

어떤 사람은 다리가 불편해도 두 팔로 기어서 산의 정상을 정복할 수 있다. 어떤 사람은 다리가 불편하여 다른 사람의 등에 업혀서 정상을 오를 수 있다. 또 누구는 돈이 많아서 다리가 불편하지만 헬리콥터로 정상보다 더 높은 곳을 볼 수 있을지도 모른다. 다리가 멀쩡한 사람들도 이 세가지 사례를 마찬가지로 적용할 수 있다. 그럼 그중 누군가는 ‘진짜 기쁨’을 느끼지 못한채 가짜로 기뻐하고 있는 것일까?

누구든지 ‘이렇게 해보니까 더 즐겁더라’라고 제안할 수는 있다. 자신이 소중하게 생각하는 어떤 사람에게 그렇게 꼭 해보라고 강력하게 추천할 수도 있다. 자신의 명령을 들어야 하는 하급자에게 굳이 그렇게 해야만 한다고 지시할 수도 있다. 하지만 그렇게 했을 때 너는 반드시 기뻐해야 한다고 강요할 수 없고, 그렇게 하지 않았을 때 진짜 즐거움을 모를 거라고 강요할 수도 없다. 강요된 감정이야말로 진짜 즐거움과 거리가 멀다.

과학을 연구하다가 등장한 질문은 단 하나도 버릴 질문이 없다.

음…

다 논문거리다. 이미 논문이 나와 있거나, 아직 논문이 나오지 않았거나 둘 중 하나일 뿐.

http://imnews.imbc.com/replay/nwdesk/article/3329790_5780.html

도난 경보기에 걸리지 않는 특수가방을 이용해서 도둑질을 하던 부부가 붙잡혔다.

이 특수가방의 원리는 ‘패러데이의 새장’이라고 부른다. 그런 원리가 있다.

http://en.wikipedia.org/wiki/Faraday_cage

도체 안에 완전히 둘러싸여있는 물체는 외부의 전자기장의 영향을 받지 않게 된다는 원리인데, 실생활에서 여기저기 사용되거나 체험할 수 있는 원리이다.

엘리베이터에서 무선전화가 자주 끊기는 이유라든가, 전자레인지 옆에 있어도 전자파의 영향을 그다지 받지 않는다거나, 전자파 차단 필름도 같은 원리를 사용한다.

위의 범행에 사용한 특수 가방은 그다지 특별할 것 없이, 은박지 여러겹을 겹쳐서 가방 안에 잘 넣은 것이다. 만들기가 어렵지도 않고, 누구나 간단하게 만들 수 있다. 하지만 이걸 실제로 믿고 사용해도 된다는 사실을 알 정도로 물리학을 열심히 공부한 사람이 이런 나쁜 짓을 하다니 안타깝기 짝이 없는 사건이다.

도난방지 장치의 원리는 간단한 무선전력송신과 무선통신이다. 매장 입구와 출구 양 옆에 설치된 기둥을 본 적이 있을 것이다. 이 기둥에서는 전자기파가 흘러나오고 있다. 이 전자기파는 사람 몸에는 많이 흡수되지 않겠지만, 적당히 잘 설계된 안테나가 그 근처를 지나가면 그 안테나에 확 빨려들어가서 안테나에 전류를 흐르게 할 수 있다.

안테나에 흐른 전류는 도난방지장치의 전원을 켜고 작동시킨다. 아주 잠깐동안 켜진 도난방지장치는 자신이 살아있다는 목소리를 전파를 통해서 송신하고, 이 전파는 문에 설치된 기둥에 포함된 다른 안테나에 검출되어서 도난방지장치가 문짝을 지나가고 있다는 사실을 알린다. 그럼 큰 소리가 나면서 경고를 보내고 도둑이 붙잡히는 것이다.

여기서 이 특수가방은 전자파를 차단시켰으므로 전력의 송신도 안되고 도난방지장치의 전파가 밖으로 빠져나오지도 못한다. 즉, 특수가방은 도난방지장치를 성공적으로 감춰버린 것이다.

이 문제는 어떻게 해결할 수 있을까?

간단한 해결 방법으로는 도난방지장치 자체에 배터리를 달아서 외부와 통신을 계속 하고 있다가 통신이 끊기면 그 순간 소리가 나도록 하는 것이다. 소리는 외부에서 통신하고 있던 장비가 ‘어? 조금전까지 살아있던 놈이 왜 연락이 없지?’ 싶은 순간 울리도록 하면 될 것이다.

하지만 이 방법은 배터리 수명이 있기 때문에 배터리가 다 떨어지면 훔쳐가기가 더 쉬워진다는 문제가 발생한다.

따라서 앞서 이야기한 무선전력송신을 다시 한번 사용하여 배터리를 대신할 수 있다. 배터리 대신 무선전력송신 안테나를 매장 곳곳에 설치한다. 벽에 붙여서 설치하면 옷이나 다른 물건에 가려서 잘 안보이게 되므로 미관을 해치지 않고 깔끔하게 설치할 수 있다. 그럼, 마찬가지로 도난방지장치는 항상 작동하고 있게 되고, 그러다가 특수가방에 들어가버려서 갑자기 연락이 끊기면 그 순간 경보음을 울리도록 하면 될 것이다.

이 경우, 내가 도둑이라면 공범을 사용해서 매장 안의 두 곳에서 동시에 시도할 것이다. 두곳에서 동시에 특수가방에 옷을 넣으면 경보는 한번만 울리게 되고 그럼 적어도 둘 중 하나는 성공해서 매장을 빠져나갈 수 있을 테니까.

그럼, 이런 상황을 방지하기 위한 시나리오는?

도난방지장치의 위치를 추적해야 하는데, 이때 도난방지장치의 위치는 도난방지장치가 내보내는 신호의 세기로부터 알아낼 수 있다. 무선 신호의 세기는 거리의 제곱에 반비례하므로, 도난방지장치의 출력을 미리 알고 있다면 서너군데의 안테나에서 검출된 신호의 세기를 이용하여 도난방지장치가 매장 내에서 어디에 있는지 위치를 정확히 파악할 수 있다. 그럼 그 신호가 완전히 없어진 순간, 바로 직전에 있던 위치에 경고의 불이 번쩍번쩍 하도록 설치하면, 그리고 여러군데서 동시에 울릴 수도 있도록 설치한다면 이와 같은 특수 가방을 이용한 절도를 차단할 수 있다.

문제는 비용인데, 아마 이런 설비를 구현하고 설치하려면 매장 전체에 전력 송신 안테나와 위치 수신기를 설치해야 하므로 기존 장비의 열배 정도는 비용을 들여야 할 것이다. 한 벌에 수천만원 정도 하는 명품 옷이나 가방 같은 물건을 파는 매장이라면 고려해 볼 수 있을 것 같다.

이 글은 다음 글을 조금 더 고쳐서 다시 쓴 것이다.

http://snowall.tistory.com/501

파이썬 구현 예제는 생각중이다.

왜 4번이 먼저 올라왔느냐는 중요한 문제가 아니다. 나도 모른다.

# Elementary Numerical analysis 4

# based on Python

# (C) 2013. Keehwan Nam, Dept. of physics, KAIST.

# snowall@gmail.com / snowall@kaist.ac.kr

# Fast Fourier Transform

# 푸리에 변환은 적분변환이다. 다시말해서, 적분을 잘 수행하면 푸리에 변환은 계산된다.

# 0부터 1까지 구간만 생각해 보자.

# F1(k) = Integrate f(x)cos(kx)dx from 0 to 1

# F2(k) = Integrate f(x)sin(kx)dx from 0 to 1

# 위의 두 수를 k를 정해놓고 적분하면 된다. 하지만 현실은 그렇게 만만하지 않은데, 1개의 k에 대해서만 푸리에 계수를 구하는 것이 아니라, k에 대한 함수를 구해야 하는 사태가 벌어지기 때문이다.

# 따라서, 다음과 같이 하자. 조금 바꾼다.

# F = F1+i*F2

# 여기서 i는 제곱하면 -1이 나오는 허수 단위이다.

# 이렇게 되면 F는 오일러 공식을 사용해서 다음과 같이 쓸 수 있다.

# F(k) = Integrate f(x)exp(ikx)dx from 0 to 1

# 이제, 이 함수를 잘 적분하면 되는데, 우리는 지금 수치해석을 공부하고 있으므로 수치적으로 적분하면 된다.

# 0부터 1까지를 N등분하면 dx = 1/N이고 k는 k/N로 생각할 수 있다.

# F(k) = Summation f(x)exp(ikx/N) over integer x from x = 0 to x = N

# 위의 공식을 이산 푸리에 변환Discrete Fourier Transform이라고 한다.

# N등분된 구간은 k가 N개 있으므로, N개의 k값에 대한 각각의 함수를 계산하기 위해서 N개의 값을 더해야 한다. 즉, N*N번의 계산이 필요하다.

# 이제 exp(ikx/N)을 아주 잘 관찰해야 한다.

# f(x)exp(ikx/N)를 계산할 때, 잘 관찰해보면 같은 값들이 반복되는 것을 알 수 있다. k에 대한 값을 계산한 후 k+1에 대한 값을 계산한다면, f(x)exp(i(k+1)x/N)=f(x)exp(ikx/N)exp(ix/N) 이다.

# http://snowall.tistory.com/501

# 이제, 가령 점이 4개만 있다 치고, exp(ikx/N)을 대충 퉁쳐서 w(k*x)라 치고 그냥 대입해보자.

# F(0)=f(0)w(0)+f(1)w(1*0)+f(2)w(1*0)+f(3)w(1*0)

# F(1)=f(0)w(0*1)+f(1)w(1*1)+f(2)w(2*1)+f(3)w(3*1)

# F(2)=f(0)w(0*2)+f(1)w(1*2)+f(2)w(2*2)+f(3)w(3*2)

# F(3)=f(0)w(0*3)+f(1)w(1*3)+f(2)w(2*3)+f(3)w(3*3)

# 뭔가 있어 보인다면, 계산을 완성하자.

# F(0)=f(0)w(0)+f(1)w(0)+f(2)w(0)+f(3)w(0)

# F(1)=f(0)w(0)+f(1)w(1)+f(2)w(2)+f(3)w(3)

# F(2)=f(0)w(0)+f(1)w(2)+f(2)w(0)+f(3)w(2)

# F(3)=f(0)w(0)+f(1)w(3)+f(2)w(2)+f(3)w(1)

# 잘 보면 반복되는 것들이 보인다는 사실을 알 수 있다.

# 따라서 반복되는 부분을 미리 계산하고 재활용하면 된다.

# E1=f(0)w(0)+f(2)w(0)

# E2=f(0)w(0)+f(2)w(2)

# E3=f(1)w(0)+f(3)w(0)

# E4=f(1)w(1)+f(3)w(3)

# F(0)=E1+E3

# F(1)=E2+E4

# F(2)=E1+E3*w(2)

# F(3)=E2+E4*w(2)

# 원래는 16번의 곱셈과 12번의 덧셈이 있었는데, 잘 고쳤더니 10번의 곱셈과 8번의 덧셈으로 줄어들었다.

# 만약 점의 수가 수십만개로 늘어난다면, 아마 그 효율은 매우 좋아질 것이다. 즉 N*N번의 연산이 N*logN번의 연산으로 확 줄어들게 된다.

# 이 알고리즘은 Cooley-Tukey 알고리즘의 특수한 경우이며, 가우스가 최초에 발견하고 쿨리와 튜키가 나중에 다시 발견했다.

# http://ko.wikipedia.org/wiki/고속_푸리에_변환

# 파이썬 구현은 생략한다.

# http://docs.scipy.org/doc/numpy/reference/routines.fft.html 파이썬에서는 numpy 패키지에서 제공한다.

http://www.zdnet.co.kr/news/news_view.asp?artice_id=20130902091807&type=xml

지디넷에서 별을 삼키지 않는 블랙홀을 찾았다고 보도했다. 뭐 좀 이상해 보이는 블랙홀인데, 어쨌든 그 이름을 새지태리어스 A*로 이름을 붙였다고 한다.

새지태리어스는 궁수자리다.

http://en.wikipedia.org/wiki/Sagittarius_A*

그리고 궁수자리 A*는 꽤 오래전에 발견된 블랙홀이다.

음… 즉, 이번에 새로이 발견된 블랙홀이 아니라, 예전에 발견했던 블랙홀을 계속 관찰해서 새로운 사실을 알아냈다는 뜻이다.

http://www.sci-news.com/astronomy/science-chandra-gas-supermassive-black-hole-sagittarius-a-01349.html

이번에 알아낸 것은 다른 영역에서는 꽤 밝은 편인데 엑스선 영역은 왜 어두운지 알아냈다는 것이다.

외국 언론사의 보도 내용을 갖고 오는건 좋지만, 그래도 뭘 좀 아는 기자가 제대로 보도해준다면 좀 더 흥미로운 소식이 되지 않을까.