01_infinite_series.py

Download
python 298 lines 9.1 KB
  1"""
  2Infinite Series - Convergence Tests and Taylor/Maclaurin Expansions
  3
  4This script demonstrates:
  5- Series convergence tests (ratio, root, comparison)
  6- Partial sums and convergence visualization
  7- Taylor and Maclaurin series expansions
  8- Euler summation technique
  9- Power series radius of convergence
 10"""
 11
 12import numpy as np
 13
 14try:
 15    import matplotlib.pyplot as plt
 16    HAS_MATPLOTLIB = True
 17except ImportError:
 18    HAS_MATPLOTLIB = False
 19    print("Warning: matplotlib not available, skipping visualizations")
 20
 21
 22def ratio_test(a_n_func, n_max=100):
 23    """
 24    Apply ratio test for convergence: lim |a_{n+1}/a_n|
 25    If limit < 1: converges, > 1: diverges, = 1: inconclusive
 26    """
 27    ratios = []
 28    for n in range(1, n_max):
 29        a_n = a_n_func(n)
 30        a_n1 = a_n_func(n + 1)
 31        if abs(a_n) > 1e-15:
 32            ratios.append(abs(a_n1 / a_n))
 33
 34    if ratios:
 35        limit = ratios[-1]
 36        return limit, ratios
 37    return None, []
 38
 39
 40def root_test(a_n_func, n_max=100):
 41    """
 42    Apply root test for convergence: lim |a_n|^(1/n)
 43    If limit < 1: converges, > 1: diverges, = 1: inconclusive
 44    """
 45    roots = []
 46    for n in range(1, n_max):
 47        a_n = a_n_func(n)
 48        if a_n != 0:
 49            roots.append(abs(a_n) ** (1/n))
 50
 51    if roots:
 52        limit = roots[-1]
 53        return limit, roots
 54    return None, []
 55
 56
 57def partial_sum(a_n_func, n_terms):
 58    """Compute partial sum S_n = sum_{k=1}^{n} a_k"""
 59    total = sum(a_n_func(k) for k in range(1, n_terms + 1))
 60    return total
 61
 62
 63def taylor_series_exp(x, n_terms=10):
 64    """Taylor series for e^x centered at 0"""
 65    result = 0
 66    for n in range(n_terms):
 67        result += x**n / np.math.factorial(n)
 68    return result
 69
 70
 71def taylor_series_sin(x, n_terms=10):
 72    """Taylor series for sin(x) centered at 0"""
 73    result = 0
 74    for n in range(n_terms):
 75        sign = (-1)**n
 76        result += sign * x**(2*n + 1) / np.math.factorial(2*n + 1)
 77    return result
 78
 79
 80def taylor_series_cos(x, n_terms=10):
 81    """Taylor series for cos(x) centered at 0"""
 82    result = 0
 83    for n in range(n_terms):
 84        sign = (-1)**n
 85        result += sign * x**(2*n) / np.math.factorial(2*n)
 86    return result
 87
 88
 89def euler_summation(a_n_func, n_terms=50):
 90    """
 91    Euler summation for accelerating slowly convergent series
 92    Uses Euler transform: S ≈ (1/2)[S_even + S_odd]
 93    """
 94    partial_sums = [partial_sum(a_n_func, n) for n in range(1, n_terms + 1)]
 95
 96    # Apply Euler transform iteratively
 97    transformed = partial_sums[:]
 98    for iteration in range(5):
 99        new_transform = []
100        for i in range(len(transformed) - 1):
101            new_transform.append((transformed[i] + transformed[i+1]) / 2)
102        if not new_transform:
103            break
104        transformed = new_transform
105
106    return transformed[-1] if transformed else partial_sums[-1]
107
108
109def radius_of_convergence_ratio(c_n_func, n_test=100):
110    """
111    Compute radius of convergence using ratio test
112    R = lim |c_n / c_{n+1}|
113    """
114    ratios = []
115    for n in range(1, n_test):
116        c_n = c_n_func(n)
117        c_n1 = c_n_func(n + 1)
118        if abs(c_n1) > 1e-15:
119            ratios.append(abs(c_n / c_n1))
120
121    if ratios:
122        return ratios[-1]
123    return None
124
125
126# ============================================================================
127# MAIN DEMONSTRATIONS
128# ============================================================================
129
130print("=" * 70)
131print("INFINITE SERIES - CONVERGENCE TESTS AND TAYLOR EXPANSIONS")
132print("=" * 70)
133
134# Test 1: Ratio test on geometric series
135print("\n1. RATIO TEST - Geometric series a_n = (1/2)^n")
136print("-" * 70)
137geometric = lambda n: (1/2)**n
138limit, ratios = ratio_test(geometric, 50)
139print(f"Ratio test limit: {limit:.6f}")
140print(f"Expected: 0.5 (converges since < 1)")
141print(f"Partial sum (50 terms): {partial_sum(geometric, 50):.6f}")
142print(f"Exact sum: {1 / (1 - 0.5):.6f}")
143
144# Test 2: Ratio test on divergent series
145print("\n2. RATIO TEST - Divergent series a_n = n!")
146print("-" * 70)
147factorial_series = lambda n: np.math.factorial(n)
148# This grows too fast, test on smaller range
149limit_approx = factorial_series(11) / factorial_series(10)
150print(f"Ratio |a_11 / a_10| = {limit_approx:.1f}")
151print(f"Expected: grows to infinity (diverges)")
152
153# Test 3: Root test
154print("\n3. ROOT TEST - Series a_n = (1/3)^n")
155print("-" * 70)
156series_root = lambda n: (1/3)**n
157limit, roots = root_test(series_root, 50)
158print(f"Root test limit: {limit:.6f}")
159print(f"Expected: 0.333... (converges since < 1)")
160
161# Test 4: Alternating harmonic series
162print("\n4. ALTERNATING HARMONIC SERIES - sum (-1)^{n+1}/n")
163print("-" * 70)
164alt_harmonic = lambda n: ((-1)**(n+1)) / n
165s_100 = partial_sum(alt_harmonic, 100)
166s_1000 = partial_sum(alt_harmonic, 1000)
167print(f"Partial sum (100 terms): {s_100:.6f}")
168print(f"Partial sum (1000 terms): {s_1000:.6f}")
169print(f"Converges to ln(2) = {np.log(2):.6f}")
170
171# Test 5: Euler summation for ln(2)
172print("\n5. EULER SUMMATION - Accelerating convergence to ln(2)")
173print("-" * 70)
174euler_result = euler_summation(alt_harmonic, 50)
175print(f"Standard partial sum (50 terms): {partial_sum(alt_harmonic, 50):.6f}")
176print(f"Euler summation (50 terms): {euler_result:.6f}")
177print(f"Exact value ln(2): {np.log(2):.6f}")
178
179# Test 6: Taylor series for e^x
180print("\n6. TAYLOR SERIES - e^x at x=1")
181print("-" * 70)
182x_val = 1.0
183for n in [5, 10, 20]:
184    approx = taylor_series_exp(x_val, n)
185    error = abs(approx - np.exp(x_val))
186    print(f"n={n:2d} terms: {approx:.10f}, error: {error:.2e}")
187print(f"Exact e^1: {np.exp(x_val):.10f}")
188
189# Test 7: Taylor series for sin(x)
190print("\n7. TAYLOR SERIES - sin(x) at x=π/4")
191print("-" * 70)
192x_val = np.pi / 4
193for n in [3, 5, 10]:
194    approx = taylor_series_sin(x_val, n)
195    error = abs(approx - np.sin(x_val))
196    print(f"n={n:2d} terms: {approx:.10f}, error: {error:.2e}")
197print(f"Exact sin(π/4): {np.sin(x_val):.10f}")
198
199# Test 8: Taylor series for cos(x)
200print("\n8. TAYLOR SERIES - cos(x) at x=π/6")
201print("-" * 70)
202x_val = np.pi / 6
203for n in [3, 5, 10]:
204    approx = taylor_series_cos(x_val, n)
205    error = abs(approx - np.cos(x_val))
206    print(f"n={n:2d} terms: {approx:.10f}, error: {error:.2e}")
207print(f"Exact cos(π/6): {np.cos(x_val):.10f}")
208
209# Test 9: Radius of convergence
210print("\n9. RADIUS OF CONVERGENCE - Power series sum c_n x^n")
211print("-" * 70)
212# Series: sum (1/n^2) x^n
213c_n = lambda n: 1 / n**2
214R = radius_of_convergence_ratio(c_n, 100)
215print(f"Series: sum (1/n^2) x^n")
216print(f"Radius of convergence R = {R:.6f}")
217print(f"Expected: R = 1")
218
219# Series: sum (n!) x^n
220c_n_fact = lambda n: np.math.factorial(n) if n <= 20 else np.inf
221print(f"\nSeries: sum (n!) x^n")
222print(f"Radius of convergence R = 0 (series diverges for all x ≠ 0)")
223
224# Series: sum (1/n!) x^n (exponential)
225c_n_exp = lambda n: 1 / np.math.factorial(n)
226R_exp = radius_of_convergence_ratio(c_n_exp, 50)
227print(f"\nSeries: sum (1/n!) x^n")
228print(f"Radius of convergence R → ∞ (ratio → 0)")
229
230# Visualization
231if HAS_MATPLOTLIB:
232    print("\n" + "=" * 70)
233    print("VISUALIZATIONS")
234    print("=" * 70)
235
236    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
237
238    # Plot 1: Convergence of geometric series
239    ax = axes[0, 0]
240    n_vals = np.arange(1, 51)
241    partial_sums = [partial_sum(geometric, n) for n in n_vals]
242    ax.plot(n_vals, partial_sums, 'b-', linewidth=2, label='Partial sums')
243    ax.axhline(y=1.0, color='r', linestyle='--', label='Limit = 1')
244    ax.set_xlabel('Number of terms')
245    ax.set_ylabel('Partial sum')
246    ax.set_title('Convergence: Geometric series (1/2)^n')
247    ax.legend()
248    ax.grid(True, alpha=0.3)
249
250    # Plot 2: Taylor series for e^x
251    ax = axes[0, 1]
252    x_range = np.linspace(-2, 2, 100)
253    ax.plot(x_range, np.exp(x_range), 'k-', linewidth=2, label='e^x (exact)')
254    for n in [3, 5, 10]:
255        y_approx = [taylor_series_exp(x, n) for x in x_range]
256        ax.plot(x_range, y_approx, '--', label=f'n={n} terms')
257    ax.set_xlabel('x')
258    ax.set_ylabel('y')
259    ax.set_title('Taylor Series: e^x')
260    ax.legend()
261    ax.grid(True, alpha=0.3)
262    ax.set_ylim([-1, 8])
263
264    # Plot 3: Taylor series for sin(x)
265    ax = axes[1, 0]
266    x_range = np.linspace(-2*np.pi, 2*np.pi, 200)
267    ax.plot(x_range, np.sin(x_range), 'k-', linewidth=2, label='sin(x) exact')
268    for n in [2, 4, 6]:
269        y_approx = [taylor_series_sin(x, n) for x in x_range]
270        ax.plot(x_range, y_approx, '--', label=f'n={n} terms')
271    ax.set_xlabel('x')
272    ax.set_ylabel('y')
273    ax.set_title('Taylor Series: sin(x)')
274    ax.legend()
275    ax.grid(True, alpha=0.3)
276    ax.set_ylim([-2, 2])
277
278    # Plot 4: Convergence of alternating harmonic series
279    ax = axes[1, 1]
280    n_vals = np.arange(1, 101)
281    partial_sums = [partial_sum(alt_harmonic, n) for n in n_vals]
282    ax.plot(n_vals, partial_sums, 'b-', linewidth=1.5, label='Partial sums')
283    ax.axhline(y=np.log(2), color='r', linestyle='--', label='Limit = ln(2)')
284    ax.set_xlabel('Number of terms')
285    ax.set_ylabel('Partial sum')
286    ax.set_title('Convergence: Alternating harmonic series')
287    ax.legend()
288    ax.grid(True, alpha=0.3)
289
290    plt.tight_layout()
291    plt.savefig('/opt/projects/01_Personal/03_Study/examples/Mathematical_Methods/01_series_convergence.png', dpi=150)
292    print("Saved visualization: 01_series_convergence.png")
293    plt.close()
294
295print("\n" + "=" * 70)
296print("DEMONSTRATION COMPLETE")
297print("=" * 70)