02_complex_numbers.py

Download
python 326 lines 9.2 KB
  1"""
  2Complex Numbers - Arithmetic, Polar Form, and Conformal Mappings
  3
  4This script demonstrates:
  5- Complex arithmetic operations
  6- Polar and exponential forms
  7- Euler's formula
  8- Roots of unity
  9- Conformal mappings (z^2, 1/z, exp(z))
 10- Simple Mandelbrot set visualization
 11"""
 12
 13import numpy as np
 14
 15try:
 16    import matplotlib.pyplot as plt
 17    HAS_MATPLOTLIB = True
 18except ImportError:
 19    HAS_MATPLOTLIB = False
 20    print("Warning: matplotlib not available, skipping visualizations")
 21
 22
 23def complex_to_polar(z):
 24    """Convert complex number to polar form (r, theta)"""
 25    r = abs(z)
 26    theta = np.angle(z)
 27    return r, theta
 28
 29
 30def polar_to_complex(r, theta):
 31    """Convert polar form to complex number"""
 32    return r * np.exp(1j * theta)
 33
 34
 35def nth_roots_of_unity(n):
 36    """Compute n-th roots of unity: exp(2πik/n) for k=0,1,...,n-1"""
 37    roots = []
 38    for k in range(n):
 39        theta = 2 * np.pi * k / n
 40        root = np.exp(1j * theta)
 41        roots.append(root)
 42    return roots
 43
 44
 45def conformal_map_z_squared(z):
 46    """Conformal mapping w = z^2"""
 47    return z**2
 48
 49
 50def conformal_map_inverse(z):
 51    """Conformal mapping w = 1/z"""
 52    if z == 0:
 53        return np.inf
 54    return 1 / z
 55
 56
 57def conformal_map_exp(z):
 58    """Conformal mapping w = exp(z)"""
 59    return np.exp(z)
 60
 61
 62def mandelbrot_iteration(c, max_iter=100):
 63    """
 64    Mandelbrot iteration: z_{n+1} = z_n^2 + c
 65    Returns number of iterations before |z| > 2
 66    """
 67    z = 0
 68    for n in range(max_iter):
 69        if abs(z) > 2:
 70            return n
 71        z = z**2 + c
 72    return max_iter
 73
 74
 75# ============================================================================
 76# MAIN DEMONSTRATIONS
 77# ============================================================================
 78
 79print("=" * 70)
 80print("COMPLEX NUMBERS - ARITHMETIC, POLAR FORM, AND MAPPINGS")
 81print("=" * 70)
 82
 83# Test 1: Complex arithmetic
 84print("\n1. COMPLEX ARITHMETIC")
 85print("-" * 70)
 86z1 = 3 + 4j
 87z2 = 1 - 2j
 88print(f"z1 = {z1}")
 89print(f"z2 = {z2}")
 90print(f"z1 + z2 = {z1 + z2}")
 91print(f"z1 - z2 = {z1 - z2}")
 92print(f"z1 * z2 = {z1 * z2}")
 93print(f"z1 / z2 = {z1 / z2:.4f}")
 94print(f"z1 conjugate = {np.conj(z1)}")
 95print(f"|z1| = {abs(z1):.4f}")
 96
 97# Test 2: Polar form
 98print("\n2. POLAR FORM AND EULER'S FORMULA")
 99print("-" * 70)
100z = 1 + 1j
101r, theta = complex_to_polar(z)
102print(f"z = {z}")
103print(f"Polar form: r = {r:.4f}, θ = {theta:.4f} rad ({np.degrees(theta):.2f}°)")
104print(f"Verification: r*e^(iθ) = {polar_to_complex(r, theta):.4f}")
105print(f"\nEuler's formula: e^(iπ) + 1 = {np.exp(1j * np.pi) + 1:.10f}")
106print(f"e^(iπ/2) = {np.exp(1j * np.pi / 2):.4f} (should be i)")
107
108# Test 3: Roots of unity
109print("\n3. ROOTS OF UNITY")
110print("-" * 70)
111for n in [3, 4, 5]:
112    print(f"\n{n}-th roots of unity:")
113    roots = nth_roots_of_unity(n)
114    for k, root in enumerate(roots):
115        r, theta = complex_to_polar(root)
116        print(f"  ω_{k} = {root:.4f}, |ω| = {r:.4f}, arg = {np.degrees(theta):6.2f}°")
117    # Verify they are roots
118    verification = sum(root**n for root in roots)
119    print(f"  Verification: sum(ω^{n}) = {verification:.6f} (should be {n})")
120
121# Test 4: De Moivre's theorem
122print("\n4. DE MOIVRE'S THEOREM")
123print("-" * 70)
124z = np.exp(1j * np.pi / 6)  # e^(iπ/6)
125n = 3
126z_n_direct = z**n
127r, theta = complex_to_polar(z)
128z_n_demoivre = polar_to_complex(r**n, n * theta)
129print(f"z = e^(iπ/6) = {z:.4f}")
130print(f"n = {n}")
131print(f"z^n (direct) = {z_n_direct:.4f}")
132print(f"z^n (De Moivre) = {z_n_demoivre:.4f}")
133print(f"Expected: e^(iπ/2) = {np.exp(1j * np.pi / 2):.4f}")
134
135# Test 5: Complex logarithm
136print("\n5. COMPLEX LOGARITHM")
137print("-" * 70)
138z = 1 + 1j
139log_z = np.log(z)
140print(f"z = {z}")
141print(f"ln(z) = {log_z:.4f}")
142r, theta = complex_to_polar(z)
143print(f"ln(z) = ln(r) + iθ = {np.log(r):.4f} + {theta:.4f}i")
144print(f"Verification: e^(ln(z)) = {np.exp(log_z):.4f}")
145
146# Test 6: Conformal mappings
147print("\n6. CONFORMAL MAPPINGS")
148print("-" * 70)
149
150# w = z^2
151z_test = 1 + 1j
152w = conformal_map_z_squared(z_test)
153print(f"\nMapping w = z^2:")
154print(f"z = {z_test} → w = {w:.4f}")
155r, theta = complex_to_polar(z_test)
156r_w, theta_w = complex_to_polar(w)
157print(f"|z| = {r:.4f}, arg(z) = {np.degrees(theta):.2f}°")
158print(f"|w| = {r_w:.4f}, arg(w) = {np.degrees(theta_w):.2f}°")
159print(f"Note: |w| = |z|^2, arg(w) = 2*arg(z)")
160
161# w = 1/z
162z_test = 2 + 0j
163w = conformal_map_inverse(z_test)
164print(f"\nMapping w = 1/z:")
165print(f"z = {z_test} → w = {w:.4f}")
166z_test = 1j
167w = conformal_map_inverse(z_test)
168print(f"z = {z_test} → w = {w:.4f}")
169
170# w = exp(z)
171z_test = 1 + 1j * np.pi
172w = conformal_map_exp(z_test)
173print(f"\nMapping w = exp(z):")
174print(f"z = 1 + iπ → w = {w:.4f}")
175print(f"Expected: e^1 * e^(iπ) = e * (-1) = {np.e * (-1):.4f}")
176
177# Test 7: Mandelbrot set sample points
178print("\n7. MANDELBROT SET - SAMPLE POINTS")
179print("-" * 70)
180test_points = [
181    (0 + 0j, "Origin"),
182    (-1 + 0j, "(-1, 0)"),
183    (0.25 + 0j, "(0.25, 0)"),
184    (-0.5 + 0.5j, "(-0.5, 0.5)"),
185    (1 + 0j, "(1, 0) - outside"),
186]
187
188for c, description in test_points:
189    iterations = mandelbrot_iteration(c, max_iter=100)
190    in_set = "IN SET" if iterations == 100 else "ESCAPES"
191    print(f"c = {c:>10} ({description:20s}): {iterations:3d} iterations - {in_set}")
192
193# Visualization
194if HAS_MATPLOTLIB:
195    print("\n" + "=" * 70)
196    print("VISUALIZATIONS")
197    print("=" * 70)
198
199    fig = plt.figure(figsize=(14, 10))
200
201    # Plot 1: Roots of unity
202    ax1 = plt.subplot(2, 3, 1, projection='polar')
203    for n, color in [(3, 'r'), (4, 'b'), (5, 'g')]:
204        roots = nth_roots_of_unity(n)
205        theta_vals = [np.angle(r) for r in roots]
206        r_vals = [abs(r) for r in roots]
207        ax1.plot(theta_vals + [theta_vals[0]], r_vals + [r_vals[0]],
208                'o-', label=f'n={n}', markersize=8, color=color)
209    ax1.set_title('Roots of Unity')
210    ax1.legend()
211    ax1.grid(True)
212
213    # Plot 2: Conformal mapping z^2
214    ax2 = plt.subplot(2, 3, 2)
215    x = np.linspace(-2, 2, 20)
216    y = np.linspace(-2, 2, 20)
217
218    # Draw grid in z-plane
219    for xi in x[::4]:
220        z_line = xi + 1j * y
221        w_line = conformal_map_z_squared(z_line)
222        ax2.plot(w_line.real, w_line.imag, 'b-', alpha=0.5, linewidth=0.5)
223    for yi in y[::4]:
224        z_line = x + 1j * yi
225        w_line = conformal_map_z_squared(z_line)
226        ax2.plot(w_line.real, w_line.imag, 'r-', alpha=0.5, linewidth=0.5)
227
228    ax2.set_xlabel('Re(w)')
229    ax2.set_ylabel('Im(w)')
230    ax2.set_title('Conformal Map: w = z²')
231    ax2.grid(True, alpha=0.3)
232    ax2.axis('equal')
233    ax2.set_xlim(-5, 5)
234    ax2.set_ylim(-5, 5)
235
236    # Plot 3: Conformal mapping 1/z
237    ax3 = plt.subplot(2, 3, 3)
238    x = np.linspace(-2, 2, 20)
239    y = np.linspace(-2, 2, 20)
240
241    for xi in x[::4]:
242        if abs(xi) > 0.1:
243            z_line = xi + 1j * y
244            w_line = 1 / z_line
245            w_line = w_line[np.abs(w_line) < 10]
246            ax3.plot(w_line.real, w_line.imag, 'b-', alpha=0.5, linewidth=0.5)
247
248    for yi in y[::4]:
249        if abs(yi) > 0.1:
250            z_line = x + 1j * yi
251            w_line = 1 / z_line
252            w_line = w_line[np.abs(w_line) < 10]
253            ax3.plot(w_line.real, w_line.imag, 'r-', alpha=0.5, linewidth=0.5)
254
255    ax3.set_xlabel('Re(w)')
256    ax3.set_ylabel('Im(w)')
257    ax3.set_title('Conformal Map: w = 1/z')
258    ax3.grid(True, alpha=0.3)
259    ax3.axis('equal')
260    ax3.set_xlim(-3, 3)
261    ax3.set_ylim(-3, 3)
262
263    # Plot 4: Mandelbrot set
264    ax4 = plt.subplot(2, 3, 4)
265    xmin, xmax, ymin, ymax = -2.5, 1.0, -1.25, 1.25
266    width, height = 400, 300
267
268    x_vals = np.linspace(xmin, xmax, width)
269    y_vals = np.linspace(ymin, ymax, height)
270    mandelbrot_set = np.zeros((height, width))
271
272    for i, y in enumerate(y_vals):
273        for j, x in enumerate(x_vals):
274            c = x + 1j * y
275            mandelbrot_set[i, j] = mandelbrot_iteration(c, max_iter=50)
276
277    ax4.imshow(mandelbrot_set, extent=[xmin, xmax, ymin, ymax],
278               cmap='hot', origin='lower', interpolation='bilinear')
279    ax4.set_xlabel('Re(c)')
280    ax4.set_ylabel('Im(c)')
281    ax4.set_title('Mandelbrot Set')
282
283    # Plot 5: Complex exponential spiral
284    ax5 = plt.subplot(2, 3, 5)
285    t = np.linspace(0, 4*np.pi, 500)
286    z = (0.1 + 0.1j) * t * np.exp(1j * t)
287    ax5.plot(z.real, z.imag, 'b-', linewidth=1.5)
288    ax5.set_xlabel('Re(z)')
289    ax5.set_ylabel('Im(z)')
290    ax5.set_title('Complex Exponential Spiral')
291    ax5.grid(True, alpha=0.3)
292    ax5.axis('equal')
293
294    # Plot 6: Euler's identity visualization
295    ax6 = plt.subplot(2, 3, 6)
296    theta = np.linspace(0, 2*np.pi, 100)
297    z_circle = np.exp(1j * theta)
298    ax6.plot(z_circle.real, z_circle.imag, 'b-', linewidth=2, label='|z|=1')
299
300    # Mark special points
301    special_points = [
302        (0, 1, '0'),
303        (np.pi/2, 1j, 'π/2'),
304        (np.pi, -1, 'π'),
305        (3*np.pi/2, -1j, '3π/2'),
306    ]
307    for angle, z, label in special_points:
308        ax6.plot(z.real, z.imag, 'ro', markersize=8)
309        ax6.text(z.real*1.2, z.imag*1.2, f'e^(i{label})', ha='center')
310
311    ax6.set_xlabel('Real')
312    ax6.set_ylabel('Imaginary')
313    ax6.set_title("Euler's Formula: e^(iθ)")
314    ax6.grid(True, alpha=0.3)
315    ax6.axis('equal')
316    ax6.legend()
317
318    plt.tight_layout()
319    plt.savefig('/opt/projects/01_Personal/03_Study/examples/Mathematical_Methods/02_complex_numbers.png', dpi=150)
320    print("Saved visualization: 02_complex_numbers.png")
321    plt.close()
322
323print("\n" + "=" * 70)
324print("DEMONSTRATION COMPLETE")
325print("=" * 70)