러스트로 만드는 수학 및 과학 계산 라이브러리 시리즈 - 6편: 편미분방정식(PDE) 기초와 유한차분법(Finite Difference Method)

앞선 글에서 우리는 선형대수, FFT, 적분·미분방정식(ODE) 풀이까지 구현하며 러스트 기반 수학/과학 계산 라이브러리를 점차 확장해 왔습니다. 이제는 한 단계 더 나아가, 편미분방정식(PDE, Partial Differential Equation)의 기초 해법에 도전해보겠습니다. PDE는 열 방정식, 파동 방정식, 라플라스 방정식 등 물리, 공학, 금융수학까지 다양한 분야에 등장하며, 이를 수치적으로 해결하기 위해 유한차분법(FDM, Finite Difference Method) 같은 기초적인 방법을 활용할 수 있습니다.

이번 글에서는 다음과 같은 주제를 다룹니다.

1. 1차원 열방정식(Heat Equation) 예제를 통해 PDE 개념 소개
2. 유한차분법(FDM): 격자점에 대해 미분 연산을 근사하여 PDE를 시간 스텝별로 풀어가는 기본 아이디어
3. 간단한 예: 1D 열방정식 u_t = α u_xx를 Euler 전진과 중앙이차차분 방식으로 근사하고, 시간 발전을 계산해보는 예제를 구현

PDE 해석 개요

1D 열방정식 예:

u_t = α u_xx,   0 < x < L, t > 0
u(x,0) = f(x)   초기조건
u(0,t)=u(L,t)=0  경계조건(Dirichlet) 가정

여기서 u_t는 시간에 대한 편미분, u_xx는 공간에 대한 2차 편미분.

유한차분법의 기본 아이디어:

• 공간을 N개의 격자로 나누어 x_i = iΔx (i=0..N)
• 시간도 스텝 단위로 t_n = nΔt 분할
• u(x_i, t_n)를 u_i^n이라 표기
• 2차 미분 근사: u_xx ≈ (u_{i+1}^n - 2u_i^n + u_{i-1}^n) / (Δx²)
• 시간 전진: u_i^{n+1} = u_i^n + Δt * α * ( (u_{i+1}^n - 2u_i^n + u_{i-1}^n)/Δx² )

이 간단한 전진 시간, 중앙 공간 차분(FTCS) 기법은 안정성 조건(αΔt/Δx² ≤ 1/2) 등이 필요하지만, 개념적 입문 예제로 충분합니다.

구현 예제

// src/pde.rs

/// 1D 열방정식: u_t = alpha * u_xx
/// 경계조건: u(0,t)=0, u(L,t)=0
/// 초기조건: u(x,0)=f(x)
///
/// 파라미터:
/// f: 초기 조건 함수
/// L: 공간 영역 길이
/// N: 공간 격자 수 (내부점: 1..N-1, 0과 N은 경계)
/// dt: 시간 스텝
/// steps: 시간 스텝 수
/// alpha: 열전달계수
pub fn heat_equation_1d<F>(f: F, L: f64, N: usize, dt: f64, steps: usize, alpha: f64) -> Vec<Vec<f64>>
where
    F: Fn(f64)->f64
{
    let dx = L/(N as f64);
    let r = alpha*dt/(dx*dx);
    assert!(r <= 0.5, "안정성 조건 위반: r=alpha*dt/dx² <=0.5 필요");

    // 초기 조건
    let mut u = vec![0.0; N+1];
    for i in 0..=N {
        let x = i as f64 * dx;
        u[i] = f(x);
    }
    // 경계조건 u(0,t)=u(L,t)=0이므로 u[0]=u[N]=0 유지

    let mut results = Vec::with_capacity(steps+1);
    results.push(u.clone());

    let mut u_new = u.clone();
    for _ in 0..steps {
        for i in 1..N {
            u_new[i] = u[i] + r*(u[i+1]-2.0*u[i]+u[i-1]);
        }
        // 경계조건
        u_new[0]=0.0;
        u_new[N]=0.0;

        u.copy_from_slice(&u_new);
        results.push(u.clone());
    }

    results
}

예제: 초기 조건과 시각화

예를 들어 f(x)=sin(πx/L) 형태로 초기 온도 분포를 두고, 시간이 지남에 따라 열이 퍼져가는 상황을 시뮬레이션할 수 있습니다.

fn main() {
    let L=1.0;
    let N=50;
    let dt=0.0005;
    let steps=2000;
    let alpha=0.01;

    let results = heat_equation_1d(|x| (std::f64::consts::PI*x/L).sin(), L,N,dt,steps,alpha);

    // results[k]: 시간 k*dt후의 온도 분포 u(x_i)
    // 여기서 간단히 결과를 출력하거나 파일로 저장해 Python 등으로 시각화 가능
    println!("t={}, u(x=0.5)={}", steps as f64 * dt, results.last().unwrap()[N/2]);
}

시간이 지날수록 중앙점 (x=0.5) 의 온도분포가 0에 가까워지는(열이 퍼져 균일해지는) 경향을 관찰할 수 있습니다.

향후 개선점

• 고차 정확도 스킴: Centered differencing 외에 Crank-Nicolson, implicit scheme, 더 정교한 PDE 솔버를 추가 가능
• 2D/3D 확장: 2차원, 3차원 격자로 확장, 라플라스 방정식, 파동 방정식 등 다양한 PDE 해결 가능
• FEM(Finite Element Method), FVM(Finite Volume Method): 더 정교한 PDE 솔버 방법론 적용
• 경계조건 다양화: Neumann, Robin 경계조건 적용
• 수치적 안정성, 오차 추정, 어댑티브 격자: 실전 PDE 솔버 성능 향상

결론

이번 글에서는 PDE 솔빙의 기초로서 1D 열방정식을 유한차분법으로 근사 풀이하는 과정을 구현했습니다. 이는 PDE 영역에 대한 첫 발걸음일 뿐이며, 실제 응용에서는 더 안정적이고 정교한 스킴, 다차원 문제, 복잡한 경계/초기 조건 등을 처리해야 합니다. 그럼에도 러스트의 장점을 통해 모듈성, 안전성, 확장성을 갖춘 PDE 솔버 개발이 가능함을 확인할 수 있습니다.

유용한 링크와 리소스

• Finite Difference Method(FDM): LeVeque "Finite Difference Methods for Ordinary and Partial Differential Equations"
• PDE Numerical Methods: Strikwerda "Finite Difference Schemes and Partial Differential Equations"
• Rust PDE 예제: 현재 많지 않지만, Rust+ndarray+rayon 등 조합으로 구현 사례 탐색 가능