from numpy import * from matplotlib.pyplot import * from scipy.interpolate import interp1d def SciPySpline(x, y): return interp1d(x,y,'cubic',bounds_error=False) class SplineInterpolation: def __init__(self, x, y, leftddx = 0, rightddx = 0): n = len(x) if n < 3: raise Exception("need at least three support points to compute a spline") self.x = x.copy() self.y = y.copy() self.n = n # compute lambda, mu and S mu = zeros(n-1) mu[1:]=(x[1:-1]-x[:-2])/(x[2:]-x[:-2]) l = zeros(n-1) l[1:]=(x[2:]-x[1:-1])/(x[2:]-x[:-2]) S = zeros(n) S[1:-1]=6*(((y[2:] -y[1:-1])/(x[2:] -x[1:-1]))- \ ((y[1:-1]-y[:-2]) /(x[1:-1]-x[:-2]))) / (x[2:]-x[:-2]) # second moment boundary conditions S[0] = leftddx S[n-1] = rightddx # construct coeff matrix tmp = zeros((n, n)) tmp[0,0] = 1 for i in range(1, n-1): tmp[i,i] = 2 tmp[i,i-1] = mu[i] tmp[i,i+1] = l[i] tmp[n-1,n-1] = 1 # compute coeff self.M = linalg.solve(tmp, S) # compute the dependent coeffs self.alpha = (self.M[1:]-self.M[:-1])/(self.x[1:]-self.x[:-1]) self.m = (y[1:]-y[:-1])/(x[1:]-x[:-1]) - \ 1./6.*(x[1:]-x[:-1])*(2*self.M[:-1]+self.M[1:]) def __call__(self, x): y = zeros_like(x) for i in range(0, len(x)): index=0 while index < len(self.x) - 1 and self.x[index] <= x[i]: index += 1 if index > 0: index -= 1 y[i] = self.y[index] + \ self.m[index]*(x[i] - self.x[index]) + \ 0.5*self.M[index]*(x[i] - self.x[index])**2 + \ 1./6.*self.alpha[index]*(x[i] - self.x[index])**3 return array(y) # smoothest spline from SciPy spline = SciPySpline # natural spline spline = SplineInterpolation def f(x): return sin(x) x5 = linspace(0, 2.*pi, 5) fp5 = spline(x5, f(x5)) x10 = linspace(0, 2.*pi, 10) fp10 = spline(x10, f(x10)) x15 = linspace(0, 2.*pi, 15) fp15 = spline(x15, f(x15)) x = linspace(0, 2.*pi, 1000) figure() plot(x, f(x), label='f(x)') plot(x, fp5(x), label='5 points') plot(x, fp10(x), label='10 points') plot(x, fp15(x), label='15 points') legend() def g(x): return 1/(1+x*x) x5 = linspace(-5., 5., 5) gp5 = spline(x5, g(x5)) x10 = linspace(-5., 5., 10) gp10 = spline(x10, g(x10)) x15 = linspace(-5., 5., 15) gp15 = spline(x15, g(x15)) x = linspace(-5, 5, 1000) figure() plot(x, g(x), label='g(x)') plot(x, gp5(x), label='5 points') plot(x, gp10(x), label='10 points') plot(x, gp15(x), label='15 points') legend() def h(x): return x**-12 - x**-6 x5 = linspace(1, 5., 5) hp5 = spline(x5, h(x5)) x10 = linspace(1, 5, 10) hp10 = spline(x10, h(x10)) x15 = linspace(1, 5., 15) hp15 = spline(x15, h(x15)) x = linspace(1, 5, 1000) figure() plot(x, h(x), label='h(x)') plot(x, hp5(x), label='5 points') plot(x, hp10(x), label='10 points') plot(x, hp15(x), label='15 points') legend() show()