Search code examples
matplotlibinterpolationastronomycubic-spline

Cubic spline interpolation drops out halfway

I am trying to make a cubic spline interpolation and for some reason, the interpolation drops off in the middle of it. It's very mysterious and I can't find any mention of similar occurrences anywhere online.

This is for my dissertation so I have excluded some labels etc. to keep it obscure intentionally, but all the relevant code is as follows. For context, this is an astronomy related plot.

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

W = np.array([0.435,0.606,0.814,1.05,1.25,1.40,1.60])
sum_all = np.array([sum435,sum606,sum814,sum105,sum125,sum140,sum160])
sum_can = np.array([sumc435,sumc606,sumc814,sumc105,sumc125,sumc140,sumc160])

fall = CubicSpline(W,sum_all)
newallx=np.arange(0.435,1.6,0.001)
newally=fall(newallx)

fcan = CubicSpline(W,sum_can)
newcanx=np.arange(0.435,1.6,0.001)
newcany=fcan(newcanx)

#----plot

plt.plot(newallx,newally)
plt.plot(newcanx,newcany)
plt.plot(W,sum_all,marker='o',color='r',linestyle='')
plt.plot(W,sum_can,marker='o',color='b',linestyle='')
plt.yscale("log")
plt.ylabel("Flux S$_v$ [erg s$^-$$^1$ cm$^-$$^2$ Hz$^-$$^1$]")
plt.xlabel("Wavelength [n$\lambda$]")
plt.show()

The plot that I get from that comes out like this, with a clear gap in the interpolation:

img

And in case you are wondering, these are the values in the sum_all and sum_can arrays (I assume it doesn't matter, but just in case you want the numbers to plot it yourself):

sum_all:
[  3.87282732e+32   8.79993191e+32   1.74866333e+33   1.59946687e+33
   9.08556547e+33   6.70458731e+33   9.84832359e+33]
can_all:
[  2.98381061e+28   1.26194810e+28   3.30328780e+28   2.90254609e+29
   3.65117723e+29   3.46256846e+29   3.64483736e+29]

The gap happens between [0.606,1.26194810e+28] and [0.814,3.30328780e+28]. If I change the intervals from 0.001 to something higher, it's obvious that the plot doesn't actually break off but merely dips below 0 on the y-axis (but the plot is continuous). So why does it do that? Surely that's not a correct interpolation? Just looking with our eyes, that's clearly not a well-interpolated connection between those two points.

Any tips or comments would be extremely appreciated. Thank you so much in advance!

Solution

  • The reason for the breakdown can be better observed on a linear scale.

    enter image description here

    We see that the spline actually passes below 0, which is undefined on a log scale.

    So I would suggest to first take the logarithm of the data, perform the spline interpolation on the logarithmically scaled data, and then scale back by the 10th power.

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

W = np.array([0.435,0.606,0.814,1.05,1.25,1.40,1.60])
sum_all = np.array([  3.87282732e+32,   8.79993191e+32,   1.74866333e+33,   1.59946687e+33,
   9.08556547e+33,   6.70458731e+33,   9.84832359e+33])
sum_can = np.array([  2.98381061e+28,   1.26194810e+28,   3.30328780e+28,   2.90254609e+29,
   3.65117723e+29,   3.46256846e+29,   3.64483736e+29])

fall = CubicSpline(W,np.log10(sum_all))
newallx=np.arange(0.435,1.6,0.001)
newally=fall(newallx)

fcan = CubicSpline(W,np.log10(sum_can))
newcanx=np.arange(0.435,1.6,0.01)
newcany=fcan(newcanx)


plt.plot(newallx,10**newally)
plt.plot(newcanx,10**newcany)
plt.plot(W,sum_all,marker='o',color='r',linestyle='')
plt.plot(W,sum_can,marker='o',color='b',linestyle='')
plt.yscale("log")

plt.ylabel("Flux S$_v$ [erg s$^-$$^1$ cm$^-$$^2$ Hz$^-$$^1$]")
plt.xlabel("Wavelength [n$\lambda$]")
plt.show()

    enter image description here