행렬과 벡터 연산 (1) - 확산 방정식 예제

이번 시리즈 물에서는 순수 파이썬 구현과 numpy 구현의 속도 차이가 얼마만큼 나고 이러한 차이가 왜 발생하는지 알아보려 합니다. 순수 파이썬으로 문제를 구현한 이후 지난 포스트에서 살펴본 line profiler 를 이용해 연산 속도를 점진적으로 개선시킵니다. 이후 벡터 연산을 지원하는 numpy 배열을 이용해 구현하고 얼마만큼의 속도 개선을 이룰 수 있는지 보려고 합니다. 이를 살펴보기 위한 예제로는 메모리와 연산량을 매우 많이 잡아먹는 적절한 예제가 필요하겠죠. 바로 확산 방정식입니다.

확산 방정식

확산은 유체가 공간에서 퍼지는 양상을 나타낸 수식입니다. 예를 들어 실온의 물에 물감을 몇 방울 떨어뜨리면 물감이 천천히 퍼지면서 물과 완전히 섞일 때까지 움직이는데, 물을 젓는 등의 외부 유인이 없다면 두 액체가 섞이는 과정은 확산으로 이루어집니다. 따라서 시간과 공간에 따른 유체의 양으로 표현할 수 있는데, 1차원 상에서는 다음 수식과 같습니다.

$\frac{d}{dt}u(x,t)=D\cdot \frac{d^2}{dx^2}u(x,t)$

이 공식에서 $u$는 확산하는 유체의 양을 나타내고 이 값이 0이라면 물만 존재하고 1이라면 물감만 존재합니다. $D$는 확산 계수로 값이 클수록 쉽게 확산되는데, 여기에서는 1로 고정합니다. 이 수식은 유체의 시간에 따른 변화량 (왼쪽)과 유체의 공간에 따른 2차 변화량을 확산 계수로 이어주는데, 오일러 방법을 이용해 도함수를 차로 기술하여 불연속적인 공간과 시간의 근삿값으로 유도할 수 있습니다.

$\frac{d}{dt}u(x,t) \approx \frac{u(x, t+dt) - u(x, t)}{dt}$

$dt$는 시간의 단위로서 시간의 정밀도를 나타냅니다. 당연히 $dt$가 줄어들수록 원래의 확산 방정식을 매우 정확하게 근사할 수 있을 겁니다. 오른쪽의 공간 확산 방정식을 오일러 방법으로 풀어쓰고 $u(x, t)$가 주어졌을 때 $u(x, t+dt)$를 계산할 수 있도록 고쳐쓰면 다음과 같습니다. $dx$는 공간 상의 해상도를 나타내고 이 값이 작아질수록 행렬의 셀 크기도 점점 작아집니다. 이 값 또한 1로 고정합니다.

$u(x, t+dt)=u(x,t)+dt\cdot D \cdot \frac{u(x+dx, t)+u(x-dx, t) - 2\cdot u(x,t)}{dx^2}$

위 식을 이용해 문제를 풀 때 먼저 $x$ 값에 따른 경계 조건을 설정해주어야 합니다. $x$는 행렬의 색인 (공간 상에서의 위치) 이므로 이 값이 valid 한 행렬 색인 아닐 시에는 이에 맞는 조건을 걸어주어야 합니다. 여기서는 행렬 한 차원의 크기가 $N$일 때 그 차원에서 색인이 -1 인 값은 $N-1$ 값이고 $N$ 위치의 값은 색인이 0인 주기적 경계 조건을 사용합니다. 또한, 시간의 변화에 따른 $u$ 값 또한 계산하려는 시점마다 현재 상태와 다음 상태를 나타내는 최소 두 개가 필요합니다.

 

Psuedo code

위 수식을 바탕으로 한 의사 코드는 다음과 같습니다.

## 초기 조건 생성
u = vector of length N
for i in range(N):
    u = 0  # water
    
D = 1
t = 0
dt = 0.001
while True:
    unew = vector of length N
    
    ## 행렬 값 갱신 based on 확산 방정식
    for i in range(N):
        unew[i] = u[i] + D * dt * (u[(i+1)%N] + u[(i-1)%N) - 2*u[i])
    ## 새로운 u
    u = unew

 

2차원으로의 확장

2차원 / 3차원으로의 확산 방정식 확장도 간단합니다. 앞서 $u(x, t)$에서 $x$가 1차원 행렬을 나타내는 색인이었다면 2차원을 나타내는 $y$ 매개 변수를 이런 식으로 $u(x, y, t)$ 추가하면 되겠죠. 2차원 확산 방정식은 다음과 같습니다.

$\frac{d}{dt}u(x, y, t)=D \cdot (\frac{d^2}{dx^2}u(x, y, t) + \frac{d^2}{dy^2}u(x, y, t))$

이를 코드로 나타내면 다음과 같습니다.

for i in range(N):
    for j in range(M):
        unew[i][j] = u[i][j] + dt * D * (
            (u[(i+1)%N][j] + u[(i-1)%N][j] - 2*u[i][j]) + 
            (u[i][(j+1)%N] + u[i][(j-1)%N] - 2*u[i][j])
        )

 

---

행렬과 벡터 연산 (2) - 확산 방정식 순수 파이썬 구현

행렬과 벡터 연산 (3) - 파이썬 리스트와 numpy 연산 속도 차이

행렬과 벡터 연산 (4) - numpy 배열을 이용한 확산 방정식

행렬과 벡터 연산 (5) - numpy 배열 메모리 최적화

행렬과 벡터 연산 (6) - numexpr 모듈 이용하기