공업수학 완전 정복 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)$$
진동하면서 감쇠합니다. 통신 회로, 진동 시스템에서 흔히 나타납니다.
**물리적 의미 요약**:
| 조건 | 명칭 | 응답 특성 |
| ----------------- | -------- | ------------------------ |
| $R > 2\sqrt{L/C}$ | 과감쇠 | 진동 없이 느리게 감쇠 |
| $R = 2\sqrt{L/C}$ | 임계감쇠 | 진동 없이 가장 빠른 감쇠 |
| $R < 2\sqrt{L/C}$ | 부족감쇠 | 진동하며 감쇠 |
| $R = 0$ | 무감쇠 | 지속 진동 |
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$의 형태를 추측합니다.
| $f(x)$의 형태 | $y_p$의 추측 형태 |
| ------------------------------------ | --------------------------------------- |
| $ke^{ax}$ | $Ae^{ax}$ |
| $kx^n$ | $A_n x^n + \cdots + A_1 x + A_0$ |
| $k\cos\omega x$ 또는 $k\sin\omega x$ | $A\cos\omega x + B\sin\omega x$ |
| $ke^{ax}\cos\omega x$ | $e^{ax}(A\cos\omega x + B\sin\omega x)$ |
**주의**: $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 주요 변환쌍
| $f(t)$ | $F(s) = \mathcal{L}\{f(t)\}$ |
| -------------------- | ---------------------------------- |
| $1$ (단위 계단) | $\dfrac{1}{s}$ |
| $t$ | $\dfrac{1}{s^2}$ |
| $t^n$ | $\dfrac{n!}{s^{n+1}}$ |
| $e^{at}$ | $\dfrac{1}{s-a}$ |
| $\sin\omega t$ | $\dfrac{\omega}{s^2+\omega^2}$ |
| $\cos\omega t$ | $\dfrac{s}{s^2+\omega^2}$ |
| $e^{at}\sin\omega t$ | $\dfrac{\omega}{(s-a)^2+\omega^2}$ |
| $e^{at}\cos\omega t$ | $\dfrac{s-a}{(s-a)^2+\omega^2}$ |
| $\delta(t)$ (임펄스) | $1$ |
| $u(t-a)$ (지연 계단) | $\dfrac{e^{-as}}{s}$ |
**유도 예시**: $\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 활용
from scipy.integrate import solve_ivp
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로 해석해 구하기
기호 정의
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 다양한 적분기 비교
오일러법 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으로 자동화
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 솔버 공식 문서](https://docs.scipy.org/doc/scipy/reference/integrate.html)
- [SymPy 라플라스 변환 문서](https://docs.sympy.org/latest/modules/integrals/integrals.html)
현재 단락 (1/426)
공학도라면 피해갈 수 없는 과목이 바로 **공업수학**입니다. 특히 상미분방정식(ODE, Ordinary Differential Equation)은 전자회로 분석, 제어 시스템 설...