이전 글까지 우리는 선형대수(행렬 분해, 고유값 문제), 복소수 및 FFT를 통한 다항식 연산 최적화 등 다양한 주제를 다루며 러스트 기반 수학/과학 계산 라이브러리의 기초를 닦아 왔습니다. 이제 범위를 한 단계 넓혀, 수치적 적분(Numerical Integration)과 상미분방정식(ODE, Ordinary Differential Equation) 풀이 방법을 다루어보겠습니다.
수치 적분과 ODE 풀이는 물리, 화학, 생물, 경제학, 머신러닝, 신호처리 등 광범위한 분야에서 필수적인 역할을 합니다. 컴퓨터로 적분값을 근사하거나, 시간에 따른 동적 시스템의 변화를 예측하기 위해 Euler, RK4(Runge-Kutta 4차) 같은 고전적이고 믿을 수 있는 알고리즘이 자주 사용됩니다.
목표
- 수치 적분:
- 단순 1차원 적분: 사다리꼴 적분법, Simpson 적분법을 구현해 f(x) 함수를 정해진 구간에서 근사적 적분값을 계산해봅니다.
- 다중 적분 확장 가능성도 언급(시간관계상 개념만).
- ODE 풀이 기초:
- Euler 방법: 가장 단순한 1계 ODE y'(t) = f(t,y)의 근사적 해를 시간 스텝 단위로 구함.
- RK4 방법: 더 정확한 고전적 방법으로 ODE를 풀어봅니다.
이들은 모두 라이브러리 사용자가 f(t), f(t,y) 같은 콜백 함수를 전달하고, 우리가 구현한 수치 알고리즘을 통해 근사적 결과를 얻는 식의 인터페이스를 가지게 될 것입니다.
수치 적분 예시: Simpson 법칙
∫_a^b f(x) dx를 근사하기 위해 Simpson 적분법을 구현해보겠습니다. Simpson 법칙은 구간을 짝수 개로 나누고, 홀수 인덱스에서의 함수값에 가중치를 더해 정밀한 근사값을 얻습니다.
// src/integration.rs
pub fn simpson_integration<F>(f: F, a: f64, b: f64, n: usize) -> f64
where
F: Fn(f64)->f64
{
assert!(n > 0 && n % 2 == 0, "n은 짝수여야 함");
let h = (b - a) / (n as f64);
let mut s = f(a) + f(b);
for i in 1..n {
let x = a + i as f64 * h;
if i % 2 == 0 {
s += 2.0 * f(x);
} else {
s += 4.0 * f(x);
}
}
s * h / 3.0
}
// 예: f(x)=x²를 0~1구간에서 적분 -> 실제값 1/3≈0.3333
// simpson_integration(|x| x*x, 0.0,1.0,1000) -> 근사값
ODE 풀이: Euler 방법
Euler 방법은 ODE y'(t)=f(t,y), y(t0)=y0에 대해 한 스텝 h만큼 전진하여 y(t+h)=y(t)+h*f(t,y)를 통해 근사합니다. 단순하지만 정확도가 제한적입니다.
// src/ode.rs
pub fn euler<F>(f: F, t0: f64, y0: f64, h: f64, steps: usize) -> (Vec<f64>, Vec<f64>)
where
F: Fn(f64, f64)->f64
{
let mut t = t0;
let mut y = y0;
let mut ts = Vec::with_capacity(steps+1);
let mut ys = Vec::with_capacity(steps+1);
ts.push(t);
ys.push(y);
for _ in 0..steps {
let dydt = f(t,y);
y = y + h*dydt;
t = t + h;
ts.push(t);
ys.push(y);
}
(ts, ys)
}
// 예: y'(t)=y, y(0)=1 일 때 해는 y(t)=e^t
// euler(|t,y| y, 0.0, 1.0, 0.1, 100) -> t=10일 때 근사값 y(10)는 e^10≈22026.4657 근사
ODE 풀이: Runge-Kutta 4차(RK4)
RK4는 4번의 중간 평가를 통해 더 정확한 1스텝 전진을 구현합니다. 공식:
k1 = f(t, y)
k2 = f(t+h/2, y+h*k1/2)
k3 = f(t+h/2, y+h*k2/2)
k4 = f(t+h, y+h*k3)
y(t+h) = y(t) + (h/6)*(k1+2k2+2k3+k4)
// src/ode.rs (추가)
pub fn rk4<F>(f: F, t0: f64, y0: f64, h: f64, steps: usize) -> (Vec<f64>, Vec<f64>)
where
F: Fn(f64, f64)->f64
{
let mut t = t0;
let mut y = y0;
let mut ts = Vec::with_capacity(steps+1);
let mut ys = Vec::with_capacity(steps+1);
ts.push(t);
ys.push(y);
for _ in 0..steps {
let k1 = f(t, y);
let k2 = f(t + h*0.5, y + h*k1*0.5);
let k3 = f(t + h*0.5, y + h*k2*0.5);
let k4 = f(t + h, y + h*k3);
y = y + (h/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4);
t = t + h;
ts.push(t);
ys.push(y);
}
(ts, ys)
}
// 동일한 y'(t)=y 예제로 rk4 적용 시, 훨씬 정확한 근사값 획득 가능.
예제 및 검증
fn main() {
// 적분 예제: ∫0² x dx = 2²/2=2
let approx = simpson_integration(|x| x, 0.0, 2.0, 1000);
println!("Simpson 적분 근사: {}", approx); //≈2.0 근접
// ODE 예제: y'(t)=y, y(0)=1, 정해해 e^t
let (ts_euler, ys_euler) = euler(|t,y| y, 0.0,1.0,0.1,100);
let (ts_rk4, ys_rk4) = rk4(|t,y| y, 0.0,1.0,0.1,100);
let true_val = (ts_rk4.last().unwrap()).exp();
println!("Euler at t=10: {}, error={}", ys_euler.last().unwrap(), (ys_euler.last().unwrap()-true_val).abs());
println!("RK4 at t=10: {}, error={}", ys_rk4.last().unwrap(), (ys_rk4.last().unwrap()-true_val).abs());
// RK4가 훨씬 작은 오차를 가짐
}
향후 과제
5편에서는 수치적 적분과 ODE 해석의 기초를 다루었습니다. 앞으로 다음과 같은 방향으로 라이브러리를 확장할 수 있습니다.
- 고차 적분법, 어댑티브 스텝: 스텝 크기 h를 가변적으로 조정해 오차를 일정 수준으로 관리
- 고차 ODE 시스템, stiff 방정식 해결: 다변수 ODE, 매우 빠른 변화 스케일을 가진 stiff 문제를 위한 BDF나 Rosenbrock-W methods
- 편미분방정식(PDE) 솔버: 유한차분(FDM), 유한요소(FEM) 기법을 통해 PDE 해석
- 최적화 알고리즘, 확률/통계 함수 구현: 머신러닝, 데이터 분석에 유용한 기능 추가
결론
이번 글에서 적분 및 ODE 풀이 기초를 다룸으로써, 러스트 기반 수학 라이브러리에서 선형대수, FFT, 고급 분해 기법뿐 아니라 동적 시스템, 해석적 계산의 한 영역까지 진입할 수 있음을 확인했습니다. 이러한 기능들은 물리 시뮬레이션, 화학 반응 속도 해석, 기계학습 모델 훈련 전단계 시뮬레이션 등 다양한 실제 문제에 응용될 수 있습니다.
유용한 링크와 리소스
- ODE 풀이 알고리즘: Gerald & Wheatley "Applied Numerical Analysis", Hairer & Wanner "Solving Ordinary Differential Equations"
- 수치적 적분 기법: Numerical Recipes, Wikipedia "Numerical Integration"
- Rust ODE/적분 참고 라이브러리: ODESolver Crate, SciRust 등 참고