10. Higher-Order ODE and Systems
10. Higher-Order ODE and Systems¶
Learning Objectives¶
- Formulate the characteristic equation for nth-order constant-coefficient linear ODEs and find the general solution
- Find particular solutions for non-homogeneous ODEs using the method of undetermined coefficients and variation of parameters
- Express systems of ODEs in vector-matrix form and solve using eigenvalues/eigenvectors and matrix exponentials
- Classify equilibrium points in the phase plane and determine stability
- Apply linearization techniques to nonlinear systems and analyze representative physical/biological models
- Find normal modes of coupled oscillators and understand their connection to Lagrangian mechanics
1. Higher-Order Linear ODEs¶
1.1 nth-Order Constant-Coefficient ODEs¶
The general form of an nth-order constant-coefficient linear ODE is:
$$a_n y^{(n)} + a_{n-1} y^{(n-1)} + \cdots + a_1 y' + a_0 y = f(x)$$
where $a_0, a_1, \ldots, a_n$ are constants. If $f(x) = 0$, it's homogeneous, and if $f(x) \neq 0$, it's non-homogeneous.
Key Principle - Superposition: Solutions of homogeneous equations form a linear space, so the general solution with $n$ linearly independent solutions $y_1, y_2, \ldots, y_n$ is:
$$y_h = c_1 y_1 + c_2 y_2 + \cdots + c_n y_n$$
The general solution of a non-homogeneous equation is:
$$y = y_h + y_p$$
where $y_p$ is one particular solution.
1.2 Characteristic Equation and General Solution¶
Assuming a solution of the form $y = e^{rx}$ for the homogeneous equation gives the characteristic equation:
$$a_n r^n + a_{n-1} r^{n-1} + \cdots + a_1 r + a_0 = 0$$
The general solution is determined by the types of characteristic roots:
| Type of Characteristic Roots | Solution Form |
|---|---|
| Distinct real roots $r_1, r_2, \ldots, r_n$ | $c_1 e^{r_1 x} + c_2 e^{r_2 x} + \cdots$ |
| Repeated root $r$ (multiplicity $m$) | $(c_1 + c_2 x + \cdots + c_m x^{m-1}) e^{rx}$ |
| Complex roots $\alpha \pm i\beta$ | $e^{\alpha x}(c_1 \cos\beta x + c_2 \sin\beta x)$ |
| Repeated complex roots ($\alpha \pm i\beta$, multiplicity $m$) | $e^{\alpha x}\sum_{k=0}^{m-1} x^k (a_k \cos\beta x + b_k \sin\beta x)$ |
Example: 4th-order ODE $y^{(4)} - 5y'' + 4y = 0$
import numpy as np
import sympy as sp
# --- SymPy๋ก ํน์ฑ๋ฐฉ์ ์ ํ๊ธฐ ---
r = sp.Symbol('r')
char_eq = r**4 - 5*r**2 + 4 # ํน์ฑ๋ฐฉ์ ์
roots = sp.solve(char_eq, r)
print(f"ํน์ฑ๊ทผ: {roots}")
# ์ถ๋ ฅ: ํน์ฑ๊ทผ: [-2, -1, 1, 2]
# --- ์ผ๋ฐํด ๊ตฌํ๊ธฐ ---
x = sp.Symbol('x')
y = sp.Function('y')
ode = sp.Eq(y(x).diff(x, 4) - 5*y(x).diff(x, 2) + 4*y(x), 0)
general_sol = sp.dsolve(ode, y(x))
print(f"์ผ๋ฐํด: {general_sol}")
# ์ถ๋ ฅ: y(x) = C1*exp(-2*x) + C2*exp(-x) + C3*exp(x) + C4*exp(2*x)
Example with repeated roots: $y''' - 3y'' + 3y' - y = 0$
Characteristic equation $r^3 - 3r^2 + 3r - 1 = (r-1)^3 = 0$, so $r = 1$ (multiplicity 3)
$$y = (c_1 + c_2 x + c_3 x^2) e^x$$
# ์ค๋ณต๊ทผ์ด ์๋ 3์ฐจ ODE
ode2 = sp.Eq(y(x).diff(x, 3) - 3*y(x).diff(x, 2) + 3*y(x).diff(x) - y(x), 0)
sol2 = sp.dsolve(ode2, y(x))
print(f"์ค๋ณต๊ทผ ์ผ๋ฐํด: {sol2}")
# ์ถ๋ ฅ: y(x) = (C1 + C2*x + C3*x**2)*exp(x)
1.3 Particular Solutions for Non-Homogeneous Problems¶
Method of Undetermined Coefficients¶
Applicable when $f(x)$ consists of polynomials, exponentials, trigonometric functions, or their products.
| Form of $f(x)$ | Assumed form of $y_p$ |
|---|---|
| $P_n(x)$ (nth-degree polynomial) | $A_n x^n + A_{n-1} x^{n-1} + \cdots + A_0$ |
| $e^{\alpha x}$ | $A e^{\alpha x}$ |
| $\cos\beta x$ or $\sin\beta x$ | $A \cos\beta x + B \sin\beta x$ |
| $e^{\alpha x} P_n(x)$ | $e^{\alpha x}(A_n x^n + \cdots + A_0)$ |
If $\alpha$ is a characteristic root, multiply by $x^s$ ($s$ is the multiplicity of $\alpha$).
Variation of Parameters¶
A general method applicable for any form of $f(x)$. For second-order ODE $y'' + p(x)y' + q(x)y = f(x)$:
$$y_p = -y_1 \int \frac{y_2 f}{W} dx + y_2 \int \frac{y_1 f}{W} dx$$
where $W = y_1 y_2' - y_2 y_1'$ is the Wronskian.
# --- ๋น์ ์ฐจ ODE: y'' + y = sec(x) ---
# ๋ฏธ์ ๊ณ์๋ฒ์ผ๋ก๋ ํ ์ ์์ โ ๋งค๊ฐ๋ณ์ ๋ณํ๋ฒ ์ฌ์ฉ
x = sp.Symbol('x')
y = sp.Function('y')
ode_nh = sp.Eq(y(x).diff(x, 2) + y(x), sp.sec(x))
sol_nh = sp.dsolve(ode_nh, y(x))
print(f"๋งค๊ฐ๋ณ์ ๋ณํ๋ฒ ๊ฒฐ๊ณผ: {sol_nh}")
# ๋ก ์คํค์ ์ง์ ๊ณ์ฐ
y1 = sp.cos(x)
y2 = sp.sin(x)
W = y1 * sp.diff(y2, x) - y2 * sp.diff(y1, x)
print(f"๋ก ์คํค์: W = {sp.simplify(W)}")
# ์ถ๋ ฅ: W = 1
# ํน์ํด ๊ณ์ฐ
f_x = sp.sec(x)
yp = -y1 * sp.integrate(y2 * f_x / W, x) + y2 * sp.integrate(y1 * f_x / W, x)
yp_simplified = sp.simplify(yp)
print(f"ํน์ํด: y_p = {yp_simplified}")
2. Systems of ODEs¶
2.1 Vector-Matrix Representation¶
A system of first-order ODEs is expressed concisely in vector-matrix form:
$$\frac{d\mathbf{x}}{dt} = A\mathbf{x} + \mathbf{f}(t)$$
where $\mathbf{x} = \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix}$, and $A$ is an $n \times n$ coefficient matrix.
Important: Any nth-order ODE can be converted to a system of first-order ODEs. For example, $y'' + 3y' + 2y = 0$ becomes:
$$x_1 = y, \quad x_2 = y'$$
$$\frac{d}{dt}\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -2 & -3 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}$$
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# --- 2์ฐจ ODE โ ์ฐ๋ฆฝ 1์ฐจ ODE ๋ณํ ๋ฐ ์์น ํ์ด ---
# y'' + 3y' + 2y = 0, y(0) = 1, y'(0) = 0
def system(t, x):
"""x[0] = y, x[1] = y'"""
return [x[1], -2*x[0] - 3*x[1]]
t_span = (0, 5)
x0 = [1.0, 0.0]
sol = solve_ivp(system, t_span, x0, t_eval=np.linspace(0, 5, 200), method='RK45')
# ํด์ํด์ ๋น๊ต
t_exact = np.linspace(0, 5, 200)
# ํน์ฑ๊ทผ: r = -1, -2 โ y = 2e^{-t} - e^{-2t}
y_exact = 2*np.exp(-t_exact) - np.exp(-2*t_exact)
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].plot(sol.t, sol.y[0], 'b-', label='์์นํด (RK45)')
axes[0].plot(t_exact, y_exact, 'r--', label='ํด์ํด')
axes[0].set_xlabel('t')
axes[0].set_ylabel('y(t)')
axes[0].set_title('ํด ๋น๊ต')
axes[0].legend()
axes[1].plot(sol.y[0], sol.y[1], 'b-')
axes[1].set_xlabel('y')
axes[1].set_ylabel("y'")
axes[1].set_title('์์ ํ๋ฉด ๊ถค์ ')
axes[1].plot(x0[0], x0[1], 'ro', markersize=8, label='์ด๊ธฐ๊ฐ')
axes[1].legend()
plt.tight_layout()
plt.savefig('second_order_ode_solution.png', dpi=150, bbox_inches='tight')
plt.show()
2.2 Solution Using Eigenvalues/Eigenvectors¶
For the homogeneous system $\mathbf{x}' = A\mathbf{x}$, assuming a solution $\mathbf{x} = \mathbf{v} e^{\lambda t}$ gives:
$$A\mathbf{v} = \lambda \mathbf{v}$$
This reduces to finding the eigenvalues $\lambda$ and eigenvectors $\mathbf{v}$ of $A$.
Case 1: Distinct real eigenvalues $\lambda_1, \lambda_2$
$$\mathbf{x}(t) = c_1 \mathbf{v}_1 e^{\lambda_1 t} + c_2 \mathbf{v}_2 e^{\lambda_2 t}$$
Case 2: Complex eigenvalues $\lambda = \alpha \pm i\beta$
$$\mathbf{x}(t) = e^{\alpha t}\left[c_1(\mathbf{a}\cos\beta t - \mathbf{b}\sin\beta t) + c_2(\mathbf{a}\sin\beta t + \mathbf{b}\cos\beta t)\right]$$
where $\mathbf{v} = \mathbf{a} + i\mathbf{b}$.
Case 3: Repeated eigenvalue ($\lambda$, multiplicity 2, 1 eigenvector) - requires generalized eigenvector $\mathbf{w}$:
$$(A - \lambda I)\mathbf{w} = \mathbf{v}$$
$$\mathbf{x}(t) = c_1 \mathbf{v} e^{\lambda t} + c_2 (\mathbf{v} t + \mathbf{w}) e^{\lambda t}$$
import numpy as np
from scipy.linalg import eig
# --- ์ฐ๋ฆฝ ODE์ ๊ณ ์ ๊ฐ/๊ณ ์ ๋ฒกํฐ ํ์ด ---
# x' = Ax, A = [[1, 1], [4, -2]]
A = np.array([[1, 1], [4, -2]])
eigenvalues, eigenvectors = eig(A)
print("๊ณ ์ ๊ฐ:", eigenvalues)
print("๊ณ ์ ๋ฒกํฐ (์ด๋ฒกํฐ):")
print(eigenvectors)
# ๊ณ ์ ๊ฐ: [2, -3]
# ๊ณ ์ ๋ฒกํฐ: v1 = [1, 1], v2 = [1, -4] (์ ๊ทํ๋จ)
# ์ผ๋ฐํด ๊ตฌ์ฑ
t = np.linspace(0, 3, 200)
# ์ด๊ธฐ์กฐ๊ฑด x(0) = [3, 2]์์ c1, c2 ๊ฒฐ์
x0 = np.array([3, 2])
# c1*v1 + c2*v2 = x0 โ [c1, c2] = V^{-1} x0
V = eigenvectors
c = np.linalg.solve(V, x0)
print(f"๊ณ์: c1 = {c[0]:.4f}, c2 = {c[1]:.4f}")
# ํด์ํด ๊ณ์ฐ
x1_sol = c[0] * V[0, 0] * np.exp(eigenvalues[0].real * t) + \
c[1] * V[0, 1] * np.exp(eigenvalues[1].real * t)
x2_sol = c[0] * V[1, 0] * np.exp(eigenvalues[0].real * t) + \
c[1] * V[1, 1] * np.exp(eigenvalues[1].real * t)
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# ์๊ฐ์ ๋ฐ๋ฅธ ํด
axes[0].plot(t, x1_sol.real, 'b-', label='$x_1(t)$')
axes[0].plot(t, x2_sol.real, 'r-', label='$x_2(t)$')
axes[0].set_xlabel('t')
axes[0].set_title('์ฐ๋ฆฝ ODE์ ํด')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# ์์ ํ๋ฉด
axes[1].plot(x1_sol.real, x2_sol.real, 'b-', linewidth=2)
axes[1].plot(x0[0], x0[1], 'ro', markersize=8, label='์ด๊ธฐ๊ฐ')
axes[1].set_xlabel('$x_1$')
axes[1].set_ylabel('$x_2$')
axes[1].set_title('์์ ํ๋ฉด ๊ถค์ ')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('system_ode_eigenvalue.png', dpi=150, bbox_inches='tight')
plt.show()
2.3 Matrix Exponential¶
The solution of the system $\mathbf{x}' = A\mathbf{x}$, $\mathbf{x}(0) = \mathbf{x}_0$ is expressed using the matrix exponential:
$$\mathbf{x}(t) = e^{At} \mathbf{x}_0$$
The matrix exponential extends the series definition of the scalar exponential function to matrices:
$$e^{At} = I + At + \frac{(At)^2}{2!} + \frac{(At)^3}{3!} + \cdots = \sum_{k=0}^{\infty} \frac{(At)^k}{k!}$$
Properties: - $e^{0} = I$ (identity matrix) - $\frac{d}{dt} e^{At} = A e^{At}$ - If $A$ is diagonalizable: $A = PDP^{-1}$ โ $e^{At} = P e^{Dt} P^{-1}$
from scipy.linalg import expm
# --- ํ๋ ฌ ์ง์๋ฅผ ์ด์ฉํ ์ฐ๋ฆฝ ODE ํ์ด ---
A = np.array([[1, 1], [4, -2]])
x0 = np.array([3, 2])
t_vals = np.linspace(0, 3, 200)
solutions = np.array([expm(A * t) @ x0 for t in t_vals])
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(t_vals, solutions[:, 0], 'b-', linewidth=2, label='$x_1(t)$')
ax.plot(t_vals, solutions[:, 1], 'r-', linewidth=2, label='$x_2(t)$')
ax.set_xlabel('t', fontsize=12)
ax.set_ylabel('x(t)', fontsize=12)
ax.set_title('ํ๋ ฌ ์ง์๋ฅผ ์ด์ฉํ ํด', fontsize=14)
ax.legend(fontsize=12)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('matrix_exponential.png', dpi=150, bbox_inches='tight')
plt.show()
# --- SymPy๋ก ์ฌ๋ณผ๋ฆญ ํ๋ ฌ ์ง์ ๊ณ์ฐ ---
t_sym = sp.Symbol('t')
A_sym = sp.Matrix([[1, 1], [4, -2]])
exp_At = sp.exp(A_sym * t_sym) # ์ฌ๋ณผ๋ฆญ ํ๋ ฌ ์ง์
print("e^{At} =")
sp.pprint(exp_At)
Non-homogeneous system $\mathbf{x}' = A\mathbf{x} + \mathbf{f}(t)$ has solution:
$$\mathbf{x}(t) = e^{At}\mathbf{x}_0 + \int_0^t e^{A(t-s)} \mathbf{f}(s) \, ds$$
This is called Duhamel's integral or the variation of constants formula.
3. Phase Plane Analysis¶
For a 2D autonomous system $\mathbf{x}' = A\mathbf{x}$, the eigenvalues of $A$ determine the qualitative behavior of the system.
3.1 Classification of Equilibrium Points (Nodes, Saddles, Spirals, Centers)¶
An equilibrium point $\mathbf{x}^*$ satisfies $A\mathbf{x}^* = \mathbf{0}$. If $A$ is invertible, the origin $\mathbf{x}^* = \mathbf{0}$ is the unique equilibrium.
Classification by eigenvalues $\lambda_1, \lambda_2$ of $A$:
| Type of Eigenvalues | Equilibrium Type | Stability |
|---|---|---|
| $\lambda_1 < \lambda_2 < 0$ (real, negative) | Stable node | Asymptotically stable |
| $\lambda_1 > \lambda_2 > 0$ (real, positive) | Unstable node | Unstable |
| $\lambda_1 < 0 < \lambda_2$ (real, opposite signs) | Saddle point | Unstable |
| $\alpha \pm i\beta$, $\alpha < 0$ | Stable spiral | Asymptotically stable |
| $\alpha \pm i\beta$, $\alpha > 0$ | Unstable spiral | Unstable |
| $\pm i\beta$ (purely imaginary) | Center | Stable (not asymptotically) |
Trace-determinant plane: Using $\tau = \text{tr}(A) = \lambda_1 + \lambda_2$, $\Delta = \det(A) = \lambda_1 \lambda_2$: - $\Delta < 0$: Saddle point - $\Delta > 0$, $\tau < 0$: Stable (node or spiral) - $\Delta > 0$, $\tau > 0$: Unstable (node or spiral) - $\Delta > 0$, $\tau = 0$: Center - $\tau^2 - 4\Delta > 0$: Node, $\tau^2 - 4\Delta < 0$: Spiral
3.2 Stability Criteria¶
Definition: An equilibrium point $\mathbf{x}^*$ is - Stable: For all $\epsilon > 0$, if $\|\mathbf{x}(0) - \mathbf{x}^*\| < \delta$, then $\|\mathbf{x}(t) - \mathbf{x}^*\| < \epsilon$ for all $t > 0$ - Asymptotically stable: Stable and $\mathbf{x}(t) \to \mathbf{x}^*$ as $t \to \infty$ - Unstable: Not stable
Stability criteria for linear systems: - All eigenvalues have negative real parts โ Asymptotically stable - At least one eigenvalue has positive real part โ Unstable - All eigenvalues have real parts โค 0, with at least one zero โ Additional analysis needed
3.3 Drawing Phase Portraits¶
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
def plot_phase_portrait(A, title, ax, xlim=(-3, 3), ylim=(-3, 3)):
"""2D ์ ํ ์์คํ
์ ์์ ์ด์ํ๋ฅผ ๊ทธ๋ฆฐ๋ค."""
# ๋ฒกํฐ์ฅ (streamplot)
x1 = np.linspace(xlim[0], xlim[1], 20)
x2 = np.linspace(ylim[0], ylim[1], 20)
X1, X2 = np.meshgrid(x1, x2)
U = A[0, 0] * X1 + A[0, 1] * X2
V = A[1, 0] * X1 + A[1, 1] * X2
ax.streamplot(X1, X2, U, V, color='steelblue', density=1.5, linewidth=0.8,
arrowsize=1.2)
# ์ฌ๋ฌ ์ด๊ธฐ์กฐ๊ฑด์์ ๊ถค์ ์ถ๊ฐ
for angle in np.linspace(0, 2*np.pi, 12, endpoint=False):
r0 = 2.5
ic = [r0 * np.cos(angle), r0 * np.sin(angle)]
def rhs(t, x):
return A @ x
sol = solve_ivp(rhs, [0, 10], ic, t_eval=np.linspace(0, 10, 500),
method='RK45')
ax.plot(sol.y[0], sol.y[1], 'b-', alpha=0.4, linewidth=0.8)
# ๊ณ ์ ๊ฐ/๊ณ ์ ๋ฒกํฐ
eigenvalues, eigenvectors = np.linalg.eig(A)
eigval_str = ", ".join([f"{ev:.2f}" for ev in eigenvalues])
ax.set_title(f"{title}\n$\\lambda = {eigval_str}$", fontsize=11)
# ๊ณ ์ ๋ฒกํฐ ๋ฐฉํฅ ํ์ (์ค์์ธ ๊ฒฝ์ฐ)
for i in range(2):
if np.isreal(eigenvalues[i]):
v = eigenvectors[:, i].real
ax.arrow(0, 0, v[0]*1.5, v[1]*1.5, head_width=0.15,
head_length=0.1, fc='red', ec='red', linewidth=1.5)
ax.arrow(0, 0, -v[0]*1.5, -v[1]*1.5, head_width=0.15,
head_length=0.1, fc='red', ec='red', linewidth=1.5)
ax.plot(0, 0, 'ko', markersize=6)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_xlim(xlim)
ax.set_ylim(ylim)
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
# --- 4๊ฐ์ง ํํ์ ์ ํ์ ์์ ์ด์ํ ---
fig, axes = plt.subplots(2, 2, figsize=(12, 12))
# (a) ์์ ๋
ธ๋: ๊ณ ์ ๊ฐ ๋ชจ๋ ์์
A_stable_node = np.array([[-2, 0], [0, -1]])
plot_phase_portrait(A_stable_node, '์์ ๋
ธ๋ (Stable Node)', axes[0, 0])
# (b) ์์ฅ์ : ๊ณ ์ ๊ฐ์ด ์ด๋ถํธ
A_saddle = np.array([[1, 0], [0, -2]])
plot_phase_portrait(A_saddle, '์์ฅ์ (Saddle Point)', axes[0, 1])
# (c) ์์ ์์ฉ๋์ด: ๋ณต์ ๊ณ ์ ๊ฐ, ์ค์๋ถ ์์
A_spiral = np.array([[-0.5, 2], [-2, -0.5]])
plot_phase_portrait(A_spiral, '์์ ์์ฉ๋์ด (Stable Spiral)', axes[1, 0])
# (d) ์ค์ฌ: ์ํ์ ๊ณ ์ ๊ฐ
A_center = np.array([[0, 1], [-4, 0]])
plot_phase_portrait(A_center, '์ค์ฌ (Center)', axes[1, 1])
plt.tight_layout()
plt.savefig('phase_portraits.png', dpi=150, bbox_inches='tight')
plt.show()
4. Introduction to Nonlinear Systems¶
4.1 Linearization¶
The behavior near an equilibrium point $\mathbf{x}^*$ of a nonlinear autonomous system $\mathbf{x}' = \mathbf{F}(\mathbf{x})$ is analyzed through linearization using the Jacobian matrix:
$$\mathbf{x}' \approx J(\mathbf{x}^*) (\mathbf{x} - \mathbf{x}^*)$$
where the Jacobian matrix is:
$$J = \begin{pmatrix} \frac{\partial F_1}{\partial x_1} & \frac{\partial F_1}{\partial x_2} \\ \frac{\partial F_2}{\partial x_1} & \frac{\partial F_2}{\partial x_2} \end{pmatrix}_{\mathbf{x} = \mathbf{x}^*}$$
Hartman-Grobman theorem: If an equilibrium point is hyperbolic (i.e., all eigenvalues of the Jacobian have nonzero real parts), the topological behavior of the nonlinear system is topologically equivalent to the linearized system.
Caution: Centers are not hyperbolic, so linearization alone cannot determine the exact behavior of the nonlinear system.
4.2 Lotka-Volterra Equations¶
A classical nonlinear system modeling predator-prey interaction:
$$\frac{dx}{dt} = \alpha x - \beta x y \quad \text{(prey)}$$
$$\frac{dy}{dt} = -\gamma y + \delta x y \quad \text{(predator)}$$
Equilibrium points: 1. $(0, 0)$ - Extinction (trivial) 2. $(\gamma/\delta, \alpha/\beta)$ - Coexistence
Jacobian at the coexistence equilibrium:
$$J = \begin{pmatrix} 0 & -\beta\gamma/\delta \\ \delta\alpha/\beta & 0 \end{pmatrix}$$
Eigenvalues are $\lambda = \pm i\sqrt{\alpha\gamma}$ (purely imaginary), so linear analysis gives a center. Nonlinear analysis also confirms closed orbits (periodic motion) due to conserved quantity $V = \delta x - \gamma \ln x + \beta y - \alpha \ln y$.
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# --- ๋กํธ์นด-๋ณผํ
๋ผ ์๋ฎฌ๋ ์ด์
---
alpha, beta, gamma, delta = 1.0, 0.5, 0.75, 0.25
def lotka_volterra(t, z):
x, y = z
dxdt = alpha * x - beta * x * y
dydt = -gamma * y + delta * x * y
return [dxdt, dydt]
# ์ฌ๋ฌ ์ด๊ธฐ์กฐ๊ฑด
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
initial_conditions = [[2, 1], [4, 2], [1, 3], [6, 1]]
colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728']
for ic, color in zip(initial_conditions, colors):
sol = solve_ivp(lotka_volterra, [0, 30], ic,
t_eval=np.linspace(0, 30, 1000), method='RK45',
rtol=1e-10, atol=1e-12)
# ์๊ฐ ์์ญ
axes[0].plot(sol.t, sol.y[0], '-', color=color, alpha=0.8,
label=f'ํผ์์ ({ic[0]},{ic[1]})')
axes[0].plot(sol.t, sol.y[1], '--', color=color, alpha=0.8,
label=f'ํฌ์์ ({ic[0]},{ic[1]})')
# ์์ ํ๋ฉด
axes[1].plot(sol.y[0], sol.y[1], '-', color=color, linewidth=1.5,
label=f'IC=({ic[0]},{ic[1]})')
axes[1].plot(ic[0], ic[1], 'o', color=color, markersize=6)
# ํํ์ ํ์
x_eq, y_eq = gamma / delta, alpha / beta
axes[1].plot(x_eq, y_eq, 'k*', markersize=15, zorder=5, label='ํํ์ ')
axes[0].set_xlabel('์๊ฐ t')
axes[0].set_ylabel('๊ฐ์ฒด์')
axes[0].set_title('๋กํธ์นด-๋ณผํ
๋ผ: ์๊ฐ ์์ญ')
axes[0].legend(fontsize=7, ncol=2)
axes[0].grid(True, alpha=0.3)
axes[1].set_xlabel('ํผ์์ x')
axes[1].set_ylabel('ํฌ์์ y')
axes[1].set_title('๋กํธ์นด-๋ณผํ
๋ผ: ์์ ํ๋ฉด')
axes[1].legend(fontsize=8)
axes[1].grid(True, alpha=0.3)
# ๋ฒกํฐ์ฅ (streamplot)
x_range = np.linspace(0.2, 8, 20)
y_range = np.linspace(0.2, 5, 20)
X, Y = np.meshgrid(x_range, y_range)
U = alpha * X - beta * X * Y
V = -gamma * Y + delta * X * Y
speed = np.sqrt(U**2 + V**2)
axes[2].streamplot(X, Y, U, V, color=speed, cmap='coolwarm', density=1.5,
linewidth=0.8)
axes[2].plot(x_eq, y_eq, 'k*', markersize=15, label='ํํ์ ')
axes[2].set_xlabel('ํผ์์ x')
axes[2].set_ylabel('ํฌ์์ y')
axes[2].set_title('๋กํธ์นด-๋ณผํ
๋ผ: ๋ฒกํฐ์ฅ')
axes[2].legend()
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('lotka_volterra.png', dpi=150, bbox_inches='tight')
plt.show()
4.3 Van der Pol Oscillator¶
An oscillator with nonlinear damping, widely used in electronics and biological rhythm modeling:
$$\ddot{x} - \mu(1 - x^2)\dot{x} + x = 0$$
- When $|x| < 1$: Negative damping (energy supply) โ Amplitude increases
- When $|x| > 1$: Positive damping (energy dissipation) โ Amplitude decreases
This competition leads to a limit cycle when $\mu > 0$ (around $|x| \approx 2$).
System form: $x_1 = x$, $x_2 = \dot{x}$
$$\dot{x}_1 = x_2, \quad \dot{x}_2 = \mu(1 - x_1^2) x_2 - x_1$$
Jacobian at the origin $(0, 0)$:
$$J = \begin{pmatrix} 0 & 1 \\ -1 & \mu \end{pmatrix}$$
If $\mu > 0$, then $\text{tr}(J) = \mu > 0$, so the origin is unstable.
# --- ๋ฐ๋ฐ๋ฅดํด ์ง๋์ ์๋ฎฌ๋ ์ด์
---
def van_der_pol(t, z, mu):
x, xdot = z
return [xdot, mu * (1 - x**2) * xdot - x]
fig, axes = plt.subplots(2, 3, figsize=(16, 10))
mu_values = [0.1, 1.0, 5.0]
for idx, mu in enumerate(mu_values):
# ์ฌ๋ฌ ์ด๊ธฐ์กฐ๊ฑด
ics = [[0.1, 0], [4, 0], [0, 5], [2, -3]]
for ic in ics:
sol = solve_ivp(van_der_pol, [0, 50], ic,
args=(mu,), t_eval=np.linspace(0, 50, 2000),
method='RK45', rtol=1e-10, atol=1e-12)
# ์๊ฐ ์์ญ
axes[0, idx].plot(sol.t, sol.y[0], linewidth=0.8, alpha=0.8)
# ์์ ํ๋ฉด
axes[1, idx].plot(sol.y[0], sol.y[1], linewidth=0.8, alpha=0.8)
axes[0, idx].set_xlabel('t')
axes[0, idx].set_ylabel('x(t)')
axes[0, idx].set_title(f'Van der Pol ($\\mu$ = {mu}): ์๊ฐ ์์ญ')
axes[0, idx].grid(True, alpha=0.3)
axes[1, idx].set_xlabel('x')
axes[1, idx].set_ylabel('$\\dot{x}$')
axes[1, idx].set_title(f'Van der Pol ($\\mu$ = {mu}): ์์ ํ๋ฉด')
axes[1, idx].plot(0, 0, 'ro', markersize=5)
axes[1, idx].grid(True, alpha=0.3)
axes[1, idx].set_aspect('equal')
plt.tight_layout()
plt.savefig('van_der_pol.png', dpi=150, bbox_inches='tight')
plt.show()
Observations: - $\mu \to 0$: Nearly sinusoidal oscillation (weak nonlinearity) - $\mu = 1$: Clear convergence to limit cycle - $\mu \gg 1$: Relaxation oscillation - Alternating slow/fast segments
5. Physics Applications¶
5.1 Coupled Oscillators¶
Consider the motion of two masses connected by springs:
๋ฒฝ โโโ k โโโ [mโ] โโโ k_c โโโ [mโ] โโโ k โโโ ๋ฒฝ
Newton's equations of motion:
$$m_1 \ddot{x}_1 = -k x_1 - k_c (x_1 - x_2)$$
$$m_2 \ddot{x}_2 = -k x_2 - k_c (x_2 - x_1)$$
For the symmetric case $m_1 = m_2 = m$, in matrix form:
$$m \begin{pmatrix} \ddot{x}_1 \\ \ddot{x}_2 \end{pmatrix} = -\begin{pmatrix} k + k_c & -k_c \\ -k_c & k + k_c \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}$$
5.2 Normal Modes¶
Normal modes are special motion patterns where all particles oscillate at the same frequency.
Substituting $\mathbf{x}(t) = \mathbf{u} e^{i\omega t}$ gives an eigenvalue problem:
$$K\mathbf{u} = \omega^2 M\mathbf{u}$$
or $M^{-1}K\mathbf{u} = \omega^2 \mathbf{u}$
Normal modes for the symmetric case:
| Mode | Frequency | Pattern | Physical Meaning |
|---|---|---|---|
| Mode 1 (in-phase) | $\omega_1 = \sqrt{k/m}$ | $\mathbf{u}_1 = (1, 1)^T$ | Both masses move together |
| Mode 2 (out-of-phase) | $\omega_2 = \sqrt{(k + 2k_c)/m}$ | $\mathbf{u}_2 = (1, -1)^T$ | Masses move opposite |
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# --- ๊ฒฐํฉ ์ง๋์์ ์ ๊ท ๋ชจ๋ ---
m = 1.0 # ์ง๋
k = 1.0 # ๋ฒฝ ์คํ๋ง ์์
kc = 0.5 # ๊ฒฐํฉ ์คํ๋ง ์์
# ๊ฐ์ฑ ํ๋ ฌ๊ณผ ์ง๋ ํ๋ ฌ
K = np.array([[k + kc, -kc],
[-kc, k + kc]])
M = np.array([[m, 0],
[0, m]])
# ์ ๊ท ๋ชจ๋ (์ผ๋ฐํ๋ ๊ณ ์ ๊ฐ ๋ฌธ์ )
from scipy.linalg import eigh
omega_sq, modes = eigh(K, M)
omega = np.sqrt(omega_sq)
print("์ ๊ท ์ฃผํ์:")
for i, w in enumerate(omega):
print(f" omega_{i+1} = {w:.4f} rad/s (f = {w/(2*np.pi):.4f} Hz)")
print(f"\n์ ๊ท ๋ชจ๋ ๋ฒกํฐ:")
print(f" ๋ชจ๋ 1 (๋์์): {modes[:, 0]}")
print(f" ๋ชจ๋ 2 (์ญ์์): {modes[:, 1]}")
# ์๋ฎฌ๋ ์ด์
: ์ด๊ธฐ์ ์ง๋ 1๋ง ๋ณ์
def coupled_oscillator(t, z):
x1, x2, v1, v2 = z
a1 = (-k * x1 - kc * (x1 - x2)) / m
a2 = (-k * x2 - kc * (x2 - x1)) / m
return [v1, v2, a1, a2]
# ์ด๊ธฐ์กฐ๊ฑด: x1(0) = 1, x2(0) = 0 (๋นํธ ํ์ ๊ด์ฐฐ)
sol = solve_ivp(coupled_oscillator, [0, 60], [1, 0, 0, 0],
t_eval=np.linspace(0, 60, 2000), method='RK45',
rtol=1e-12, atol=1e-14)
fig, axes = plt.subplots(3, 1, figsize=(14, 10))
# ๊ฐ ์ง๋์ ๋ณ์
axes[0].plot(sol.t, sol.y[0], 'b-', linewidth=1, label='$x_1(t)$ (์ง๋ 1)')
axes[0].plot(sol.t, sol.y[1], 'r-', linewidth=1, label='$x_2(t)$ (์ง๋ 2)')
axes[0].set_xlabel('์๊ฐ t')
axes[0].set_ylabel('๋ณ์')
axes[0].set_title('๊ฒฐํฉ ์ง๋์: ์๋์ง ์ ๋ฌ (๋นํธ ํ์)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# ์ ๊ท ์ขํ
q1 = (sol.y[0] + sol.y[1]) / np.sqrt(2) # ๋์์ ๋ชจ๋
q2 = (sol.y[0] - sol.y[1]) / np.sqrt(2) # ์ญ์์ ๋ชจ๋
axes[1].plot(sol.t, q1, 'g-', linewidth=1, label='$q_1(t)$ (๋์์ ๋ชจ๋)')
axes[1].plot(sol.t, q2, 'm-', linewidth=1, label='$q_2(t)$ (์ญ์์ ๋ชจ๋)')
axes[1].set_xlabel('์๊ฐ t')
axes[1].set_ylabel('์ ๊ท ์ขํ')
axes[1].set_title('์ ๊ท ์ขํ: ๊ฐ ๋ชจ๋์ ๋
๋ฆฝ์ ์ง๋')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
# ์๋์ง ๊ตํ
E1 = 0.5 * m * sol.y[2]**2 + 0.5 * k * sol.y[0]**2
E2 = 0.5 * m * sol.y[3]**2 + 0.5 * k * sol.y[1]**2
E_coupling = 0.5 * kc * (sol.y[0] - sol.y[1])**2
axes[2].plot(sol.t, E1, 'b-', linewidth=1, label='์ง๋ 1 ์๋์ง')
axes[2].plot(sol.t, E2, 'r-', linewidth=1, label='์ง๋ 2 ์๋์ง')
axes[2].plot(sol.t, E1 + E2 + E_coupling, 'k--', linewidth=0.8,
label='์ด ์๋์ง', alpha=0.5)
axes[2].set_xlabel('์๊ฐ t')
axes[2].set_ylabel('์๋์ง')
axes[2].set_title('์๋์ง ๊ตํ: ๋นํธ ์ง๋์ = $|\\omega_2 - \\omega_1|/2$')
axes[2].legend()
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('coupled_oscillators.png', dpi=150, bbox_inches='tight')
plt.show()
# ๋นํธ ์ง๋์ ๊ณ์ฐ
omega_beat = abs(omega[1] - omega[0]) / 2
T_beat = 2 * np.pi / omega_beat if omega_beat > 0 else float('inf')
print(f"\n๋นํธ ์ง๋์: {omega_beat:.4f} rad/s")
print(f"๋นํธ ์ฃผ๊ธฐ: {T_beat:.2f} s")
5.3 Systems of ODEs from Lagrangian Mechanics¶
The double pendulum is a representative nonlinear system from Lagrangian mechanics.
From the Lagrangian $L = T - V$, the Euler-Lagrange equations are:
$$\frac{d}{dt} \frac{\partial L}{\partial \dot{q}_i} - \frac{\partial L}{\partial q_i} = 0$$
Small angle approximation linearizes the double pendulum equations:
$$\begin{pmatrix} (m_1 + m_2) l_1 & m_2 l_2 \\ l_1 & l_2 \end{pmatrix} \begin{pmatrix} \ddot{\theta}_1 \\ \ddot{\theta}_2 \end{pmatrix} = -g \begin{pmatrix} (m_1 + m_2) \theta_1 \\ \theta_2 \end{pmatrix}$$
This reduces to a generalized eigenvalue problem $K\mathbf{u} = \omega^2 M\mathbf{u}$.
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# --- ์ด์ค ์ง์ (๋น์ ํ ์ ์ฒด ๋ฐฉ์ ์) ---
g = 9.81
m1, m2 = 1.0, 1.0
l1, l2 = 1.0, 1.0
def double_pendulum(t, z):
th1, th2, w1, w2 = z
delta = th1 - th2
den1 = (m1 + m2) * l1 - m2 * l1 * np.cos(delta)**2
den2 = l2 / l1 * den1
dw1 = (-m2 * l1 * w1**2 * np.sin(delta) * np.cos(delta)
- m2 * l2 * w2**2 * np.sin(delta)
- (m1 + m2) * g * np.sin(th1)
+ m2 * g * np.sin(th2) * np.cos(delta)) / den1
dw2 = (m2 * l2 * w2**2 * np.sin(delta) * np.cos(delta)
+ (m1 + m2) * l1 * w1**2 * np.sin(delta)
+ (m1 + m2) * g * np.sin(th1) * np.cos(delta)
- (m1 + m2) * g * np.sin(th2)) / den2
return [w1, w2, dw1, dw2]
# ๋ ๊ฐ์ง ์ฝ๊ฐ ๋ค๋ฅธ ์ด๊ธฐ์กฐ๊ฑด โ ์นด์ค์ค ๋ฏผ๊ฐ์ฑ
t_span = (0, 20)
t_eval = np.linspace(0, 20, 5000)
ic1 = [np.pi/2, np.pi/2, 0, 0]
ic2 = [np.pi/2 + 0.001, np.pi/2, 0, 0] # theta1์ 0.001 rad ์ฐจ์ด
sol1 = solve_ivp(double_pendulum, t_span, ic1, t_eval=t_eval,
method='RK45', rtol=1e-12, atol=1e-14)
sol2 = solve_ivp(double_pendulum, t_span, ic2, t_eval=t_eval,
method='RK45', rtol=1e-12, atol=1e-14)
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# theta1 ๋น๊ต
axes[0, 0].plot(sol1.t, sol1.y[0], 'b-', linewidth=0.8, label='IC 1')
axes[0, 0].plot(sol2.t, sol2.y[0], 'r--', linewidth=0.8, label='IC 2')
axes[0, 0].set_xlabel('t')
axes[0, 0].set_ylabel('$\\theta_1$')
axes[0, 0].set_title('$\\theta_1(t)$: ์ด๊ธฐ์กฐ๊ฑด ๋ฏผ๊ฐ์ฑ (์นด์ค์ค)')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)
# theta2 ๋น๊ต
axes[0, 1].plot(sol1.t, sol1.y[1], 'b-', linewidth=0.8, label='IC 1')
axes[0, 1].plot(sol2.t, sol2.y[1], 'r--', linewidth=0.8, label='IC 2')
axes[0, 1].set_xlabel('t')
axes[0, 1].set_ylabel('$\\theta_2$')
axes[0, 1].set_title('$\\theta_2(t)$: ์ด๊ธฐ์กฐ๊ฑด ๋ฏผ๊ฐ์ฑ (์นด์ค์ค)')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
# ์์ ๊ณต๊ฐ (theta1 vs omega1)
axes[1, 0].plot(sol1.y[0], sol1.y[2], 'b-', linewidth=0.3, alpha=0.7)
axes[1, 0].set_xlabel('$\\theta_1$')
axes[1, 0].set_ylabel('$\\dot{\\theta}_1$')
axes[1, 0].set_title('์์ ๊ณต๊ฐ: $(\\theta_1, \\dot{\\theta}_1)$')
axes[1, 0].grid(True, alpha=0.3)
# ์ง์ ๋ ๊ถค์
x1 = l1 * np.sin(sol1.y[0])
y1 = -l1 * np.cos(sol1.y[0])
x2 = x1 + l2 * np.sin(sol1.y[1])
y2 = y1 - l2 * np.cos(sol1.y[1])
axes[1, 1].plot(x2, y2, 'b-', linewidth=0.2, alpha=0.5)
axes[1, 1].set_xlabel('x')
axes[1, 1].set_ylabel('y')
axes[1, 1].set_title('์ด์ค ์ง์ ๋์ ๊ถค์ ')
axes[1, 1].set_aspect('equal')
axes[1, 1].grid(True, alpha=0.3)
plt.suptitle('์ด์ค ์ง์ (Double Pendulum) - ์นด์ค์ค์ ์ญํ', fontsize=14, y=1.02)
plt.tight_layout()
plt.savefig('double_pendulum.png', dpi=150, bbox_inches='tight')
plt.show()
# --- ์์งํญ ์ ๊ท ๋ชจ๋ ๋ถ์ ---
M_mat = np.array([[(m1 + m2) * l1, m2 * l2],
[l1, l2]])
K_mat = np.array([[(m1 + m2) * g, 0],
[0, g]])
# ์ผ๋ฐํ๋ ๊ณ ์ ๊ฐ ๋ฌธ์ : K u = omega^2 M u
from scipy.linalg import eig
eigenvalues, eigenvectors = eig(K_mat, M_mat)
omega_normal = np.sqrt(eigenvalues.real)
omega_normal = np.sort(omega_normal)
print("\n์ด์ค ์ง์ (์์งํญ) ์ ๊ท ์ฃผํ์:")
print(f" omega_1 = {omega_normal[0]:.4f} rad/s (๋์์ ๋ชจ๋)")
print(f" omega_2 = {omega_normal[1]:.4f} rad/s (์ญ์์ ๋ชจ๋)")
print(f" ๋น์จ omega_2/omega_1 = {omega_normal[1]/omega_normal[0]:.4f}")
Exercises¶
Basic Problems¶
Problem 1: Find the general solution of the following ODEs. - (a) $y^{(4)} + 4y'' = 0$ - (b) $y''' - y = 0$ - (c) $y'' + 4y' + 13y = 0$
Solution Hint
(a) Characteristic equation: $r^4 + 4r^2 = r^2(r^2 + 4) = 0$ โ $r = 0, 0, \pm 2i$ General solution: $y = c_1 + c_2 x + c_3 \cos 2x + c_4 \sin 2x$ (b) $r^3 - 1 = (r-1)(r^2+r+1) = 0$ โ $r = 1, -\frac{1}{2} \pm \frac{\sqrt{3}}{2}i$ (c) $r^2 + 4r + 13 = 0$ โ $r = -2 \pm 3i$# ๊ฒ์ฆ
import sympy as sp
x = sp.Symbol('x')
y = sp.Function('y')
# (a)
sol_a = sp.dsolve(y(x).diff(x, 4) + 4*y(x).diff(x, 2), y(x))
print(f"(a): {sol_a}")
# (b)
sol_b = sp.dsolve(y(x).diff(x, 3) - y(x), y(x))
print(f"(b): {sol_b}")
# (c)
sol_c = sp.dsolve(y(x).diff(x, 2) + 4*y(x).diff(x) + 13*y(x), y(x))
print(f"(c): {sol_c}")
Problem 2: Find the general solution of the following system using eigenvalues/eigenvectors.
$$\frac{d}{dt}\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 3 & -2 \\ 2 & -2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}$$
Solution Hint
Characteristic equation: $\lambda^2 - \lambda - 2 = 0$ โ $\lambda_1 = 2$, $\lambda_2 = -1$ Finding eigenvectors for each eigenvalue and constructing the general solution: $$\mathbf{x}(t) = c_1 \begin{pmatrix} 2 \\ 1 \end{pmatrix} e^{2t} + c_2 \begin{pmatrix} 1 \\ 2 \end{pmatrix} e^{-t}$$Applied Problems¶
Problem 3 (Coupled oscillators): Three masses $m$ are connected by identical springs with spring constant $k$ (both ends fixed to walls). Find the normal frequencies and normal modes.
Solution Hint
Stiffness matrix: $$K = k \begin{pmatrix} 2 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 2 \end{pmatrix}$$ Eigenvalues: $\omega_n^2 = \frac{k}{m}(2 - 2\cos\frac{n\pi}{4})$, $n = 1, 2, 3$import numpy as np
from scipy.linalg import eigh
k, m_val = 1.0, 1.0
K = k * np.array([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])
M = m_val * np.eye(3)
omega_sq, modes = eigh(K, M)
print("์ ๊ท ์ฃผํ์:", np.sqrt(omega_sq))
print("์ ๊ท ๋ชจ๋:\n", modes)
Problem 4 (Nonlinear analysis): Find all equilibrium points of the following nonlinear system and classify the stability at each point.
$$\dot{x} = x(3 - x - 2y), \quad \dot{y} = y(2 - x - y)$$
Solution Hint
Equilibrium points: $(0,0)$, $(3,0)$, $(0,2)$, $(1,1)$ โ points where $\dot{x}=0$, $\dot{y}=0$ simultaneously. Jacobian at each equilibrium: $$J = \begin{pmatrix} 3 - 2x - 2y & -2x \\ -y & 2 - x - 2y \end{pmatrix}$$import sympy as sp
x, y = sp.symbols('x y')
F1 = x * (3 - x - 2*y)
F2 = y * (2 - x - y)
# ํํ์
eq_pts = sp.solve([F1, F2], [x, y])
print("ํํ์ :", eq_pts)
# ์ผ์ฝ๋น ํ๋ ฌ
J = sp.Matrix([[sp.diff(F1, x), sp.diff(F1, y)],
[sp.diff(F2, x), sp.diff(F2, y)]])
for pt in eq_pts:
J_at = J.subs([(x, pt[0]), (y, pt[1])])
eigenvals = J_at.eigenvals()
print(f"\nํํ์ {pt}: ๊ณ ์ ๊ฐ = {eigenvals}")
Problem 5 (Phase portraits): Analyze how the phase portrait of the following system changes with parameter $\mu$ (Hopf bifurcation):
$$\dot{x} = \mu x - y - x(x^2 + y^2), \quad \dot{y} = x + \mu y - y(x^2 + y^2)$$
Solution Hint
Eigenvalues of the Jacobian at the origin: $\lambda = \mu \pm i$ - $\mu < 0$: Stable spiral - $\mu = 0$: Bifurcation point (non-hyperbolic) - $\mu > 0$: Unstable spiral + stable limit cycle (radius $r = \sqrt{\mu}$) In polar coordinates $(r, \theta)$: $\dot{r} = r(\mu - r^2)$, $\dot{\theta} = 1$import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
def hopf_system(t, z, mu):
x, y = z
r_sq = x**2 + y**2
return [mu*x - y - x*r_sq, x + mu*y - y*r_sq]
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
mu_vals = [-0.5, 0.0, 0.5]
for idx, mu in enumerate(mu_vals):
for r0 in [0.1, 0.3, 0.8, 1.2]:
for theta0 in [0, np.pi/2, np.pi]:
ic = [r0*np.cos(theta0), r0*np.sin(theta0)]
sol = solve_ivp(hopf_system, [0, 30], ic, args=(mu,),
t_eval=np.linspace(0, 30, 2000),
method='RK45', rtol=1e-10)
axes[idx].plot(sol.y[0], sol.y[1], linewidth=0.5, alpha=0.6)
axes[idx].set_title(f'$\\mu = {mu}$')
axes[idx].set_xlabel('x')
axes[idx].set_ylabel('y')
axes[idx].set_aspect('equal')
axes[idx].grid(True, alpha=0.3)
if mu > 0:
theta = np.linspace(0, 2*np.pi, 100)
axes[idx].plot(np.sqrt(mu)*np.cos(theta), np.sqrt(mu)*np.sin(theta),
'r-', linewidth=2, label=f'ํ๊ณ ์ฃผ๊ธฐ $r=\\sqrt{{{mu}}}$')
axes[idx].legend()
plt.suptitle('ํธํ ๋ถ๊ธฐ (Hopf Bifurcation)', fontsize=14)
plt.tight_layout()
plt.savefig('hopf_bifurcation.png', dpi=150, bbox_inches='tight')
plt.show()
References¶
Textbooks¶
- Boas, M. L. (2005). Mathematical Methods in the Physical Sciences, 3rd ed., Chapter 8. Wiley.
- Strogatz, S. H. (2018). Nonlinear Dynamics and Chaos, 2nd ed. CRC Press.
- Standard textbook for nonlinear dynamics and phase plane analysis
- Hirsch, M. W., Smale, S., & Devaney, R. L. (2013). Differential Equations, Dynamical Systems, and an Introduction to Chaos, 3rd ed. Academic Press.
- Arfken, G. B. et al. (2012). Mathematical Methods for Physicists, 7th ed., Chapter 10. Academic Press.
Online Resources¶
- MIT OCW 18.03: Differential Equations (Arthur Mattuck)
- 3Blue1Brown: Differential equations, studying the unsolvable
- Steve Brunton (YouTube): Data-Driven Dynamical Systems series
Core Library Documentation¶
- SciPy
solve_ivp: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html - SymPy
dsolve: https://docs.sympy.org/latest/modules/solvers/ode.html - Matplotlib
streamplot: https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.streamplot.html
Next Lesson¶
In 09. Series Solutions and Special Functions, we'll cover series solutions to ODEs using the Frobenius method, and study the properties and applications of the most important special functions in physics (Bessel functions, Legendre polynomials, Hermite functions, Laguerre functions).