This question is a generalization of the previous one for multiple dimensions. In the answer to that question an implementation for the clamped spline interpolation for 1D case and arbitrary degree of spline is given. How can it be extended for multidimensional case?
1 Answer
I still haven't figured out how to write a routine for arbitrary dimension, but I'm posting my (incomplete!) solution in case people might have ideas on extending what I have.
Bivariate interpolant
Here is a random bivariate polynomial, which we'll use for generating test data:
f[x_, y_] := -2 + 4 x^2 + 4 x^3 - 3 y - 5 x^2 y + 5 y^2 + 5 x y^2 + y^3
Here's some test data from f[x, y]
, sampled at non-equispaced points:
da = Flatten[Table[{{x, y}, f[x, y]},
{x, {-2, -4/3, -1/5, 1/9, 3/4, 1}},
{y, {0, 1/6, 3/8, 9/5, 2}}], 1];
Here's the reference interpolating function:
{p, q} = {2, 3}; (* spline degrees in the two variables *)
ipf = Interpolation[da, InterpolationOrder -> {p, q}, Method -> "Spline"];
Some preliminary processing to separate out independent and dependent variables:
{pts, vals} = Transpose[SplitBy[SortBy[da, First], #[[1, 1]] &], {3, 2, 1}];
Make the knot sequence for each independent variable:
makeKnots[list_?VectorQ, deg_Integer?NonNegative] :=
With[{n = Length[list]},
Join[ConstantArray[list[[1]], deg + 1],
If[deg + 2 <= n, MovingAverage[ArrayPad[list, -1], deg], {}],
ConstantArray[list[[-1]], deg + 1]]]
{u, v} = {pts[[1, All, 1]], pts[[All, 1, 2]]};
{uk, vk} = MapThread[makeKnots, {{u, v}, {p, q}}];
Build the control points:
{m, n} = {Length[u], Length[v]};
usol = LinearSolve[Outer[BSplineBasis[{p, uk}, #2, #1] &,
u, Range[0, m - 1], 1]];
vsol = LinearSolve[Outer[BSplineBasis[{q, vk}, #2, #1] &,
v, Range[0, n - 1], 1]];
cpts = vsol /@ Transpose[usol /@ vals];
Finally, the bivariate interpolating spline:
spf[x_, y_] = Fold[Dot, cpts, {Table[BSplineBasis[{q, vk}, k - 1, y], {k, n}],
Table[BSplineBasis[{p, uk}, k - 1, x], {k, m}]}];
Plot the two interpolants and the original function:
MapThread[Plot3D[#1[x, y], {x, -2, 1}, {y, 0, 2}, PlotLabel -> #2] &,
{{f, ipf, spf}, {"True", "InterpolatingFunction", "B-spline"}}]
// GraphicsRow
Evaluate ipf[]
and spf[]
at the same argument:
{ipf[-1, 1], spf[-1, 1]}
{-8.59478, -552429212/64275003}
Note that only the second function gave exact output.
The difference between ipf[]
and spf[]
, showing good agreement:
Plot3D[ipf[x, y] - spf[x, y], {x, -2, 1}, {y, 0, 2}, PlotRange -> All]
If you want further confirmation, you can try recovering the control points and knots from ipf[]
, using a procedure similar to the one in my previous answer, and compare them with the control points and knots I generated here.
Trivariate interpolant
Hopefully, you can see the similarities and differences between the previous example and this one:
(* random polynomial *)
mp[x_, y_, z_] := 3 - 7 x^2 + 2 x^3 - 2 y + 8 y^2 + 5 x y^2 + 8 y^3 + 4 x y^3 -
2 y^4 + 6 x z + 2 x^2 z - 6 x^3 z + 2 y z + 4 x y z -
2 x^2 y z - 3 y^2 z - 8 x y^2 z + 7 y^3 z - 5 z^2 + x z^2 +
2 x^2 z^2 + 6 y z^2 - 4 y^2 z^2 + 9 z^3 - 4 x z^3 -
3 y z^3 - 9 z^4
(* random data *)
da = Flatten[Table[{{x, y, z}, mp[x, y, z]}, {x, {-2, -3/2, -9/7, 7/9, 2, 3}},
{y, {1, 5/3, 5/2, 9/2, 5}}, {z, {-1, 5/7, 11/10, 2}}], 2];
{p, q, r} = {4, 3, 2}; (* B-spline degree for each variable *)
ipf = Interpolation[da, InterpolationOrder -> {p, q, r}, Method -> "Spline"];
{pts, vals} = Transpose[GatherBy[da, {#[[1, 1]] &, #[[1, 2]] &}], {4, 3, 2, 1}];
(* make knots *)
{u, v, w} = {pts[[1, 1, All, 1]], pts[[1, All, 1, 2]], pts[[All, 1, 1, 3]]};
{uk, vk, wk} = MapThread[makeKnots, {{u, v, w}, {p, q, r}}];
(* make control points *)
{l, m, n} = Length /@ {u, v, w};
usol = LinearSolve[Outer[BSplineBasis[{p, uk}, #2, #1] &, u, Range[0, l - 1], 1]];
vsol = LinearSolve[Outer[BSplineBasis[{q, vk}, #2, #1] &, v, Range[0, m - 1], 1]];
wsol = LinearSolve[Outer[BSplineBasis[{r, wk}, #2, #1] &, w, Range[0, n - 1], 1]];
cpts = Map[wsol, Transpose[Map[vsol,
Transpose[Map[usol, vals, {2}], {2, 3, 1}], {2}], {1, 3, 2}], {2}];
(* B-spline interpolant *)
spf[x_, y_, z_] = Fold[Dot, cpts, {Table[BSplineBasis[{r, wk}, k - 1, z], {k, n}],
Table[BSplineBasis[{q, vk}, k - 1, y], {k, m}],
Table[BSplineBasis[{p, uk}, k - 1, x], {k, l}]}];
Tests:
{ipf[2, 3, 1], spf[2, 3, 1]}
{383.531, 10447533501473/27240371550}
ContourPlot3D[#[x, y, z], {x, -2, 3}, {y, 1, 5}, {z, -1, 2},
BoxRatios -> Automatic, MaxRecursion -> 0] & /@
{ipf, spf} // GraphicsRow
-
$\begingroup$ The program in Trivariate interpolant was very time-consuming. $\endgroup$– xyzOct 19, 2015 at 7:40
-
$\begingroup$ Agreed. For some reason, it seemed faster in version 8 than in version 10… $\endgroup$ Oct 19, 2015 at 7:51
-
$\begingroup$ Although the precision of
usol
isInfinity
, when you useContourPlot3D
, the precison will change to theMachinePrecison
. In addition, you use the property of B-spline basis. Namely, In generally, when $u_0 \in \left[u _i ,u _{i+1} \right)$, the B-spline function basis $N_{i-p,p},\cdots ,N_{i,p}$ are not equal to 0. $\endgroup$– xyzOct 20, 2015 at 2:07 -
$\begingroup$ That's correct. The only reason for showing the plot is to visually demonstrate that the explicitly-constructed B-spline is equivalent to the output of
Interpolation[]
. If you want to force the use of higher precision, the plotting functions all accept aWorkingPrecision
setting. The explicitly constructed B-spline will keep up; theInterpolatingFunction[]
will still be stuck atMachinePrecision
. $\endgroup$ Oct 20, 2015 at 2:19 -
1$\begingroup$ @Shutao, to repeat what I said: "look at the low order cases, try to determine a pattern"; sometimes, one just starts by guessing. :) As a simpler example: if you keep differentiating $x^k$ $n$ times for various small values of $k$ and $n$, you'll eventually notice the pattern $\frac{k!}{(k-n)!}x^{k-n}$, which you then proceed to prove inductively. A lot of work in combinatorics actually proceeds in this way. $\endgroup$ Oct 20, 2015 at 3:25
arbitrary precision
mean? Could you give the difinition or a demo to showarbitrary precision interpolation
? Or I would like to know which case did you need toarbitrary precision interpolation
? $\endgroup$MachinePrecsion
is enough in my work. $\endgroup$