블로그에 재미난 글을 써보려고 고민을 좀 하다가 문득 제가 요즘 즐기는 게임인 PLAYERUNKNOWN'S BATTLEGROUND(이하 배그)에서 나오는 여러 가지 현상을 물리학적으로 나름 엄밀하게(기준은 최소 고등학교 물I, 물II) 분석해보면 재미있을 거 같아서 새로운 주제의 글을 쓰게 되었습니다. 물I, 물II에서 재미없게 배운 물리학 지식을 이런 곳에서 적용해보면 나름 좋을 거 같아서 말입니다. 그래서 오늘 쓰는 글은 배그에서 가장 처음으로 만날 수 있는 물리학적 현상인 '종단 속도'입니다.

- 배그에서 낙하 속도는 최대 234km/h가 나온다. -

 

  물II에서 포물선 운동을 배울 때 자유낙하라는 것을 배우긴 하지만, 이때는 공기 저항을 고려하지 않기 때문에 단순히 연직방향으로 중력가속도만 작용해서 등가속도 직선 운동을 합니다. 하지만 현실 세계에는 공기 저항이 작용하여 어느 이상의 속력을 내지 못하도록 만들죠. 만약 물II에서 배운 것처럼 공기 저항 없는 자유 낙하가 현실 세계에서 성립했다면, 빗방울 맞고도 골로 가겠죠. 물II에서 배우진 않았지만, 현실 세계에서 공기 저항이 존재한다는 것은 누구나 아는 사실입니다. 그런데 과연 배그에 나오는 종단 속도 234km/h는 정말 현실 세계에서도 똑같이 나올까요?


  먼저 이를 물리적으로 엄밀하게 분석하려면 자유 낙하하는 신체에는 어떤 힘이 작용하는지 알아봐야 합니다. 먼저, 중력이 작용할테고, 물I 유체역학 부분에서 배운 부력도 작용하고, 공기 저항으로 인해 발생하는 항력(drag force)이 있을 것입니다. 각 힘을 $F_g, F_b, F_d$ 라 하면 아래와 같이 표현할 수 있습니다.

 

이제 수식으로 정리해서 종단 속도를 유도해내면 됩니다. 하지만 고등학교 과정에서 항력은 배우지 않았기 때문에 좀 찾아봤습니다. 항력 공식은 아래와 같습니다.

$$F_d  = -\frac{1}{2}\rho v^2 A C_d \hat{\mathrm{v}} $$

설명을 하자면, $\rho$는 유체의 밀도, $v$는 유체에 대한 물체의 상대 속도, $A$는 운동 방향에 수직인 평면에 대한 정사영 면적, $C_d$는 항력 계수, $\hat{\mathrm{v}}$는 속도의 방향을 나타낸 단위 벡터입니다.

 

이제 이들 공식을 이용해서 운동 방정식을 세워보면 다음과 같습니다.

$$F_d + F_b + F_g = ma$$

$$-\frac{1}{2}\rho v^2 A C_d \hat{\mathrm{v}} - \rho g V + mg = ma$$

이때 $a = 0$ 일 때 힘의 평형을 이루므로 종단속도를 가지게 됩니다. 따라서 식을 정리해보면 아래와 같습니다.

$$v = \sqrt{\frac{2\left(mg - \rho g V  \right)}{\rho A C_d}}$$

 

[2] 자료를 참고해서 상수를 아래와 같이 결정해보고 구해보겠습니다.

$m$: 70kg

$g$: 9.8m/s2

$C_d$: 0.7

$\rho$: 1 kg/m3

$A$:  0.18 m2(머리 방향으로 낙하할 때)

$V$: 0.075 m3 ([3]에서 인체의 평균 밀도 참고)

 

$$v = \sqrt{\frac{2\left( 70 \times 9.8 - 1 \times 9.8 \times 0.075\right)}{1 \times 0.18 \times 0.7}}  = 약\ 104.29 m/s = 375.444 km/h$$

* 사실 부력은 되게 작아서 없다고 보고 계산해도 됩니다.

 

흠 생각보다 더 빠르게 나오네요. [2]에 따르면 배를 아래로 했을 경우에는 200 km/h, 머리를 향했을 때 240~290 km/h가 나오고 항력을 최소한으로 했을 때가 480 km/h 나온다네요. 인체의 몸에 맞게 $C_d$를 좀 더 수정하면 실제 값에 가까워질 거 같네요.

 

  다음으론 단순히 종단 속도만 구하긴 아쉬우니까 어떻게 가속이 되어 종단 속도까지 도달하는지 분석해봅시다. 위에 있는 운동 방정식에서 a = 0으로 두지 않고 미분 방정식을 풀면 됩니다. 1계 비선형 미분방정식인 거 같네요. 어떻게 풀질 몰라서 WolframAlpha를 참고하면서 풀어봤습니다.


일단 아래와 같이 상수를 좀 줄여서 깔끔하게 한 상태로 방정식을 정리하고 풀어봤습니다.

$$ a = -\frac{1}{2}\rho A C_d,\ b = - \rho g V + mg,\  c = m \\ a v^2 + b = c \frac{dv}{dt}\\ 1 = \frac{c}{av^2 + b} \frac{dv}{dt}\\ \int\, dt = \frac{c}{b} \int \frac{1}{\frac{a}{b} v^2 + 1}\, dv \\ \sqrt{\frac{a}{b}} v = \tan u,\ \sqrt{\frac{a}{b}}\ dv = \sec^2 u\ du \\ \begin{align} \frac{c}{b} \int \frac{1}{\frac{a}{b} v^2 + 1}\ dv &= \frac{c}{b} \times \sqrt{\frac{b}{a}} \int \frac{\sec^2 u}{\tan^2 u + 1} \, du \\ &= \frac{c}{b} \times \sqrt{\frac{b}{a}} \int\, du \\ &= \frac{c}{b} \times \sqrt{\frac{b}{a}} \tan^{-1}\left(\sqrt{\frac{a}{b}} v\right) + C \end{align} \\\\  t = \frac{c}{b} \times \sqrt{\frac{b}{a}} \tan^{-1}\left(\sqrt{\frac{a}{b}} v\right) + C \\ \therefore v = \sqrt{\frac{b}{a}} \tan \left(\frac{\sqrt{ab}}{c}\left(t + C\right)\right)$$


이때, $t = 0$일 때 $v = 0$ 이므로, 속도는 아래와 같습니다.

$$v = \sqrt{\frac{b}{a}} \tan \left(\frac{\sqrt{ab}}{c}t \right)$$

(어째 tan 함수인게 불안한데...)


이제 위의 상수값을 대입하여 a, b, c의 값을 계산하고 정리하여 속도 방정식을 완성해보겠습니다.

$$a = -\frac{1}{2}\rho A C_d = -\frac{1}{2} \times 1 \times 0.18 \times 0.7 = -0.063 \\ b = - \rho g V + mg = - 1 \times 9.8 \times 0.075 + 70 \times 9.8 = 685.265 \\ c = m = 70 \\ \therefore v = -104.29 i \tan \left(0.09i\ t \right)$$

* 104.29i 가 아니라 -104.29i인 이유는 $a < 0, b > 0$에서 $\sqrt{\frac{b}{a}} = \sqrt{\frac{b}{-a}} \times \frac{1}{i} = - i \sqrt{\frac{b}{-a}}$ 이기 때문입니다. 이거 덕분에 함수가 음수 나와서 20분 동안 고민했네요;;; 역시 수1 같은 기초를 열심히 해둬야...


아무튼 이렇게 상수까지 다 때려박아보니 속도가 한없이 증가할 거 같이 생긴 tan 함수에 고2 이후론 거의 보지도 못한 i가 함수 안쪽과 바깥쪽에 등장하고 있네요. 뭔가 싶어서 WolframAlpha에 때려 넣어보니 tanh 함수로 바꿔줍니다. 따라서 아래와 같이 됩니다.

$$v = 104.29 \tanh \left(0.09\ t \right)$$


그려보면 아래와 같습니다. (x축이 시간(s)축, y축이 속도(m/s)축)


이렇게 해서 종단 속도도 구해보고, 시간에 따라 실제 낙하 속도가 어떻게 변하는지 확인해봤습니다. 다음 번에는 뭘 해볼지 고민 해봐야 겠습니다.

(정작 종단 속도 구하는 것보다 미방 푸는데 시간을 훨씬 많이 쓴 듯)


참고 자료

[1] https://ko.wikipedia.org/wiki/%ED%95%AD%EB%A0%A5

[2] https://en.wikipedia.org/wiki/Speed_skydiving

[3] https://en.wikipedia.org/wiki/Orders_of_magnitude_(density)

  1. Minjae Isaac Kwon 2018.04.01 14:00 신고

    F_r 이 v^2 = (dx/dt)^2 에 Dependent 하도록 풀면 (2계도) 비선형 미방이라 ♩♩♬맞습니다.
    그냥 근사 쳐서 v=dx/dt 에 맞게 풀면 일반해도 있고 리니어텀이라서 공간차원 추가해도 먹고 살기 편합니다.
    (2계도 2차 비선형 미방은 2차원으로 가는순간 수치해석밖에 노답)

  지난 번 글에서는 간단하게 왜 fractal 그리는 작업에 GPU가 필요한지에 대해 간단히 이야기를 해봤습니다. 이번에는 그 fractal을 그리는데 사용되는 CUDA에 대해 설명을 해보도록 하겠습니다.


  CUDA(Compute Unified Device Architecture)는 C언어(정확히는 CUDA C)를 사용해서 GPU에서 알고리즘을 작성할 수 있는 GPGPU 기술입니다. 정확히는 오직 Nvidia GPU에서만 지원하는 기술이죠. AMD GPU나 기타 다른 GPU에서도 비슷한 역할을 하는 OpenCL(Nvidia GPU에서도 돌아갑니다)이라는 것도 있지만, 제가 알기로는 CUDA가 더 좋다고 알고 있습니다. 요즘엔 딥 러닝 때문에 기본적으로 연산을 위해 GPU를 쓰는 경우가 있는데, 대체적으로 딥 러닝 프레임워크가 CUDA를 지원해서 딥 러닝 할 때는 어쩔 수 없이 Nvidia GPU를 고를 수 밖에 없는 일이 발생하기도 했습니다.


  원래 GPU는 CPU와 다르기 때문에 기본적으로 CUDA를 이용해 병렬 프로그래밍을 짤 때에도 기존 프로그래밍과는 조금 다른 방식으로 프로그래밍을 해야 합니다. 무엇보다 RAM을 동시에 같이 쓰지 않기 때문에 CUDA에서는 컴퓨터(Host)와 GPU(Device) 사이의 메모리 복사를 통해 작업물과 결과물을 주고 받습니다. 그리고 Kernel이라는 GPU에서의 function을 짜고 이를 GPU에서 돌려서 연산을 수행합니다. 즉, 기본적인 알고리즘은 아래와 같이 이뤄지게 됩니다.


1. CPU에서 GPU로 데이터 전송(Host to Device 메모리 복사)

2. GPU에서 kernel을 통해 데이터 연산을 수행(kernel 연산)

3. 연산 결과물을 GPU에서 CPU로 데이터 전송(Device to Host 메모리 복사)


이 과정을 원래는 CUDA의 기본 지원 언어인 C/C++에서 짤 때 일일히 구현해야 하는데, 제가 fractal 그리는 프로그램을 짤 때는 PyCUDA를 사용했기 때문에 저 과정을 간단하게 작성할 수 있었습니다. 제 프로그램에서 CUDA 부분만 가져와보면 아래와 같습니다.


kernel 부분

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
__global__ void calcMandelbrot(bool *isJulia, double *c, double *matrix, int *result) {
    int idx_x = blockIdx.x * blockDim.x + threadIdx.x;
    int idx_y = blockIdx.y * blockDim.y + threadIdx.y;
    double a = matrix[0+ (float)idx_x * matrix[2];
    double b = matrix[1+ (float)idx_y * matrix[2];
 
    double an = a, bn = b;
    double aan = an * an, bbn = bn * bn;
 
    for (int i = 0; i < 400; i+=5) {
        // iteration 1
        if (isJulia[0]) { // Julia Set
            bn = 2 * an * bn + c[1];
            an = aan - bbn + c[0];
            aan = an * an;
            bbn = bn * bn;
        } else { // Mandelbrot Set
            bn = 2 * an * bn + b;
            an = aan - bbn + a;
            aan = an * an;
            bbn = bn * bn;
        }
 
        if (an * an + bn * bn > 4.0f) {
            result[idx_x * gridDim.y * blockDim.y + idx_y] = i;
            break;
        }
 
        ... 생략 ...
 
    }
}
cs


main 부분

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
def calcMandelbrot(type, start, size, unit, c=0 + 0j):
    """ Calculating and drawing mandelbrot set and julia set
    :param type: fractal type
    :param start: center point in plane
    :param size: image size(width, height)
    :param unit: gap among pixels
    :param c: julia constant
    :return: image 2D array (format: RGBA)
    """
 
    s = time()
 
    # change size to multiple of 16(bigger or equal to prior size). Because of threads per block is 16
    size = ((size + 255/ 256).astype(np.uint32) * 256
 
    result = np.empty((size[0], size[1]), np.int32)
 
    func = cu.get_function("calcMandelbrot")
    func(cuda.In(np.array([type == Fractal.JULIA], np.bool)),
         cuda.In(np.array([c.real, c.imag], np.float64)),
         cuda.In(np.array([start[0- size[0/ 2 * unit,
                           start[1- size[1/ 2 * unit, unit], np.float64)),
         cuda.Out(result),
         block=(16161), grid=(int(size[0/ 16), int(size[1/ 16)))
 
    # Because in image symmetric transformation occurs between x axis and y axis
    result = np.transpose(result)
 
    return array2imgarray(result, cc.m_cyclic_wrwbw_40_90_c42_s25)
cs


main부분을 보시면 16번째 줄에서 결과물을 복사해서 담을 CPU 메모리를 할당하고 있고, 18번째 줄에서 function을 가져와서 19번째 줄에서 호출을 합니다. 원래 C/C++ 같았으면 19번째 줄에서 파라미터 부분(cuda.In, cuda.Out)을 싸그리 할당하고 복사하는 복잡한 과정을 거쳐야 호출할 수 있는데, PyCUDA에서는 간단하게 GPU에 입력만 될 파라미터는 cuda.In으로, GPU에서 출력만 될 파라미터는 cuda.Out으로 처리만 하면 나머지 할당하고 복사하는 과정을 알아서 해줍니다. 만약에 GPU에 입력된 후 다시 출력될 파라미터 같은 경우는 cuda.InOut으로 해주면 간단하게 처리됩니다.


그런데 잘 보시면 func 호출 시 단순히 파라미터만 있는 것이 아니라 block이나 grid와 같은 것이 보입니다. 이것은 GPU에 얼마 만큼의 병렬 처리 단위를 생성해서 계산을 할 것인지를 정하는 것입니다. 제가 겪었던 첫번째 시행착오가 바로 이 부분이었죠.


CUDA에서 병렬 처리 단위는 block과 thread로 나눠집니다. 기본적으로 여러 개의 block이 있고, 이 block 안에 thread가 있는 것입니다. 그리고 각 thread는 연산을 할 때 자신의 위치를 지표로 삼아 데이터를 불러와서 연산을 합니다. 예를 들어 fractal 이미지 같은 경우 픽셀의 (x, y)값마다 thread가 있을테니 각 thread는 자신의 (x, y)에 대한 값을 연산합니다. CUDA에서는 이 (x, y)와 같은 thread index가 3차원까지 되어 있으며 각 차원 별로, 그리고 각 단위(block, thread) 별로 제한값이 있습니다. 이것은 GPU마다 다른데 특히 block 당 thread의 최대 개수는 최근 device 경우에는 1024개로 제한이 되어 있습니다. 이 제한 값은 Linux에서 CUDA 설치 시 제공되는 Sample 중 deviceQuery라는 것을 컴파일해서 실행해보면 알 수 있습니다. 제 경우(GTX 970)엔 아래와 같이 확인할 수 있습니다.


GTX 970 Device Query


이렇게 각 block과 thread 개수를 정해서 GPU에서 kernel을 실행하면 개수에 맞게 thread가 생성이 되어 kernel 함수가 연산을 수행합니다. 이때 자기 thread가 어디 위치에 있는 것인지 알아야 하기 때문에 CUDA에서는 기본적으로 gridDim, blockIdx, blockDim, threadIdx라는 상수를 제공합니다. 제가 이걸 프로그래밍하면서 당연하게 'blockIdx는 block의 위치값, blockDim은 block의 차원일테니, blockDim으로 block의 총 개수를 보면 되겠다!'라는 생각을 했었는데, 이게 아니더라고요. block의 총 개수: gridDim, thread의 총 개수: blockDim 입니다. blockIdx는 block의 index, threadIdx는 thread의 index이고요. 이 내용을 간단하게 그림으로 나타내면 아래와 같습니다. 3D는 그리기 귀찮아서, 2D로 그려봤습니다.


  암튼 이렇게 CUDA의 기본 구조를 파악하고 나서 간단하게 fractal 그리기 프로그램을 만들었습니다. 초기 버전은 단순하게 block 당 thread를 1개로만 해서 만들었습니다. 이때는 뭐 거의 800x800 그리는데 0.3초 정도로 매우 느렸습니다. 그래서 어떻게 해야 하나 싶어서 CUDA sample을 뒤져서 Mandelbrot 예제(옛날에 CUDA sample 뒤지다가 발견했었거든요 ㅇㅇ)를 찾아내서 코드를 비교해봤습니다. 전체적인 구조는 비슷하더라고요. 차이점을 찾아보니 일단 저는 단일 thread만 사용했고, 이상하게 CUDA sample은 loop를 일부 풀어놨습니다(다음 글에서 소개할 기술입니다). 아무래도 thread가 하나인 게 문제인 거 같아서 thread를 여러 개로 늘려봤더니, 속도가 100배 정도(...) 빨라졌습니다. 위에서 말했다시피 thread는 개수가 겨우 1024개밖에 안 되기 때문에 몇 개일 때 최적의 값을 가질지 한 번 검색을 해봤더니, 알고리즘마다 case by case인 것 같더라고요. 제 알고리즘에서 이를 계산해보니 아래와 같은 그래프가 나왔습니다. 네, 정말, thread 한 개 쓰는 건 무식한 짓이라는 것이 너무 자명하게 보이죠. (지금 다시 그래프를 그려보니 100배까진 차이가 안 나네요, 처음 알고리즘을 도대체 어떻게 짰었길래...)

사실상 42 thread 정도 돌리면 이후는 거의 엇비슷한 속도가 나옵니다. 그 이후는 돌릴 때마다 누가 빠른지 조금씩 다르더라고요. 그래서 저는 fractal 그리는 프로그램은 162 thread로 작성했습니다.


이렇게 간단히 fractal 그리는 프로그램을 만들면서 기본적인 최적화를 해봤습니다. 그리고 CUDA sample에 나와있는 의미심장한 loop를 풀어 둔 것을 검색하면서 또 다른 것을 배웠죠. 다음 글에서 다뤄볼 예정입니다.


오랜만에 글을 써 봅니다. 근래에 학업 공부를 하고 그러다 보니까 블로그에 공부한 걸 잘 안 올렸네요.(나름 자료를 찾아보고 틈틈히 보긴 했는데 쩝)


  2018년 겨울방학에 들어서고, 꽤 전부터 생각해오던 CUDA를 이용한 프렉탈 그리기를 직접 프로그래밍했습니다. 중학교 때 프렉탈을 처음 알게 되고 나서, 그 신비함과 아름다움에 빠져서 언젠가는 꼭 프렉탈을 그려내는 프로그래밍을 하겠다고 다짐을 했는데(그 옛날에 julia set의 상수마다 이미지를 만들어서 1TB 하드를 꽉 채우겠다는 생각도 해본 거 같습니다. 뭐 정작 여러 가지 그려보다 보니 이쁜 모양이 나오는 julia set 상수는 그렇게 많진 않았지만요)겨울방학 들어서고 하게 되었네요.


  요즘엔 안드로이드 프로그래밍을 거의 안 하고(앱 유지보수 정도...) 딥 러닝 등을 Python으로 하고 있습니다. 그래서 이번 프렉탈 그리기는 제가 옛날부터 주로 쓰던 Java AWT라던가 그런 걸로 만들지 않았고, Python 3에다가 PyQt5와 PyCUDA를 활용해서 만들어 봤습니다. 아래 이미지들이 제가 지금 제작한(그리고 제작하고 있는) 앱의 스크린샷입니다.

 

(GitHub 링크: https://github.com/AMakeApp/FractalViewer)


  이러한 프렉탈을 그리는데는 사실 CPU를 써도 크게 상관은 없습니다. 단지 GPU를 쓰는 이유는 프렉탈을 그리는데 들어가는 알고리즘이 상당한 연산을 필요로 하기 때문입니다. 이 연산 과정을 알려면 위의 프렉탈들이 어떻게 만들어지는지 알아볼 필요가 있습니다.


  제가 현재 앱에 구현해 둔 프렉탈은 위 두 이미지인 Mandelbrot Set과 Julia Set입니다. 사실 이 둘은 모양이 상이함에도 상당한 상관관계가 있습니다. 왜냐하면 이 둘의 집합을 계산하는 과정이 상당히 유사하기 때문이죠. 그렇기에 좀 더 많이 알려져 있는 Mandelbrot Set에 대해 알아봅시다. Mandelbrot Set은 이름에서 보이듯이 Mandelbrot가 고안한 집합입니다. 사실 많이 알려져 있듯이 저렇게 알록달록 해야지 Mandelbrot Set은 아니고 수학적으로는 아래 이미지처럼 복소평면 상에서 검정색 점들의 집합이 Mandelbrot Set입니다.

Mandelset hires


이 집합은 일반적으로 우리가 집합을 정의할 때와는 조금 다르게 특이하게 정의되어 있는데, 다음 점화식에서 수열이 발산하지 않게 하는 복소수 $c$의 집합으로 정의되어 있습니다.

$$z_0 = 0, z_{n + 1} = z_n^2 + c$$


즉, 이 집합은 무한 집합일테고 복소평면 상의 어느 한 점에 대해 이 점이 Mandelbrot Set에 속하는지 계산하려면 각 복소수에 대해 수많은 반복 연산을 통해 수열의 발산 여부를 확인해야 합니다. 다행히도 $| z_n | > 4$라면 발산한다고 볼 수 있기 때문에 일정한 반복 연산 후 $z_n$의 절대값을 비교하면 좀 더 수월하게 계산 할 수 있습니다.(아마 이 부분이 수학적으로 증명되어 있는 거 같습니다. 이 문서를 참고하세요) 하지만 그렇다고 해도 수많은 반복 연산이 필요한 것은 사실입니다. 초중고등학교에서 배우는 것처럼 단순히 함수꼴로 구할 수 있는 것도 아니고, 각 점에 대해 일일히 노가다를 해야하는 것과 마찬가지니까요. 각 점에 대해서는 제가 대충 해본 결과 200번 정도 이상하면 눈으로 보기에는 프렉탈이 충분히 이쁘게 나옵니다. 커다란 프렉탈 이미지를 생성하려고 2000 x 2000의 프렉탈 이미지를 생성한다고 합시다. 대략 이정도면 200 * 2000 * 2000 = 8억 번 정도의 연산이 필요로 합니다. 일반적으로 한 알고리즘에 대해서 1억 번 연산에 대한 소요시간을 1초로 잡을텐데(요즘 CPU는 더 빨라서 더 많이 연산을 할 수는 있습니다) 프렉탈 이미지를 계산하기 위해서는 8초나 걸린다는 겁니다. 요즘 CPU는 중사양 CPU도 4코어 8쓰레드 정도 가질텐데, 멀티 코어 프로그래밍을 한다고 해도 최소 1초는 걸린다는 이야기가 됩니다.


  이렇게 CPU는 프렉탈 이미지를 계산하는 데는 상당히 오래 걸리기 때문에 GPU의 수많은 코어를 사용해서 더 빠르게 계산을 할 수 있습니다. 비유를 하자면 CPU는 고급 인력 4~8명이 일하는 것이라면 GPU는 초보자 수 천 명이 일을 하는 것과 비슷합니다. 프렉탈 연산이 많은 반복 연산이 필요로 하지만 정작 그 연산 자체가 그렇게 복잡하지 않기 때문에(특히 분기 구문 같은 것이 많지 않기 때문에) 이러한 경우에는 GPU에서 빠르게 계산할 수 있습니다. 이러한 이유로 프렉탈을 빠르게 그려내기 위해서 GPU를 활용했습니다.


  GPU로 그려내기로 결정하고 나서, 제 컴퓨터에 장착되어 있는 GTX970에서 GPGPU 프로그래밍을 하기 위해서 Nvidia에서 지원하는 CUDA를 활용했습니다. 여기서 CUDA는 일반적으로 CUDA C를 사용해서 C/C++과 함께 사용되는데, GUI로 이미지를 보여주면서 편하게 프로그래밍하고 싶었기 때문에 Python 3를 기반으로 CUDA를 사용할 수 있게 만들어 둔 PyCUDA를 사용했습니다. 그리고 유려한 GUI를 보여주는 Qt를 Python에 접목시킨 PyQt5를 사용했고요.


다음 글부터는 이 Fractal을 그리는데 어떻게 CUDA 프로그래밍을 했는지 제 시행착오와 함께 CUDA에 대해 이야기해보겠습니다.

 음...딥러닝쪽 게시글을 오래 전부터 안 올리고 있었는데, 사실 좀 귀찮았습니다(ㅋㅋㅋㅋㅋㅋ)

  뭐 정리는 해야할 거 같지만, 데이터만 쌓아두고 여러 모델 구축해서 돌리면서 재미있게 놀았죠. 최근에 Mathematica를 구매해서 내장 딥러닝 모듈을 통해 CIFAR-10을 학습하여 76%인가?(30분 학습했는데, 이정도면 대단한 거겠죠? ㅋㅋㅋ) 나오기도 했고, 방학 때는 소논문으로 날씨 예측을 주제로 삼아 새롭게 Keras를 공부하여 패턴 예측을 시도해보기도 하였습니다. 뭐, 나중에 간단하게 모델들 업로드하면서 글로 정리하기로 하고... 오늘은 제가 학수고대한 NCS가 어제 DHL로 도착해서 간단히 리뷰를 해보려고 합니다!


단촐하지만 위엄있는 USB 스틱!


  NCS(Neural Compute Stick)은 아마 인텔의 자회사?로 추정되는 Movidius에서 개발한 딥러닝 inference용(즉, 학습은 불가능합니다) 컴퓨팅 스틱으로 저전력으로 딥러닝을 할 수 있도록 만들어진 것입니다. 79달러이며 저는 Mouser에서 9만원 정도로 구매하였습니다. 다른 리셀러는 모르겠는데, Mouser는 60달러 이상 무료 배송이고 어차피 79달러에서 할인을 하진 않을테니 그냥 여기서 구매했습니다. 아마 여기가 가장 낫지 않을까 싶네요. 지금은 주문이 많이 밀려서 그런지 꽤 걸리는 거 같고, 제가 주문하고 오기까지 한 1달 정도(배송은 1주일도 안 걸렸을 거예요) 걸렸네요.

  Myriad 2를 기반으로 하며 성능은 대략 100Gflops 정도이며, 전력을 겨우 1W만 먹는다고 합니다(!) 그래서인지 Movidius에서 지금 RPI에서 돌아가도록 SDK 배포판을 만들어서 제공하고 있으며, Movidius에선 RPI3를 예시로 보여주지만 제가 가지고 있는 RPI B(1버전)에서도 잘 돌아가더라고요. 아쉽게도 현재 지원 중인 딥러닝 프레임워크는 Caffe 밖에 없습니다. 예제로 들어있는 AlexNet, GoogLeNet, Gender, SqueezeNet을 돌려보면 각각 283ms, 569ms, 237ms, 304ms 정도 걸립니다(Inference Time 기준). 이미지 크기가 크지만 않다면 대역폭 문제가 크게 작용하지 않아서 그런지, 제 컴퓨터에서 USB 3.0으로 할 때나, RPI에서 USB 2.0(정확히는 2.0보다 안 좋습니다. RPI는 젠더로 USB 대역폭을 제공해서...)으로 할 때나 걸리는 시간은 거의 같습니다. 생각보다 딜레이 시간이 커도 각 딜레이 시간에 맞게 카메라에서 frame을 받아서 inference하면 뭐 그래도 쓸만한 거 같네요(기본 제공 예제에 있습니다). 그리고 생각보다(?) 79달러라는 싼 가격 때문인지 많이 구매해서 쓸 수 있게(?) 여러 개의 USB 스틱을 사용할 수 있도록 지원하는 거 같네요.


  지금 고민 중인 RPI Zero W를 구매하게 된다면, 카메라도 같이 구매하고, Movidius랑 같이 엮어서 Multiple Object Detection을 해보려고 생각하고 있습니다. 잘만 응용하면 시각 장애인에게 음성으로 자신 앞에 어떤 물체가 있는지 알려줄 수 있기 때문에(실제품이 나온다고 들었긴 했는데, 한 번 제 손으로 만들어보고 싶네요 ㅎㅎ) 한 번 만들어볼 생각을 하고 있습니다. 크기도 작고, 전력도 작게 먹는데다가, 무려 RPI를 지원하기 때문에 응용할 곳은 무궁무진할 거 같네요. 나중에 TensorFlow랑 Keras만 지원되면 정말 좋을 거 같습니다 ㅎㅎ

  중3때 삼각비라는 부분에서 sin, cos, tan을 배웁니다. 특수각인 0°, 30°, 45°, 60°, 90°에서의 sin, cos, tan값을 배우죠. 고2에 오면 이과 과정 미적분 II에서 삼각함수 단원에서 좀 더 자세히 배웁니다. 60분법 대신에 호도법을 사용하며 정의역도 실수 전체로 확장하죠. 삼각함수 덧셈, 뺄셈정리를 통해 15°, 75° 등의 각을 삼각함수의 특수각 값을 이용하여 구하죠. 하지만 현실에서 각도는 특수각의 합, 차, 1/2만 존재하지 않습니다. 6°, 34.5° 등의 일반적인 각에서 우리는 삼각함수의 값을 구할 필요가 있을 때가 있죠. 그러나 우리가 배운 내용으로는 이러한 값을 구할 수 없습니다. 그러나 컴퓨터는 이러한 값들을 아주 손쉽게 계산해냅니다. 어떻게 계산해내는 걸까요?


  그 과정을 알아보기 위하여 C언어의 math.h에서 사용하는 삼각함수 라이브러리(glibc)를 뜯어보았습니다. C언어에서는 math.h를 include하여 sin, cos, tan 등의 삼각함수를 사용할 수 있습니다. 이 링크에서 아래 코드의 원본을 보실 수 있습니다.

/* k_sinf.c -- float version of k_sin.c
 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
 */

/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================
 */

#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: k_sinf.c,v 1.4 1995/05/10 20:46:33 jtc Exp $";
#endif

#include <math.h>
#include <math_private.h>

static const float 
half =  5.0000000000e-01,/* 0x3f000000 */
S1  = -1.6666667163e-01, /* 0xbe2aaaab */
S2  =  8.3333337680e-03, /* 0x3c088889 */
S3  = -1.9841270114e-04, /* 0xb9500d01 */
S4  =  2.7557314297e-06, /* 0x3638ef1b */
S5  = -2.5050759689e-08, /* 0xb2d72f34 */
S6  =  1.5896910177e-10; /* 0x2f2ec9d3 */

float __kernel_sinf(float x, float y, int iy)
{
    float z,r,v;
    int32_t ix;
    GET_FLOAT_WORD(ix,x);
    ix &= 0x7fffffff;           /* high word of x */
    if(ix<0x32000000)           /* |x| < 2**-27 */
       {if((int)x==0) return x;}        /* generate inexact */
    z   =  x*x;
    v   =  z*x;
    r   =  S2+z*(S3+z*(S4+z*(S5+z*S6)));
    if(iy==0) return x+v*(S1+z*r);
    else      return x-((z*(half*y-v*r)-y)-v*S1);
}

  이 kernel_sin을 사용하여 실제의 sin함수를 구현하는 부분인 s_sinf.c 코드를 보면 $-{\pi \over 4} < x < {\pi \over 4}$ 에서는 간단히 iy = 0으로 정의되어 있을 것을 볼 수 있으니 iy = 0일 때 알고리즘을 살펴봅시다. 크기가 크지 않으니 z, v, r을 계산하여 x + v * (S1 + z * r)을 계산하여 넘기는 것으로 볼 수 있습니다. 이를 식으로 정리해보면 아래와 같습니다.

$$ \sin x = x + S1 x^3 + S2 x^5 + S3 x^7 + S4 x^9 + S5 x^{11} + S6 x^{13} $$


  이제 위에 하드코딩된 값을 대입하여 이 그래프를 그려보겠습니다. Legend에 나와있는 것처럼 f(x)가 위의 함수고, sin(x)는 Mathmatica 기본 내장 함수입니다.


사실상 일치해서 sin(x) 값이 다 덮어씌운 것만 보입니다. 이제 정의역을 f(x)에 정의된 것보다 좀 더 늘려볼까요?

정의역을 $-2\pi < x < 2\pi$ 까지 늘려봤습니다. 이제는 어느정도 차이가 두드러지게 나타나네요! 실제 sin 함수를 구현하는 s_sinf.c 코드를 보시면 $-{\pi \over 4} < x < {\pi \over 4}$ 이외에는 오류처리거나 reduction을 통해 표현하는 것을 볼 수 있습니다. 왜냐면 sin은 $0 \le x \le {\pi \over 4}$ 만 알면 정의역이 모든 실수일 때 정의할 수 있기 때문입니다.


  위와 같이 다항식으로 해석함수(고등학교 과정에서는 무한히 미분 가능한 함수 정도? 자세한 건 맨 아래에 있습니다)를 표현하는 것을 테일러 급수(Taylor series)라고 합니다. 이제 테일러 급수가 어떻게 미분 가능한 함수를 표현하는지 알아봅시다. 일단 미분 가능한 함수에 대한 n차 근사식을 $p_n(x)$ 라고 합시다. 그리고 우리는 한 점 $\left(x_0,\ f(x_0)\right)$ 과 그 점에 대한 미분계수( $f'(x)$ )가 주어졌을 때 그 점에서의 접선을 나타낼 수 있습니다 공식은 아래와 같죠.

$$ p_1(x) = f(x_0) + f'(x_0)(x - x_0)$$


이 공식에서 우리는 아래와 같은 성질을 찾을 수 있습니다.

$$ p_1(x_0) = f(x_0),\ p_1'(x_0) = f'(x_0) $$


따라서 이와 같은 성질을 이용하여 이번엔 2차 다항식에서 근사식을 만들어봅시다. 위의 성질을 사용하면 아래와 같은 조건을 만족하는 $p_2(x)$ 를 구해야겠죠.

$$ p_2(x_0) = f(x_0),\ p_2'(x_0) = f'(x_0),\ p_2''(x_0) = f''(x_0) $$


이에 대해 만족하는 식을 구하면 아래와 같습니다.

$$ p_2(x) = f(x_0) + f'(x_0)(x - x_0) + {f''(x_0) \over 2}(x - x_0)^2 $$


이와 같은 방식을 계속 반복하면 n차 미분한 식을 $f^n(x)$ 라고 했을 때 아래와 같은 n차 근사식을 구할 수 있습니다.

$$ p_n(x) = \sum^n_{k=0}{f^k(x - x_0)^k \over k!} $$


이때 아래와 같은 표현으로 $f(x)$ 를 나타내는데, 이 표현을 테일러 정리라고 하며 $R_n$ 을 나머지 항이라고 부릅니다.

$$ \begin{align*} f(x) &= p_n(x) + R_n \\&= f(x_0) + f'(x_0)(x - x_0) + {f''(x_0) \over 2}(x - x_0)^2 + \cdots + {f^{n + 1}(x_0) \over (n + 1)!}(x - x_0)^{n + 1} \end{align*} \\ R_n(x) = {f^{n + 1}(x_0) \over (n + 1)!}(x - x_0)^{n + 1}$$


나머지항은 일종의 테일러 정리의 식과 실제 $f(x)$ 의 식 간의 차이(오류)를 보여주므로, $R_n(x)$ 가 $n = \infty$ 일 때 0으로 수렴하면 테일러 정리의 식과 실제 식과 같아지므로 이때 우리는 위 무한급수를 테일러 급수라고 하며, 특별히 $x_0 = 0$ 인 경우에는 매클로린 급수라고 합니다.


이러한 방식으로 테일러 급수를 유도해낼 수 있는데, 수학적으로 엄밀하게 증명해보려면 미분법과 적분법을 사용하는 방법이 있습니다. 굳이 쓸 필요는 없어 보여서 이 링크를 참고하시기 바랍니다. 검색해도 잘 나와요.


자, 이제 대략적으로 테일러 급수에 대한 내용을 정리해보았습니다. 이제 그러면 적용을 해봅시다. 위에서 예시를 든 것처럼 $\sin x$ 에 대해서 테일러 급수를 써봅시다. 아래 그래프는 $\sin x$ 의 그래프(검정)와, $x_0 = 0$ 에서 테일러 급수의 1차(빨강), 3차(주황), 5차(노랑), 7차(초록), 9차(파랑), 11차(보라) 항까지 합한 것의 그래프입니다.

보이는 것처럼 합하는 항수가 커질수록 실제 $\sin x$ 에 다가가는 것을 볼 수 있습니다.


이제 다음으로 파데 근사(Padé approximant)에 대해 써보려고 합니다. 어떤 함수를 다항식으로 표현한 것이 테일러 급수라면, 어떤 함수를 유리식으로 표현한 것이 파데 근사입니다. 테일러 급수의 일반화이며 상위호환이라고 하는 거 같네요. 음이 아닌 정수 $m,\ n$에 대하여 파데 근사의 꼴은 아래와 같습니다. 이때 $n = 0$ 이 되면 매클로린 급수와 같게 됩니다.

$$ R(x)_{m/n} = {\sum_{j=0}^m{a_jx^j} \over 1 + \sum_{i=1}^n{b_kx^k}} $$


...파데 근사 계산을 한 번 해서 정리해보고 싶었는데 내용 이해가 잘 안 되는 관계로(쿨럭) 링크는 남겨둘테니 궁금하신 분들은 여기서 살펴보세요. 계산법을 제대로 이해하긴 힘들거 같고, 대신에 Mathematica에서 PadeApproximant를 함수로 제공하기 때문에 이것으로 $x_0=0$ 일 때 계산한 $\sin x$ 를 아래 표로 정리해보았습니다.

 m\n

1(빨강)

3(주황)

5(초록)

7(파랑)

 0(0)

$x$

$x-\frac{x^3}{6}$

 $x-\frac{x^3}{6}+\frac{x^5}{120}$

$x-\frac{x^3}{6}+\frac{x^5}{120}-\frac{x^7}{5040}$

 2(0.3)

 $\dfrac{x}{1 + \frac{x^2}{6}}$

 $\dfrac{x - \frac{7x^3}{60}}{1 + \frac{x^2}{20}}$

$\dfrac{x - \frac{x^3}{7} + \frac{11x^5}{2520}}{1 + \frac{x^2}{42}}$

 $\dfrac{x - \frac{11x^3}{72} + \frac{13x^5}{2160} - \frac{x^7}{12096}}{1 + \frac{x^2}{72}}$

 4(0.7)

  $\dfrac{x}{1 + \frac{x^2}{6} + \frac{7x^4}{360}}$

$\dfrac{x - \frac{31x^3}{294}}{1 + \frac{3x^2}{49}+\frac{11x^4}{5880}}$

 $\dfrac{x - \frac{53x^3}{396} + \frac{551x^5}{166320}}{1 + \frac{13x^2}{396} + \frac{5x^4}{11088}}$

 $\dfrac{x - \frac{241x^3}{1650} + \frac{601x^5}{118800} - \frac{121x^7}{2268000}}{1 + \frac{17x^2}{825} + \frac{19x^4}{118800}}$


그래프로 나타내면 아래와 같습니다. 조금 굵은 검정 선은 $\sin x$ 그래프이고, n 옆에 써 있는 색상에 m 옆에 써 있는 어두움을 가진 선이 그 (m, n)의 파데 근사 그래프입니다.

(무슨 그래픽 기교하는 줄)

이렇게 컴퓨터는 함수를 근사하여 계산합니다. 나중에 기회가 되면 파데 근사를 좀 더 조사해서 올려보도록 하겠습니다.


p.s 참고 자료: Introduction to Numerical Analysis - Alastair Wood, 위키백과-테일러 급수, 네이버 지식백과, 네이버캐스트, 위키백과-파데 근사, WolframMathWorld-Padé Approximant

  학교에서 본 미적분 영상에서 최단강하곡선 문제에 대한 이야기를 봐서 이번 글에선 이 내용을 다뤄보려고 합니다. 최단강하곡선은 2차원 상에 주어진 두 점을 잇는 곡선 중에서 두 점 중 높은 곳에 물체를 올려뒀을 때 가장 빨리 낮은 곳으로 내려가는 곡선입니다. 역사적으로 보면 1689년에 요한 베르누이가 유럽의 최고 수학자들에게 보냈다는 문제로 유명하죠. 사실 이 문제는 베르누이가 뉴턴을 시험해보려던 문제였는데, 뉴턴은 퇴근 후 몇 시간 끄적거려서 풀어내고 익명으로 풀이를 보냈죠. 이름 없는 풀이를 본 베르누이는 뉴턴인 것을 알아봤다고 합니다.



(이거 만드느라 2시간 쓴 건 함정) (검정: 선분, 파란색: 사이클로이드, 빨간색: 호)


  위 그림이 최단강하곡선 문제를 시각적으로 나타낸 것입니다. $(0, 1)$에서 $(1, 0)$까지 이은 여러 곡선(곡선의 정의에 따라 선분도 포함됩니다) 중에서 간단히 호, 선분, 그리고 이 글의 주인공인 사이클로이드를 그려봤습니다. 언핏 봤을 땐 직선으로 내려가는게 가장 빠르다고 생각할 수도 있습니다. 하지만 이 직관이라는 것은 생각보다 신뢰도가 떨어지고, 무엇보다 수학은 직관으로 시작할 수는 있어도 직관으로만 결과를 도출해내는 것은 대단히 위험하기 때문에 이제 과연 어떤 곡선이 최단강하곡선인지 유도해낼 것입니다.


  일단 이 최단강하곡선을 유도해내기 전에 이 문제를 통해 새롭게 탄생한 변분법이라는 것을 알아야합니다. 변분법이란 범함수(함수 집합을 정의역으로 하는 함수)의 최대 또는 최소를 구해내는 방법입니다. 이 변분법을 쓰는 예로는 여기서 소개하는 사이클로이드의 유도, 평면에 있는 두 점을 잇는 곡선 유도 등이 있습니다. 이제 저는 변분법에서 쓰이는 오일러-라그랑주 방정식을 통해 이 문제를 해결할 겁니다.




  먼저 중력의 방향을 위 이미지처럼 $+y$ 방향이라고 잡고 어떤 두 점 $A(0, 0),\ B(x_0, y_0)\ (x_0, y_0 > 0)$ (계산의 편의를 위해 간단히 잡았습니다)에서 물체를 놨을 때 최단시간으로 떨어지는 곡선을 구해야하기 때문에 결과값이 시간인 범함수를 만들어야 합니다. 이때 중학교 때부터 알고 있던 간단한 공식을 이용합니다.  $ t = {s \over v}$ 이죠. 이때 우리는 곡선을 다루기 때문에 범함수 $L$에 대해 아래와 같이 공식을 세울 수 있습니다.

$$ L = \int_0^{x_0} {ds \over v} $$


이때 이 곡선에서는 역학적 에너지가 보존되기 때문에 아래와 같이 $v$ 를 구할 수 있습니다.

$$ mgy = {1 \over 2}mv^2 \\[10pt] \therefore v = \sqrt{2gy} $$


$ds$ 는 고등학교 기하와 벡터 과정에서 구하는 곡선의 길이 공식을 통해 유도해낼 수 있습니다.

$$ ds = \sqrt{ dx^2 + dy^2 } = \sqrt{1 + x'^2 }\ dy $$


즉, 범함수 $L$은 아래와 같이 나타낼 수 있습니다.

$$ L = \int_0^{x_0} {\sqrt{1 + x'^2 \over 2gy}\ dy} = {1 \over \sqrt{2g}} \int_0^{x_0} {\sqrt{1 + x'^2 \over y}\ dy} $$


이제 오일러-라그랑주 방정식을 봅시다. 범함수 $J$에 대해서 다음과 같이 나옵니다. 이 방정식의 유도과정은 위키백과에 검색하면 잘 나오니 궁금하시면 검색하시면 됩니다.

$$ J = \int_a^b F(x, f(x), f'(x))\ dx,\ {\partial F \over \partial f} - {d \over dx} {\partial F \over \partial f'} = 0 $$


우리의 범함수 $L$에 대해 변형하면 아래와 같이 나옵니다.

$$ F(y, x, x') = \sqrt{1 + x'^2 \over y},\ {\partial F \over \partial x} - {d \over dy} {\partial F \over \partial x'} = 0 $$


이제 쭉 풀어봅시다.

$$ {\partial F \over \partial x} - {d \over dy} {\partial F \over \partial x'} = 0 \\[15pt] {d \over dy} {\partial F \over \partial x'} = 0 \left(\because {\partial F \over \partial x} = 0\right) \\[15pt] {d \over dy} \left( \sqrt{y \over 1 + x'^2} \cdot {x' \over y} \right) = 0 \\[15pt] \sqrt{y \over 1 + x'^2} \cdot {x' \over y} = C \\[15pt] x' = \pm\sqrt{C^2 y \over 1 - C^2 y} $$


$x'$ 꼴로 나왔으니 이제 적분을 해주면 됩니다. 하지만 루트 안에 $y$ 가 있기에 그냥 적분하기엔 불가능할 정도이므로 치환적분을 써줍니다. 중간에 나오는 $\sin^2 t$는 삼각함수 반각공식(삼각함수 덧셈정리 응용)을 통해 풀어줍니다.

$$ y = {1 \over C^2} \sin^2 t,\ dy = {2 \over C^2} \sin t \cos t\ dt\\[15pt] \begin{align*} x &= \pm \int \sqrt{C^2 y \over 1 - C^2 y}\ dy \\[15pt] &= \pm {2 \over C^2} \int \sin^2 t \ dt \\[15pt] &= \pm {1 \over C^2} \int 1 - \cos 2t \ dt \\[15pt] &= \pm {1 \over C^2}\left(t - {1 \over 2}\sin 2t\right) + C'\end{align*}$$


이때 위에 점 $A$ 를 대입해보면 $C' = 0$ 을 얻을 수 있으며 위에 $y$ 의 $\pm$ 은 양수이므로 사라집니다. 좀 더 깔끔하게 정리하기 위해서 $x$ , $y$ 를 정리해보면 아래와 같습니다.

$$ y = {1 \over C^2} \sin^2 t = {1 \over 2C^2} \left(1 - \cos 2t\right) \\[15pt] 2t\text{를 } \theta\text{로 변환시키면} \\[15pt] \begin{cases} \displaystyle x = {1 \over 2C^2} \left(\theta - \sin \theta\right) \\ \displaystyle y = {1 \over 2C^2} \left(1 - \cos \theta \right)\end{cases}$$


여기까지가 일반적으로 알려져있는 사이클로이드 곡선의 식이며 더 나아가 $A$ , $B$ 를 잇는 사이클로이드 곡선의 식을 유도해보겠습니다. 중력의 방향이 $+y$ 방향으로 정해져있기 때문에 $x$축 대칭을 하고 원래 좌표평면에 맞게 좌표를 바꾼 점 $B(a, b)$를 대입해서 정리하면 됩니다.

$$ \begin{cases} \displaystyle x = {1 \over 2C^2} \left(\theta - \sin \theta\right) \\ \displaystyle y = -{1 \over 2C^2} \left(1 - \cos \theta \right)\end{cases} \\[15pt] x\text{축 대칭한 점 } B(a, b) \text{를 대입하고 이때 성립하는 } \theta \text{를 } \theta_0 \text{라 두면} \\[15pt] \begin{cases} \displaystyle a = {1 \over 2C^2} \left(\theta_0 - \sin \theta_0\right) \\ \displaystyle b = -{1 \over 2C^2} \left(1 - \cos \theta _0 \right)\end{cases} \\[15pt] \text{이때의 } \theta_0\text{을 구하는 방정식은 두 식을 나눠서 정리하면 아래와 같다} \\[15pt] {\theta_0 - \sin \theta_0 \over 1 - \cos \theta_0} = -{a \over b} \quad \left(\text{단, }0 < \theta_0 < 2\pi \right) \\[15pt] \text{따라서 이를 정리하면} \\[15pt] \begin{cases} \displaystyle x = -{b \over 1 - \cos \theta_0} \left(\theta - \sin \theta\right) \\ \displaystyle y = {b \over 1 - \cos \theta_0} \left(1 - \cos \theta \right)\quad \left(\text{단, }0 \leq \theta \leq \theta_0 \right) \end{cases}$$


$\theta_0$  같은 경우엔 일반적 풀이방법이 없기 때문에 수치적으로 풀면 됩니다. 저 같은 경우엔 Mathematica 11의 FindRoot 기능을 통해 구한 후 맨 위 이미지를 그려냈습니다.


... 힘드네요 ㅠㅠ(장장 7시간 정도 걸쳐서 쓴....)


  아무튼 이렇게 사이클로이드 곡선이 유도가 됩니다. 당연히 구한대로 실제로 실험해봐도 사이클로이드 곡선에서 물체를 굴릴 때가 가장 빠르게 도착지점에 도착합니다.


참고 사이트

https://ko.wikipedia.org/wiki/%EB%B3%80%EB%B6%84%EB%B2%95

http://mathworld.wolfram.com/BrachistochroneProblem.html

http://zetablog.tistory.com/33

http://mathsci.kaist.ac.kr/~nipl/am621/lecturenotes/Euler-Lagrange_equation.pdf

'수학, 과학 > 수학 이모저모' 카테고리의 다른 글

최단강하곡선 문제(Brachistochrone Problem)  (0) 2017.06.18
피보나치 수열 일반항 유도  (0) 2017.06.04
회전변환  (0) 2017.06.03
MathJax 사용하기  (0) 2017.05.28

  어느날 갑자기 피보나치 수열의 일반항을 구하는 것이 궁금해서 찾아두고, 여기에 제가 살을 몇 개 붙인 내용입니다. 일단 피보나치 수열의 점화식은 아래와 같습니다.

$$ \Large a_{n + 2} = a_{n + 1} + a_n,\ a_1 = 1,\ a_2 = 1 $$


  여기서 두 가지의 방법이 있습니다. 하나는 수열의 특성 방정식을 이용한 풀이와, 또 다른 하나는 고등학교 과정에서(09년 개정 교육과정 이전 점화식 내용) 나오는 트릭을 이용하는 풀이가 있습니다. 순서대로 두 풀이 모두 써보겠습니다.


  일단 먼저 특성 방정식을 사용하는 풀이를 보여드리겠습니다. 이상하게 이 특성 방정식을 이용한 풀이는 특성 방정식이 정확히 무엇인지를 설명하진 않더라고요...ㅇㅂㅇ... 나름 제대로 찾아보려고 해도 잘 안 나와서 위키백과 및 주워들은 것으로 설명을 해보겠습니다. 예를 들어 $a_{n + 2} = pa_{n + 1} + qa_{n}$라는 피보나치 수열과 형태가 같은 어떤 점화식이 있다고 봅시다. 이 점화식은 특성 방정식 $x^2 = px + q$의 해 $\alpha, \beta$를 통해 $a_n = P\alpha^n + Q\beta^n$의 형태로 나타낼 수 있다고 합니다. 한 번 증명을 어느 블로그(클릭!)를 참고해서 써보자면 아래와 같습니다.


일단 특성 방정식의 근과 계수와의 관계를 이용해 아래와 같이 나타낼 수 있습니다.

$$ p = \left(\alpha + \beta\right),\ q = -\alpha\beta $$


따라서 원래의 점화식을 아래와 같이 정리할 수 있습니다.


$$\displaystyle{ i)\ \alpha - \beta \ne 0 \\[10pt] a_{n + 2} = \left(\alpha + \beta\right)a_{n + 1} - \alpha\beta a_n \\[10pt] a_{n + 2} - \alpha a_{n + 1} = \beta\left(a_{n + 1} - \alpha a_n\right),\ a_{n + 2} - \beta a_{n + 1} = \alpha\left(a_{n + 1} - \beta a_n\right) \\[10pt] a_{n + 2} - \alpha a_{n + 1} = \left(a_2 - \alpha a_1\right) \beta^{n - 1}\ \cdots\ A,\ a_{n + 2} - \beta a_{n + 1} = \left(a_2 - \beta a_1\right) \alpha^{n - 1}\ \cdots\ B \\[10pt] B - A\text{ 를 하면},\ \left(\alpha - \beta\right)a_{n + 1} = \left(a_2 - \beta a_1\right) \alpha^{n - 1} - \left(a_2 - \alpha a_1\right) \beta^{n - 1} \\[10pt] \therefore a_n = { a_2 - \beta a_1 \over \alpha - \beta } \alpha^{n - 1} - { a_2 - \alpha a_1 \over \alpha - \beta }\beta^{n - 1} \\[10pt]~\\[10pt] ii)\ \alpha - \beta = 0 \\[10pt] a_{n + 2} = 2 \alpha a_{n + 1} - \alpha^2a_n \\[10pt] a_{n + 2} - \alpha a_{n + 1} = \alpha\left(a_{n + 1} - \alpha a_{n}\right) \\[10pt] a_{n + 1} = \alpha a_n + \left(a_2 - \alpha a_1\right) \alpha^{n - 1} \\[10pt] \therefore a_n = a_1 \alpha^{n - 1} + n\left(a_2 - \alpha a_1\right) \alpha^{n - 2}} $$


이렇게 증명할 수 있고, 이제 이것을 이용해 단번에 피보나치 수열의 점화식을 일반항으로 나타내봅시다.


$$ \text{특성 방정식을 구해보면},\ x^2 - x - 1 = 0 \\[10pt] \therefore \alpha = {1 + \sqrt{5} \over 2},\ \beta = { 1 - \sqrt{5} \over 2} \\[10pt] \alpha \ne \beta \text{ 이므로} \\[10pt] \begin{align*} a_n &= { a_2 - \beta a_1 \over \alpha - \beta } \alpha^{n - 1} - { a_2 - \alpha a_1 \over \alpha - \beta }\beta^{n - 1} \\[10pt] &= { 1 - \beta \over \alpha - \beta } \alpha^{n - 1} - { 1 - \alpha \over \alpha - \beta }\beta^{n - 1} \\[10pt] &= {1 \over \sqrt{5}} \left(\left(1 + \sqrt{5} \over 2\right)^n - \left(1 - \sqrt{5} \over 2\right)^n\right) \end{align*}$$


  그 다음으로 고등학교 과정에 나오는 트릭을 이용해 일반항을 유도해보겠습니다. (사실 특성 방정식 풀이랑 거의 비슷하게 나옵니다) 위의 피보나지 수열 점화식에서 $a_{n + 1}$의 일부를 좌측으로 이항하여 $b_n$ 꼴을 만들어 낼 수 있도록 합니다. 그러면 이제 아래와 같이 결국 특성 방정식에서 나오는 꼴의 점화식이 되기 때문에 이후 과정은 같게 나오고 일반항이 유도가 됩니다.


$$ a_{n + 2} = a_{n + 1} + a_n \\[10pt] a_{n + 2} - \alpha a_{n + 1} = \beta \left(a_{n + 1} + {1 \over \beta} a_n\right)\ (\alpha + \beta = 1) \\[10pt] -\alpha = {1 \over \beta} \\[10pt] \therefore a_{n + 2} = \left(\alpha + \beta\right)a_{n + 1} - \alpha\beta a_n$$


아무튼 이렇게 피보나치 수열의 일반항을 유도해 낼 수 있습니다.


p.s 고등학교 과정 풀이를 유도하다가 아예 $\alpha$ 만 쓰고 풀이를 쭉 전개해봤는데 막 $\alpha^4$이 나오고 그러더라고요 허허.... 위에 풀이 방식대로 하는게 정신건강에 좋습니다.


'수학, 과학 > 수학 이모저모' 카테고리의 다른 글

최단강하곡선 문제(Brachistochrone Problem)  (0) 2017.06.18
피보나치 수열 일반항 유도  (0) 2017.06.04
회전변환  (0) 2017.06.03
MathJax 사용하기  (0) 2017.05.28

  최근에 일품 기하와 백터를 풀다가 삼각함수로 일일히 좌표를 구해서 벡터를 직접 구해서 푸는 문제가 나오더라고요. 이때 선생님께서 회전변환을 알려주셔서 간단하게 정리해보려고 합니다.



  위와 같이 $x^2 + y^2 = r^2$ 인 원이 있고 원 위의 한 점을 $B$ 라고 잡습니다. 그리고 점 $B$ 에서 $x$ 축으로 수선의 발을 내린 것을 점 $C$ 라고 잡고, $\overline{BC} // \overline{AD}$, $\overline{BC} = \overline{AD}$ 인 점 $D$ 를 잡습니다. 이때 $\Box ACBD$ 가 직사각형이 되기 때문에 $\overrightarrow{AB} = \overrightarrow{AC} + \overrightarrow{AD}$ 가 됩니다. $B = \left(x, y\right)$ 라 하면, $C = \left(x, 0\right), D = \left(0, y\right)$ 입니다.


  자 이제, 점 $A$ 를 기준으로 점 $B$ , 점 $C$ , 점 $D$ 를 $\theta$ 만큼 시계 반대 방향으로 회전시킵니다. 회전시킨 후 점들은 각각 점 $B'$ , 점 $C'$ , 점 $D'$ 라고 합시다. 이때 $C' = \left(x \cos \theta, x \sin \theta \right), D' = \left(y \cos \left({\pi \over 2} + \theta \right), y \sin \left({\pi \over 2} + \theta \right) \right)$ 입니다. $\overrightarrow{AC'} + \overrightarrow{AD'} = \overrightarrow{AB'}$ 이므로

$$ \begin{align*} B' &= \left(x \cos \theta + y \cos \left({\pi \over 2} + \theta \right), x \sin \theta + y \sin \left({\pi \over 2} + \theta \right) \right) \\ &= \left( x \cos \theta - y \sin \theta, x \sin \theta + y \cos \theta \right) \end{align*} $$

이 됩니다.


  따라서 어떤 점의 좌표를 $\left(x, y\right)$ , 이를 $\theta$ 만큼 회전변환시킨 좌표를 $\left(x', y'\right)$ 라 하고, 이를 행렬로 정리해보면,

$$ \begin{pmatrix} x' \\ y' \end{pmatrix} = \begin{pmatrix} \cos \theta & - \sin \theta \\ \sin \theta & \cos \theta \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} $$

가 됩니다.


... 2차원에서 이렇게 증명(?)한 것을 토대로, 3차원으로 회전변환을 확대시키고 싶긴 한데 제가 기하를 못해서... 당장 머릿속으로는 유도하기가 쉽지가 않네요. 나중에 가능하면 해보겠습니다.


'수학, 과학 > 수학 이모저모' 카테고리의 다른 글

최단강하곡선 문제(Brachistochrone Problem)  (0) 2017.06.18
피보나치 수열 일반항 유도  (0) 2017.06.04
회전변환  (0) 2017.06.03
MathJax 사용하기  (0) 2017.05.28

원래는 이거 말고 다른 주제로 글을 쓸려고 했는데...MathJax가 저만 말썽인지, 다른 블로그에서 소개한 게 틀린 건지 제대로 안 되는 경우가 있어서 MathJax 사용법을 간단히 올려보려고 합니다.


일단 Tistroy 경우에는 수식을 입력하기 위해서 아래와 같이 수식입력 탭이 있습니다.


실제로 이를 이용하여 수식을 작성하면 이미지로 들어가기 때문에 아래와 같이 조금 못나게 수식이 보입니다. 게다가 이미지이기 때문에 확대하면 계단현상을 볼 수 있죠.


하지만 MathJax를 사용하여 이쁘게 수식을 랜더링하면 아래와 같이 이쁘게 나옵니다. 게다가 맘에 드는 랜더링 설정을 아래 수식을 우클릭해서 고를 수도 있습니다.


$$ \zeta \left( s \right) = \sum_{n=1}^{\infty} {{1 \over {n^s}}} $$


또한 이미지로 글자 사이에 수식을 넣을 때는 아래와 같이 별로 이쁘지 않습니다. 위치 조정하기도 귀찮을 뿐더러, 위치 조정을 해도 다른 글자들과 어울리지 않습니다. (저 같은 경우엔 이미지여서 뭔가 이질감이 느껴집니다.)


이 3 이상의 정수일 때, 을 만족하는 양의 정수 는 존재하지 않는다.


하지만 이것도 MathJax로 작성한다면?!?!


$n$이 3 이상의 정수일 때, $ x^n + y^n = z^n $을 만족하는 양의 정수 $x, y, z$는 존재하지 않는다.


크, 아릅답네요.


이렇게 좋은 MathJax를 쓰려면 아래 코드를 Tistory에서 관리-HTML/CSS 편집에 들어가셔서 HTML 코드에서 <head> 태그 밑에 넣어주시면 됩니다.


<script type="text/javascript" async
  src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-MML-AM_CHTML">
    MathJax.Hub.Config({tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]}});
</script>


그리고 수식은 $$ TeX문법 $$(한 줄 전체를 차지하는 수식) 또는 $ TeX문법 $(다른 글자들과 같이 사용할 때)로 넣으시면 됩니다.



'수학, 과학 > 수학 이모저모' 카테고리의 다른 글

최단강하곡선 문제(Brachistochrone Problem)  (0) 2017.06.18
피보나치 수열 일반항 유도  (0) 2017.06.04
회전변환  (0) 2017.06.03
MathJax 사용하기  (0) 2017.05.28

2017년 3월 22일에 나온 아직도 따끈따끈한 논문입니다. GitHub링크는 아래와 같습니다.

https://github.com/luanfujun/deep-photo-styletransfer

간단히 논문에 대해 요약을 하자면 아래와 같습니다.


기존까지 있었던 image style transfer(유명한 neural-style 등)은 이미지에 스타일을 입히는 게 얼추 가능은 했지만, style faithfulness가 약하거나 photorealism이 약했습니다. photorealism이 약한 이유는 사진 내부에 있는 직선 성분과 택스쳐 성분이 왜곡되기 때문이고, style faithfulness가 약한 이유는 예를 들어서 sky style은 sky에 매칭이 되야 하고, building style은 building에 매칭되어야 하는데 이게 잘 안돼서 연관없는 것이 매칭되었기 때문입니다.  따라서 Deep Photo Style Transfer는 이것들을 해결하고자 뭐 여러가지 적용했습니다. 그래서 이전의 알고리즘과는 달리 위 두 가지를 충족시키는 좋은 알고리즘이 되었다 뭐 이정도입니다. 그래서 실제로 이걸 돌려보면 사진의 전체적인 분위기(시각, 계절, 날씨, 택스쳐 재배치 등)이 제대로 됩니다.


아래가 그 예시입니다. (원본을 GitHub에서 링크 따온 것입니다.) 왼쪽부터 순서대로 원본, 스타일, 결과물입니다.



코드가 GtiHub에 있기 때문에 돌려볼려고 clone도 받고 Octave도 설치해서 라플라시안 돌리고 그러려고 했는데... GPU 램 부족하다고 안 되네요 쩝... 나중에 좋은 장비 구축하면 그때가서 돌려보겠습니다.


p.s. 논문은 대충 읽어봤는데 처음해보는 리뷰라 그런지 잘 안 되네요 쿨럭...  나중에 시간이 되면 더 꼼꼼히 읽어서 리뷰해보겠습니다

  1. 스터디중 2017.10.24 12:52 신고

    안녕하세요. 리뷰 잘 읽었습니다.
    한가지 질문을 드려도 될까요?
    segmentation mask에 맞게 transfer 되는건 알겠는데, 그럼 다른 mask 는 transfer에 전혀 관여를 안하는건가요?

    • makeapp 2017.10.28 00:12 신고

      간단하게 논문 조금 훑어보고 예제를 살펴본 정도라서, 저도 정확하게 답변드릴순 없을 거 같습니다. 직접 논문 읽어보시고, 프로그램 실행해보시는게 좋을 거 같습니다.

+ Recent posts

티스토리 툴바