1"""
2Complex Analysis - Cauchy Integral, Residues, Laurent Series, Analytic Continuation
3
4This script demonstrates:
5- Cauchy integral formula (numerical)
6- Residue theorem and computation
7- Laurent series expansion
8- Poles and essential singularities
9- Analytic continuation concepts
10- Conformal mapping properties
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 cauchy_integral_formula(f, z0, radius=1.0, n_points=1000):
24 """
25 Cauchy integral formula: f(z0) = (1/2πi) ∮ f(z)/(z-z0) dz
26 Integrate around circle |z - z0| = radius
27 """
28 theta = np.linspace(0, 2*np.pi, n_points)
29 z = z0 + radius * np.exp(1j * theta)
30 dz = 1j * radius * np.exp(1j * theta)
31
32 integrand = f(z) / (z - z0) * dz
33 integral = np.sum(integrand) * (2*np.pi / n_points)
34
35 # Divide by 2πi
36 result = integral / (2j * np.pi)
37 return result
38
39
40def cauchy_derivative_formula(f, z0, n=1, radius=0.5, n_points=1000):
41 """
42 Cauchy formula for derivatives:
43 f^(n)(z0) = (n!/(2πi)) ∮ f(z)/(z-z0)^(n+1) dz
44 """
45 theta = np.linspace(0, 2*np.pi, n_points)
46 z = z0 + radius * np.exp(1j * theta)
47 dz = 1j * radius * np.exp(1j * theta)
48
49 integrand = f(z) / (z - z0)**(n+1) * dz
50 integral = np.sum(integrand) * (2*np.pi / n_points)
51
52 # Multiply by n!/(2πi)
53 result = np.math.factorial(n) * integral / (2j * np.pi)
54 return result
55
56
57def residue_at_pole(f, z0, order=1, h=1e-4):
58 """
59 Compute residue at pole z0 of order n
60 Res(f, z0) = lim_{z→z0} (1/(n-1)!) d^(n-1)/dz^(n-1) [(z-z0)^n f(z)]
61 """
62 if order == 1:
63 # Simple pole: Res = lim_{z→z0} (z-z0)f(z)
64 z_near = z0 + h
65 residue = (z_near - z0) * f(z_near)
66 return residue
67 else:
68 # Higher order pole: use numerical differentiation
69 g = lambda z: (z - z0)**order * f(z)
70
71 # Numerical (n-1)th derivative
72 derivative = g(z0)
73 for k in range(1, order):
74 # Use central difference
75 derivative = (g(z0 + h) - g(z0 - h)) / (2 * h)
76
77 residue = derivative / np.math.factorial(order - 1)
78 return residue
79
80
81def contour_integral(f, path_func, t_range, n_points=1000):
82 """
83 Compute contour integral ∫_C f(z) dz
84 path_func: parametric path z(t)
85 """
86 t = np.linspace(t_range[0], t_range[1], n_points)
87 dt = (t_range[1] - t_range[0]) / (n_points - 1)
88
89 z = path_func(t)
90 # Compute dz/dt numerically
91 dz_dt = np.gradient(z, dt)
92
93 integrand = f(z) * dz_dt
94 integral = np.trapz(integrand, t)
95
96 return integral
97
98
99def laurent_series_coefficients(f, z0, r_inner, r_outer, n_terms=10):
100 """
101 Compute Laurent series coefficients:
102 f(z) = ∑_{n=-∞}^{∞} a_n (z-z0)^n
103 a_n = (1/2πi) ∮ f(z)/(z-z0)^(n+1) dz
104 """
105 coeffs = {}
106
107 # Positive powers (analytic part)
108 radius = (r_inner + r_outer) / 2
109 n_points = 1000
110 theta = np.linspace(0, 2*np.pi, n_points)
111
112 for n in range(-n_terms, n_terms + 1):
113 z = z0 + radius * np.exp(1j * theta)
114 dz = 1j * radius * np.exp(1j * theta)
115
116 integrand = f(z) / (z - z0)**(n+1) * dz
117 integral = np.sum(integrand) * (2*np.pi / n_points)
118
119 coeffs[n] = integral / (2j * np.pi)
120
121 return coeffs
122
123
124# ============================================================================
125# MAIN DEMONSTRATIONS
126# ============================================================================
127
128print("=" * 70)
129print("COMPLEX ANALYSIS - CAUCHY, RESIDUES, LAURENT SERIES")
130print("=" * 70)
131
132# Test 1: Cauchy integral formula
133print("\n1. CAUCHY INTEGRAL FORMULA")
134print("-" * 70)
135
136# f(z) = z^2
137f_poly = lambda z: z**2
138z0 = 1 + 1j
139
140result = cauchy_integral_formula(f_poly, z0, radius=0.5)
141exact = f_poly(z0)
142
143print(f"f(z) = z²")
144print(f"z₀ = {z0}")
145print(f"Cauchy formula: f(z₀) = {result:.6f}")
146print(f"Direct evaluation: {exact:.6f}")
147print(f"Error: {abs(result - exact):.2e}")
148
149# f(z) = e^z
150f_exp = lambda z: np.exp(z)
151z0 = 0 + 0j
152
153result = cauchy_integral_formula(f_exp, z0, radius=1.0)
154exact = f_exp(z0)
155
156print(f"\nf(z) = e^z")
157print(f"z₀ = {z0}")
158print(f"Cauchy formula: f(z₀) = {result:.6f}")
159print(f"Direct evaluation: {exact:.6f}")
160print(f"Error: {abs(result - exact):.2e}")
161
162# Test 2: Cauchy derivative formula
163print("\n2. CAUCHY DERIVATIVE FORMULA")
164print("-" * 70)
165
166f_exp = lambda z: np.exp(z)
167z0 = 0 + 0j
168
169for n in range(1, 4):
170 result = cauchy_derivative_formula(f_exp, z0, n=n, radius=0.5)
171 exact = np.exp(z0) # All derivatives of e^z equal e^z
172
173 print(f"f^({n})(0) for f(z)=e^z:")
174 print(f" Cauchy formula: {result:.6f}")
175 print(f" Exact: {exact:.6f}")
176 print(f" Error: {abs(result - exact):.2e}")
177
178# Test 3: Residue theorem
179print("\n3. RESIDUE THEOREM")
180print("-" * 70)
181
182# f(z) = 1/z (simple pole at z=0)
183print("f(z) = 1/z, simple pole at z=0")
184f_simple_pole = lambda z: 1 / z if abs(z) > 1e-10 else np.inf
185residue = residue_at_pole(f_simple_pole, 0 + 0j, order=1, h=0.01)
186print(f"Residue at z=0: {residue:.6f}")
187print(f"Expected: 1.0")
188
189# Verify with contour integral around circle
190circle_path = lambda t: 1.0 * np.exp(1j * t)
191integral = contour_integral(f_simple_pole, circle_path, (0, 2*np.pi))
192print(f"∮ f(z)dz around |z|=1: {integral:.6f}")
193print(f"2πi × Residue = {2j * np.pi * residue:.6f}")
194
195# f(z) = 1/(z-1)(z-2) with poles at z=1 and z=2
196print("\nf(z) = 1/[(z-1)(z-2)]")
197f_two_poles = lambda z: 1 / ((z - 1) * (z - 2)) if abs((z-1)*(z-2)) > 1e-10 else np.inf
198
199# Residues at each pole
200res_1 = residue_at_pole(lambda z: 1/(z-2), 1 + 0j, order=1, h=0.01)
201res_2 = residue_at_pole(lambda z: 1/(z-1), 2 + 0j, order=1, h=0.01)
202
203print(f"Residue at z=1: {res_1:.6f} (expected: -1.0)")
204print(f"Residue at z=2: {res_2:.6f} (expected: 1.0)")
205
206# Contour enclosing both poles
207circle_large = lambda t: 3.0 * np.exp(1j * t)
208integral_large = contour_integral(f_two_poles, circle_large, (0, 2*np.pi))
209print(f"∮ f(z)dz around |z|=3: {integral_large:.6f}")
210print(f"2πi × (Res₁ + Res₂) = {2j * np.pi * (res_1 + res_2):.6f}")
211
212# Test 4: Laurent series
213print("\n4. LAURENT SERIES")
214print("-" * 70)
215
216# f(z) = 1/z(z-1) expanded around z=0
217print("f(z) = 1/[z(z-1)], expanded around z=0")
218print("Region: 0 < |z| < 1")
219
220f_laurent = lambda z: 1 / (z * (z - 1)) if abs(z * (z-1)) > 1e-10 else np.inf
221
222# For this function: 1/[z(z-1)] = -1/z - 1/(z-1) = -1/z + 1/(1-z)
223# In region 0 < |z| < 1: = -1/z + ∑(z^n) = -1/z + 1 + z + z² + ...
224
225coeffs = laurent_series_coefficients(f_laurent, 0 + 0j, 0.1, 0.9, n_terms=5)
226
227print("\nLaurent coefficients a_n:")
228for n in range(-5, 6):
229 if n in coeffs:
230 coeff = coeffs[n]
231 if abs(coeff) > 1e-6:
232 print(f" a_{n:2d} = {coeff.real:8.4f} + {coeff.imag:8.4f}i")
233
234# Test 5: Essential singularity
235print("\n5. ESSENTIAL SINGULARITY")
236print("-" * 70)
237
238# f(z) = e^(1/z) has essential singularity at z=0
239print("f(z) = e^(1/z), essential singularity at z=0")
240
241f_essential = lambda z: np.exp(1/z) if abs(z) > 1e-10 else np.inf
242
243# Laurent series: e^(1/z) = ∑_{n=0}^∞ (1/z)^n / n! = ∑_{n=-∞}^0 z^n / (-n)!
244coeffs_essential = laurent_series_coefficients(f_essential, 0 + 0j, 0.2, 0.8, n_terms=5)
245
246print("\nLaurent coefficients (principal part is infinite):")
247for n in range(-5, 1):
248 if n in coeffs_essential:
249 coeff = coeffs_essential[n]
250 expected = 1 / np.math.factorial(-n) if n <= 0 else 0
251 print(f" a_{n:2d} = {coeff.real:10.6f} (expected: {expected:.6f})")
252
253# Residue (a_{-1})
254residue_essential = coeffs_essential.get(-1, 0)
255print(f"\nResidue (a_-1): {residue_essential.real:.6f}")
256print(f"Expected: 1.0")
257
258# Test 6: Conformal mapping properties
259print("\n6. CONFORMAL MAPPING PROPERTIES")
260print("-" * 70)
261
262# Test that analytic functions preserve angles
263# f(z) = z²
264f_map = lambda z: z**2
265
266z0 = 1 + 0.5j
267h = 0.01
268
269# Two directions from z0
270direction1 = h * (1 + 0j) # Horizontal
271direction2 = h * (0 + 1j) # Vertical
272
273# Map to w-plane
274w0 = f_map(z0)
275w1 = f_map(z0 + direction1)
276w2 = f_map(z0 + direction2)
277
278# Angle in z-plane
279angle_z = np.angle(direction2 / direction1)
280
281# Angle in w-plane
282dw1 = w1 - w0
283dw2 = w2 - w0
284angle_w = np.angle(dw2 / dw1)
285
286print(f"Mapping: w = z²")
287print(f"Point: z₀ = {z0}")
288print(f"\nAngles between directions:")
289print(f" In z-plane: {np.degrees(angle_z):.2f}°")
290print(f" In w-plane: {np.degrees(angle_w):.2f}°")
291print(f" Difference: {np.degrees(abs(angle_z - angle_w)):.2f}° (should be ≈0)")
292
293# Test 7: Analytic continuation
294print("\n7. ANALYTIC CONTINUATION CONCEPT")
295print("-" * 70)
296
297# Example: f(z) = ∑ z^n (geometric series)
298# Converges for |z| < 1, equals 1/(1-z)
299# Can be continued to entire complex plane except z=1
300
301print("f(z) = ∑_{n=0}^∞ z^n, |z| < 1")
302print("Analytic continuation: f(z) = 1/(1-z), z ≠ 1")
303
304z_inside = 0.5 + 0j
305z_outside = 1.5 + 0j
306
307# Series approximation (works only inside disk)
308def geometric_series(z, n_terms=50):
309 if abs(z) >= 1:
310 return np.nan
311 return sum(z**n for n in range(n_terms))
312
313series_val = geometric_series(z_inside)
314continued_val = 1 / (1 - z_inside)
315
316print(f"\nInside disk |z| < 1:")
317print(f" z = {z_inside}")
318print(f" Series sum: {series_val:.6f}")
319print(f" Continuation: {continued_val:.6f}")
320
321print(f"\nOutside disk |z| > 1:")
322print(f" z = {z_outside}")
323print(f" Series: diverges")
324print(f" Continuation: 1/(1-z) = {1/(1-z_outside):.6f}")
325
326# Visualization
327if HAS_MATPLOTLIB:
328 print("\n" + "=" * 70)
329 print("VISUALIZATIONS")
330 print("=" * 70)
331
332 fig = plt.figure(figsize=(15, 10))
333
334 # Plot 1: Contour for Cauchy integral
335 ax1 = plt.subplot(2, 3, 1)
336 z0 = 1 + 0.5j
337 radius = 0.7
338 theta = np.linspace(0, 2*np.pi, 100)
339 z_circle = z0 + radius * np.exp(1j * theta)
340
341 ax1.plot(z_circle.real, z_circle.imag, 'b-', linewidth=2, label='Contour')
342 ax1.plot(z0.real, z0.imag, 'ro', markersize=10, label='z₀')
343 ax1.arrow(z_circle.real[25], z_circle.imag[25],
344 z_circle.real[26]-z_circle.real[25],
345 z_circle.imag[26]-z_circle.imag[25],
346 head_width=0.1, head_length=0.05, fc='blue', ec='blue')
347 ax1.set_xlabel('Re(z)')
348 ax1.set_ylabel('Im(z)')
349 ax1.set_title('Cauchy Integral Formula Contour')
350 ax1.legend()
351 ax1.grid(True, alpha=0.3)
352 ax1.set_aspect('equal')
353
354 # Plot 2: Poles and residues
355 ax2 = plt.subplot(2, 3, 2)
356 poles = [1 + 0j, 2 + 0j]
357 circle_large = 3.0 * np.exp(1j * theta)
358
359 ax2.plot(circle_large.real, circle_large.imag, 'b-', linewidth=2, label='Contour')
360 for i, pole in enumerate(poles):
361 ax2.plot(pole.real, pole.imag, 'rx', markersize=15, markeredgewidth=3)
362 ax2.text(pole.real, pole.imag + 0.3, f'z={pole.real:.0f}', ha='center')
363
364 ax2.set_xlabel('Re(z)')
365 ax2.set_ylabel('Im(z)')
366 ax2.set_title('Residue Theorem: f(z)=1/[(z-1)(z-2)]')
367 ax2.grid(True, alpha=0.3)
368 ax2.set_aspect('equal')
369 ax2.set_xlim(-3.5, 3.5)
370 ax2.set_ylim(-3.5, 3.5)
371
372 # Plot 3: Laurent series annulus
373 ax3 = plt.subplot(2, 3, 3)
374 r_inner = 0.3
375 r_outer = 0.9
376 circle_inner = r_inner * np.exp(1j * theta)
377 circle_outer = r_outer * np.exp(1j * theta)
378
379 ax3.fill_between(circle_outer.real, circle_outer.imag,
380 alpha=0.3, label='Annulus')
381 ax3.plot(circle_inner.real, circle_inner.imag, 'r--', linewidth=2)
382 ax3.plot(circle_outer.real, circle_outer.imag, 'b--', linewidth=2)
383 ax3.plot(0, 0, 'ko', markersize=8, label='Singularity')
384
385 ax3.set_xlabel('Re(z)')
386 ax3.set_ylabel('Im(z)')
387 ax3.set_title('Laurent Series: Region of Convergence')
388 ax3.legend()
389 ax3.grid(True, alpha=0.3)
390 ax3.set_aspect('equal')
391
392 # Plot 4: Conformal mapping w = z²
393 ax4 = plt.subplot(2, 3, 4)
394
395 # Grid in z-plane
396 x_lines = np.linspace(-1.5, 1.5, 7)
397 y_lines = np.linspace(-1.5, 1.5, 7)
398
399 for x in x_lines:
400 y = np.linspace(-1.5, 1.5, 100)
401 z = x + 1j * y
402 w = z**2
403 ax4.plot(w.real, w.imag, 'b-', alpha=0.5, linewidth=0.8)
404
405 for y in y_lines:
406 x = np.linspace(-1.5, 1.5, 100)
407 z = x + 1j * y
408 w = z**2
409 ax4.plot(w.real, w.imag, 'r-', alpha=0.5, linewidth=0.8)
410
411 ax4.set_xlabel('Re(w)')
412 ax4.set_ylabel('Im(w)')
413 ax4.set_title('Conformal Map: w = z²')
414 ax4.grid(True, alpha=0.3)
415 ax4.set_aspect('equal')
416
417 # Plot 5: Analytic continuation
418 ax5 = plt.subplot(2, 3, 5)
419
420 # Unit circle (convergence of series)
421 circle_unit = np.exp(1j * theta)
422 ax5.plot(circle_unit.real, circle_unit.imag, 'b-', linewidth=2,
423 label='|z| = 1 (series boundary)')
424 ax5.fill(circle_unit.real, circle_unit.imag, alpha=0.2, color='blue')
425
426 # Pole at z=1
427 ax5.plot(1, 0, 'rx', markersize=15, markeredgewidth=3, label='Pole at z=1')
428
429 # Sample points
430 z_in = 0.5 + 0j
431 z_out = 1.5 + 0j
432 ax5.plot(z_in.real, z_in.imag, 'go', markersize=10, label='Series works')
433 ax5.plot(z_out.real, z_out.imag, 'mo', markersize=10, label='Need continuation')
434
435 ax5.set_xlabel('Re(z)')
436 ax5.set_ylabel('Im(z)')
437 ax5.set_title('Analytic Continuation: ∑z^n = 1/(1-z)')
438 ax5.legend()
439 ax5.grid(True, alpha=0.3)
440 ax5.set_aspect('equal')
441 ax5.set_xlim(-1.5, 2)
442 ax5.set_ylim(-1.5, 1.5)
443
444 # Plot 6: Branch cut (example: log(z))
445 ax6 = plt.subplot(2, 3, 6)
446
447 # Draw branch cut along negative real axis
448 ax6.plot([-2, 0], [0, 0], 'r-', linewidth=4, label='Branch cut')
449 ax6.plot(0, 0, 'ko', markersize=8, label='Branch point')
450
451 # Contour avoiding branch cut
452 theta_avoid = np.linspace(0.1, 2*np.pi - 0.1, 100)
453 z_avoid = 1.5 * np.exp(1j * theta_avoid)
454 ax6.plot(z_avoid.real, z_avoid.imag, 'b--', linewidth=2,
455 label='Valid contour')
456
457 ax6.set_xlabel('Re(z)')
458 ax6.set_ylabel('Im(z)')
459 ax6.set_title('Branch Cut: log(z)')
460 ax6.legend()
461 ax6.grid(True, alpha=0.3)
462 ax6.set_aspect('equal')
463 ax6.set_xlim(-2.5, 2.5)
464 ax6.set_ylim(-2.5, 2.5)
465
466 plt.tight_layout()
467 plt.savefig('/opt/projects/01_Personal/03_Study/examples/Mathematical_Methods/11_complex_analysis.png', dpi=150)
468 print("Saved visualization: 11_complex_analysis.png")
469 plt.close()
470
471print("\n" + "=" * 70)
472print("DEMONSTRATION COMPLETE")
473print("=" * 70)