01. Infinite Series and Convergence

Previous Mathematical Methods in the Physical Sciences - Overview
๐Ÿ“‹ Topic list
Next 02. Complex Numbers

01. Infinite Series and Convergence

Learning Objectives

  • Define convergence and divergence of series and determine convergence using partial sums
  • Apply major convergence tests: comparison test, ratio test, root test, integral test, alternating series test
  • Find the radius of convergence for power series and perform term-by-term differentiation/integration
  • Approximate functions using Taylor series and Maclaurin series and apply them to physics problems
  • Understand the concepts of asymptotic series and Stirling's approximation and utilize them in physical applications

1. Basic Concepts of Series

1.1 Sequences and Series

A sequence is a function from natural numbers to real (or complex) numbers:

$$a_1, a_2, a_3, \ldots, a_n, \ldots$$

A sequence $\{a_n\}$ converges to a specific value $L$ if:

$$\lim_{n \to \infty} a_n = L$$

This means that for sufficiently large $n$, $a_n$ becomes arbitrarily close to $L$.

A series is the sum of the terms of a sequence:

$$S = \sum_{n=1}^{\infty} a_n = a_1 + a_2 + a_3 + \cdots$$

Necessary condition for series convergence: $\lim_{n \to \infty} a_n = 0$

Note: This is only a necessary condition, not sufficient. The harmonic series $\sum \frac{1}{n}$ has $a_n \to 0$ but diverges.

1.2 Partial Sums and Convergence

A partial sum $S_N$ is the sum of the first $N$ terms of the series:

$$S_N = \sum_{n=1}^{N} a_n = a_1 + a_2 + \cdots + a_N$$

Convergence of a series: If the sequence of partial sums $\{S_N\}$ converges to a finite value $S$, we say the series $\sum a_n$ converges to $S$.

import numpy as np
import matplotlib.pyplot as plt

# ์˜ˆ์ œ: ๊ธฐํ•˜๊ธ‰์ˆ˜ (geometric series)์˜ ๋ถ€๋ถ„ํ•ฉ
# sum_{n=0}^{inf} r^n = 1/(1-r)  (|r| < 1)

def geometric_partial_sums(r, N_max):
    """๊ธฐํ•˜๊ธ‰์ˆ˜์˜ ๋ถ€๋ถ„ํ•ฉ์„ ๊ณ„์‚ฐํ•ฉ๋‹ˆ๋‹ค."""
    n_values = np.arange(N_max + 1)
    terms = r ** n_values
    partial_sums = np.cumsum(terms)
    return n_values, partial_sums

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# ์ˆ˜๋ ดํ•˜๋Š” ๊ฒฝ์šฐ: |r| < 1
r_conv = 0.5
n_vals, S_n = geometric_partial_sums(r_conv, 20)
exact = 1 / (1 - r_conv)
axes[0].plot(n_vals, S_n, 'bo-', markersize=4, label=f'S_N (r={r_conv})')
axes[0].axhline(y=exact, color='r', linestyle='--', label=f'S = {exact:.4f}')
axes[0].set_xlabel('N')
axes[0].set_ylabel('S_N')
axes[0].set_title(f'์ˆ˜๋ ดํ•˜๋Š” ๊ธฐํ•˜๊ธ‰์ˆ˜ (r = {r_conv})')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# ์กฐํ™”๊ธ‰์ˆ˜ vs p-๊ธ‰์ˆ˜ ๋น„๊ต
N_max = 100
n = np.arange(1, N_max + 1)
harmonic = np.cumsum(1.0 / n)           # p = 1 (๋ฐœ์‚ฐ)
p2_series = np.cumsum(1.0 / n**2)       # p = 2 (์ˆ˜๋ ด)
p3_series = np.cumsum(1.0 / n**3)       # p = 3 (์ˆ˜๋ ด)

axes[1].plot(n, harmonic, 'r-', label='p=1 (์กฐํ™”๊ธ‰์ˆ˜, ๋ฐœ์‚ฐ)')
axes[1].plot(n, p2_series, 'b-', label=f'p=2 (์ˆ˜๋ ด, pi^2/6 = {np.pi**2/6:.4f})')
axes[1].plot(n, p3_series, 'g-', label='p=3 (์ˆ˜๋ ด)')
axes[1].axhline(y=np.pi**2/6, color='b', linestyle='--', alpha=0.5)
axes[1].set_xlabel('N')
axes[1].set_ylabel('S_N')
axes[1].set_title('p-๊ธ‰์ˆ˜ ๋น„๊ต: sum(1/n^p)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('series_convergence.png', dpi=150, bbox_inches='tight')
plt.show()

Key series results:

Series Convergence condition Sum (when convergent)
Geometric series $\sum r^n$ $|r| < 1$ $\frac{1}{1-r}$
p-series $\sum \frac{1}{n^p}$ $p > 1$ $\zeta(p)$
Harmonic series $\sum \frac{1}{n}$ Diverges -

2. Convergence Tests

2.1 Comparison Test

If $0 \leq a_n \leq b_n$ (for sufficiently large $n$), then:

  • $\sum b_n$ converges $\Rightarrow$ $\sum a_n$ converges
  • $\sum a_n$ diverges $\Rightarrow$ $\sum b_n$ diverges

Limit Comparison Test: If $a_n > 0$, $b_n > 0$ and

$$\lim_{n \to \infty} \frac{a_n}{b_n} = L \quad (0 < L < \infty)$$

then $\sum a_n$ and $\sum b_n$ either both converge or both diverge.

2.2 Ratio Test

$$\rho = \lim_{n \to \infty} \left|\frac{a_{n+1}}{a_n}\right|$$

  • $\rho < 1$ implies absolutely convergent
  • $\rho > 1$ implies divergent
  • $\rho = 1$ implies inconclusive (use another test)

Physics Tip: The ratio test is particularly useful for series containing factorials or exponential functions.

2.3 Root Test

$$\rho = \lim_{n \to \infty} |a_n|^{1/n}$$

  • $\rho < 1$ implies absolutely convergent
  • $\rho > 1$ implies divergent
  • $\rho = 1$ implies inconclusive

2.4 Integral Test

If $f(x)$ is positive, continuous, and decreasing on $[1, \infty)$ with $f(n) = a_n$, then:

$$\sum_{n=1}^{\infty} a_n \quad \text{and} \quad \int_1^{\infty} f(x) \, dx \quad \text{both converge or both diverge}$$

Example: Deriving the convergence condition for p-series

$$\int_1^{\infty} \frac{1}{x^p} \, dx = \left[\frac{x^{1-p}}{1-p}\right]_1^{\infty}$$

If $p > 1$, the integral converges, so $\sum \frac{1}{n^p}$ also converges.

2.5 Alternating Series Test

An alternating series $\sum (-1)^n b_n$ ($b_n > 0$) converges if:

  1. $b_{n+1} \leq b_n$ (monotone decreasing)
  2. $\lim_{n \to \infty} b_n = 0$

Error bound for alternating series: $|S - S_N| \leq b_{N+1}$ (less than or equal to the absolute value of the next term)

import numpy as np
import sympy as sp

# ์ˆ˜๋ ด ํŒ์ •๋ฒ• ์ž๋™ ์ ์šฉ ๋„๊ตฌ
def test_convergence(a_n_func, name="๊ธ‰์ˆ˜"):
    """์ฃผ์–ด์ง„ ๊ธ‰์ˆ˜์— ๋Œ€ํ•ด ๋น„์œจ ํŒ์ •๋ฒ•๊ณผ ๊ทผ ํŒ์ •๋ฒ•์„ ์ ์šฉํ•ฉ๋‹ˆ๋‹ค."""
    n = sp.Symbol('n', positive=True, integer=True)
    a_n = a_n_func(n)

    print(f"=== {name}: sum a_n, a_n = {a_n} ===\n")

    # ๋น„์œจ ํŒ์ •๋ฒ•
    ratio = sp.simplify(a_n_func(n + 1) / a_n)
    rho_ratio = sp.limit(sp.Abs(ratio), n, sp.oo)
    print(f"๋น„์œจ ํŒ์ •๋ฒ•: |a_(n+1)/a_n| -> {rho_ratio}")
    if rho_ratio < 1:
        print("  => ์ ˆ๋Œ€ ์ˆ˜๋ ด\n")
    elif rho_ratio > 1:
        print("  => ๋ฐœ์‚ฐ\n")
    else:
        print("  => ํŒ์ • ๋ถˆ๋Šฅ\n")

    # ๊ทผ ํŒ์ •๋ฒ•
    root = sp.Abs(a_n) ** (sp.Rational(1, n))
    rho_root = sp.limit(root, n, sp.oo)
    print(f"๊ทผ ํŒ์ •๋ฒ•: |a_n|^(1/n) -> {rho_root}")
    if rho_root < 1:
        print("  => ์ ˆ๋Œ€ ์ˆ˜๋ ด\n")
    elif rho_root > 1:
        print("  => ๋ฐœ์‚ฐ\n")
    else:
        print("  => ํŒ์ • ๋ถˆ๋Šฅ\n")

# ํ…Œ์ŠคํŠธ ์˜ˆ์ œ๋“ค
n = sp.Symbol('n', positive=True, integer=True)

# 1) sum n! / n^n  (์ˆ˜๋ ด)
test_convergence(lambda n: sp.factorial(n) / n**n, "n!/n^n")

# 2) sum 1/n^2  (์ˆ˜๋ ด)
test_convergence(lambda n: 1 / n**2, "1/n^2")

# 3) sum n^2 / 2^n  (์ˆ˜๋ ด)
test_convergence(lambda n: n**2 / 2**n, "n^2/2^n")

import numpy as np
import matplotlib.pyplot as plt

# ๊ต๋Œ€๊ธ‰์ˆ˜ ์‹œ๊ฐํ™”: ln(2) = sum_{n=1}^{inf} (-1)^{n+1} / n
N = 30
n = np.arange(1, N + 1)
terms = (-1.0) ** (n + 1) / n
partial_sums = np.cumsum(terms)

plt.figure(figsize=(10, 5))
plt.plot(n, partial_sums, 'bo-', markersize=5, label='๋ถ€๋ถ„ํ•ฉ S_N')
plt.axhline(y=np.log(2), color='r', linestyle='--',
            label=f'ln(2) = {np.log(2):.6f}')

# ์˜ค์ฐจ ๋ฒ”์œ„ ํ‘œ์‹œ
for i in range(len(n)):
    error_bound = 1.0 / (n[i] + 1)
    plt.plot([n[i], n[i]],
             [partial_sums[i] - error_bound, partial_sums[i] + error_bound],
             'g-', alpha=0.3, linewidth=2)

plt.xlabel('N')
plt.ylabel('S_N')
plt.title('๊ต๋Œ€๊ธ‰์ˆ˜: ln(2) = 1 - 1/2 + 1/3 - 1/4 + ...')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('alternating_series.png', dpi=150, bbox_inches='tight')
plt.show()

3. Power Series

3.1 Radius of Convergence

A power series has the form:

$$f(x) = \sum_{n=0}^{\infty} c_n (x - a)^n = c_0 + c_1(x-a) + c_2(x-a)^2 + \cdots$$

where $a$ is the center and $c_n$ are coefficients.

Radius of convergence $R$:

$$R = \lim_{n \to \infty} \left|\frac{c_n}{c_{n+1}}\right| \quad \text{or} \quad R = \frac{1}{\lim_{n \to \infty} |c_n|^{1/n}}$$

  • $|x - a| < R$ implies absolute convergence
  • $|x - a| > R$ implies divergence
  • $|x - a| = R$ requires separate investigation

Interval of convergence:

      Diverges    Converges    Diverges
  -------|========|========|--------->  x
        a-R       a       a+R
         <---- R ---->

3.2 Term-by-Term Differentiation and Integration

Within the interval of convergence, power series can be differentiated and integrated term-by-term:

$$f'(x) = \sum_{n=1}^{\infty} n c_n (x-a)^{n-1}$$

$$\int f(x) \, dx = C + \sum_{n=0}^{\infty} \frac{c_n (x-a)^{n+1}}{n+1}$$

The radius of convergence remains the same after differentiation or integration (though convergence at endpoints may change).

import sympy as sp

x = sp.Symbol('x')

# ๋ฉฑ๊ธ‰์ˆ˜์˜ ์ˆ˜๋ ด ๋ฐ˜๊ฒฝ ๊ณ„์‚ฐ ์˜ˆ์ œ
print("=== ๋ฉฑ๊ธ‰์ˆ˜์˜ ์ˆ˜๋ ด ๋ฐ˜๊ฒฝ ===\n")

# 1) sum x^n / n!  (์ง€์ˆ˜ํ•จ์ˆ˜)
n = sp.Symbol('n', positive=True, integer=True)
c_n = 1 / sp.factorial(n)
c_n1 = 1 / sp.factorial(n + 1)
R1 = sp.limit(sp.Abs(c_n / c_n1), n, sp.oo)
print(f"sum x^n/n!: R = {R1}  (๋ชจ๋“  x์—์„œ ์ˆ˜๋ ด)")

# 2) sum n * x^n  (R = 1)
c_n = n
c_n1 = n + 1
R2 = sp.limit(sp.Abs(c_n / c_n1), n, sp.oo)
print(f"sum n*x^n: R = {R2}")

# 3) sum n! * x^n  (R = 0, x=0์—์„œ๋งŒ ์ˆ˜๋ ด)
c_n = sp.factorial(n)
c_n1 = sp.factorial(n + 1)
R3 = sp.limit(sp.Abs(c_n / c_n1), n, sp.oo)
print(f"sum n!*x^n: R = {R3}  (x=0์—์„œ๋งŒ ์ˆ˜๋ ด)")

# ํ•ญ๋ณ„ ๋ฏธ๋ถ„ ์‹œ๊ฐํ™”
print("\n=== ํ•ญ๋ณ„ ๋ฏธ๋ถ„: d/dx [1/(1-x)] = 1/(1-x)^2 ===")
# 1/(1-x) = sum x^n  =>  ๋ฏธ๋ถ„ํ•˜๋ฉด  1/(1-x)^2 = sum n*x^{n-1}
x_vals = np.linspace(-0.9, 0.9, 200)

import numpy as np
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# ์›๋ž˜ ํ•จ์ˆ˜์™€ ๊ธ‰์ˆ˜
for N in [3, 5, 10, 20]:
    series_sum = sum(x_vals**k for k in range(N))
    axes[0].plot(x_vals, series_sum, alpha=0.7, label=f'N={N}')

axes[0].plot(x_vals, 1/(1-x_vals), 'k--', linewidth=2, label='1/(1-x)')
axes[0].set_ylim(-5, 15)
axes[0].set_title('1/(1-x) = sum x^n')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# ๋ฏธ๋ถ„ํ•œ ํ•จ์ˆ˜์™€ ๊ธ‰์ˆ˜
for N in [3, 5, 10, 20]:
    deriv_sum = sum(k * x_vals**(k-1) for k in range(1, N))
    axes[1].plot(x_vals, deriv_sum, alpha=0.7, label=f'N={N}')

axes[1].plot(x_vals, 1/(1-x_vals)**2, 'k--', linewidth=2, label='1/(1-x)^2')
axes[1].set_ylim(-5, 30)
axes[1].set_title('1/(1-x)^2 = sum n*x^{n-1}')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('power_series_differentiation.png', dpi=150, bbox_inches='tight')
plt.show()

4. Taylor Series and Maclaurin Series

4.1 Taylor Expansion

If a function $f(x)$ is infinitely differentiable near $x = a$, it can be expanded as a Taylor series:

$$f(x) = \sum_{n=0}^{\infty} \frac{f^{(n)}(a)}{n!} (x - a)^n$$

$$= f(a) + f'(a)(x-a) + \frac{f''(a)}{2!} (x-a)^2 + \frac{f'''(a)}{3!} (x-a)^3 + \cdots$$

When $a = 0$, this is called a Maclaurin series.

Taylor remainder $R_N$: the error after expanding to $N$ terms

$$R_N(x) = \frac{f^{(N+1)}(c)}{(N+1)!} (x-a)^{N+1} \quad \text{(Lagrange remainder, } a < c < x\text{)}$$

4.2 Series Representations of Common Functions

Most frequently used series expansions in physics:

$$e^x = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \cdots \quad \text{(all } x\text{)}$$

$$\sin(x) = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \cdots \quad \text{(all } x\text{)}$$

$$\cos(x) = 1 - \frac{x^2}{2!} + \frac{x^4}{4!} - \frac{x^6}{6!} + \cdots \quad \text{(all } x\text{)}$$

$$\ln(1+x) = x - \frac{x^2}{2} + \frac{x^3}{3} - \frac{x^4}{4} + \cdots \quad (-1 < x \leq 1)$$

$$(1+x)^p = 1 + px + \frac{p(p-1)}{2!} x^2 + \cdots \quad (|x| < 1, \text{ binomial series})$$

$$\frac{1}{1-x} = 1 + x + x^2 + x^3 + \cdots \quad (|x| < 1)$$

$$\arctan(x) = x - \frac{x^3}{3} + \frac{x^5}{5} - \frac{x^7}{7} + \cdots \quad (|x| \leq 1)$$

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt

# SymPy๋ฅผ ์ด์šฉํ•œ ํ…Œ์ผ๋Ÿฌ ์ „๊ฐœ
x = sp.Symbol('x')

functions = {
    'e^x': sp.exp(x),
    'sin(x)': sp.sin(x),
    'cos(x)': sp.cos(x),
    'ln(1+x)': sp.ln(1 + x),
}

print("=== ์ฃผ์š” ํ•จ์ˆ˜์˜ ํ…Œ์ผ๋Ÿฌ ๊ธ‰์ˆ˜ (a=0, 6์ฐจ๊นŒ์ง€) ===\n")
for name, func in functions.items():
    taylor = sp.series(func, x, 0, n=7)
    print(f"{name} = {taylor}\n")

# ํ…Œ์ผ๋Ÿฌ ๊ธ‰์ˆ˜์˜ ์ˆ˜๋ ด ์‹œ๊ฐํ™”
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

x_vals = np.linspace(-2 * np.pi, 2 * np.pi, 500)

# sin(x)์˜ ํ…Œ์ผ๋Ÿฌ ๊ทผ์‚ฌ
ax = axes[0, 0]
ax.plot(x_vals, np.sin(x_vals), 'k-', linewidth=2, label='sin(x)')
for N in [1, 3, 5, 7, 9]:
    approx = sum((-1)**k * x_vals**(2*k+1) / np.math.factorial(2*k+1)
                 for k in range(N))
    ax.plot(x_vals, approx, alpha=0.7, label=f'{2*N-1}์ฐจ')
ax.set_ylim(-2, 2)
ax.set_title('sin(x)์˜ ํ…Œ์ผ๋Ÿฌ ๊ทผ์‚ฌ')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# cos(x)์˜ ํ…Œ์ผ๋Ÿฌ ๊ทผ์‚ฌ
ax = axes[0, 1]
ax.plot(x_vals, np.cos(x_vals), 'k-', linewidth=2, label='cos(x)')
for N in [1, 2, 4, 6, 8]:
    approx = sum((-1)**k * x_vals**(2*k) / np.math.factorial(2*k)
                 for k in range(N))
    ax.plot(x_vals, approx, alpha=0.7, label=f'{2*N-2}์ฐจ')
ax.set_ylim(-2, 2)
ax.set_title('cos(x)์˜ ํ…Œ์ผ๋Ÿฌ ๊ทผ์‚ฌ')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# e^x์˜ ํ…Œ์ผ๋Ÿฌ ๊ทผ์‚ฌ
x_exp = np.linspace(-3, 3, 300)
ax = axes[1, 0]
ax.plot(x_exp, np.exp(x_exp), 'k-', linewidth=2, label='e^x')
for N in [1, 2, 4, 6, 10]:
    approx = sum(x_exp**k / np.math.factorial(k) for k in range(N + 1))
    ax.plot(x_exp, approx, alpha=0.7, label=f'{N}์ฐจ')
ax.set_ylim(-5, 20)
ax.set_title('e^x์˜ ํ…Œ์ผ๋Ÿฌ ๊ทผ์‚ฌ')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# ln(1+x)์˜ ํ…Œ์ผ๋Ÿฌ ๊ทผ์‚ฌ
x_ln = np.linspace(-0.95, 2, 300)
ax = axes[1, 1]
ax.plot(x_ln, np.log(1 + x_ln), 'k-', linewidth=2, label='ln(1+x)')
for N in [1, 3, 5, 10, 20]:
    x_series = np.linspace(-0.95, 0.95, 300)  # ์ˆ˜๋ ด ๊ตฌ๊ฐ„ ๋‚ด
    approx = sum((-1)**(k+1) * x_series**k / k for k in range(1, N + 1))
    ax.plot(x_series, approx, alpha=0.7, label=f'{N}์ฐจ')
ax.axvline(x=1, color='gray', linestyle=':', alpha=0.5, label='x=1 (๊ฒฝ๊ณ„)')
ax.axvline(x=-1, color='gray', linestyle=':', alpha=0.5, label='x=-1 (๊ฒฝ๊ณ„)')
ax.set_ylim(-4, 2)
ax.set_title('ln(1+x)์˜ ํ…Œ์ผ๋Ÿฌ ๊ทผ์‚ฌ')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('taylor_series_visualization.png', dpi=150, bbox_inches='tight')
plt.show()

4.3 Approximation Using Series

The most important application of series expansion in physics is small-quantity approximation.

When $\varepsilon \ll 1$:

$$(1 + \varepsilon)^p \approx 1 + p\varepsilon + \frac{p(p-1)}{2} \varepsilon^2 + \cdots$$

$$e^\varepsilon \approx 1 + \varepsilon + \frac{\varepsilon^2}{2} + \cdots$$

$$\sin(\varepsilon) \approx \varepsilon - \frac{\varepsilon^3}{6} + \cdots$$

$$\cos(\varepsilon) \approx 1 - \frac{\varepsilon^2}{2} + \cdots$$

$$\tan(\varepsilon) \approx \varepsilon + \frac{\varepsilon^3}{3} + \cdots$$

$$\frac{1}{1+\varepsilon} \approx 1 - \varepsilon + \varepsilon^2 - \cdots$$

import numpy as np
import matplotlib.pyplot as plt

# ์†Œ๋Ÿ‰ ๊ทผ์‚ฌ์˜ ์ •ํ™•๋„ ๋น„๊ต
epsilon = np.linspace(0, 1, 100)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# sqrt(1 + epsilon) = (1+epsilon)^{1/2} ์˜ ๊ทผ์‚ฌ
exact = np.sqrt(1 + epsilon)
order1 = 1 + epsilon / 2                       # 1์ฐจ ๊ทผ์‚ฌ
order2 = 1 + epsilon / 2 - epsilon**2 / 8      # 2์ฐจ ๊ทผ์‚ฌ
order3 = 1 + epsilon/2 - epsilon**2/8 + epsilon**3/16  # 3์ฐจ

axes[0].plot(epsilon, exact, 'k-', linewidth=2, label='์ •ํ™•๊ฐ’')
axes[0].plot(epsilon, order1, 'b--', label='1์ฐจ: 1 + e/2')
axes[0].plot(epsilon, order2, 'r--', label='2์ฐจ: 1 + e/2 - e^2/8')
axes[0].plot(epsilon, order3, 'g--', label='3์ฐจ')
axes[0].set_xlabel('epsilon')
axes[0].set_ylabel('sqrt(1 + epsilon)')
axes[0].set_title('์ดํ•ญ๊ธ‰์ˆ˜ ๊ทผ์‚ฌ: (1+e)^{1/2}')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# ์ƒ๋Œ€ ์˜ค์ฐจ
rel_err1 = np.abs(order1 - exact) / exact * 100
rel_err2 = np.abs(order2 - exact) / exact * 100
rel_err3 = np.abs(order3 - exact) / exact * 100

axes[1].semilogy(epsilon, rel_err1, 'b-', label='1์ฐจ ๊ทผ์‚ฌ')
axes[1].semilogy(epsilon, rel_err2, 'r-', label='2์ฐจ ๊ทผ์‚ฌ')
axes[1].semilogy(epsilon, rel_err3, 'g-', label='3์ฐจ ๊ทผ์‚ฌ')
axes[1].set_xlabel('epsilon')
axes[1].set_ylabel('์ƒ๋Œ€ ์˜ค์ฐจ (%)')
axes[1].set_title('๊ทผ์‚ฌ ์ฐจ์ˆ˜์— ๋”ฐ๋ฅธ ์ƒ๋Œ€ ์˜ค์ฐจ')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('approximation_error.png', dpi=150, bbox_inches='tight')
plt.show()

5. Asymptotic Series

5.1 Divergent but Useful Series

An asymptotic series diverges when summed to infinity, but finite partial sums provide excellent approximations to a function.

Asymptotic expansion of a function $f(x)$:

$$f(x) \sim \sum_{n=0}^{\infty} \frac{a_n}{x^n} \quad (x \to \infty)$$

This means that for the partial sum of $N$ terms, the following holds:

$$\lim_{x \to \infty} x^N \left[f(x) - \sum_{n=0}^{N} \frac{a_n}{x^n}\right] = 0 \quad \text{(for each } N\text{)}$$

Key Point: The entire series diverges, but truncating at the optimal number of terms yields an excellent approximation. Generally, the optimal truncation point is where the terms stop decreasing and start increasing again.

5.2 Stirling's Approximation

Stirling's approximation is an asymptotic approximation for factorials:

$$n! \sim \sqrt{2\pi n} \left(\frac{n}{e}\right)^n \quad (n \to \infty)$$

$$\ln(n!) \sim n\ln(n) - n + \frac{1}{2}\ln(2\pi n)$$

This approximation is crucial in statistical mechanics for calculating Boltzmann entropy, combinatorial problems, etc.

More precise Stirling series:

$$n! \sim \sqrt{2\pi n} \left(\frac{n}{e}\right)^n \left(1 + \frac{1}{12n} + \frac{1}{288n^2} - \frac{139}{51840n^3} + \cdots\right)$$

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma, factorial

# ์Šคํ„ธ๋ง ๊ทผ์‚ฌ์˜ ์ •ํ™•๋„
n_vals = np.arange(1, 51)

# ์ •ํ™•ํ•œ ln(n!)
exact_lnfact = np.array([np.sum(np.log(np.arange(1, n+1))) for n in n_vals])

# ์Šคํ„ธ๋ง ๊ทผ์‚ฌ (์—ฌ๋Ÿฌ ์ฐจ์ˆ˜)
stirling_0 = n_vals * np.log(n_vals) - n_vals  # ๊ฐ€์žฅ ๊ฐ„๋‹จํ•œ ํ˜•ํƒœ
stirling_1 = stirling_0 + 0.5 * np.log(2 * np.pi * n_vals)  # 1์ฐจ ๋ณด์ •
stirling_2 = stirling_1 + 1 / (12 * n_vals)  # 2์ฐจ ๋ณด์ •

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# ์ ˆ๋Œ€๊ฐ’ ๋น„๊ต
axes[0].plot(n_vals, exact_lnfact, 'k-', linewidth=2, label='ln(n!) ์ •ํ™•๊ฐ’')
axes[0].plot(n_vals, stirling_0, 'b--', label='n*ln(n) - n')
axes[0].plot(n_vals, stirling_1, 'r--', label='+ (1/2)ln(2*pi*n)')
axes[0].set_xlabel('n')
axes[0].set_ylabel('ln(n!)')
axes[0].set_title('์Šคํ„ธ๋ง ๊ทผ์‚ฌ')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# ์ƒ๋Œ€ ์˜ค์ฐจ
rel_err0 = np.abs(stirling_0 - exact_lnfact) / exact_lnfact * 100
rel_err1 = np.abs(stirling_1 - exact_lnfact) / exact_lnfact * 100
rel_err2 = np.abs(stirling_2 - exact_lnfact) / exact_lnfact * 100

axes[1].semilogy(n_vals, rel_err0, 'b-', label='0์ฐจ: n*ln(n)-n')
axes[1].semilogy(n_vals, rel_err1, 'r-', label='1์ฐจ: + ln(sqrt(2*pi*n))')
axes[1].semilogy(n_vals, rel_err2, 'g-', label='2์ฐจ: + 1/(12n)')
axes[1].set_xlabel('n')
axes[1].set_ylabel('์ƒ๋Œ€ ์˜ค์ฐจ (%)')
axes[1].set_title('์Šคํ„ธ๋ง ๊ทผ์‚ฌ์˜ ์ƒ๋Œ€ ์˜ค์ฐจ')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('stirling_approximation.png', dpi=150, bbox_inches='tight')
plt.show()

# ํ†ต๊ณ„์—ญํ•™ ์‘์šฉ: ๋ณผ์ธ ๋งŒ ์—”ํŠธ๋กœํ”ผ
print("=== ํ†ต๊ณ„์—ญํ•™ ์‘์šฉ: ์ด์ƒ๊ธฐ์ฒด ์—”ํŠธ๋กœํ”ผ ===\n")
print("N๊ฐœ ์ž…์ž๊ฐ€ W๊ฐœ ๋ฏธ์‹œ์ƒํƒœ๋ฅผ ๊ฐ€์งˆ ๋•Œ:")
print("S = k_B * ln(W)")
print("์Šคํ„ธ๋ง ๊ทผ์‚ฌ๋ฅผ ์‚ฌ์šฉํ•˜๋ฉด:")
print("ln(N!) ~ N*ln(N) - N")
print()
print("์˜ˆ: N = 10^23 (์•„๋ณด๊ฐ€๋“œ๋กœ ์ˆ˜ ๊ทœ๋ชจ)")
N = 1e23
print(f"ln(N!) ~ N*ln(N) - N = {N * np.log(N) - N:.4e}")
print("์ง์ ‘ ๊ณ„์‚ฐ์€ ๋ถˆ๊ฐ€๋Šฅํ•˜์ง€๋งŒ, ์Šคํ„ธ๋ง ๊ทผ์‚ฌ๋กœ ์‰ฝ๊ฒŒ ๊ณ„์‚ฐ ๊ฐ€๋Šฅ!")

6. Applications in Physics

6.1 Period of a Pendulum (Series Expansion)

The exact period of a simple pendulum is given by a complete elliptic integral:

$$T = 4\sqrt{\frac{L}{g}} K\left(\sin\frac{\theta_0}{2}\right)$$

where $K(k)$ is the complete elliptic integral of the first kind. Expanding as a series:

$$T = 2\pi\sqrt{\frac{L}{g}} \left[1 + \frac{1}{4}\sin^2\frac{\theta_0}{2} + \frac{9}{64}\sin^4\frac{\theta_0}{2} + \cdots\right]$$

$$\approx T_0 \left[1 + \frac{\theta_0^2}{16} + \frac{11\theta_0^4}{3072} + \cdots\right]$$

where $T_0 = 2\pi\sqrt{L/g}$ is the period in the small-angle approximation.

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import ellipk

# ๋‹จ์ง„์ž์˜ ์ฃผ๊ธฐ: ์ •ํ™•๊ฐ’ vs ๊ธ‰์ˆ˜ ๊ทผ์‚ฌ
g = 9.81  # m/s^2
L = 1.0   # m
T0 = 2 * np.pi * np.sqrt(L / g)  # ์†Œ๊ฐ ๊ทผ์‚ฌ ์ฃผ๊ธฐ

theta0 = np.linspace(0.01, np.pi * 0.95, 200)  # ์ดˆ๊ธฐ ๊ฐ๋„ (rad)
k = np.sin(theta0 / 2)

# ์ •ํ™•ํ•œ ์ฃผ๊ธฐ (์™„์ „ ํƒ€์› ์ ๋ถ„)
T_exact = 4 * np.sqrt(L / g) * ellipk(k**2)

# ๊ธ‰์ˆ˜ ๊ทผ์‚ฌ (์†Œ๊ฐ ๊ทผ์‚ฌ๋ถ€ํ„ฐ ๊ณ ์ฐจํ•ญ๊นŒ์ง€)
T_approx0 = T0 * np.ones_like(theta0)                          # 0์ฐจ (์†Œ๊ฐ)
T_approx2 = T0 * (1 + (1/4) * k**2)                           # 2์ฐจ
T_approx4 = T0 * (1 + (1/4)*k**2 + (9/64)*k**4)               # 4์ฐจ
T_approx6 = T0 * (1 + (1/4)*k**2 + (9/64)*k**4 + (25/256)*k**6)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# ์ฃผ๊ธฐ ๋น„๊ต
theta_deg = np.degrees(theta0)
axes[0].plot(theta_deg, T_exact / T0, 'k-', linewidth=2, label='์ •ํ™•๊ฐ’')
axes[0].plot(theta_deg, T_approx0 / T0, 'b--', label='0์ฐจ (์†Œ๊ฐ๊ทผ์‚ฌ)')
axes[0].plot(theta_deg, T_approx2 / T0, 'r--', label='2์ฐจ ๋ณด์ •')
axes[0].plot(theta_deg, T_approx4 / T0, 'g--', label='4์ฐจ ๋ณด์ •')
axes[0].plot(theta_deg, T_approx6 / T0, 'm--', label='6์ฐจ ๋ณด์ •')
axes[0].set_xlabel('์ดˆ๊ธฐ ๊ฐ๋„ (๋„)')
axes[0].set_ylabel('T / T_0')
axes[0].set_title('๋‹จ์ง„์ž ์ฃผ๊ธฐ: ๊ธ‰์ˆ˜ ์ „๊ฐœ ๊ทผ์‚ฌ')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# ์ƒ๋Œ€ ์˜ค์ฐจ
axes[1].semilogy(theta_deg, np.abs(T_approx0 - T_exact)/T_exact * 100,
                 'b-', label='0์ฐจ (์†Œ๊ฐ)')
axes[1].semilogy(theta_deg, np.abs(T_approx2 - T_exact)/T_exact * 100,
                 'r-', label='2์ฐจ')
axes[1].semilogy(theta_deg, np.abs(T_approx4 - T_exact)/T_exact * 100,
                 'g-', label='4์ฐจ')
axes[1].semilogy(theta_deg, np.abs(T_approx6 - T_exact)/T_exact * 100,
                 'm-', label='6์ฐจ')
axes[1].set_xlabel('์ดˆ๊ธฐ ๊ฐ๋„ (๋„)')
axes[1].set_ylabel('์ƒ๋Œ€ ์˜ค์ฐจ (%)')
axes[1].set_title('๊ธ‰์ˆ˜ ๊ทผ์‚ฌ์˜ ์ƒ๋Œ€ ์˜ค์ฐจ')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('pendulum_period.png', dpi=150, bbox_inches='tight')
plt.show()

# ์ˆ˜์น˜ ๊ฒฐ๊ณผ
print("=== ๋‹จ์ง„์ž ์ฃผ๊ธฐ ๋น„๊ต (L=1m) ===\n")
print(f"{'๊ฐ๋„':>6s}  {'์ •ํ™•๊ฐ’(s)':>10s}  {'์†Œ๊ฐ๊ทผ์‚ฌ(s)':>12s}  {'์˜ค์ฐจ(%)':>8s}")
print("-" * 45)
for angle_deg in [5, 10, 30, 45, 60, 90, 120, 150]:
    angle_rad = np.radians(angle_deg)
    k_val = np.sin(angle_rad / 2)
    T_ex = 4 * np.sqrt(L / g) * ellipk(k_val**2)
    err = (T0 - T_ex) / T_ex * 100
    print(f"{angle_deg:>5d}   {T_ex:>10.6f}  {T0:>12.6f}  {err:>8.3f}")

6.2 Non-relativistic Approximation of Relativistic Energy

Einstein's relativistic energy-momentum relation:

$$E = \gamma mc^2 = \frac{mc^2}{\sqrt{1 - v^2/c^2}}$$

Expanding the kinetic energy $K = E - mc^2$ as a series (when $\beta = v/c \ll 1$):

$$\gamma = \frac{1}{\sqrt{1 - \beta^2}} = 1 + \frac{1}{2}\beta^2 + \frac{3}{8}\beta^4 + \frac{5}{16}\beta^6 + \cdots$$

$$K = mc^2 (\gamma - 1)$$

$$= \underbrace{\frac{1}{2}mv^2}_{\text{Newtonian}} + \underbrace{\frac{3}{8}m\frac{v^4}{c^2} + \frac{5}{16}m\frac{v^6}{c^4} + \cdots}_{\text{Relativistic correction}}$$

import numpy as np
import matplotlib.pyplot as plt

# ์ƒ๋Œ€๋ก ์  ์šด๋™ ์—๋„ˆ์ง€์˜ ๊ธ‰์ˆ˜ ๊ทผ์‚ฌ
beta = np.linspace(0.001, 0.99, 500)  # v/c

# ์ •ํ™•ํ•œ ์ƒ๋Œ€๋ก ์  ์šด๋™ ์—๋„ˆ์ง€: K/(mc^2) = gamma - 1
gamma = 1 / np.sqrt(1 - beta**2)
K_exact = gamma - 1  # mc^2 ๋‹จ์œ„

# ๊ธ‰์ˆ˜ ๊ทผ์‚ฌ (์ดํ•ญ๊ธ‰์ˆ˜ ์ „๊ฐœ)
K_order2 = 0.5 * beta**2                                         # ๋‰ดํ„ด ์—ญํ•™
K_order4 = 0.5 * beta**2 + (3/8) * beta**4                      # 1์ฐจ ๋ณด์ •
K_order6 = 0.5 * beta**2 + (3/8) * beta**4 + (5/16) * beta**6  # 2์ฐจ ๋ณด์ •

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# ์šด๋™ ์—๋„ˆ์ง€ ๋น„๊ต
axes[0].plot(beta, K_exact, 'k-', linewidth=2, label='์ƒ๋Œ€๋ก ์  (์ •ํ™•)')
axes[0].plot(beta, K_order2, 'b--', label='๋‰ดํ„ด: (1/2)mv^2')
axes[0].plot(beta, K_order4, 'r--', label='1์ฐจ ๋ณด์ •')
axes[0].plot(beta, K_order6, 'g--', label='2์ฐจ ๋ณด์ •')
axes[0].set_xlabel('beta = v/c')
axes[0].set_ylabel('K / (mc^2)')
axes[0].set_title('์šด๋™ ์—๋„ˆ์ง€: ์ƒ๋Œ€๋ก  vs ๋‰ดํ„ด')
axes[0].set_ylim(0, 5)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# ๋‰ดํ„ด ์—ญํ•™์˜ ์ƒ๋Œ€ ์˜ค์ฐจ
beta_low = np.linspace(0.001, 0.5, 300)
gamma_low = 1 / np.sqrt(1 - beta_low**2)
K_exact_low = gamma_low - 1
K_newton = 0.5 * beta_low**2
K_corr1 = 0.5 * beta_low**2 + (3/8) * beta_low**4

rel_err_newton = np.abs(K_newton - K_exact_low) / K_exact_low * 100
rel_err_corr1 = np.abs(K_corr1 - K_exact_low) / K_exact_low * 100

axes[1].semilogy(beta_low, rel_err_newton, 'b-', label='๋‰ดํ„ด ์—ญํ•™ ์˜ค์ฐจ')
axes[1].semilogy(beta_low, rel_err_corr1, 'r-', label='1์ฐจ ๋ณด์ • ํ›„ ์˜ค์ฐจ')
axes[1].axhline(y=1, color='gray', linestyle=':', alpha=0.5, label='1% ์˜ค์ฐจ์„ ')
axes[1].set_xlabel('beta = v/c')
axes[1].set_ylabel('์ƒ๋Œ€ ์˜ค์ฐจ (%)')
axes[1].set_title('๋‰ดํ„ด ์—ญํ•™์˜ ์œ ํšจ ๋ฒ”์œ„')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('relativistic_energy.png', dpi=150, bbox_inches='tight')
plt.show()

# ๊ฒฐ๊ณผ ์ถœ๋ ฅ
print("=== ๋‰ดํ„ด ์—ญํ•™์ด ์œ ํšจํ•œ ์†๋„ ๋ฒ”์œ„ ===\n")
print(f"{'v/c':>6s}  {'์ƒ๋Œ€๋ก ์  K':>12s}  {'๋‰ดํ„ด K':>10s}  {'์˜ค์ฐจ(%)':>8s}")
print("-" * 42)
for b in [0.01, 0.05, 0.1, 0.2, 0.3, 0.5, 0.7, 0.9]:
    g = 1 / np.sqrt(1 - b**2)
    K_rel = g - 1
    K_new = 0.5 * b**2
    err = abs(K_new - K_rel) / K_rel * 100
    print(f"{b:>6.2f}  {K_rel:>12.6f}  {K_new:>10.6f}  {err:>8.3f}")

6.3 Electric Multipole Expansion

When a point charge $q$ is located at distance $d$ from the origin, the potential at a point far from the origin ($r \gg d$) can be expanded as a series.

    Observation point P(r, theta)
    *
    |\
    | \  r
    |  \
    |   \
    |theta\
    |------*--- Charge q (distance d, on z-axis)
    |
    Origin O

Distance to observation point:

$$\frac{1}{|\mathbf{r} - \mathbf{d}|} = \frac{1}{r} \sum_{l=0}^{\infty} \left(\frac{d}{r}\right)^l P_l(\cos\theta)$$

where $P_l$ are Legendre polynomials.

For an electric dipole ($+q$ and $-q$ separated by distance $d$):

$$V(r, \theta) \approx \frac{1}{4\pi\varepsilon_0} \frac{p\cos\theta}{r^2} \quad (r \gg d)$$

where $p = qd$ is the dipole moment.

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import legendre

# ๋‹ค์ค‘๊ทน ์ „๊ฐœ (Multipole Expansion) ์‹œ๊ฐํ™”
# ์ „ํ•˜ q๊ฐ€ z์ถ• ์œ„ d = 1์— ์žˆ์„ ๋•Œ, 1/|r - d*z_hat|์˜ ๊ธ‰์ˆ˜ ์ „๊ฐœ

def exact_potential_ratio(r, theta, d=1.0):
    """์ •ํ™•ํ•œ ์ „์œ„ ๋น„์œจ: r / |r - d*z_hat| (1/r ์ œ์™ธํ•œ ๋ถ€๋ถ„)"""
    # |r - d*z_hat|^2 = r^2 - 2*r*d*cos(theta) + d^2
    dist = np.sqrt(r**2 - 2*r*d*np.cos(theta) + d**2)
    return r / dist

def multipole_approx(r, theta, d=1.0, L_max=5):
    """๋‹ค์ค‘๊ทน ์ „๊ฐœ: sum_{l=0}^{L_max} (d/r)^l * P_l(cos(theta))"""
    result = np.zeros_like(r)
    cos_theta = np.cos(theta)
    for l in range(L_max + 1):
        Pl = legendre(l)
        result += (d / r)**l * Pl(cos_theta)
    return result

# r/d ๋น„์œจ์— ๋”ฐ๋ฅธ ๋‹ค์ค‘๊ทน ์ „๊ฐœ์˜ ์ •ํ™•๋„
r_values = np.linspace(1.5, 10, 200)
theta = np.pi / 4  # 45๋„

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

exact = exact_potential_ratio(r_values, theta)

for L_max in [0, 1, 2, 5, 10]:
    approx = multipole_approx(r_values, theta, L_max=L_max)
    axes[0].plot(r_values, approx, label=f'L_max = {L_max}')

axes[0].plot(r_values, exact, 'k--', linewidth=2, label='์ •ํ™•๊ฐ’')
axes[0].set_xlabel('r/d')
axes[0].set_ylabel('r * V(r) / (kq)')
axes[0].set_title(f'๋‹ค์ค‘๊ทน ์ „๊ฐœ (theta = 45๋„)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# ์Œ๊ทน์ž ์ „์œ„ ํŒจํ„ด (2D)
r_grid = np.linspace(0.5, 5, 200)
theta_grid = np.linspace(0, 2*np.pi, 200)
R, Theta = np.meshgrid(r_grid, theta_grid)
X = R * np.sin(Theta)
Z = R * np.cos(Theta)

# ์Œ๊ทน์ž ์ „์œ„: V ~ cos(theta)/r^2 (๋‹จ์œ„ ์ƒ๋žต)
V_dipole = np.cos(Theta) / R**2

axes[1].contourf(X, Z, V_dipole, levels=np.linspace(-2, 2, 41),
                 cmap='RdBu_r', extend='both')
axes[1].set_xlabel('x')
axes[1].set_ylabel('z')
axes[1].set_title('์ „๊ธฐ ์Œ๊ทน์ž ์ „์œ„ ํŒจํ„ด')
axes[1].set_aspect('equal')
axes[1].set_xlim(-3, 3)
axes[1].set_ylim(-3, 3)
# ์ „ํ•˜ ์œ„์น˜ ํ‘œ์‹œ
axes[1].plot(0, 0.1, 'r+', markersize=15, markeredgewidth=2)
axes[1].plot(0, -0.1, 'b_', markersize=15, markeredgewidth=2)

plt.tight_layout()
plt.savefig('multipole_expansion.png', dpi=150, bbox_inches='tight')
plt.show()

# ๋ฅด์žฅ๋“œ๋ฅด ๋‹คํ•ญ์‹์˜ ์ฒ˜์Œ ๋ช‡ ๊ฐœ
print("=== ๋ฅด์žฅ๋“œ๋ฅด ๋‹คํ•ญ์‹ P_l(x) ===\n")
import sympy as sp
x_sym = sp.Symbol('x')
for l in range(6):
    Pl = sp.legendre(l, x_sym)
    print(f"P_{l}(x) = {sp.expand(Pl)}")

Practice Problems

Problem 1. Convergence Test (Basic)

Determine the convergence/divergence of the following series. State the test used.

(a) $\sum_{n=1}^{\infty} \frac{n^2}{3^n}$

(b) $\sum_{n=2}^{\infty} \frac{1}{n \ln^2(n)}$

(c) $\sum_{n=1}^{\infty} \frac{(-1)^n n}{n^2 + 1}$

Hint: (a) ratio test, (b) integral test, (c) alternating series test

Problem 2. Radius of Convergence (Basic)

Find the radius of convergence $R$ for the following power series.

(a) $\sum_{n=0}^{\infty} \frac{n! x^n}{(2n)!}$

(b) $\sum_{n=1}^{\infty} \frac{(-1)^{n+1} x^n}{n \cdot 3^n}$

Problem 3. Taylor Series Approximation (Intermediate)

Expand $f(x) = \sqrt{1 + x}$ as a Taylor series about $x = 0$ up to 3rd order, and:

(a) Find the approximate value of $\sqrt{1.1}$ (substitute $x = 0.1$)

(b) Calculate the relative error by comparing with the exact value

(c) Use Python to plot a comparison of errors for 1st, 2nd, and 3rd order approximations

Problem 4. Physics Application (Intermediate)

For a pendulum with maximum angle $\theta_0 = 30$ degrees:

(a) Calculate the exact period (elliptic integral)

(b) Compare with the small-angle approximation period and find the error

(c) Find the error when including the 2nd order correction term

Problem 5. Stirling's Approximation (Intermediate)

(a) Use Stirling's approximation to find $\log_{10}(100!)$

(b) Compare with Python's exact calculation

(c) For an ideal gas with $N$ molecules having $W$ microstates $W = N! / (n_1! n_2! \cdots n_k!)$, explain why Stirling's approximation is essential

Problem 6. Comprehensive Application (Advanced)

When an electron has velocity $v = 0.1c$:

(a) Find the exact relativistic kinetic energy in $mc^2$ units

(b) Compare with Newtonian kinetic energy and find the relative error

(c) Calculate how much the error reduces when including the 1st order relativistic correction term $\frac{3}{8}m\frac{v^4}{c^2}$

(d) Find the maximum $v/c$ for which Newtonian mechanics has an error less than 1%

# ์—ฐ์Šต ๋ฌธ์ œ ํ’€์ด ๋„์šฐ๋ฏธ
import numpy as np
import sympy as sp

# ๋ฌธ์ œ 1(a) ๊ฒ€์ฆ
n = sp.Symbol('n', positive=True, integer=True)
a_n = n**2 / 3**n
ratio = sp.simplify(sp.Abs((n+1)**2 / 3**(n+1)) / (n**2 / 3**n))
rho = sp.limit(ratio, n, sp.oo)
print(f"๋ฌธ์ œ 1(a): ๋น„์œจ = {rho} < 1 => ์ˆ˜๋ ด")

# ๋ฌธ์ œ 4 ๊ฒ€์ฆ
from scipy.special import ellipk
L, g = 1.0, 9.81
theta0 = np.radians(30)
T0 = 2 * np.pi * np.sqrt(L / g)
k = np.sin(theta0 / 2)
T_exact = 4 * np.sqrt(L / g) * ellipk(k**2)
print(f"\n๋ฌธ์ œ 4: T_exact = {T_exact:.6f}s, T_0 = {T0:.6f}s")
print(f"์†Œ๊ฐ ๊ทผ์‚ฌ ์˜ค์ฐจ = {abs(T0 - T_exact)/T_exact * 100:.4f}%")

T_corrected = T0 * (1 + (1/4) * k**2)
print(f"2์ฐจ ๋ณด์ • ํ›„ ์˜ค์ฐจ = {abs(T_corrected - T_exact)/T_exact * 100:.6f}%")

# ๋ฌธ์ œ 5 ๊ฒ€์ฆ
import math
stirling_log10_100 = (100 * np.log(100) - 100 + 0.5 * np.log(200 * np.pi)) / np.log(10)
exact_log10_100 = np.log10(float(math.factorial(100)))
print(f"\n๋ฌธ์ œ 5: log10(100!) ์Šคํ„ธ๋ง = {stirling_log10_100:.4f}")
print(f"         log10(100!) ์ •ํ™•๊ฐ’ = {exact_log10_100:.4f}")
print(f"         ์˜ค์ฐจ = {abs(stirling_log10_100 - exact_log10_100):.6f}")

References

Textbooks

  1. Boas, M. L. (2005). Mathematical Methods in the Physical Sciences, 3rd ed., Chapter 1. Wiley.
  2. Arfken, G. B., Weber, H. J., & Harris, F. E. (2012). Mathematical Methods for Physicists, 7th ed., Chapter 5. Academic Press.
  3. Riley, K. F. et al. (2006). Mathematical Methods for Physics and Engineering, Chapter 4. Cambridge University Press.

Online Resources

  1. MIT OCW 18.01SC: Single Variable Calculus - Sequences and Series
  2. Paul's Online Math Notes: Series & Sequences (https://tutorial.math.lamar.edu/)
  3. 3Blue1Brown: Taylor series visualization

Python Tools

  • sympy.series(): symbolic Taylor series calculation
  • scipy.special.ellipk(): complete elliptic integral
  • scipy.special.factorial(): factorial calculation
  • numpy.cumsum(): partial sum calculation

Next Lesson

02. Complex Numbers covers algebraic operations with complex numbers, polar and exponential representations, De Moivre's theorem, and Euler's formula which is essential in physics. We will use Taylor expansion learned in series to derive $e^{ix} = \cos(x) + i\sin(x)$.

Previous Mathematical Methods in the Physical Sciences - Overview
๐Ÿ“‹ Topic list
Next 02. Complex Numbers
← → to navigate between lessons