확산-반응 방정식을 로그 공간으로 변환시키는 과정과 실제 적용

# 확산-반응 방정식 로그 공간(pH) 전개 과정

### 1. 기본 확산-반응 방정식 (농도 기준)
\[
\frac{\partial [H^+]}{\partial t} = D \nabla^2 [H^+] - R([H^+])
\]
- $[H^+] = [H^+](\mathbf{x}, t)$ : $H^+$ 농도  
- $D$ : 확산계수
- $R([H^+])$ : 반응항

### 2. 로그 변수 정의
\[
pH = -\log_{10} [H^+] \implies [H^+] = 10^{-pH}
\]

### 3. 시간 미분 변환
\[
\frac{\partial [H^+]}{\partial t} = \frac{\partial}{\partial t} 10^{-pH} = -10^{-pH} \ln10 \cdot \frac{\partial pH}{\partial t} = -[H^+]\ln10 \cdot \frac{\partial pH}{\partial t}
\]
따라서,
\[
\frac{\partial pH}{\partial t} = -\frac{1}{[H^+]\ln10} \frac{\partial [H^+]}{\partial t}
\]

### 4. 공간 미분 변환 (라플라시안)
\[
\nabla [H^+] = \nabla 10^{-pH} = -10^{-pH} \ln10 \cdot \nabla pH = -[H^+] \ln 10 \nabla pH
\]
\[
\nabla^2 [H^+] = \nabla \cdot \nabla [H^+] = \nabla \left( -[H^+] \ln 10 \nabla pH \right) =  [H^+]\ln10 \left[ \ln10(\nabla pH)^2 - \nabla^2 pH \right]
\]
여기서,
\[
(\nabla pH)^2 = \left(\frac{\partial pH}{\partial x}\right)^2 + \left(\frac{\partial pH}{\partial y}\right)^2
\]

### 5. 원래 방정식에 대입
\[
\frac{\partial [H^+]}{\partial t} = D \nabla^2 [H^+] - R([H^+])
\]
\[
-[H^+] \ln10 \frac{\partial pH}{\partial t} = D [H^+] \ln10\left[ \ln10 (\nabla pH)^2 - \nabla^2 pH \right] - R([H^+])
\]

### 6. 정리
\[
\boxed{
\frac{\partial pH}{\partial t} = D \left[ \nabla^2 pH - \ln10 (\nabla pH)^2 \right] + \frac{R([H^+])}{[H^+] \ln10}
}
\]
또는
\[
\boxed{
\frac{\partial pH}{\partial t} = D \left[ \nabla^2 pH - \ln10 (\nabla pH)^2 \right] + \frac{R(10^{-pH})}{10^{-pH} \ln10}
}
\]

### 7. 해석적 의미 및 구현 포인트
- $pH$의 시간 변화는 라플라시안 $\nabla^2 pH$와 그라디언트 제곱 $(\nabla pH)^2$에 의존  
- 반응항은 농도 기반 $R([H^+])$ 형태로 계산 후 변환  
- 수치 미분 시 $\nabla pH, \nabla^2 pH$를 각각 구해 사용

 

이상 여기까지가 GPT를 기반으로 쓴 확산-반응 방정식의 로그 공간에서의 전개다.

그러면 이걸 어디다가 쓰느냐?

시뮬레이션을 제작할 때 우리는 몰농도를 1e-7, 1e-8, 1e-9 등등으로 두고는 한다. 0.0000001 0.00000001 0.000000001 이 녀석들 말이다.

이제는 모두가 아는 상식, 컴퓨터는 실수를 부동소수점으로 저장하기 때문에 소수점이 많아지면 많아질 수록 오차율이 굉장히 크다.

순서대로 test1.py test2.py test3.py pH_diffusion.gif

그렇지 않고서야 내가 이렇게 많은 시행 착오를 거칠 일이 없었겠지. 저주한다 부동소수점!!!!!!111

특히 이 녀석이 문제다.

# 반응-확산 식 계산
dC_H = (D_H * L_H - R) * dt

C_H는 $[H^+]$, D_H는 $H^+$의 확산 계수 $D$, L_H는 $\nabla [H^+]$로 생각을 하면 되겠다. 즉 저기서 $\partial [H^+]$를 계산해주고 있는 셈이다. 근데 문제는 저기서 저녀석이 0이 되어버리면 Kw / C_H로 구하는 C_OH는 무한대로 발산해버려서 전체 용액이 그냥 강강강산이 되어버린다는 거지 아하하하하하

왜 0이 되는가????????당연히 부동소수점이다. 24개의 비트로 소수부분을 저장하는 것은 미친짓이다!

그러면 여기서 어떤 해법을 취해줄 수 있는가? 다음과 같이 정리하는 것이다.

 

어떠한 실수 $R$을 $R=r\cdot10^{-n}$과 같이 나타내고, $R$을 저장하는 대신 $r$과 $n$을 저장한다.

 

그렇게 하면 뭐가 낫느냐? 우선 저장하는 것이 매우 작은 실수가 아니라 적당한 크기의 실수 $r$, 정수 $n$을 저장하는 것이기 때문에 부동소수점 오차가 거의 없어진다. 특히 수소이온 농도같은 경우는 $pH=-\log_{10}[H^+]$의 관계에 있고, 이를 다시 풀어쓰면 $[H^+]=10^{-pH}$에 있기 때문에, 지수인 pH를 저장해 이에 대해 연산을 진행하면 된다.

그리고 위에서 $[H^+]$로 나타내어진 방정식을 $pH$에 대해 굳이굳이 다시 나타낸 이유 또한 그것이다.

 

그래서 실제 적용하면 어떻게 바뀌게 되는가?

원래는 용액의 pH를 1e-8로만 해도 뻑이 가는 녀석을 r=1, e=11로 설정해도 계산에 아무 지장이 없다!
dtstep=10만 되어도 zerodivide 에러가 뜨던 게 완벽히 고쳐진다!

확산이 쥐꼬리만큼 되던 것이 아주 잘 퍼진다! 이게 확산이지!

최근 코딩에 대한 지식 없이도 AI에게 맡겨서 모든 코딩을 해버리는 걸 바이브 코딩이라고 부르는 유행이 있는 모양이다

하지만 모두들 간단한 기술적 지식만 응용해도 세상은 아름다워진다! 모두들 기술을 공부하자!

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from tqdm import tqdm
from numba import njit

nr, nz = 100, 200
dr, dz = 0.001, 0.001
dt = 1e-4  # 시간 간격

# 확산계수 근사값, pH 단위 간 확산 정도 (실제로 확산은 농도기준임)
D_pH = 1e-2  

# 반응 속도 상수 (pH 간 차이 반응)
k = 1e-2  

Kw_pH = 14.0  # 물의 중성 pH 기준

# 초기 pH (중성)
n = np.full((nr, nz), 11.0)

# 경계 조건 (중성 고정)
n[0, :] = 7.0
n[-1, :] = 7.0
n[:, 0] = 7.0
n[:, -1] = 7.0

# 산성 영역 설정 (pH 낮게 설정)
n[nr//2-10:nr//2+10, nz-120:nz] = 2.0

@njit
def laplacian_pH(pH, dr, dz, nr, nz):
    L = np.zeros_like(pH)
    for i in range(1, nr-1):
        for j in range(1, nz-1):
            r = i * dr
            if r == 0:
                r = dr / 2
            L[i, j] = (
                (pH[i+1, j] - 2*pH[i, j] + pH[i-1, j]) / dr**2 +
                (1/r) * (pH[i+1, j] - pH[i-1, j]) / (2*dr) +
                (pH[i, j+1] - 2*pH[i, j] + pH[i, j-1]) / dz**2
            )
    return L

frames = []
steps = 20000
max_delta_pH = 0.02  # 한 스텝에 pH 변화 제한

for step in tqdm(range(steps)):
    L_n = laplacian_pH(n, dr, dz, nr, nz)
    
    # 반응 항: pH가 물의 중성 pH에서 얼마나 벗어났는지에 비례
    # 산성(pH < 7)은 중성으로 회복하려는 경향, 염기성(pH > 7)도 마찬가지
    R = -k * (n - 7.0)  # 단순 회복 반응
    
    dn_dt = D_pH * L_n + R
    n_new = n + dn_dt * dt
    
    # 변화 제한
    delta = n_new - n
    delta = np.clip(delta, -max_delta_pH, max_delta_pH)
    n += delta
    
    # 경계 조건 유지
    n[0, :] = 7.0
    n[-1, :] = 7.0
    n[:, 0] = 7.0
    n[:, -1] = 7.0
    
    if step % 10 == 0:
        frames.append(n.copy())

# 애니메이션 그리기
fig, ax = plt.subplots()
im = ax.imshow(frames[0].T, origin='lower', extent=[0, nr*dr, 0, nz*dz], vmin=0, vmax=14, cmap='plasma')
plt.colorbar(im, label='pH')
ax.set_title("pH Distribution over Time")
ax.set_xlabel("r (m)")
ax.set_ylabel("z (m)")

def update(frame_idx):
    im.set_array(frames[frame_idx].T)
    ax.set_title(f"pH Distribution at t = {frame_idx*dt*10:.5f} s")
    return [im]

ani = FuncAnimation(fig, update, frames=len(frames), interval=50)
plt.show()