Search code examples
pythonscipyinterpolationradius

scipy interpolate gives varying results

I'm trying to get the same answer when reading of an interpolated function in Python using scipy.interpolate.interp1d but when I change the size of the x linspace I get different results.

Below is a simplified case where I feed the interpolated functions different radii and they return drastically different results. I can't work out why this is happening so any help would be greatly appreciated. 

from scipy.interpolate import interp1d
import numpy as np
import matplotlib.pyplot as plt

plt.close('all')

M_centre = 2e30
G = 1.67e-11
m_test = 6e24
radius = np.linspace(5,1e2,1000)
radius2 = np.linspace(5,1e21,1000)
V_circ = np.sqrt(G*M_centre/radius)
V_circ2 = np.sqrt(G*M_centre/radius2)

velocities_circ = interp1d(radius,V_circ)
test_r = velocities_circ(50)
print(test_r)

velocities_circ2 = interp1d(radius2,V_circ2)
test_r2 = velocities_circ2(50)
print(test_r2)


Out: 
817312853.7629617
2584569596.664017

I've thought about maybe the step size of the linspace is causing the varying reading on the interpolated function but it surely can't vary by an order of magnitude can it?

Edit: I have also tried this method using numpy.interp but it gives the same results as above.

Solution

  • Just a quick illustration of the problem with reduced numbers for contrast:

    from scipy.interpolate import interp1d
import numpy as np
import matplotlib.pyplot as plt

M_centre = 2e30
G = 1.67e-11
m_test = 6e24

radius1 = np.linspace(5,1e3,10)
radius2 = np.linspace(5,1e2,10)

V_circ1 = np.sqrt(G*M_centre/radius1)
V_circ2 = np.sqrt(G*M_centre/radius2)

velocities_circ1 = interp1d(radius1,V_circ1)
test_r1 = velocities_circ1(50)
print(test_r1)

velocities_circ2 = interp1d(radius2,V_circ2)
test_r2 = velocities_circ2(50)
print(test_r2)

plt.plot(radius1, V_circ1, "ro", label = "radius1")
plt.plot(radius2, V_circ2, "bx", label = "radius2")

plt.plot(radius1, velocities_circ1(radius1), "r")
plt.plot(radius2, velocities_circ2(radius2), "b")
plt.legend()
plt.xlim(0, 400)
plt.show()

    I think the reason for the different output is obvious. enter image description here

    And the equivalent diagram for same range (5, 1e2), but different number of points (3 vs 10):

    enter image description here