할 일

파이썬에서 라플라시안을 포함한 미분 방정식을 시뮬레이션하고 싶었지만, 마지막 날에 라플라시안의 차이를 for 문 이외로 처리하는 방법을 생각해 냈습니다 (바퀴의 재발명을 한) 때문에 확산 방정식의 애니메이션을 썼습니다. 했다.

생각나고 나서 다시 조사해 보니 몇 년 전에 자신이 한 일의 상위 호환을 하고 있는 쪽이 계셨습니다. 그 기사는 상당히 도움이되었습니다. NumPy · SciPy를 이용한 수치 계산의 고속화 : 응용 1

풀고있는 미분 방정식은 다음 식으로 표현되는 1 차원 확산 방정식입니다.
$$
\frac{\partial }{\partial t}\psi(x,t)= D\frac{\partial^2}{\partial x^2}\psi(x,t)
$$
초기 조건은 가우시안으로 주었고, 시간 발전은 우선 오일러법으로 했습니다.

작성한 코드

Diffusion_1D_ani.ipynb

%matplotlib nbagg
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import ArtistAnimation

# calculate Laplacian with periodic boundary condition
def Lap(L, M):
    # expand vector
    VL = [L[0]]
    VR = [L[-1]]
    Le = np.hstack([VR ,L, VL])
    Lc = np.dot(M, Le) # central diff
    return Lc[1:N+1]

# generate 1-D lattice and initialize
N = 64
x = np.linspace(-10, 10, N)
psi = np.exp(-x**2/7)

# time and diffusion const
dt = 0.05
T = np.arange(0.0, 37, dt)
D = 3

# make (N+2)x(N+2) matrix to calc Laplacian
I = np.eye(N+2,N+2)
S = np.vstack((I[1:],np.array([0]*(N+2))))
M = S + S.T -2*I 

# make animation
fig = plt.figure()
ims = []
for i in range(len(T)):
    psi += dt*(D*Lap(psi, M))
    line = plt.plot(x, psi, 'b')
    ims.append(line)

ani = ArtistAnimation(fig, ims, interval=5, blit=True, repeat=True)

plt.show()

실행 예

궁리한 점

라플라시안 계산을 할 때 원래의 $N$차원 배열 $\psi$의 우단과 좌단의 값을 각각 좌단과 우단에 np.hstack로 붙여 $N+2$차원 배열로 했습니다. 이를 통해 미리 준비된 $ (N + 2)\times (N + 2) $ 행렬 $ M $과 내적을 취하면 자동으로주기 경계 조건을 포함하는 확산 항을 계산할 수있었습니다. 경계 조건용으로 추가한 부분은 불필요하므로 함수의 끝에 버리고 있습니다.

총괄

과제점 : 우선 보이는 형태로 하고 싶었으므로, 수치적인 안정성이나 오차는 1mm도 고려하고 있지 않습니다. NumPy에 빠르고 안전하게 수치 적분 해 주는 무언가가 있다고는 생각하기 때문에 조사해 어떻게든 하려고 합니다. 이러한 과제에 대해 앞으로 어떤 진척이 태어나면 추기합니다.

향후 : 이번 1차원 버전을 적당히 확장하면 2차원 버전도 만들 수 있을 것 같기 때문에, 반응 확산계의 시뮬레이션을 해보고 싶다고 생각하고 있습니다. (어쩐지 에그인 비선형항을 포함한 방정식이 아니면 어떻게든 되겠지...)

추가:

(11/23) 전진 차분으로 타협하면, 3행으로 비선형항을 넣을 수 있었습니다. 일단 1차원의 Navier-Stokes 방정식도 할 수 있을 것 같습니다.

def nl(L, N):
    Lf = np.hstack((L,L[0]))
    return L * np.diff(Lf)

전진 차분이라면 레이놀즈 수뿐만 아니라 초기 조건이나 시간의 너비도 신중하게 선택하지 않으면 즉시 계산이 발산했습니다. 자, 어떻게 할까요?