Re: interpolation

*To*: mathgroup at smc.vnet.net*Subject*: [mg47046] Re: interpolation*From*: "Peltio" <peltio at twilight.zone>*Date*: Mon, 22 Mar 2004 05:19:05 -0500 (EST)*References*: <c3bgj0$7ra$1@smc.vnet.net>*Reply-to*: "Peltio" <peltioNOSP at Mdespammed.com.invalid>*Sender*: owner-wri-mathgroup at wolfram.com

"Amanda Hattaway" wrote >Is there any way to display this function f that interpolates my data >set, as a function of an independent variable, rather than just >have Mathematica generate values of f given a picked numerical domain value? If you need the polynomial of order N that interpolates N+1 points, you might want to check out the function InterpolatingPolynomial. InterpolatingPolynomial[data, x] But if' you're talking about interpolation in the real world : ), well, what you get with Interpolation is not *a* single function but rather a collection of functions that make up a piecewise cubic (if the InterpolationOrder is left to its default value of 3) polynomial. I don't know how to access the internally stored information, but I can guess it is stored in the form of a table of coefficients of the constituent polynomials over each subinterval. You could split your interval in suitable subintervals and apply InterpolationgPolynomial to each of them to retrieve the several constituents your function. By properly wrapping the interval-polynomial pairs in a Which statement you'd have _a_ single function. This looks like the easiest way to do it. For say, a 100 point interpolation it'd produce a Which statement with 99 intervals and 99 polynomials. Why do you want to look at that? : ) Perhaps you will be better off in writing your own interpolating procedure. The code for piecewise interpolation is pretty straightforward and can be found in any numerical analysis textbook and it can be pretty fast even if you choose to implement it in procedural style. You need to implement: a data structure to hold the information a procedure to evaluate the interpolating function a procedure to produce the information that make up the interpolating function I have not the code at hand for the piecewise polynomial interpolation, but as an example of what I mean, here's some code I used for toying with simple cubic splines: This is the template for storing the internal data: (* SplineData[{a, b, h}, {c[-1],c[0],..., c[N], c[N + 1]}, type] *) 'a','b' are the endpoints of the interval divided in subintervals of amplitude 'h'. The 'c's are the coefficients of the linear combination of the (translated and dilated copies of the) fundamental spline function BSpline[x_] := Which[ x <= -2, 0, -2 <= x <= -1, (2 + x)^3, -1 <= x <= 0, 1 + 3*(1 + x) + 3*(1 + x)^2 - 3*(1 + x)^3, 0 <= x <= 1, 1 + 3*(1 - x) + 3*(1 - x)^2 - 3*(1 - x)^3, 1 <= x <= 2, (2 - x)^3, x >= 2, 0] 'type' is the head of such function, and is used to identify the type of spline used (BSpline in this case). The following function computes the coefficients of the (complete) spline interpolating the points in the form of a list of coordinates (with equispaced abscissas). It requires the loading of the package Tridiagonal to work. SplineInterpolation[xylist_, yp_] := Module[ {x, y, n, h, ld, d, up, f}, {x, y} = Transpose[xylist]; n = Length[x]; h = x[[2]] - x[[1]]; ld = Table[1, {n - 1}]; ld[[-1]] = 2; ud = Table[1, {n - 1}]; ud[[1]] = 2; d = Table[4, {n}]; f = y; f[[1]] = y[[1]] + (yp[[1]]/3)*h; f[[-1]] = y[[-1]] - (yp[[-1]]/3)*h; sol = TridiagonalSolve[ld, d, ud, f]; SplineData[{x[[1]], x[[-1]], h}, Flatten[{sol[[2]] - (1/3)*h*yp[[1]], sol, sol[[-2]] + (1/3)*h*yp[[-1]]}], BSpline] ] The procedure returns the spline function in our conventional form. The following function, given a value x, determines which elemetary splines should contribute to the function at x (only three of them contribute to the linear combination at any point) splineEval[SplineData[{a_, b_, h_}, coeffs_, sp_]][x_] := Module[ {k, cs}, k = Floor[(x - a)/h] + 1; cs = Take[coeffs, {k, k + 3}]; xs = Table[a + j*h, {j, k - 2, k + 1}]; cs . (sp[(x - #1)/h] & ) /@ xs ] The whole spline function is a piecewise composition of Floor[(b- a)/h] + 1 cubic polynomials obtained by adding three adiacent 'splinelets'. In this case, to get what you want you should write a procedure that returns one such polynomial for every interval. Then feed the stuff into a Which wrapper. I did not feel at the time any urge to do so for splines, and thus I have not the code to do that. I'd do it by rigging a procedure that adopts a simplified version of splineEval to every subinterval, with a symbolic x. Obviously you'd have to do it for a polynomial interpolation. cheers, Peltio Invalid mail in reply-to. Crafty demunging required to mail me.