공업수학 완전 정복 1편: 상미분방정식(ODE)

공학도라면 피해갈 수 없는 과목이 바로 공업수학입니다. 특히 상미분방정식(ODE, Ordinary Differential Equation)은 전자회로 분석, 제어 시스템 설계, 기계 진동 해석 등 거의 모든 공학 분야의 수학적 기반이 됩니다. 이 글에서는 ODE의 기초 개념부터 라플라스 변환까지, 공학 응용 예제를 중심으로 체계적으로 설명합니다.

1. 미분방정식 기초

미분방정식이란?

미분방정식은 미지함수와 그 도함수들 사이의 관계를 나타내는 방정식입니다. 예를 들어:

$\frac{dy}{dx} = 2x, \quad y'' + 3y' + 2y = 0, \quad \frac{\partial u}{\partial t} = k\frac{\partial^2 u}{\partial x^2}$

이처럼 미분방정식은 자연 현상과 공학 시스템을 수학적으로 모델링하는 핵심 도구입니다.

분류: ODE vs PDE

상미분방정식(ODE): 하나의 독립변수에 대한 도함수만 포함합니다.

$\frac{d^2x}{dt^2} + 2\frac{dx}{dt} + x = 0$

편미분방정식(PDE): 두 개 이상의 독립변수에 대한 편미분을 포함합니다.

$\frac{\partial u}{\partial t} = c^2\frac{\partial^2 u}{\partial x^2}$

이번 글에서는 ODE에 집중합니다.

차수(Order)와 급(Degree)

차수는 방정식에 등장하는 가장 높은 도함수의 차수입니다.

• 1차 ODE: $y' + y = x$ (최고 도함수가 $y'$)
• 2차 ODE: $y'' + 2y' + y = 0$ (최고 도함수가 $y''$)

**급(degree)**은 최고 차수 도함수의 지수입니다. 대부분의 공학 문제에서 급은 1입니다.

선형 vs 비선형

선형 ODE: 미지함수와 그 도함수가 1차 항으로만 나타납니다.

$a_n(x)y^{(n)} + a_{n-1}(x)y^{(n-1)} + \cdots + a_1(x)y' + a_0(x)y = f(x)$

비선형 ODE: $y^2$, $y \cdot y'$, $\sin(y)$ 등 비선형 항이 포함됩니다.

$y'' + \sin(y) = 0 \quad \text{(진자 방정식)}$

일반해, 특수해, 특이해

• 일반해: 임의의 상수를 포함하는 완전한 해
• 특수해: 초기조건을 적용하여 상수를 결정한 해
• 특이해: 일반해로부터 얻을 수 없는 특별한 해

공학에서 ODE의 등장

공학 시스템은 대부분 ODE로 기술됩니다.

RC 충전 회로: $RC\frac{dv}{dt} + v = V_s$

스프링-질량-감쇠기 시스템: $m\ddot{x} + c\dot{x} + kx = F(t)$

인구 성장 모델: $\frac{dP}{dt} = kP$

RL 회로: $L\frac{di}{dt} + Ri = E(t)$

이 모든 방정식이 ODE이며, 풀이 방법을 배우면 공학 시스템의 동적 거동을 완전히 분석할 수 있습니다.

2. 1차 상미분방정식

2.1 변수분리법 (Separation of Variables)

형태: $\frac{dy}{dx} = f(x)g(y)$

풀이 과정:

$\frac{dy}{g(y)} = f(x)dx$

양변을 적분합니다:

$\int \frac{dy}{g(y)} = \int f(x)dx + C$

예제: RC 충전 회로

회로 방정식: $\frac{dv}{dt} = \frac{V_s - v}{RC}$

변수분리: $\frac{dv}{V_s - v} = \frac{dt}{RC}$

적분: $-\ln|V_s - v| = \frac{t}{RC} + C_1$

초기조건 $v(0) = 0$ 적용:

$V_s - v = V_s e^{-t/(RC)}$

$\boxed{v(t) = V_s\left(1 - e^{-t/(RC)}\right)}$

이 결과는 RC 회로의 **과도 응답(transient response)**을 나타냅니다. 시정수 $\tau = RC$로 결정되는 속도로 전압이 $V_s$로 수렴합니다.

물리적 의미:

• $t = \tau$일 때: $v \approx 0.632 V_s$ (63.2% 충전)
• $t = 5\tau$일 때: $v \approx 0.993 V_s$ (거의 완전 충전)

2.2 적분인자법 (Integrating Factor Method)

선형 1차 ODE의 표준형: $\frac{dy}{dx} + P(x)y = Q(x)$

적분인자: $\mu(x) = e^{\int P(x)dx}$

유도 과정:

양변에 $\mu(x)$를 곱합니다: $\mu(x)\frac{dy}{dx} + \mu(x)P(x)y = \mu(x)Q(x)$

좌변은 곱의 미분 규칙에 의해: $\frac{d}{dx}[\mu(x)y] = \mu(x)Q(x)$

적분하면: $\mu(x)y = \int \mu(x)Q(x)dx + C$

$\boxed{y = \frac{1}{\mu(x)}\left[\int \mu(x)Q(x)dx + C\right]}$

예제: RL 회로

$L\frac{di}{dt} + Ri = E$

표준형으로 변환: $\frac{di}{dt} + \frac{R}{L}i = \frac{E}{L}$

적분인자: $\mu(t) = e^{(R/L)t}$

$\frac{d}{dt}\left[e^{(R/L)t} \cdot i\right] = \frac{E}{L}e^{(R/L)t}$

적분: $e^{(R/L)t} \cdot i = \frac{E}{R}e^{(R/L)t} + C$

초기조건 $i(0) = 0$ 적용:

$\boxed{i(t) = \frac{E}{R}\left(1 - e^{-(R/L)t}\right)}$

시정수 $\tau = L/R$로 전류가 정상 상태 $E/R$에 수렴합니다.

2.3 완전 미분방정식 (Exact ODE)

형태: $M(x,y)dx + N(x,y)dy = 0$

완전 조건: $\frac{\partial M}{\partial y} = \frac{\partial N}{\partial x}$

이 조건을 만족하면 $F(x,y) = C$인 퍼텐셜 함수 $F$가 존재합니다.

풀이법:

1. $\frac{\partial F}{\partial x} = M$에서 $F$를 $x$에 대해 적분
2. $\frac{\partial F}{\partial y} = N$ 조건으로 적분 상수 결정
3. 일반해: $F(x,y) = C$

예제:

$(2xy + y^2)dx + (x^2 + 2xy)dy = 0$

완전 여부 확인: $\frac{\partial M}{\partial y} = 2x + 2y = \frac{\partial N}{\partial x}$

완전 방정식이므로: $F = \int M\,dx = \int (2xy + y^2)dx = x^2y + xy^2 + g(y)$

$\frac{\partial F}{\partial y} = x^2 + 2xy + g'(y) = N = x^2 + 2xy$이므로 $g'(y) = 0$, $g(y) = C$

일반해: $x^2y + xy^2 = C$

2.4 베르누이 방정식 (Bernoulli Equation)

형태: $\frac{dy}{dx} + P(x)y = Q(x)y^n \quad (n \neq 0, 1)$

치환: $v = y^{1-n}$으로 놓으면 선형 ODE로 변환됩니다.

$v' = (1-n)y^{-n}y'$이므로 원래 방정식을 $y^n$으로 나누면:

$y^{-n}\frac{dy}{dx} + P(x)y^{1-n} = Q(x)$

$\frac{1}{1-n}\frac{dv}{dx} + P(x)v = Q(x)$

$\frac{dv}{dx} + (1-n)P(x)v = (1-n)Q(x)$

이는 $v$에 대한 선형 1차 ODE로, 적분인자법으로 풀 수 있습니다.

응용: 물류 성장 모델 $\frac{dP}{dt} = rP\left(1 - \frac{P}{K}\right)$는 $n=2$인 베르누이 방정식입니다.

3. 2차 선형 상미분방정식

3.1 제차 방정식 (Homogeneous)

표준형: $ay'' + by' + cy = 0$

특성방정식 유도: $y = e^{rx}$를 대입합니다.

$a r^2 e^{rx} + b r e^{rx} + c e^{rx} = 0$

$\boxed{ar^2 + br + c = 0}$

판별식 $\Delta = b^2 - 4ac$의 부호에 따라 세 가지 경우가 생깁니다.

경우 1: 서로 다른 두 실근 ($\Delta > 0$)

$r_1 \neq r_2$이면: $y = C_1 e^{r_1 x} + C_2 e^{r_2 x}$

경우 2: 중복 실근 ($\Delta = 0$)

$r_1 = r_2 = r = -b/(2a)$이면: $y = (C_1 + C_2 x)e^{rx}$

$C_2 x e^{rx}$를 추가해야 하는 이유는 **론스키안(Wronskian)**이 0이 되지 않도록 독립적인 두 번째 해를 구성해야 하기 때문입니다.

경우 3: 켤레 복소근 ($\Delta < 0$)

$r = \alpha \pm i\beta$ (단, $\alpha = -b/(2a)$, $\beta = \sqrt{4ac - b^2}/(2a)$)이면: $y = e^{\alpha x}(C_1\cos\beta x + C_2\sin\beta x)$

오일러 공식 $e^{i\beta x} = \cos\beta x + i\sin\beta x$를 이용해 실수해로 변환한 형태입니다.

3.2 RLC 회로 완전 분석

회로 방정식 (전하 $q(t)$ 기준): $L\frac{d^2q}{dt^2} + R\frac{dq}{dt} + \frac{q}{C} = 0$

표준형으로 정리: $q'' + \frac{R}{L}q' + \frac{1}{LC}q = 0$

특성방정식: $r^2 + \frac{R}{L}r + \frac{1}{LC} = 0$

$r = -\frac{R}{2L} \pm \sqrt{\left(\frac{R}{2L}\right)^2 - \frac{1}{LC}}$

감쇠 계수 $\alpha = R/(2L)$, 공진 주파수 $\omega_0 = 1/\sqrt{LC}$로 정의하면:

$r = -\alpha \pm \sqrt{\alpha^2 - \omega_0^2}$

과감쇠 (Overdamped): $\alpha > \omega_0$ ($R > 2\sqrt{L/C}$)

$q(t) = C_1 e^{r_1 t} + C_2 e^{r_2 t} \quad (r_1, r_2 < 0)$

두 지수항의 합으로, 진동 없이 단조 감소합니다.

임계감쇠 (Critically Damped): $\alpha = \omega_0$ ($R = 2\sqrt{L/C}$)

$q(t) = (C_1 + C_2 t)e^{-\alpha t}$

진동 없이 가장 빠르게 평형 상태로 복귀합니다. 제어 시스템 설계에서 매우 중요합니다.

부족감쇠 (Underdamped): $\alpha < \omega_0$ ($R < 2\sqrt{L/C}$)

감쇠 진동 주파수 $\omega_d = \sqrt{\omega_0^2 - \alpha^2}$로 정의하면:

$q(t) = e^{-\alpha t}(C_1\cos\omega_d t + C_2\sin\omega_d t)$

진동하면서 감쇠합니다. 통신 회로, 진동 시스템에서 흔히 나타납니다.

물리적 의미 요약:

3.3 비제차 방정식 (Non-homogeneous)

형태: $ay'' + by' + cy = f(x)$

완전 해: $y = y_h + y_p$

• $y_h$: 제차 방정식의 일반해 (여함수, complementary function)
• $y_p$: 특수해 (particular solution)

미결정계수법 (Method of Undetermined Coefficients)

$f(x)$의 형태에 따라 $y_p$의 형태를 추측합니다.

주의: $y_p$의 형태가 $y_h$의 항과 겹치면 $x$를 곱해서 수정합니다.

예제: $y'' + 4y = 3\cos 2x$

제차해: $y_h = C_1\cos 2x + C_2\sin 2x$

$y_p$의 추측형이 $y_h$와 겹치므로 수정: $y_p = x(A\cos 2x + B\sin 2x)$

대입 후 계산하면: $A = 0$, $B = 3/4$

$y_p = \frac{3}{4}x\sin 2x$

이는 공진(resonance) 현상으로, 진폭이 시간에 비례하여 증가합니다.

강제 진동 응답:

$m\ddot{x} + c\dot{x} + kx = F_0\cos\omega t$

정상 상태 특수해: $x_p(t) = \frac{F_0/k}{\sqrt{(1-r^2)^2 + (2\zeta r)^2}}\cos(\omega t - \phi)$

여기서 진동수비 $r = \omega/\omega_n$, 감쇠비 $\zeta = c/(2\sqrt{mk})$.

$r = 1$ (공진 조건)에서 진폭이 최대가 됩니다.

매개변수 변환법 (Variation of Parameters)

미결정계수법을 적용할 수 없는 일반적인 $f(x)$에 사용합니다.

$y_h = C_1 y_1 + C_2 y_2$가 주어질 때:

$y_p = -y_1 \int \frac{y_2 f}{W}dx + y_2 \int \frac{y_1 f}{W}dx$

론스키안(Wronskian): $W = y_1 y_2' - y_1' y_2$

이 공식은 모든 연속 함수 $f(x)$에 대해 적용 가능합니다.

4. 연립 미분방정식

4.1 행렬 표현

$n$개의 1차 ODE로 이루어진 연립 방정식:

$\frac{d\mathbf{X}}{dt} = A\mathbf{X}$

여기서 $\mathbf{X} = [x_1, x_2, \ldots, x_n]^T$, $A$는 $n \times n$ 계수 행렬.

4.2 고유값/고유벡터 풀이

해의 형태: $\mathbf{X} = \mathbf{v}e^{\lambda t}$

대입하면: $\lambda \mathbf{v} = A\mathbf{v}$

즉, $\lambda$는 $A$의 고유값, $\mathbf{v}$는 대응하는 고유벡터입니다.

특성방정식: $\det(A - \lambda I) = 0$

예제: 두 개의 연결된 탱크 시스템

$\frac{d}{dt}\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} -2 & 1 \\ 1 & -2 \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \end{pmatrix}$

특성방정식: $(\lambda+2)^2 - 1 = 0$

$\lambda_1 = -1, \quad \lambda_2 = -3$

고유벡터: $\mathbf{v}_1 = [1, 1]^T$, $\mathbf{v}_2 = [1, -1]^T$

일반해: $\mathbf{X} = C_1 \begin{pmatrix} 1 \\ 1 \end{pmatrix} e^{-t} + C_2 \begin{pmatrix} 1 \\ -1 \end{pmatrix} e^{-3t}$

4.3 회로망 해석 적용

두 루프가 있는 회로 (키르히호프 전압 법칙):

$L_1\frac{di_1}{dt} + R_1 i_1 - M\frac{di_2}{dt} = E_1(t)$ $L_2\frac{di_2}{dt} + R_2 i_2 - M\frac{di_1}{dt} = E_2(t)$

여기서 $M$은 상호 인덕턴스. 이를 행렬 형태로 정리하여 고유값법 또는 라플라스 변환으로 풀 수 있습니다.

5. 라플라스 변환 (Laplace Transform)

5.1 정의와 수렴

정의: $\mathcal{L}\{f(t)\} = F(s) = \int_0^\infty f(t)e^{-st}dt$

$s = \sigma + j\omega$는 복소 변수이며, 적분이 수렴하려면 $\text{Re}(s) > \sigma_c$ (수렴 좌평면)를 만족해야 합니다.

왜 라플라스 변환을 쓰는가?

• 미분 연산이 대수 연산으로 변환됩니다
• 초기조건을 자동으로 처리합니다
• 복잡한 ODE를 쉬운 대수 방정식으로 풀 수 있습니다
• 전달 함수(Transfer Function) 개념의 기반입니다

5.2 주요 변환쌍

유도 예시: $\mathcal{L}\{e^{at}\}$

$\mathcal{L}\{e^{at}\} = \int_0^\infty e^{at} e^{-st} dt = \int_0^\infty e^{-(s-a)t} dt = \left[\frac{e^{-(s-a)t}}{-(s-a)}\right]_0^\infty = \frac{1}{s-a}$

$s > a$일 때 수렴합니다.

5.3 주요 성질

선형성: $\mathcal{L}\{af(t) + bg(t)\} = aF(s) + bG(s)$

시간 이동 (제1 이동 정리): $\mathcal{L}\{e^{at}f(t)\} = F(s-a)$

주파수 이동 (제2 이동 정리): $\mathcal{L}\{f(t-a)u(t-a)\} = e^{-as}F(s)$

미분 변환 (가장 중요한 성질): $\mathcal{L}\{f'(t)\} = sF(s) - f(0)$ $\mathcal{L}\{f''(t)\} = s^2F(s) - sf(0) - f'(0)$ $\mathcal{L}\{f^{(n)}(t)\} = s^nF(s) - s^{n-1}f(0) - s^{n-2}f'(0) - \cdots - f^{(n-1)}(0)$

초기조건이 $s$-도메인에서 자동으로 포함됩니다.

적분 변환: $\mathcal{L}\left\{\int_0^t f(\tau)d\tau\right\} = \frac{F(s)}{s}$

합성곱 정리(Convolution Theorem): $\mathcal{L}\{(f*g)(t)\} = F(s)G(s)$

여기서 $(f*g)(t) = \int_0^t f(\tau)g(t-\tau)d\tau$

이 정리는 시스템의 임펄스 응답과 입력 신호의 합성곱이 출력임을 수학적으로 증명합니다.

최종값 정리: $\lim_{t\to\infty} f(t) = \lim_{s\to 0} sF(s)$

정상 상태 값을 $t$-도메인 풀이 없이 구할 수 있습니다.

초기값 정리: $\lim_{t\to 0^+} f(t) = \lim_{s\to\infty} sF(s)$

5.4 역 라플라스 변환

부분 분수 분해(Partial Fraction Decomposition):

$F(s) = P(s)/Q(s)$에서 분모를 인수분해하여 표 변환쌍의 합으로 분해합니다.

경우 1: 서로 다른 실수 극점

$F(s) = \frac{N(s)}{(s-p_1)(s-p_2)\cdots(s-p_n)} = \frac{A_1}{s-p_1} + \frac{A_2}{s-p_2} + \cdots$

$A_k = \lim_{s\to p_k}(s-p_k)F(s)$

경우 2: 중복 극점 ($s = p$가 $m$중근)

$F(s) = \cdots + \frac{A_m}{(s-p)^m} + \frac{A_{m-1}}{(s-p)^{m-1}} + \cdots + \frac{A_1}{s-p}$

$A_k = \frac{1}{(m-k)!}\left[\frac{d^{m-k}}{ds^{m-k}}(s-p)^m F(s)\right]_{s=p}$

경우 3: 켤레 복소 극점

$\frac{As + B}{s^2 + 2\alpha s + (\alpha^2 + \beta^2)} \Rightarrow e^{-\alpha t}(C_1\cos\beta t + C_2\sin\beta t)$

헤비사이드 전개 정리(Heaviside Expansion):

분모가 서로 다른 1차 인수의 곱일 때: $F(s) = \sum_{k=1}^n \frac{N(p_k)}{Q'(p_k)} \cdot \frac{1}{s - p_k}$

5.5 라플라스 변환으로 ODE 풀기

예제: RLC 직렬 회로 과도 응답

$L\frac{d^2i}{dt^2} + R\frac{di}{dt} + \frac{i}{C} = E_0\delta(t)$

초기조건: $i(0) = 0$, $i'(0) = E_0/L$

풀이:

양변에 라플라스 변환 적용:

$L[s^2I(s) - si(0) - i'(0)] + R[sI(s) - i(0)] + \frac{I(s)}{C} = E_0$

초기조건 대입:

$Ls^2 I(s) - LE_0/L + RsI(s) + \frac{I(s)}{C} = E_0$

$I(s)\left(Ls^2 + Rs + \frac{1}{C}\right) = 2E_0 - E_0 = E_0$

잠깐, 초기조건 수정: $i(0) = 0$, $i'(0) = 0$ (캐패시터에 전하 없음)이면 임펄스 입력으로 인해:

$I(s) = \frac{E_0}{Ls^2 + Rs + 1/C} = \frac{E_0/L}{s^2 + (R/L)s + 1/(LC)}$

$\alpha = R/(2L)$, $\omega_0 = 1/\sqrt{LC}$, $\omega_d = \sqrt{\omega_0^2 - \alpha^2}$ (부족감쇠 경우)로 정의하면:

$I(s) = \frac{E_0/L}{(s+\alpha)^2 + \omega_d^2}$

역 라플라스 변환 (표에서 $\mathcal{L}\{e^{-\alpha t}\sin\omega_d t\} = \omega_d/[(s+\alpha)^2+\omega_d^2]$):

$i(t) = \frac{E_0}{L\omega_d}e^{-\alpha t}\sin\omega_d t$

이것이 RLC 회로의 임펄스 응답이며, 동시에 그린 함수이기도 합니다.

6. Python으로 ODE 수치 풀기

6.1 scipy.integrate.solve_ivp 활용

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# RLC 회로 시뮬레이션 (2차 ODE -> 연립 1차 ODE)
# L * q'' + R * q' + q/C = Vs (계단 입력)
# 상태 변수: y[0] = q (전하), y[1] = i = dq/dt (전류)

def rlc_circuit(t, y, R, L, C, Vs):
    q, dq_dt = y
    # q'' = (Vs - R*i - q/C) / L
    d2q_dt2 = (Vs - R * dq_dt - q / C) / L
    return [dq_dt, d2q_dt2]

# 회로 파라미터
R = 100      # 저항 (Ohm)
L = 0.1      # 인덕턴스 (H)
C = 1e-6     # 캐패시턴스 (F)
Vs = 10.0    # 전원 전압 (V)

# 공진 주파수와 감쇠비 계산
omega_0 = 1.0 / np.sqrt(L * C)
alpha = R / (2 * L)
zeta = alpha / omega_0
print(f"공진 주파수: {omega_0/(2*np.pi):.1f} Hz")
print(f"감쇠비 zeta: {zeta:.3f}")
if zeta > 1:
    print("-> 과감쇠")
elif zeta == 1:
    print("-> 임계감쇠")
else:
    print("-> 부족감쇠")

# 수치 풀이
y0 = [0.0, 0.0]  # 초기조건: q(0) = 0, i(0) = 0
t_span = (0, 0.01)
t_eval = np.linspace(0, 0.01, 2000)

sol = solve_ivp(
    rlc_circuit, t_span, y0,
    t_eval=t_eval,
    args=(R, L, C, Vs),
    method='RK45',
    rtol=1e-8,
    atol=1e-10
)

# 전압 계산: v_C = q/C
v_C = sol.y[0] / C
i = sol.y[1]

# 시각화
fig, axes = plt.subplots(2, 1, figsize=(10, 8))

axes[0].plot(sol.t * 1e3, v_C, 'b-', linewidth=2, label='커패시터 전압')
axes[0].axhline(y=Vs, color='r', linestyle='--', label=f'정상 상태 {Vs}V')
axes[0].set_xlabel('시간 (ms)')
axes[0].set_ylabel('전압 (V)')
axes[0].set_title('RLC 직렬 회로 - 계단 응답')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(sol.t * 1e3, i * 1e3, 'g-', linewidth=2, label='전류')
axes[1].set_xlabel('시간 (ms)')
axes[1].set_ylabel('전류 (mA)')
axes[1].set_title('RLC 회로 전류')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('rlc_response.png', dpi=150)
plt.show()

6.2 sympy로 해석해 구하기

import sympy as sp

# 기호 정의
t, s = sp.symbols('t s')
R_sym, L_sym, C_sym = sp.symbols('R L C', positive=True)
Vs_sym = sp.Symbol('Vs')

# RC 회로 해석해 (sympy로 미분방정식 풀기)
v = sp.Function('v')
tau = sp.Symbol('tau', positive=True)

# dv/dt = (Vs - v) / (R*C)
ode = sp.Eq(v(t).diff(t), (Vs_sym - v(t)) / (R_sym * C_sym))

# 초기조건 v(0) = 0
solution = sp.dsolve(ode, v(t), ics={v(0): 0})
print("RC 회로 해:")
print(solution)
sp.pprint(solution)

# 라플라스 변환 계산
f = sp.exp(-t) * sp.sin(2*t)
F = sp.laplace_transform(f, t, s, noconds=True)
print(f"\nL{{e^(-t) sin(2t)}} = {F}")

# 역 라플라스 변환
F_inv = sp.inverse_laplace_transform(
    1 / (s**2 + 2*s + 5), s, t
)
print(f"L^(-1){{1/(s^2+2s+5)}} = {F_inv}")

6.3 다양한 적분기 비교

import numpy as np
import matplotlib.pyplot as plt

# 오일러법 vs RK4 비교
# 예제: y' = -y, y(0) = 1, 정확해 y = e^(-t)

def f(t, y):
    return -y

def euler_method(f, t0, y0, h, n):
    t = np.zeros(n+1)
    y = np.zeros(n+1)
    t[0], y[0] = t0, y0
    for i in range(n):
        y[i+1] = y[i] + h * f(t[i], y[i])
        t[i+1] = t[i] + h
    return t, y

def rk4_method(f, t0, y0, h, n):
    t = np.zeros(n+1)
    y = np.zeros(n+1)
    t[0], y[0] = t0, y0
    for i in range(n):
        k1 = h * f(t[i], y[i])
        k2 = h * f(t[i] + h/2, y[i] + k1/2)
        k3 = h * f(t[i] + h/2, y[i] + k2/2)
        k4 = h * f(t[i] + h, y[i] + k3)
        y[i+1] = y[i] + (k1 + 2*k2 + 2*k3 + k4) / 6
        t[i+1] = t[i] + h
    return t, y

h = 0.5
n = 10
t0, y0 = 0.0, 1.0

t_euler, y_euler = euler_method(f, t0, y0, h, n)
t_rk4, y_rk4 = rk4_method(f, t0, y0, h, n)
t_exact = np.linspace(0, n*h, 200)
y_exact = np.exp(-t_exact)

plt.figure(figsize=(10, 5))
plt.plot(t_exact, y_exact, 'k-', linewidth=2, label='정확해')
plt.plot(t_euler, y_euler, 'r--o', label=f'오일러법 (h={h})')
plt.plot(t_rk4, y_rk4, 'b--s', label=f'RK4 (h={h})')
plt.xlabel('t')
plt.ylabel("y(t)")
plt.title("오일러법 vs RK4 비교: y' = -y")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# 오차 비교
print(f"t=5에서 오일러법 오차: {abs(y_euler[-1] - np.exp(-5)):.6f}")
print(f"t=5에서 RK4 오차: {abs(y_rk4[-1] - np.exp(-5)):.6f}")

6.4 라플라스 변환을 Python으로 자동화

import sympy as sp

s, t = sp.symbols('s t', real=True)
a, omega = sp.symbols('a omega', positive=True)

# 라플라스 변환 표 자동 생성
functions = {
    '1': sp.Integer(1),
    't': t,
    't^2': t**2,
    't^n': t**3,
    'e^(at)': sp.exp(a*t),
    'sin(wt)': sp.sin(omega*t),
    'cos(wt)': sp.cos(omega*t),
    't*e^(at)': t * sp.exp(a*t),
    'e^(at)*sin(wt)': sp.exp(a*t) * sp.sin(omega*t),
}

print("라플라스 변환표:")
print("-" * 50)
for name, func in functions.items():
    F = sp.laplace_transform(func, t, s, noconds=True)
    print(f"L{{{name}}} = {sp.simplify(F)}")

7. 공학 응용 종합 예제

7.1 전원이 있는 RLC 회로 과도 분석

문제: $R = 50\,\Omega$, $L = 0.2\,\text{H}$, $C = 100\,\mu\text{F}$인 직렬 RLC 회로에 $t=0$에서 $V_s = 100\,\text{V}$ DC 전원이 연결된다. 초기조건은 $v_C(0) = 0$, $i(0) = 0$이다.

라플라스 변환으로 풀기:

키르히호프 전압 법칙: $L\frac{di}{dt} + Ri + \frac{1}{C}\int i\,dt = V_s u(t)$

$i = C dv_C/dt$로 치환하여 전하 $q$로 정리: $L\ddot{q} + R\dot{q} + \frac{q}{C} = V_s$

라플라스 변환: $\left(Ls^2 + Rs + \frac{1}{C}\right)Q(s) = \frac{V_s}{s}$

$Q(s) = \frac{V_s}{s\left(Ls^2 + Rs + \frac{1}{C}\right)} = \frac{V_s/L}{s\left(s^2 + \frac{R}{L}s + \frac{1}{LC}\right)}$

수치 대입: $\alpha = R/(2L) = 125\,\text{s}^{-1}$, $\omega_0 = 1/\sqrt{LC} \approx 223.6\,\text{s}^{-1}$

$\omega_0 > \alpha$이므로 부족감쇠입니다.

$\omega_d = \sqrt{\omega_0^2 - \alpha^2} \approx 185.0\,\text{s}^{-1}$

부분 분수 분해 후 역 라플라스 변환:

$v_C(t) = V_s\left[1 - e^{-\alpha t}\left(\cos\omega_d t + \frac{\alpha}{\omega_d}\sin\omega_d t\right)\right]$

$= 100\left[1 - e^{-125t}\left(\cos 185t + 0.676\sin 185t\right)\right]\,\text{V}$

7.2 스프링-질량-감쇠기 공진 분석

문제: $m = 1\,\text{kg}$, $k = 100\,\mathrm{N/m}$, $c = 4\,\mathrm{N \cdot s/m}$, $F(t) = 10\cos\omega t\,\mathrm{N}$

고유 진동수: $\omega_n = \sqrt{k/m} = 10\,\text{rad/s}$

감쇠비: $\zeta = c/(2\sqrt{mk}) = 0.2$ (부족감쇠)

정상 상태 응답 진폭: $X = \frac{F_0/k}{\sqrt{(1-r^2)^2 + (2\zeta r)^2}}$

$r = \omega/\omega_n = 1$ (공진)일 때: $X_{\mathrm{res}} = \frac{F_0/k}{2\zeta} = \frac{10/100}{2 \times 0.2} = 0.25\,\text{m}$

공진 없을 때 정적 처짐 $F_0/k = 0.1\,\text{m}$의 2.5배입니다. 이것이 공진의 위험성입니다.

정리 및 다음 단계

이번 글에서 다룬 내용:

1. 미분방정식 기초: 분류, 차수, 선형성, 해의 종류
2. 1차 ODE 풀이법: 변수분리법, 적분인자법, 완전 ODE, 베르누이 방정식
3. 2차 선형 ODE: 특성방정식, 세 가지 경우, RLC 회로 완전 분석
4. 비제차 방정식: 미결정계수법, 매개변수 변환법, 강제 진동
5. 연립 ODE: 행렬 표현, 고유값/고유벡터 풀이
6. 라플라스 변환: 정의, 변환쌍, 성질, 역변환, ODE 풀이
7. Python 구현: scipy, sympy로 수치/기호 풀이

다음 편에서는 **푸리에 급수/변환과 편미분방정식(PDE)**을 다룰 예정입니다. 신호처리와 전자기학의 핵심 수학 도구입니다.

참고 자료

• Kreyszig, E. "Advanced Engineering Mathematics", 10th Edition, Wiley
• Boyce, W. & DiPrima, R. "Elementary Differential Equations", 11th Edition, Wiley
• Simmons, G. "Differential Equations with Applications and Historical Notes", 3rd Edition
• SciPy ODE 솔버 공식 문서
• SymPy 라플라스 변환 문서