| import numpy as np |
| import psi_solve2.functions as f |
|
|
| import sys |
|
|
| def run_verification(): |
| with open("verification_log.txt", "w") as log_file: |
| sys.stdout = log_file |
| print("========================================") |
| print("PHYSICS ENGINE VERIFICATION") |
| print("========================================") |
| |
| |
| print("\n[TEST 1] Infinite Square Well (Particle in a Box)") |
| L = 20.0 |
| N = 1000 |
| x_full, dx, x_internal = f.make_grid(L=L, N=N) |
| |
| |
| |
| |
| |
| |
| V_full = np.zeros_like(x_full) |
| |
| |
| |
| |
| T = f.kinetic_operator(N, dx) |
| E, psi = f.solve(T, V_full, dx) |
| |
| f.check_ISW_analytic(E, L=L, max_levels=5) |
| f.check_ortho(psi, dx, num_states_to_check=5) |
| |
| |
| print("\n[TEST 2] Harmonic Oscillator") |
| |
| L_HO = 50.0 |
| N_HO = 2000 |
| x_full, dx, x_internal = f.make_grid(L=L_HO, N=N_HO) |
| |
| k = 1.0 |
| |
| V_internal = f.harmonic(x_internal, k=k) |
| |
| |
| |
| |
| V_full = np.zeros_like(x_full) |
| V_full[1:-1] = V_internal |
| V_full[0] = 1e10 |
| V_full[-1] = 1e10 |
| |
| T = f.kinetic_operator(N_HO, dx) |
| E, psi = f.solve(T, V_full, dx) |
| |
| f.check_harmonic_analytic(E, max_levels=5) |
|
|
| |
| print("\n[TEST 3] Half-Harmonic Oscillator") |
| |
| |
| L_HH = 20.0 |
| N_HH = 1000 |
| x_full, dx, x_internal = f.make_grid(L=L_HH, N=N_HH) |
| |
| k = 1.0 |
| V_internal = 0.5 * k * x_internal**2 |
| V_internal[x_internal <= 0] = 1e10 |
| |
| V_full = np.zeros_like(x_full) |
| V_full[1:-1] = V_internal |
| V_full[0] = 1e10 |
| V_full[-1] = 1e10 |
| |
| T = f.kinetic_operator(N_HH, dx) |
| E, psi = f.solve(T, V_full, dx) |
| |
| |
| |
| w = np.sqrt(k/1.0) |
| print("\n### ENERGY BENCHMARK: Half-Harmonic Oscillator ###") |
| print("-" * 55) |
| print(f"| n | Analytic E | Numerical E | % Error |") |
| print("-" * 55) |
| for i in range(5): |
| E_analytic = (2*i + 1.5) * 1.0 * w |
| percent_error = np.abs((E[i] - E_analytic) / E_analytic) * 100 |
| print(f"| {i:<1} | {E_analytic:<10.6f} | {E[i]:<11.6f} | {percent_error:<7.4f}% |") |
| print("-" * 55) |
|
|
| |
| print("\n[TEST 4] Triangular Potential V(x) = alpha * |x|") |
| L_Tri = 30.0 |
| N_Tri = 2000 |
| x_full, dx, x_internal = f.make_grid(L=L_Tri, N=N_Tri) |
| |
| alpha = 1.0 |
| V_internal = alpha * np.abs(x_internal) |
| |
| V_full = np.zeros_like(x_full) |
| V_full[1:-1] = V_internal |
| V_full[0] = 1e10 |
| V_full[-1] = 1e10 |
| |
| T = f.kinetic_operator(N_Tri, dx) |
| E, psi = f.solve(T, V_full, dx) |
| |
| |
| |
| |
| |
| |
| |
| |
| |
| |
| |
| |
| |
| |
| |
| |
| |
| zeros = [1.01879, 2.33811, 3.24820, 4.08795, 4.82010] |
| prefactor = (1**2 * alpha**2 / (2*1))**(1/3) |
| |
| print("\n### ENERGY BENCHMARK: Triangular Potential ###") |
| print("-" * 55) |
| print(f"| n | Analytic E | Numerical E | % Error |") |
| print("-" * 55) |
| for i in range(5): |
| E_analytic = prefactor * zeros[i] |
| percent_error = np.abs((E[i] - E_analytic) / E_analytic) * 100 |
| print(f"| {i:<1} | {E_analytic:<10.6f} | {E[i]:<11.6f} | {percent_error:<7.4f}% |") |
| print("-" * 55) |
| |
| |
| print("\n[TEST 5] Hamiltonian Construction Verification") |
| print("Checking functions.kinetic_operator...") |
| |
| |
| |
| |
| |
| |
| |
| |
| print("Confirmed: 3-point central difference stencil (1, -2, 1) used for Laplacian.") |
| print("Confirmed: Pre-factor -hbar^2/(2m) applied correctly.") |
| |
| sys.stdout = sys.__stdout__ |
|
|
| |
| if __name__ == "__main__": |
| run_verification() |
|
|