06_finite_difference.py

  1"""
  2유한 차분법 (Finite Difference Method)
  3Finite Difference Methods for PDEs
  4
  5편미분방정식(PDE)을 수치적으로 푸는 방법입니다.
  6"""
  7
  8import numpy as np
  9import matplotlib.pyplot as plt
 10from typing import Callable, Tuple
 11
 12
 13# =============================================================================
 14# 1. 1D 열 방정식 (Heat Equation)
 15# =============================================================================
 16def heat_equation_explicit(
 17    L: float,
 18    T: float,
 19    nx: int,
 20    nt: int,
 21    alpha: float,
 22    initial_condition: Callable[[np.ndarray], np.ndarray],
 23    boundary_left: float = 0,
 24    boundary_right: float = 0
 25) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
 26    """
 27    1D 열 방정식 (명시적 방법)
 28
 29    ∂u/∂t = α ∂²u/∂x²
 30
 31    FTCS (Forward Time, Central Space):
 32    u(i,n+1) = u(i,n) + r[u(i+1,n) - 2u(i,n) + u(i-1,n)]
 33    여기서 r = α*dt/dx²
 34
 35    안정성 조건: r ≤ 0.5
 36
 37    Args:
 38        L: 공간 영역 [0, L]
 39        T: 시간 영역 [0, T]
 40        nx: 공간 격자 수
 41        nt: 시간 스텝 수
 42        alpha: 열확산 계수
 43        initial_condition: 초기 조건 함수 u(x, 0)
 44        boundary_left, boundary_right: 경계 조건
 45
 46    Returns:
 47        (x 배열, t 배열, u 배열)
 48    """
 49    dx = L / (nx - 1)
 50    dt = T / nt
 51    r = alpha * dt / dx**2
 52
 53    print(f"dx={dx:.4f}, dt={dt:.6f}, r={r:.4f}")
 54    if r > 0.5:
 55        print(f"경고: 안정성 조건 위반 (r={r:.4f} > 0.5)")
 56
 57    x = np.linspace(0, L, nx)
 58    t = np.linspace(0, T, nt + 1)
 59    u = np.zeros((nt + 1, nx))
 60
 61    # 초기 조건
 62    u[0, :] = initial_condition(x)
 63
 64    # 경계 조건
 65    u[:, 0] = boundary_left
 66    u[:, -1] = boundary_right
 67
 68    # 시간 전진
 69    for n in range(nt):
 70        for i in range(1, nx - 1):
 71            u[n + 1, i] = u[n, i] + r * (u[n, i + 1] - 2 * u[n, i] + u[n, i - 1])
 72
 73    return x, t, u
 74
 75
 76def heat_equation_implicit(
 77    L: float,
 78    T: float,
 79    nx: int,
 80    nt: int,
 81    alpha: float,
 82    initial_condition: Callable[[np.ndarray], np.ndarray],
 83    boundary_left: float = 0,
 84    boundary_right: float = 0
 85) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
 86    """
 87    1D 열 방정식 (암시적 방법 - Crank-Nicolson)
 88
 89    무조건 안정, O(dt², dx²) 정확도
 90    """
 91    dx = L / (nx - 1)
 92    dt = T / nt
 93    r = alpha * dt / (2 * dx**2)
 94
 95    x = np.linspace(0, L, nx)
 96    t = np.linspace(0, T, nt + 1)
 97    u = np.zeros((nt + 1, nx))
 98
 99    # 초기 조건
100    u[0, :] = initial_condition(x)
101    u[:, 0] = boundary_left
102    u[:, -1] = boundary_right
103
104    # 삼대각 행렬 설정
105    n_inner = nx - 2
106    A = np.zeros((n_inner, n_inner))
107    B = np.zeros((n_inner, n_inner))
108
109    for i in range(n_inner):
110        A[i, i] = 1 + 2 * r
111        B[i, i] = 1 - 2 * r
112        if i > 0:
113            A[i, i - 1] = -r
114            B[i, i - 1] = r
115        if i < n_inner - 1:
116            A[i, i + 1] = -r
117            B[i, i + 1] = r
118
119    # 시간 전진
120    for n in range(nt):
121        # 우변 계산
122        b = B @ u[n, 1:-1]
123        b[0] += r * (u[n + 1, 0] + u[n, 0])
124        b[-1] += r * (u[n + 1, -1] + u[n, -1])
125
126        # 선형 시스템 풀기
127        u[n + 1, 1:-1] = np.linalg.solve(A, b)
128
129    return x, t, u
130
131
132# =============================================================================
133# 2. 1D 파동 방정식 (Wave Equation)
134# =============================================================================
135def wave_equation(
136    L: float,
137    T: float,
138    nx: int,
139    nt: int,
140    c: float,
141    initial_displacement: Callable[[np.ndarray], np.ndarray],
142    initial_velocity: Callable[[np.ndarray], np.ndarray] = None
143) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
144    """
145    1D 파동 방정식
146
147    ∂²u/∂t² = c² ∂²u/∂x²
148
149    유한 차분:
150    u(i,n+1) = 2u(i,n) - u(i,n-1) + s²[u(i+1,n) - 2u(i,n) + u(i-1,n)]
151    여기서 s = c*dt/dx (Courant number)
152
153    안정성 조건: s ≤ 1 (CFL 조건)
154    """
155    dx = L / (nx - 1)
156    dt = T / nt
157    s = c * dt / dx
158
159    print(f"Courant number s={s:.4f}")
160    if s > 1:
161        print(f"경고: CFL 조건 위반 (s={s:.4f} > 1)")
162
163    x = np.linspace(0, L, nx)
164    t = np.linspace(0, T, nt + 1)
165    u = np.zeros((nt + 1, nx))
166
167    # 초기 조건
168    u[0, :] = initial_displacement(x)
169
170    # 첫 번째 시간 스텝 (초기 속도 사용)
171    if initial_velocity is None:
172        initial_velocity = lambda x: np.zeros_like(x)
173
174    v0 = initial_velocity(x)
175    for i in range(1, nx - 1):
176        u[1, i] = (u[0, i] + dt * v0[i] +
177                   0.5 * s**2 * (u[0, i + 1] - 2 * u[0, i] + u[0, i - 1]))
178
179    # 경계 조건 (고정)
180    u[:, 0] = 0
181    u[:, -1] = 0
182
183    # 시간 전진
184    for n in range(1, nt):
185        for i in range(1, nx - 1):
186            u[n + 1, i] = (2 * u[n, i] - u[n - 1, i] +
187                          s**2 * (u[n, i + 1] - 2 * u[n, i] + u[n, i - 1]))
188
189    return x, t, u
190
191
192# =============================================================================
193# 3. 2D 라플라스/포아송 방정식
194# =============================================================================
195def laplace_2d(
196    Lx: float,
197    Ly: float,
198    nx: int,
199    ny: int,
200    boundary_conditions: dict,
201    max_iter: int = 10000,
202    tol: float = 1e-6
203) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
204    """
205    2D 라플라스 방정식 (Jacobi 반복법)
206
207    ∇²u = 0  또는  ∂²u/∂x² + ∂²u/∂y² = 0
208
209    Jacobi 반복:
210    u(i,j)_new = 0.25 * [u(i+1,j) + u(i-1,j) + u(i,j+1) + u(i,j-1)]
211
212    Args:
213        boundary_conditions: {'top': val, 'bottom': val, 'left': val, 'right': val}
214    """
215    dx = Lx / (nx - 1)
216    dy = Ly / (ny - 1)
217
218    x = np.linspace(0, Lx, nx)
219    y = np.linspace(0, Ly, ny)
220
221    u = np.zeros((ny, nx))
222
223    # 경계 조건 설정
224    bc = boundary_conditions
225    if callable(bc.get('top')):
226        u[-1, :] = bc['top'](x)
227    else:
228        u[-1, :] = bc.get('top', 0)
229
230    if callable(bc.get('bottom')):
231        u[0, :] = bc['bottom'](x)
232    else:
233        u[0, :] = bc.get('bottom', 0)
234
235    if callable(bc.get('left')):
236        u[:, 0] = bc['left'](y)
237    else:
238        u[:, 0] = bc.get('left', 0)
239
240    if callable(bc.get('right')):
241        u[:, -1] = bc['right'](y)
242    else:
243        u[:, -1] = bc.get('right', 0)
244
245    # Jacobi 반복
246    for iteration in range(max_iter):
247        u_old = u.copy()
248
249        for i in range(1, ny - 1):
250            for j in range(1, nx - 1):
251                u[i, j] = 0.25 * (u_old[i + 1, j] + u_old[i - 1, j] +
252                                  u_old[i, j + 1] + u_old[i, j - 1])
253
254        # 수렴 체크
255        error = np.max(np.abs(u - u_old))
256        if error < tol:
257            print(f"수렴: {iteration + 1}회 반복, 오차={error:.2e}")
258            break
259
260    return x, y, u
261
262
263def laplace_2d_sor(
264    Lx: float,
265    Ly: float,
266    nx: int,
267    ny: int,
268    boundary_conditions: dict,
269    omega: float = 1.5,
270    max_iter: int = 10000,
271    tol: float = 1e-6
272) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
273    """
274    2D 라플라스 방정식 (SOR - Successive Over-Relaxation)
275
276    Jacobi보다 빠른 수렴
277
278    Args:
279        omega: 이완 인자 (1 < ω < 2 for 가속)
280    """
281    x = np.linspace(0, Lx, nx)
282    y = np.linspace(0, Ly, ny)
283    u = np.zeros((ny, nx))
284
285    bc = boundary_conditions
286    u[-1, :] = bc.get('top', 0) if not callable(bc.get('top')) else bc['top'](x)
287    u[0, :] = bc.get('bottom', 0) if not callable(bc.get('bottom')) else bc['bottom'](x)
288    u[:, 0] = bc.get('left', 0) if not callable(bc.get('left')) else bc['left'](y)
289    u[:, -1] = bc.get('right', 0) if not callable(bc.get('right')) else bc['right'](y)
290
291    for iteration in range(max_iter):
292        max_diff = 0
293
294        for i in range(1, ny - 1):
295            for j in range(1, nx - 1):
296                u_old = u[i, j]
297                u_gs = 0.25 * (u[i + 1, j] + u[i - 1, j] +
298                              u[i, j + 1] + u[i, j - 1])
299                u[i, j] = (1 - omega) * u_old + omega * u_gs
300                max_diff = max(max_diff, abs(u[i, j] - u_old))
301
302        if max_diff < tol:
303            print(f"SOR 수렴: {iteration + 1}회 반복")
304            break
305
306    return x, y, u
307
308
309# =============================================================================
310# 시각화
311# =============================================================================
312def plot_heat_equation():
313    """열 방정식 시각화"""
314    # 초기 조건: 가운데 뜨거운 부분
315    initial = lambda x: np.sin(np.pi * x)
316
317    x, t, u = heat_equation_explicit(
318        L=1.0, T=0.5, nx=51, nt=500, alpha=0.01,
319        initial_condition=initial
320    )
321
322    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
323
324    # 시공간 등고선
325    ax = axes[0]
326    X, T_mesh = np.meshgrid(x, t)
327    contour = ax.contourf(X, T_mesh, u, levels=20, cmap='hot')
328    plt.colorbar(contour, ax=ax, label='온도')
329    ax.set_xlabel('위치 x')
330    ax.set_ylabel('시간 t')
331    ax.set_title('열 방정식 해')
332
333    # 시간별 프로파일
334    ax = axes[1]
335    times_to_plot = [0, len(t)//4, len(t)//2, 3*len(t)//4, -1]
336    for idx in times_to_plot:
337        ax.plot(x, u[idx, :], label=f't={t[idx]:.3f}')
338    ax.set_xlabel('위치 x')
339    ax.set_ylabel('온도 u')
340    ax.set_title('시간별 온도 분포')
341    ax.legend()
342    ax.grid(True)
343
344    plt.tight_layout()
345    plt.savefig('/opt/projects/01_Personal/03_Study/Numerical_Simulation/examples/heat_equation.png', dpi=150)
346    plt.close()
347    print("그래프 저장: heat_equation.png")
348
349
350def plot_wave_equation():
351    """파동 방정식 시각화"""
352    # 초기 변위: 가우시안 펄스
353    initial = lambda x: np.exp(-100 * (x - 0.5)**2)
354
355    x, t, u = wave_equation(
356        L=1.0, T=2.0, nx=101, nt=400, c=1.0,
357        initial_displacement=initial
358    )
359
360    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
361
362    # 시공간 등고선
363    ax = axes[0]
364    X, T_mesh = np.meshgrid(x, t)
365    contour = ax.contourf(X, T_mesh, u, levels=20, cmap='RdBu')
366    plt.colorbar(contour, ax=ax, label='변위')
367    ax.set_xlabel('위치 x')
368    ax.set_ylabel('시간 t')
369    ax.set_title('파동 방정식 해')
370
371    # 스냅샷
372    ax = axes[1]
373    for idx in [0, 50, 100, 150, 200]:
374        ax.plot(x, u[idx, :], label=f't={t[idx]:.2f}')
375    ax.set_xlabel('위치 x')
376    ax.set_ylabel('변위 u')
377    ax.set_title('시간별 파동 형태')
378    ax.legend()
379    ax.grid(True)
380
381    plt.tight_layout()
382    plt.savefig('/opt/projects/01_Personal/03_Study/Numerical_Simulation/examples/wave_equation.png', dpi=150)
383    plt.close()
384    print("그래프 저장: wave_equation.png")
385
386
387def plot_laplace_2d():
388    """2D 라플라스 방정식 시각화"""
389    bc = {
390        'top': lambda x: np.sin(np.pi * x),
391        'bottom': 0,
392        'left': 0,
393        'right': 0
394    }
395
396    x, y, u = laplace_2d(
397        Lx=1.0, Ly=1.0, nx=51, ny=51,
398        boundary_conditions=bc,
399        tol=1e-5
400    )
401
402    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
403
404    X, Y = np.meshgrid(x, y)
405
406    # 등고선
407    ax = axes[0]
408    contour = ax.contourf(X, Y, u, levels=20, cmap='viridis')
409    plt.colorbar(contour, ax=ax)
410    ax.set_xlabel('x')
411    ax.set_ylabel('y')
412    ax.set_title('라플라스 방정식 해')
413    ax.set_aspect('equal')
414
415    # 3D 표면
416    ax = axes[1]
417    ax = fig.add_subplot(1, 2, 2, projection='3d')
418    ax.plot_surface(X, Y, u, cmap='viridis', alpha=0.8)
419    ax.set_xlabel('x')
420    ax.set_ylabel('y')
421    ax.set_zlabel('u')
422    ax.set_title('3D 표면')
423
424    plt.tight_layout()
425    plt.savefig('/opt/projects/01_Personal/03_Study/Numerical_Simulation/examples/laplace_2d.png', dpi=150)
426    plt.close()
427    print("그래프 저장: laplace_2d.png")
428
429
430# =============================================================================
431# 테스트
432# =============================================================================
433def main():
434    print("=" * 60)
435    print("유한 차분법 (Finite Difference Method)")
436    print("=" * 60)
437
438    # 1. 열 방정식
439    print("\n[1] 1D 열 방정식")
440    print("-" * 40)
441
442    initial = lambda x: np.sin(np.pi * x)
443
444    print("명시적 방법 (FTCS):")
445    x, t, u_explicit = heat_equation_explicit(
446        L=1.0, T=0.1, nx=21, nt=100, alpha=0.1,
447        initial_condition=initial
448    )
449    print(f"  t=0.1에서 중앙 온도: {u_explicit[-1, 10]:.6f}")
450
451    print("\n암시적 방법 (Crank-Nicolson):")
452    x, t, u_implicit = heat_equation_implicit(
453        L=1.0, T=0.1, nx=21, nt=100, alpha=0.1,
454        initial_condition=initial
455    )
456    print(f"  t=0.1에서 중앙 온도: {u_implicit[-1, 10]:.6f}")
457
458    # 해석해: u = exp(-π²αt) sin(πx)
459    exact = np.exp(-np.pi**2 * 0.1 * 0.1) * np.sin(np.pi * 0.5)
460    print(f"  해석해: {exact:.6f}")
461
462    # 2. 파동 방정식
463    print("\n[2] 1D 파동 방정식")
464    print("-" * 40)
465
466    initial_wave = lambda x: np.sin(np.pi * x)
467
468    x, t, u_wave = wave_equation(
469        L=1.0, T=2.0, nx=51, nt=200, c=1.0,
470        initial_displacement=initial_wave
471    )
472    print(f"주기 T = 2L/c = 2.0")
473    print(f"t=2.0에서 중앙 변위: {u_wave[-1, 25]:.6f}")
474    print(f"(초기값과 같아야 함: {initial_wave(0.5):.6f})")
475
476    # 3. 2D 라플라스 방정식
477    print("\n[3] 2D 라플라스 방정식")
478    print("-" * 40)
479
480    bc = {'top': 100, 'bottom': 0, 'left': 0, 'right': 0}
481
482    print("Jacobi 반복:")
483    x, y, u_jacobi = laplace_2d(1.0, 1.0, 31, 31, bc, tol=1e-4)
484
485    print("\nSOR (ω=1.5):")
486    x, y, u_sor = laplace_2d_sor(1.0, 1.0, 31, 31, bc, omega=1.5, tol=1e-4)
487
488    print(f"\n중심점 온도: {u_jacobi[15, 15]:.4f}")
489
490    # 시각화
491    try:
492        plot_heat_equation()
493        plot_wave_equation()
494        plot_laplace_2d()
495    except Exception as e:
496        print(f"그래프 생성 실패: {e}")
497
498    print("\n" + "=" * 60)
499    print("유한 차분법 정리")
500    print("=" * 60)
501    print("""
502    PDE 유형과 방법:
503
504    | PDE 유형    | 대표 예      | 권장 방법              |
505    |------------|-------------|----------------------|
506    | 포물선형    | 열 방정식    | FTCS, Crank-Nicolson |
507    | 쌍곡선형    | 파동 방정식  | 중심차분, Lax-Wendroff|
508    | 타원형      | 라플라스     | Jacobi, GS, SOR      |
509
510    안정성 조건:
511    - 열 방정식 (명시적): r = αΔt/Δx² ≤ 0.5
512    - 파동 방정식: CFL = cΔt/Δx ≤ 1
513
514    차분 근사:
515    - 전진 차분: ∂u/∂t ≈ [u(t+Δt) - u(t)] / Δt
516    - 후진 차분: ∂u/∂t ≈ [u(t) - u(t-Δt)] / Δt
517    - 중심 차분: ∂²u/∂x² ≈ [u(x+Δx) - 2u(x) + u(x-Δx)] / Δx²
518
519    실무:
520    - scipy.ndimage: 간단한 필터링
521    - FEniCS, FiPy: Python PDE 프레임워크
522    - OpenFOAM: CFD (유한 체적법)
523    """)
524
525
526if __name__ == "__main__":
527    main()