solve for the spline coefficients y'' If we have scipy support, use a fast tridiagonal algorithm and numpy arrays internally to compute spline coefficients. Without scipy support, use a raw python code to compute coefficients.
Definition at line 897 of file C2Functions.py. 00897 : 00898 if _has_linalg: 00899 "y2 = spline(x_vals,y_vals, yp1=None, ypn=None) returns the y2 table for the spline as needed by splint()" 00900 00901 n=len(x) 00902 u=_numpy.zeros(n,_numpy.float64) 00903 00904 x=_numpy.asarray(x, _numpy.float64) 00905 y=_numpy.asarray(y, _numpy.float64) 00906 00907 dx=x[1:]-x[:-1] 00908 dx2=(x[2:]-x[:-2]) 00909 dy=(y[1:]-y[:-1]) 00910 dydx=dy/dx 00911 00912 # u[i]=(y[i+1]-y[i])/float(x[i+1]-x[i]) - (y[i]-y[i-1])/float(x[i]-x[i-1]) 00913 u[1:-1]=dydx[1:] 00914 u[1:-1]-=dydx[:-1] #this is an incomplete rendering of u... the rest requires recursion in the loop 00915 00916 trimat=_numpy.zeros((3, n), _numpy.float64) 00917 00918 trimat[0, 1:]=dx 00919 trimat[1, 1:-1]=dx2 00920 trimat[2, :-1]=dx 00921 trimat[0] *= (1.0/6.0) 00922 trimat[1] *= (1.0/3.0) 00923 trimat[2] *= (1.0/6.0) 00924 00925 if yp1 is None: 00926 u[0]=0.0 00927 trimat[1,0]=1 00928 trimat[0,1]=0 00929 else: 00930 u[0]=dydx[0]-yp1 00931 trimat[1,0] = dx[0]/(3.0) 00932 00933 if ypn is None: 00934 u[-1]=0.0 00935 trimat[1,-1]=1 00936 trimat[2,-2]=0 00937 else: 00938 u[-1]=ypn-dydx[-1] 00939 trimat[1,-1] = dx[-1]/(3.0) 00940 00941 y2=_linalg.solve_banded((1,1), trimat, u, debug=0) 00942 return _numeric.asarray(y2, _numeric_float) 00943 else: 00944 n=len(x) 00945 u=_numeric.zeros(n,_numeric_float) 00946 y2=_numeric.zeros(n,_numeric_float) 00947 00948 x=_numeric.asarray(x, _numeric_float) 00949 y=_numeric.asarray(y, _numeric_float) 00950 00951 dx=x[1:]-x[:-1] 00952 dxi=1.0/dx 00953 dx2i=1.0/(x[2:]-x[:-2]) 00954 dy=(y[1:]-y[:-1]) 00955 siga=dx[:-1]*dx2i 00956 dydx=dy*dxi 00957 00958 # u[i]=(y[i+1]-y[i])/float(x[i+1]-x[i]) - (y[i]-y[i-1])/float(x[i]-x[i-1]) 00959 u[1:-1]=dydx[1:] 00960 u[1:-1]-=dydx[:-1] #this is an incomplete rendering of u... the rest requires recursion in the loop 00961 00962 if yp1 is None: 00963 y2[0]=u[0]=0.0 00964 else: 00965 y2[0]= -0.5 00966 u[0]=(3.0*dxi[0])*(dy[0]*dxi[0] -yp1) 00967 00968 for i in range(1,n-1): 00969 sig=siga[i-1] 00970 p=sig*y2[i-1]+2.0 00971 y2[i]=(sig-1.0)/p 00972 u[i]=(6.0*u[i]*dx2i[i-1] - sig*u[i-1])/p 00973 00974 if ypn is None: 00975 qn=un=0.0 00976 else: 00977 qn= 0.5 00978 un=(3.0*dxi[-1])*(ypn- dy[-1]*dxi[-1] ) 00979 00980 y2[-1]=(un-qn*u[-2])/(qn*y2[-2]+1.0) 00981 for k in range(n-2,-1,-1): 00982 y2[k]=y2[k]*y2[k+1]+u[k] 00983 00984 return y2 ##
|
1.5.4