Search code examples
python-3.xmatplotlibmathpolynomial-approximations

Neukum production function


I am trying to plot a function which looks like that: Equation I want to plot. Where, D range from 0.01 to 300km, N is number of craters with dia > D per sq.km per giga year, coefficients values are in code.

I get empty plot and obviously I am doing something very wrong which I am not able to understand. I am sharing my code.

'''

a_0 = -3.0876
a_1 = -3.557528
a_2 = 0.781027
a_3 = 1.021521
a_4 = -0.156012
a_5 = -0.444058
a_6 = 0.019977
a_7 = 0.086850
a_8 = -0.005874
a_9 = -0.006809
a_10 = 8.25*10**-4
a_11 = 5.54*10**-5
a_n = a_0 + a_1 + a_2 + a_3 + a_4 + a_5 + a_6 + a_7 + a_8 + a_9 + a_10 + a_11

for d in range(1, 200000):
    n = a_0 + np.multiply(a_1, np.log(d))
    + np.multiply(a_2, np.log(d)**2)
    + np.multiply(a_3, np.log(d)**3)
    + np.multiply(a_4, np.log(d)**4)
    + np.multiply(a_5, np.log(d)**5)
    + np.multiply(a_6, np.log(d)**6)
    + np.multiply(a_7, np.log(d)**7)
    + np.multiply(a_8, np.log(d)**8)
    + np.multiply(a_9, np.log(d)**9)
    + np.multiply(a_10, np.log(d)**10)
    + np.multiply(a_11, np.log(d)**11)
print(10**n)
plt.plot(10**n, color='green')
plt.show()

''' The graphical curve should look something like thatResultant plot of equation


Solution

  • I can not understand your code, but if you want to plot the function you wrote you can try this:

    import numpy as np
    import matplotlib.pyplot as plt
    a_0 = -3.0876
    a_1 = -3.557528
    a_2 = 0.781027
    a_3 = 1.021521
    a_4 = -0.156012
    a_5 = -0.444058
    a_6 = 0.019977
    a_7 = 0.086850
    a_8 = -0.005874
    a_9 = -0.006809
    a_10 = 8.25*10**(-4)
    a_11 = 5.54*10**(-5)
    
    N_POINTS = 100
    a_coeff = np.array([a_0, a_1, a_2, a_3, a_4, a_5, a_6, a_7, a_8, a_9, a_10, a_11])
    distance = np.logspace(np.log10(0.01), np.log10(200), N_POINTS)
    
    exponents = np.arange(12)
    distance_matrix = distance[:, np.newaxis]*np.ones([N_POINTS, 12])
    N = 10**np.sum(a_coeff * (np.log10(distance_matrix)**exponents), axis=1)
    
    fig, ax = plt.subplots(1, figsize=(3, 8))
    ax.scatter(distance, N)
    ax.set(xscale='log', yscale='log', ylim=(1e-7, 1e4), xlim=(0.001, 300))
    

    enter image description here