I'm trying to write a Python program that computes the half-life of 212Po by using the Quantum tunneling model and dividing the potential barrier into multiple segments. I started by setting the number of segments to 2 in order to be able to print and see how my matrix is constructed, I compute the Coulomb potential at each segment and then the values of k at each segment, then I simply construct the matrix and calculate the half-life using the relation between the transmission coefficient and the decay time (Here I set the incoming wave from the left to A=1, and the one from left to 0, here is my code: import numpy as np #Constants hbar_c = 197.327 # MeV·fm alpha = 1 / 137.036 amu_to_MeV = 931.494 c = 299792458 #m / s # Constants for 212Po Z_alpha = 2 Z_daughter = 82 A_alpha = 4 A_daughter = 208 Q = 8.954 # MeV R_atomic = 9.01 R_max = 26.9 N = 2 def reduced_mass(A1, A2): # MeV/c^2 return (A1 * A2) / (A1 + A2) * amu_to_MeV def coulomb_potential(r, Z1, Z2): #MeV return alpha * Z1 * Z2 * hbar_c / r def compute_k(mu, Q, V): V = np.append(V, (Q-0.01)) k = np.zeros_like(V, dtype=complex) for i in range(len(V)): if Q >= V: k = (1j * np.sqrt(2 * mu * (Q - V)) / hbar_c ) else: k = (np.sqrt(2 * mu * (V - Q)) / hbar_c) return k def build_system_matrix(N, x, k): size = 2 * N + 3 M = np.zeros((size, size), dtype=complex) b = np.zeros(size, dtype=complex) # Incoming wave A = 1 b[0] = 1 M[0, 0] = 1 # A #Outermost left boundary M[1, 0] = np.exp(k[0] * x[0]) M[1, 1] = np.exp(-k[0] * x[0]) M[1, 2] = -np.exp(k[1] * x[0]) M[1, 3] = -np.exp(-k[1] * x[0]) M[2, 0] = k[0] * np.exp(k[0] * x[0]) M[2, 1] = -k[0] * np.exp(-k[0] * x[0]) M[2, 2] = -k[1] * np.exp(k[1] * x[0]) M[2, 3] = k[1] * np.exp(-k[1] * x[0]) # Internal boundaries for i in range(1, N): xi = x idx1 = 2 * i + 1 idx2 = 2 * i + 2 ai = 2 * i bi = 2 * i + 1 aip1 = 2 * i + 2 bip1 = 2 * i + 3 e_l = np.exp(k * xi) e_lm = np.exp(-k * xi) e_r = np.exp(k[i+1] * xi) e_rm = np.exp(-k[i+1] * xi) M[idx1, ai] = e_l M[idx1, bi] = e_lm M[idx1, aip1] = -e_r M[idx1, bip1] = -e_rm M[idx2, ai] = k * e_l M[idx2, bi] = -k * e_lm M[idx2, aip1] = -k[i+1] * e_r M[idx2, bip1] = k[i+1] * e_rm # Last boundary xf = x[-1] M[-2, -3] = np.exp(k[-2] * xf) M[-2, -2] = np.exp(-k[-2] * xf) M[-2, -1] = -np.exp(k[-1] * xf) M[-1, -3] = k[-2] * np.exp(k[-2] * xf) M[-1, -2] = -k[-2] * np.exp(-k[-2] * xf) M[-1, -1] = -k[-1] * np.exp(k[-1] * xf) return M, b def get_v_f(m, E): #fm, s^-1 v_by_c_2 = ( 2 * E ) / m v = np.sqrt(v_by_c_2) * c v_fm = v * (10**15) f = v_fm / (2*R_atomic) return v_fm , f def compute_lifetime(T, f): tau = 1 / (T * f) t_half = tau * np.log(2) return tau, t_half def run_decay_simulation(Q_value, N, R, R_max): #ns mu = reduced_mass(A_alpha, A_daughter) x = np.linspace(R, R_max, N+1) dx = x[1] - x[0] V = coulomb_potential(x[:-1] + dx/2, Z_alpha, Z_daughter) k = compute_k(mu, Q, V) M, b = build_system_matrix(N, x, k) solution = np.linalg.solve(M, b) F = solution[-1] T = abs(F)**2 v ,f = get_v_f((A_alpha*amu_to_MeV), Q) tau, t_half = compute_lifetime(T, f) t_half_ns = t_half * 10**(9) return t_half_ns t_half = run_decay_simulation(Q, N, R_atomic, R_max) #mu = reduced_mass(A_alpha, A_daughter) #x = np.linspace(R_atomic, R_max, N+1) #dx = x[1] - x[0] #V = coulomb_potential(x[:-1] + dx/2, Z_alpha, Z_daughter) #k = compute_k(mu, Q, V) #M, b = build_system_matrix(N, x, k) #solution = np.linalg.solve(M, b) But for some reason my solution is blowing up directly after the first iteration, i.e the coefficients gets too big, here is the solution array for N=2: solution = [1.00000000e+00+4.88145863e-17j, 2.09421821e+10+4.14232046e+05j, 1.39847295e-01-2.24207875e-02j, 4.30599920e+07+1.01357042e+03j, 6.29560298e-10+7.27660554e-10j, 3.77028289e+03+3.16638841e+03j, 4.01750469e-03-1.64008154e-03j] Where did I go wrong wrong in my code? Continue reading...
