NB. Here is a set of spline interpolation routines: NB. x spline y --> D (Find 2nd derivatives needed for "natural" spline fit) NB. x splineA y --> D (2nd derivatives for spline fits for multiple y vectors) NB. (xp0;x;xpN) splineD y --> D (a more general routine for 2nd derivatives) NB. x spline_PGA y --> D (much faster for large x, but MIGHT NOT CONVERGE?) NB. (x;y;D) spltrp x0 --> y(x0) (interpolates at x0 -- needs D) NB. (x;y;D) spltrpA x0 --> y(x0) (interpolates at x0 (needs D) for multiple y's) NB. (x;y;D) spline_fit0 x0 --> y(x0) (interpolates at x0 -- doesn't extrapolate) NB. (x;y;D) spline_der x0 --> y(x0), [dy/dx](x0) (fit & 1st derivatives of fit) NB. (x;y) splint x0 --> y(x0) (combines spline & spltrp) NB. (x;y) splint_PGA x0 --> y(x0) NB. x sinteg y --> integral of y from x0 to xN NB. x sintegA y --> integral of y from x0 to xN for multiple y vectors NB. x spintg y --> integral of y from x0 to x0,x1,x2, ... xN NB. x spintgA y --> integral from x0 to x0,x1,x2, ... xN for multiple y vectors NB. =========================================================================== DEL =: }. - }: NB. the difference between adjacent elements of array NB. Find 2nd derivatives for "natural" spline fit NB. "free endpoints",i.e., assume 2nd derivatives at endpoints = 0 NB. Usage: x spline y --> D spline=: dyad define h2=. }. h=. DEL x a=. 0.5-c=. -: h2* hh=. %h2 + }:h b=. 1#~#c u=. 3*hh*DEL (DEL y)%h 0,(tridgl a;b;c;u),0 ) NB. Find 2nd derivatives for spline fit. NB. "free endpoints", i.e., assume 2nd derivatives at endpoints=0 NB. Usage: x splineA y --> D NB. Here, y may be a (n,#x) array to yeild a (n,#x) array of D values. splineA=: dyad define h2=. }. h=. DEL x a=. 0.5-c=. -: h2* hh=. % h2+ }:h b=. 1#~#c u=. 3*hh*DEL (DEL |:y)%h |: 0,(tridgl a;b;c;u),0 ) NB. Solve tridiagonal system by direct substitution. NB. Rows are a_i*x_(i-1)+b_i*x_i+c_i*x_(i+1) = u_i NB. a,b,c,u have same length, i.e., a[0]=0 and c[n]=0. NB. Usage: soln <-- tridgl a;b;c;u NB. tridgl=: monad define 'a b c u'=. y n1=. <: n=. <: #b x=. %{.b c=. (x*{.c) 0}c u=. (x*{.u) 0}u j=. 0 while. n1>:j=.>:j do. b1=. j{b a1=. j{a c0=. (<:j){c x=. %b1-a1*c0 c1=. j{c c=. (x*c1)j}c u1=. j{u u0=. (<:j){u u=. (x*u1-a1*u0)j}u end. u=. (((n{u)-(n{a)*(n1{u))%(n{b)-(n{a)*(n1{c))n}u j=. n while. 0<:j=.<:j do. u1=. j{u c1=. j{c u2=. (>:j){u u=. (u1-c1*u2)j}u end. ) NB. More general routine for spline 2nd derivatives: NB. Specify 1st derivatives of y at endpoints: NB. Let xp0=[dy/dx] @ x_0, xpN=[dy/dx] @ x_N NB. Usage: (xp0;x;xpN) splineD y --> D NB. If one or both boundaries are free (y''=0), NB. set derivative(s) xp0 and/or xpN = 'f', e.g., NB. ('f';x;0) splineD y --> D for 1st free, last y'=0. splineD=: dyad define 'xp0 x xpN'=. x h2=. }. h=. DEL x a=. 0.5-c=. -: h2* hh=. %h2 + }:h b=. 1#~#a NB. diagonal elements are unity u=. 3*hh*DEL dydx=. h%~DEL y if. xp0 = 'f' do. NB. first free if. xpN='f' do. NB. both free 0,(tridgl a;b;c;u),0 else. NB. or just first free a=. a,0.5 b=. b,1 c=. c,0 u=. u, 3*(xpN-{:dydx)%{:h 0,(tridgl a;b;c;u) end. else. NB. or first derivative set a=. 0,a b=. 1,b c=. 0.5,c u=. (_3*(xp0-{.dydx)%{.h),u if. xpN='f' do. NB. last free (tridgl a;b;c;u),0 else. NB. or both derivatives set a=. a,0.5 b=. b,1 c=. c,0 u=. u, 3*(xpN-{:dydx)%{:h tridgl a;b;c;u end. end. ) NB. =========================================================================== NB. Find 2nd derivatives for spline fit, NB. using iterative approach: NB. Usage: x spline_PGA y --> D spline_PGA=: dyad define h2=. }. h=. DEL x a=. 0.5-c=. -: h2* hh=. %h2 + }:h b=. 1#~#c u=. 3*hh*DEL (DEL y)%h 0,(triPGA a;b;c;u),0 ) NB. ============================================================================ NB. Solve tridiagonal matrix iteratively NB. using the Parallel Gauss Algorithm triPGA=: monad define 'a d c u'=. (1{>y)&(%~) &.>y NB. normalize to unit diagonal v2_pga=: a*(_1|.c) d=. F1_pga^:_ d v2_pga=: a %(_1|.d) v1_pga=: u v1_pga=: u=. d%~ F2_pga^:_ u v2_pga=: c%d F3_pga^:_ u ) NB. -------------------- F1_pga=: 3 : '1 (0)} 1- v2_pga % _1|. y' NB. -------------------- F2_pga=: monad define z=. v1_pga - v2_pga * _1|. y ({. v1_pga) 0} z ) NB. -------------------- F3_pga=: monad define z=. v1_pga - v2_pga * 1|. y ({: v1_pga) (<: #v1_pga)} z ) NB. =============================================================================== NB. Spline interpolation in x,y(x) given 2nd derivs D(x) NB. Usage: (x;y;D) spltrp x0 --> y(x0) NB. x, y & D boxed; x0 may be an array NB. ---------------------------------------------------- spltrp=: dyad define 'x yy D'=. x f0=. (}: x)%(h=. DEL x) f1=. (}. x)%h g0=. (}: yy)%h g1=. (}. yy)%h g=. h*(g0*f1)-g1*f0 ff=. <: x I. y NB. get locations of y point(s) n=. (}.}:x) I. y NB. extraploate if y beyond ends. lin=. (n{g) + y* n{(g1-g0) NB. linear interpolation xn=. y% hn=. n{h a=. (n{f1)-xn b=. xn-(n{f0) aa=. a * <: *:a bb=. b * <: *:b cub=. 6%~(*:hn)*(aa* n{ }:D)+bb* n{ }.D NB. cubic term NB. (n=ff)=0 for points beyond ends -> linear term only lin + (n=ff)*cub ) NB. ============================================================================= NB. Spline interpolation in x,y(x) given 2nd derivs D(x) NB. Usage: (x;y;D) spltrp x0 --> y(x0) NB. x, y & D boxed; x0 may be a vector NB. Also, y & D may be arrays to give array y(x0)s NB. ---------------------------------------------------- spltrpA=: dyad define 'x yy D'=. x f0=. (}: x)%(h=. DEL x) f1=. (}. x)%h g0=. (}: |:yy)%h g1=. (}. |:yy)%h g=. h*(g0*f1)-g1*f0 ff=. <: x I. y NB. get locations of y point(s) n=. (}.}:x) I. y NB. extraploate if y beyond ends. lin=. (n{g) + y* n{(g1-g0) NB. linear interpolation xn=. y% hn=. n{h a=. (n{f1)-xn b=. xn-(n{f0) aa=. a * <: *:a bb=. b * <: *:b cub=. 6%~(*:hn)*(aa* n{ }:|:D)+bb* n{ }.|:D NB. cubic term NB. (n=ff)=0 for points beyond ends -> linear term only |: lin + (n=ff)*cub ) NB. ================================================================================= NB. Spline interpolation: Usage (x;y) splint x0 NB. x0 may be an array NB. -------------------------------------------- splint=: dyad define 'xx yy'=. x D=. xx spline yy (xx;yy;D) spltrp y ) NB. ============================================================================= NB. Iterative spline interpolation: Usage (x;y) splint_PGA x0 NB. x0 may be an array NB. -------------------------------------------- splint_PGA=: dyad define 'xx yy'=. x D=. xx spline_PGA yy (xx;yy;D) spltrp y ) NB. ================================================================================== spline_fit0=: dyad define NB. Don't extrapolate; use last value if beyond ends. NB. Usage: (x;y;D) spline_fit0 x0 --> y(x0) NB. x, y & d boxed; x0 may be an array 'x yy d'=. x h=. DEL x f0=. (}: x)%h f1=. (}. x)%h g0=. (}: yy)%h g1=. (}. yy)%h g3=. g1 - g0 NB. (y_i+1 - y_i)/(x_i+1 - x_i) g4=. ((g0*f1)-(g1*f0))*h ff=. <: x I. y NB. get locations of the point(s) t=. ff=n=. (}.}:x) I. y NB. t=0 beyond ends where n not= ff lin=. (n{g4)+(y*(n{g3)) NB. linear interpolation xn=. y% hn=. n{h a=. (n{f1)-xn b=. xn-(n{f0) aa=. a* <: *:a NB. a*(a^2 - 1)= a^3 - a bb=. b* <: *:b cub=. ((aa*(n{}:d))+(bb*(n{}.d)))*(*:hn)%6 NB. cubic term nn=. n + (n integ NB. of y(x) from x(0) to x(N) NB. -------------------------------------------- sinteg=: dyad define d=. x spline y h=. DEL x g=. +/(h%2)*((}: y)+(}. y)) NB. linear part f=. +/((h**:h)%24)*((}: d)+(}. d)) NB. cubic term g - f ) NB. ================================================================================== NB. Spline integration: Usage: x sinteg y --> integ NB. of y(x) from x(0) to x(N) NB. here, y may be a (n,#x) array of y vectors NB. -------------------------------------------- sintegA=: dyad define d=. |: x splineA y h=. DEL x g=. +/(h%2)*((}: |:y)+(}. |:y)) NB. linear part f=. +/((h**:h)%24)*((}: d)+(}. d)) NB. cubic term g - f ) NB. ================================================================================== NB. Spline integration: Usage x spintg y --> integ NB. of y(x) from x(0) to x(N); also gives the NB. partial integrals from x(0) to x(0),x(1),x(2),... NB. -------------------------------------------- spintg=: dyad define d=. x spline y h=. DEL x g=. (h%2)*((}: y)+(}. y)) NB. linear part f=. ((h**:h)%24)*((}: d)+(}. d)) NB. cubic term fg=. g - f 0,(+/\ fg) ) NB.============================================================================ NB. Spline integration: Usage x spintg y --> integ NB. of y(x) from x(0) to x(N); also gives the NB. partial integrals from x(0) to x(0),x(1),x(2),... NB. here, y may be a (n,#x) array of y vectors NB. -------------------------------------------- spintgA=: dyad define d=. |: x splineA y h=. DEL x g=. (h%2)*((}: |:y)+(}. |:y)) NB. linear part f=. ((h**:h)%24)*((}: d)+(}. d)) NB. cubic term fg=. g - f |: 0,(+/\ fg) ) NB.============================================================================== NB. Get both the spline interpolate and its derivative. NB. Beyond endpoint: extrapolate with endpoint derivative. NB. Usage: (x;y;D) spline_der x0 --> y(x0);[dy/dx](x0) NB. x, y & D boxed; x0 may be an array spline_der=: dyad define 'x yy D'=. x f0=. (}: x)%(h=. DEL x) f1=. (}. x)%h g0=. (}: yy)%h g1=. (}. yy)%h NB. g1-g0=(y_i+1 - y_i)/(x_i+1 - x_i) g=. h*(g0*f1)-g1*f0 ff=. <: x I. y NB. get locations of y point(s) n=. (}.}:x) I. y NB. extraploate if y beyond ends. dy0=. n{ g1-g0 NB. derivative of linear term lin=. (n{g) + y*dy0 NB. linear interpolation xn=. y % hn=. n{h a=. (n{f1)-xn b=. xn-(n{f0) aa=. a * <: *:a NB. a*(a^2 - 1)= a^3 - a bb=. b * <: *:b cub=. 6%~(*:hn)*(aa* n{ }:D)+bb* n{ }.D NB. cubic term NB. correction to linear derivative to give endpoint slope: clow=. -hn* +/(2 1%6)*2{. D NB. -(x_1-x_0)*(D_1/3+D_2/6) chigh=. hn* +/(1 2%6)*_2{. D NB. (x_(N-1)-x_(N-2))*(D_(N-2)/6+D_(N-1)/3) ylow=. ({.yy)+(y - {.x)*(dy0+clow) NB. extrapolate low yhigh=. ({:yy)+(y - {:x)*(dy0+chigh) NB. extrapolate high NB. (n=ff)=0 for points beyond ends -> linear term only NB. (n>ff) for points below x_0; (nff)*ylow)+((nff)*clow)+((n y(x0) NB. x & y boxed, x0 may be a vector, extrapolate off ends. lintrp=: 4 : 0 'xx yy'=. x g0=. (}: yy)% d=. DEL xx g1=. (}. yy)% d u=. (g0* }.xx) - g1* }:xx n=. xx IdX y (n{u)+ y*n{g1-g0 ) NB. x IdX x0 --> index where x0 fits in x-array, e.g.: NB. 0 1 2 3 4 5 6 <- x values NB. | | | | | | | NB. <-- 0 | 1 | 2 | 3 | 4 | 5 ---> IdX values IdX=: ] I.~ [: }: [: }. [