Untitled

895 days ago by pub

Given 3 points: x=(0,1,2) with y=(1.43,1.82,1.66).

Create a cubic spline with knots at every point and f''(endpoints)=0. Find the coefficients of the resulting cubic functions $p_0(x-x_0)$ and $p_1(x-x_1)$

This uses the tridiagonal algorithms from https://bitbucket.org/sunqiang/nmiewp/src/ac36b422a590/LUdecomp3.py

This uses the cubic spline algorithms from Jaan Kiusalaas EXCEPT integer division // in the while loop for FindSegment ... 

Coefficient formulas: http://math.fullerton.edu/mathews/n2003/splines/CubicSplinesProof.pdf 

import numpy as np import numpy.linalg 
       
# python tridiagonal algorithms as sourced from bitbucket.org def LUdecomp3(c,d,e): n = d.shape[0] for k in xrange(1,n): lam = c[k-1]/d[k-1] d[k] = d[k] - lam*e[k-1] c[k-1] = lam return c,d,e def LUsolve3(c,d,e,b): n = d.shape[0] for k in xrange(1,n): b[k] = b[k] - c[k-1]*b[k-1] b[n-1] = b[n-1]/d[n-1] for k in xrange(n-2,-1,-1): b[k] = (b[k] - e[k]*b[k+1])/d[k] return b 
       
# python cubic spline algorithms as sourced from JK def curvatures(xData,yData): n = len(xData) - 1 c = np.zeros(n) d = np.ones(n+1) e = np.zeros(n) k = np.zeros(n+1) c[0:n-1]= xData[0:n-1] - xData[1:n] d[1:n] = 2.0*(xData[0:n-1] - xData[2:n+1]) e[1:n] = xData[1:n] - xData[2:n+1] k[1:n] = 6.0*(yData[0:n-1] - yData[1:n]) \ /(xData[0:n-1] - xData[1:n]) \ -6.0*(yData[1:n] - yData[2:n+1]) \ /(xData[1:n] - xData[2:n+1]) LUdecomp3(c,d,e) LUsolve3(c,d,e,k) return k 
       
# python cubic spline algorithms as sourced from JK except // to force integer result # Here we have the curvatures k and we want to evaluate the cubic splines at a particular value x # First it locates the x-interval where the particular value of x is located (using bisection). # Then it calculates the cubic y-value using the "fast cubic" formula. def evalSpline(xData,yData,k,x): # this is a bisection algorithm to find in what interval x lies def findSegment(xData,x): iLeft = 0 iRight = len(xData)- 1 while 1: if (iRight-iLeft)<=1: return iLeft i=(iLeft+iRight)//2 if x < xData[i] : iRight=i else: iLeft=i i=findSegment(xData,x) h = xData[i]-xData[i+1] y = ((x-xData[i+1])**3/h-(x-xData[i+1])*h)*k[i]/6.0 \ - ((x-xData[i])**3/h-(x-xData[i])*h)*k[i+1]/6.0 \ + (yData[i]*(x-xData[i+1]) \ - yData[i+1]*(x-xData[i]))/h return y 
       
# Here we calculate the coefficients of each of the mth cubic polynomial (around the point xm). # My code; formulas from math.fullerton.edu def findCoefficients(xData,yData,k,m) : h=xData[m+1]-xData[m] d=(yData[m+1]-yData[m])/h s=np.array([yData[m],d-h*(2*k[m]+k[m+1])/6,k[m]/2,(k[m+1]-k[m])/6/h]) return s 
       

STARTING

# Input the data. Here we are using numpy arrays so that python function "curvatures" doesn't complain xD=np.array([0,1,2], dtype=float) yD=np.array([1.43,1.82,1.66], dtype=float) print xD, yD 
       
[ 0.  1.  2.] [ 1.43  1.82  1.66]
[ 0.  1.  2.] [ 1.43  1.82  1.66]
k= curvatures(xD,yD) print k 
       
[ 0.    -0.825  0.   ]
[ 0.    -0.825  0.   ]
y= evalSpline(xD,yD,k,0.5) print y 
       
1.6765625
1.6765625
s0=findCoefficients(xD,yD,k,0) print s0 s1=findCoefficients(xD,yD,k,1) print s1 
       
[ 1.43    0.5275  0.     -0.1375]
[ 1.82    0.115  -0.4125  0.1375]
[ 1.43    0.5275  0.     -0.1375]
[ 1.82    0.115  -0.4125  0.1375]

So:

from x=0 to x=1, we have: $p_1(x)=1.43+0.5275x-0.1375x^3$  and

from x=1 to x=2, we have: $p_2(x)=1.155+1.3525(x-1)-0.825(x-1)^2+0.1375(x-1)^3$ .

var ( 'x t' ) n=len(xD) P=list_plot(zip(xD,yD), color="red", size=30) C=sum(point((t,evalSpline(xD,yD,k,t))) for t in [xD[0]..xD[n-1], step=0.1]) P1=plot(s0[0]+s0[1]*(x-xD[0])+s0[2]*(x-xD[0])^2+s0[3]*(x-xD[0])^3,(x,xD[0],xD[1])) P2=plot(s1[0]+s1[1]*(x-xD[1])+s1[2]*(x-xD[1])^2+s1[3]*(x-xD[1])^3,(x,xD[1],xD[2])) show(P+C+P1+P2)