1"""
2Special Functions - Bessel, Legendre, Hermite, Laguerre, Spherical Harmonics
3
4This script demonstrates:
5- Bessel functions J_n(x) and Y_n(x)
6- Legendre polynomials P_n(x)
7- Hermite polynomials H_n(x)
8- Laguerre polynomials L_n(x)
9- Spherical harmonics Y_l^m(θ,φ)
10- Gamma function Γ(x)
11- Orthogonality properties
12"""
13
14import numpy as np
15
16try:
17 from scipy.special import (jv, yn, eval_legendre, eval_hermite,
18 eval_laguerre, sph_harm, gamma, factorial)
19 HAS_SCIPY = True
20except ImportError:
21 HAS_SCIPY = False
22 print("Warning: scipy not available, using limited implementations")
23 # Provide basic gamma function
24 gamma = lambda x: np.exp(np.sum(np.log(np.arange(1, int(x)))))
25
26try:
27 import matplotlib.pyplot as plt
28 from mpl_toolkits.mplot3d import Axes3D
29 HAS_MATPLOTLIB = True
30except ImportError:
31 HAS_MATPLOTLIB = False
32 print("Warning: matplotlib not available, skipping visualizations")
33
34
35def legendre_poly_manual(n, x):
36 """Compute Legendre polynomial P_n(x) using recurrence relation"""
37 if n == 0:
38 return np.ones_like(x)
39 elif n == 1:
40 return x
41
42 P_prev2 = np.ones_like(x)
43 P_prev1 = x
44
45 for k in range(2, n + 1):
46 P_n = ((2*k - 1) * x * P_prev1 - (k - 1) * P_prev2) / k
47 P_prev2 = P_prev1
48 P_prev1 = P_n
49
50 return P_n
51
52
53def hermite_poly_manual(n, x):
54 """Compute Hermite polynomial H_n(x) using recurrence relation"""
55 if n == 0:
56 return np.ones_like(x)
57 elif n == 1:
58 return 2 * x
59
60 H_prev2 = np.ones_like(x)
61 H_prev1 = 2 * x
62
63 for k in range(2, n + 1):
64 H_n = 2 * x * H_prev1 - 2 * (k - 1) * H_prev2
65 H_prev2 = H_prev1
66 H_prev1 = H_n
67
68 return H_n
69
70
71def laguerre_poly_manual(n, x):
72 """Compute Laguerre polynomial L_n(x) using recurrence relation"""
73 if n == 0:
74 return np.ones_like(x)
75 elif n == 1:
76 return 1 - x
77
78 L_prev2 = np.ones_like(x)
79 L_prev1 = 1 - x
80
81 for k in range(2, n + 1):
82 L_n = ((2*k - 1 - x) * L_prev1 - (k - 1) * L_prev2) / k
83 L_prev2 = L_prev1
84 L_prev1 = L_n
85
86 return L_n
87
88
89def check_orthogonality_legendre(n, m, num_points=1000):
90 """
91 Check orthogonality of Legendre polynomials:
92 ∫_{-1}^{1} P_n(x) P_m(x) dx = 0 if n≠m, = 2/(2n+1) if n=m
93 """
94 x = np.linspace(-1, 1, num_points)
95 if HAS_SCIPY:
96 P_n = eval_legendre(n, x)
97 P_m = eval_legendre(m, x)
98 else:
99 P_n = legendre_poly_manual(n, x)
100 P_m = legendre_poly_manual(m, x)
101
102 integrand = P_n * P_m
103 integral = np.trapz(integrand, x)
104
105 return integral
106
107
108def check_orthogonality_hermite(n, m, num_points=500):
109 """
110 Check orthogonality of Hermite polynomials:
111 ∫_{-∞}^{∞} H_n(x) H_m(x) e^(-x²) dx = 0 if n≠m
112 """
113 x = np.linspace(-5, 5, num_points)
114 if HAS_SCIPY:
115 H_n = eval_hermite(n, x)
116 H_m = eval_hermite(m, x)
117 else:
118 H_n = hermite_poly_manual(n, x)
119 H_m = hermite_poly_manual(m, x)
120
121 weight = np.exp(-x**2)
122 integrand = H_n * H_m * weight
123 integral = np.trapz(integrand, x)
124
125 return integral
126
127
128# ============================================================================
129# MAIN DEMONSTRATIONS
130# ============================================================================
131
132print("=" * 70)
133print("SPECIAL FUNCTIONS - BESSEL, LEGENDRE, HERMITE, LAGUERRE")
134print("=" * 70)
135
136# Test 1: Bessel functions
137if HAS_SCIPY:
138 print("\n1. BESSEL FUNCTIONS J_n(x)")
139 print("-" * 70)
140 x_vals = [1.0, 5.0, 10.0]
141 n_vals = [0, 1, 2]
142
143 print(f"{'x':>6s}", end='')
144 for n in n_vals:
145 print(f"{'J_'+str(n)+'(x)':>12s}", end='')
146 print()
147 print("-" * 70)
148
149 for x in x_vals:
150 print(f"{x:6.1f}", end='')
151 for n in n_vals:
152 J_n = jv(n, x)
153 print(f"{J_n:12.6f}", end='')
154 print()
155
156 # Bessel function zeros (important for boundary value problems)
157 print("\nFirst 3 zeros of J_0(x):")
158 x_range = np.linspace(0.1, 15, 500)
159 J_0 = jv(0, x_range)
160
161 # Find approximate zeros
162 zeros = []
163 for i in range(1, len(x_range)):
164 if J_0[i-1] * J_0[i] < 0:
165 zeros.append(x_range[i])
166 if len(zeros) >= 3:
167 break
168
169 for i, zero in enumerate(zeros):
170 print(f" x_{i+1} ≈ {zero:.4f}")
171
172 # Test 2: Legendre polynomials
173 print("\n2. LEGENDRE POLYNOMIALS P_n(x)")
174 print("-" * 70)
175 x_vals = [-1.0, -0.5, 0.0, 0.5, 1.0]
176 n_vals = [0, 1, 2, 3]
177
178 print(f"{'x':>6s}", end='')
179 for n in n_vals:
180 print(f"{'P_'+str(n)+'(x)':>10s}", end='')
181 print()
182 print("-" * 70)
183
184 for x in x_vals:
185 print(f"{x:6.1f}", end='')
186 for n in n_vals:
187 P_n = eval_legendre(n, x)
188 print(f"{P_n:10.4f}", end='')
189 print()
190
191 # Orthogonality check
192 print("\nOrthogonality check: ∫₋₁¹ P_n(x)P_m(x)dx")
193 print("Expected: 0 if n≠m, 2/(2n+1) if n=m")
194 for n in range(3):
195 for m in range(3):
196 integral = check_orthogonality_legendre(n, m)
197 expected = 2/(2*n+1) if n == m else 0
198 print(f" ∫P_{n}P_{m} = {integral:8.4f} (expected: {expected:6.4f})")
199
200else:
201 print("\n1-2. LEGENDRE POLYNOMIALS (manual implementation)")
202 print("-" * 70)
203 x_vals = [0.0, 0.5, 1.0]
204 for x in x_vals:
205 print(f"\nx = {x}:")
206 for n in range(4):
207 P_n = legendre_poly_manual(n, np.array([x]))[0]
208 print(f" P_{n}({x}) = {P_n:.4f}")
209
210# Test 3: Hermite polynomials
211if HAS_SCIPY:
212 print("\n3. HERMITE POLYNOMIALS H_n(x)")
213 print("-" * 70)
214 x_vals = [-2.0, -1.0, 0.0, 1.0, 2.0]
215 n_vals = [0, 1, 2, 3]
216
217 print(f"{'x':>6s}", end='')
218 for n in n_vals:
219 print(f"{'H_'+str(n)+'(x)':>10s}", end='')
220 print()
221 print("-" * 70)
222
223 for x in x_vals:
224 print(f"{x:6.1f}", end='')
225 for n in n_vals:
226 H_n = eval_hermite(n, x)
227 print(f"{H_n:10.2f}", end='')
228 print()
229
230 # Connection to quantum harmonic oscillator
231 print("\nQuantum harmonic oscillator wavefunctions:")
232 print("ψ_n(x) ∝ H_n(x) exp(-x²/2)")
233 x = 0.0
234 for n in range(4):
235 H_n = eval_hermite(n, x)
236 psi_n = H_n * np.exp(-x**2 / 2)
237 print(f" ψ_{n}(0) ∝ {psi_n:.4f}")
238
239# Test 4: Laguerre polynomials
240if HAS_SCIPY:
241 print("\n4. LAGUERRE POLYNOMIALS L_n(x)")
242 print("-" * 70)
243 x_vals = [0.0, 1.0, 2.0, 3.0]
244 n_vals = [0, 1, 2, 3]
245
246 print(f"{'x':>6s}", end='')
247 for n in n_vals:
248 print(f"{'L_'+str(n)+'(x)':>10s}", end='')
249 print()
250 print("-" * 70)
251
252 for x in x_vals:
253 print(f"{x:6.1f}", end='')
254 for n in n_vals:
255 L_n = eval_laguerre(n, x)
256 print(f"{L_n:10.4f}", end='')
257 print()
258
259 # Connection to hydrogen atom
260 print("\nHydrogen atom radial wavefunctions involve Laguerre polynomials")
261 print("R_nl(r) ∝ r^l L_{n-l-1}^{2l+1}(2r/na₀) exp(-r/na₀)")
262
263# Test 5: Spherical harmonics
264if HAS_SCIPY:
265 print("\n5. SPHERICAL HARMONICS Y_l^m(θ,φ)")
266 print("-" * 70)
267 theta_vals = [0, np.pi/4, np.pi/2]
268 phi = 0
269
270 print("At φ=0:")
271 print(f"{'θ':>10s} {'Y_0^0':>15s} {'Y_1^0':>15s} {'Y_2^0':>15s}")
272 print("-" * 70)
273
274 for theta in theta_vals:
275 print(f"{theta:10.4f}", end='')
276 for l in [0, 1, 2]:
277 m = 0
278 Y_lm = sph_harm(m, l, phi, theta)
279 # Real part (imaginary part is zero for m=0)
280 print(f"{Y_lm.real:15.6f}", end='')
281 print()
282
283 # Check normalization
284 print("\nNormalization check (numerical integration):")
285 n_theta = 50
286 n_phi = 100
287 theta = np.linspace(0, np.pi, n_theta)
288 phi = np.linspace(0, 2*np.pi, n_phi)
289 Theta, Phi = np.meshgrid(theta, phi)
290
291 for l in [0, 1, 2]:
292 for m in range(-l, l+1):
293 Y_lm = sph_harm(m, l, Phi, Theta)
294 integrand = np.abs(Y_lm)**2 * np.sin(Theta)
295
296 dtheta = np.pi / (n_theta - 1)
297 dphi = 2*np.pi / (n_phi - 1)
298 integral = np.sum(integrand) * dtheta * dphi
299
300 print(f" ∫|Y_{l}^{m:2d}|² dΩ = {integral:.4f} (expected: 1.0)")
301
302# Test 6: Gamma function
303print("\n6. GAMMA FUNCTION Γ(x)")
304print("-" * 70)
305print("Γ(n+1) = n! for integer n")
306print(f"{'n':>4s} {'n!':>12s} {'Γ(n+1)':>12s}")
307print("-" * 70)
308
309for n in range(1, 8):
310 factorial_n = np.math.factorial(n)
311 gamma_n_plus_1 = gamma(n + 1)
312 print(f"{n:4d} {factorial_n:12.0f} {gamma_n_plus_1:12.6f}")
313
314print("\nΓ(x) for half-integer values:")
315print(f"Γ(1/2) = √π = {gamma(0.5):.6f} (expected: {np.sqrt(np.pi):.6f})")
316print(f"Γ(3/2) = √π/2 = {gamma(1.5):.6f} (expected: {np.sqrt(np.pi)/2:.6f})")
317print(f"Γ(5/2) = 3√π/4 = {gamma(2.5):.6f} (expected: {3*np.sqrt(np.pi)/4:.6f})")
318
319# Visualization
320if HAS_MATPLOTLIB and HAS_SCIPY:
321 print("\n" + "=" * 70)
322 print("VISUALIZATIONS")
323 print("=" * 70)
324
325 fig = plt.figure(figsize=(15, 12))
326
327 # Plot 1: Bessel functions J_n
328 ax1 = plt.subplot(3, 3, 1)
329 x = np.linspace(0, 15, 500)
330 for n in range(4):
331 J_n = jv(n, x)
332 ax1.plot(x, J_n, linewidth=2, label=f'$J_{n}(x)$')
333 ax1.set_xlabel('x')
334 ax1.set_ylabel('$J_n(x)$')
335 ax1.set_title('Bessel Functions of First Kind')
336 ax1.legend()
337 ax1.grid(True, alpha=0.3)
338 ax1.set_ylim(-0.5, 1.0)
339
340 # Plot 2: Bessel functions Y_n
341 ax2 = plt.subplot(3, 3, 2)
342 x = np.linspace(0.1, 15, 500)
343 for n in range(3):
344 Y_n = yn(n, x)
345 ax2.plot(x, Y_n, linewidth=2, label=f'$Y_{n}(x)$')
346 ax2.set_xlabel('x')
347 ax2.set_ylabel('$Y_n(x)$')
348 ax2.set_title('Bessel Functions of Second Kind')
349 ax2.legend()
350 ax2.grid(True, alpha=0.3)
351 ax2.set_ylim(-2, 1)
352
353 # Plot 3: Legendre polynomials
354 ax3 = plt.subplot(3, 3, 3)
355 x = np.linspace(-1, 1, 500)
356 for n in range(6):
357 P_n = eval_legendre(n, x)
358 ax3.plot(x, P_n, linewidth=2, label=f'$P_{n}(x)$')
359 ax3.set_xlabel('x')
360 ax3.set_ylabel('$P_n(x)$')
361 ax3.set_title('Legendre Polynomials')
362 ax3.legend()
363 ax3.grid(True, alpha=0.3)
364 ax3.set_ylim(-1.2, 1.2)
365
366 # Plot 4: Hermite polynomials
367 ax4 = plt.subplot(3, 3, 4)
368 x = np.linspace(-3, 3, 500)
369 for n in range(5):
370 H_n = eval_hermite(n, x)
371 ax4.plot(x, H_n, linewidth=2, label=f'$H_{n}(x)$')
372 ax4.set_xlabel('x')
373 ax4.set_ylabel('$H_n(x)$')
374 ax4.set_title('Hermite Polynomials')
375 ax4.legend()
376 ax4.grid(True, alpha=0.3)
377 ax4.set_ylim(-50, 50)
378
379 # Plot 5: Quantum harmonic oscillator wavefunctions
380 ax5 = plt.subplot(3, 3, 5)
381 x = np.linspace(-4, 4, 500)
382 for n in range(4):
383 H_n = eval_hermite(n, x)
384 # Normalized wavefunction
385 norm = 1 / np.sqrt(2**n * np.math.factorial(n) * np.sqrt(np.pi))
386 psi_n = norm * H_n * np.exp(-x**2 / 2)
387 ax5.plot(x, psi_n + n, linewidth=2, label=f'$\\psi_{n}(x)$')
388
389 ax5.set_xlabel('x')
390 ax5.set_ylabel('$\\psi_n(x)$ + offset')
391 ax5.set_title('Quantum Harmonic Oscillator')
392 ax5.legend()
393 ax5.grid(True, alpha=0.3)
394
395 # Plot 6: Laguerre polynomials
396 ax6 = plt.subplot(3, 3, 6)
397 x = np.linspace(0, 10, 500)
398 for n in range(5):
399 L_n = eval_laguerre(n, x)
400 ax6.plot(x, L_n, linewidth=2, label=f'$L_{n}(x)$')
401 ax6.set_xlabel('x')
402 ax6.set_ylabel('$L_n(x)$')
403 ax6.set_title('Laguerre Polynomials')
404 ax6.legend()
405 ax6.grid(True, alpha=0.3)
406 ax6.set_ylim(-5, 5)
407
408 # Plot 7: Spherical harmonic |Y_2^0|
409 ax7 = fig.add_subplot(3, 3, 7, projection='3d')
410 theta = np.linspace(0, np.pi, 50)
411 phi = np.linspace(0, 2*np.pi, 100)
412 Theta, Phi = np.meshgrid(theta, phi)
413
414 l, m = 2, 0
415 Y_lm = sph_harm(m, l, Phi, Theta)
416 r = np.abs(Y_lm)
417
418 x_sph = r * np.sin(Theta) * np.cos(Phi)
419 y_sph = r * np.sin(Theta) * np.sin(Phi)
420 z_sph = r * np.cos(Theta)
421
422 ax7.plot_surface(x_sph, y_sph, z_sph, cmap='viridis', alpha=0.8)
423 ax7.set_xlabel('x')
424 ax7.set_ylabel('y')
425 ax7.set_zlabel('z')
426 ax7.set_title('$|Y_2^0(\\theta,\\phi)|$')
427
428 # Plot 8: Gamma function
429 ax8 = plt.subplot(3, 3, 8)
430 x_pos = np.linspace(0.1, 5, 500)
431 y_gamma = gamma(x_pos)
432
433 ax8.plot(x_pos, y_gamma, 'b-', linewidth=2)
434 ax8.plot([1, 2, 3, 4], [gamma(1), gamma(2), gamma(3), gamma(4)],
435 'ro', markersize=8, label='Integer points')
436 ax8.set_xlabel('x')
437 ax8.set_ylabel('$\\Gamma(x)$')
438 ax8.set_title('Gamma Function')
439 ax8.legend()
440 ax8.grid(True, alpha=0.3)
441 ax8.set_ylim(0, 25)
442
443 # Plot 9: Associated Legendre functions (used in spherical harmonics)
444 ax9 = plt.subplot(3, 3, 9)
445 x = np.linspace(-1, 1, 500)
446 from scipy.special import lpmv
447
448 for m in range(3):
449 P_2m = lpmv(m, 2, x)
450 ax9.plot(x, P_2m, linewidth=2, label=f'$P_2^{m}(x)$')
451
452 ax9.set_xlabel('x')
453 ax9.set_ylabel('$P_l^m(x)$')
454 ax9.set_title('Associated Legendre Functions (l=2)')
455 ax9.legend()
456 ax9.grid(True, alpha=0.3)
457
458 plt.tight_layout()
459 plt.savefig('/opt/projects/01_Personal/03_Study/examples/Mathematical_Methods/08_special_functions.png', dpi=150)
460 print("Saved visualization: 08_special_functions.png")
461 plt.close()
462
463print("\n" + "=" * 70)
464print("DEMONSTRATION COMPLETE")
465print("=" * 70)