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.