I’m working on a data that involves computations with very large numbers, and I’m encountering an issue with numerical overflow and plotting. The problem arises when I increase the value of a parameter 𝑁 (at the beginning of the code).

For example, with 𝑁=5, my code runs perfectly without any errors and the curve is generated. However, as I gradually increase 𝑁, parts of the plotted curve start to disappear. By the time 𝑁 reaches 20, the curve disappears entirely. the problem is that I should use values for N=1000 or 2000.

the following error occurs
RuntimeWarning: overflow encountered in double_scalars,

Here’s the python code:

import numpy as np
import cmath
import math
from math import e
import matplotlib.pyplot as plt
import mpmath
from scipy.special import gamma, factorial
import scipy.special as sc

N = 5 # VARIABLE OF INTEREST (SHOULD BE INCREASED)



g1, g2= 1.31928e-06, 4.947323e-07

H, ms = 1/6, 6

x_axis_points = np.arange(0, 40, 1)
E1 = 10.0 ** ((x_axis_points-30) / 10.0)/1e-11
y_axis_points = []
for QS in E1:
 
    A, B = N * 0.979405**2, N * 0.07986762

    w= (B/A)**2*(A**2/B)*(A**2/B+1)
    r= (B/A)**4*(A**2/B)*(A**2/B+1)*(A**2/B+2)*(A**2/B+3)


    f= w**2/ (r-w**2)
    c= g2**2 * QS*(r-w**2)/w
  


    BB=H *QS * g1**2
    L1, L2= 1/BB, 1/c


################### THE PROBLEM APPEARS IN THE COMMING LINES
(GIVES INFINITY WHEN INCREASING N) ##############

cons1= L1**(N*ms)/ sc.gamma(N*ms)  
cons2= L2**(f)/ sc.gamma(f  )    

sup3_num= L2**(f)*sc.gamma(N*ms+1   +  f)
sup3_den= (N*ms+1)*( L2 + L1)**( N*ms+1 +f)               
sup3=mpmath.hyp2f1(1, N*ms+1 + f, N*ms+1+1, L1/(L1+L2))

sup4_num= L1**(N*ms)*sc.gamma(f+1   +  N*ms)
sup4_den= (f+1)*( L2 + L1)**( f+1 +N*ms)                   
sup4=mpmath.hyp2f1(1, f+1 + N*ms, f+1+1, L2/(L1+L2))



Z1= cons1* sup3_num/sup3_den * sup3/ sc.gamma(f)

Z2= cons2* sup4_num/sup4_den * sup4/ sc.gamma(N*ms)   

y = math.log2(1+ Z1+Z2 ) 

y_axis_points.append((y))            

plt.plot(x_axis_points, y_axis_points,    linewidth='2.5', color='blue')

plt.xlabel('x')
plt.ylabel('y')
plt.grid()

plt.xlim([0, 40])
plt.ylim([0 , 3 ])
plt.show()