1
$\begingroup$

I have the following problem

I have the following data

Tinterspike200 = {3.01026957638`, 5.314505776636686`, 
   10.494223363285943`, 16.585365853657912`};
Tinterspike400 = {2.5609756097561167`, 3.940949935815186`, 
   6.103979460847167`, 8.921694480102463`, 12.50962772785579`, 
   17.092426187419257`, 22.13093709884531`};
Tinterspike600 = {2.3748395378690628`, 3.177150192554557`, 
   4.358151476251605`, 6.059050064184852`, 8.401797175866495`, 
   11.206675224646983`, 14.80744544287548`, 18.58793324775353`, 
   22.310654685494224`, 26.78433889602054`};
Tinterspike800 = {2.2657252888318355`, 2.7856225930680356`, 
   3.4403080872913994`, 4.2105263157894735`, 5.7124518613607185`, 
   8.318356867779203`, 11.59178433889602`, 14.441591784338895`, 
   17.77920410783055`, 21.059050064184852`, 25.532734274711167`};
Tinterspike1000 = {2.1822849807445444`, 2.593068035943517`, 
   3.0680359435173297`, 3.5879332477535297`, 4.255455712451861`, 
   5.423620025673941`, 8.164313222079588`, 11.07188703465982`, 
   13.49165596919127`, 16.084724005134788`, 19.17201540436457`, 
   22.35558408215661`, 25.84724005134788`};
Nspikes200 = {1, 2, 3, 4};
Nspikes400 = {1, 2, 3, 4, 5, 6, 7};
Nspikes600 = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
Nspikes800 = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11};
Nspikes1000 = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13};
Istim = {200, 400, 600, 800, 1000};

I arrange the data in the following way

(*DATA*)
data1 = 
  Table[{Istim[[1]], Tinterspike200[[i]], Nspikes200[[i]]}, {i, 1, 4}];
data2 = Table[{Istim[[2]], Tinterspike400[[i]], Nspikes400[[i]]}, {i, 
    1, 7}];
data3 = Table[{Istim[[3]], Tinterspike600[[i]], Nspikes600[[i]]}, {i, 
    1, 10}];
data4 = Table[{Istim[[4]], Tinterspike800[[i]], Nspikes800[[i]]}, {i, 
    1, 11}];
data5 = Table[{Istim[[5]], Tinterspike1000[[i]], 
    Nspikes1000[[i]]}, {i, 1, 13}];
data = Join[data1, data2, data3, data4, data5];
lpp = ListPointPlot3D[data, PlotStyle -> {PointSize[Large], Red}];

I determine the interpolating function

(*INTERPOLATION*)
{xmin, xmax} = MinMax[data[[All, 1]]];
{ymin, ymax} = MinMax[data[[All, 2]]];
dataInterp = {Most@#, Last@#} & /@ data;
Istim3D = Interpolation[dataInterp, InterpolationOrder -> 1]
plIstim3D = 
  Plot3D[Istim3D[x, y], {x, xmin, xmax}, {y, ymin, ymax}, 
   PlotStyle -> Opacity[0.8], 
   AxesLabel -> {"Istim", 
     "Tempi interspikes", "Numero di spikes"}, PlotRange -> All, 
   ImageSize -> 800];
Show[lpp, plIstim3D, ImageSize -> 800]

I would like to determine the analytical equation of the interpolating function, Istim3D, i.e., I would like to have f(x,y,z), where f is the function representing the interpolating function. Thanks for your help.

Improve this question
$\endgroup$
2
  • $\begingroup$ Do you mean that you want a function that fits the interpolating function? The form of this function should be inferred from the context in which the data was given. If you actually want what the InterpolatingFunction gives you under the hood, why? It's just a bunch of polynomials pasted together. Can you clarify what you really want here? $\endgroup$
    – march
    Commented Nov 6, 2020 at 17:49
  • $\begingroup$ Thank you very much for your reply. Yes, I want a function that fits the interpolating function. I would like to determine the analytical form of the interpolating surface. I need the analytical equation, even though it is a bunch of polynomials pasted together, because this function is useful for me for further analysis. I do not know whether I was clear. $\endgroup$
    – VDF
    Commented Nov 6, 2020 at 18:16

2 Answers 2

Reset to default
2
$\begingroup$

You can have a coarse approximation interpolating a polynomial as follows

n = 3;
m = 3;
Clear[f]
f[x_, y_] := Total[Flatten[Table[Subscript[a, i, j] x^i y^j, {i, 0, n}, {j, 0, m}]]];
vars = Flatten[Table[Subscript[a, i, j] , {i, 0, n}, {j, 0, m}]];
obj = Sum[(f[data[[k]][[1]], data[[k]][[2]]] - data[[k]][[3]])^2, {k, 1, Length[data]}];
sol = NMinimize[obj, vars]
fxy = f[x, y] /. sol[[2]];
gr1 = Plot3D[fxy, {x, 200, 1000}, {y, 0, 30}];
Show[lpp, gr1]

enter image description here

NOTE

Follows a script extracted from a Wolfram repository with minor changes to implement the spline issue. The data was slightly modified to cope with the needed regular mesh data.

(*AUXILIARY FUNCTIONS*)

searchSpan[knots_, u0_] :=
  With[{max = Max[knots]},
    If[u0 == max, Position[knots, max][[1, 1]] - 2,
      Ordering[UnitStep[u0 - knots], 1][[1]] - 2]
  ]
  calcInterpCtrlPoints[{deg_, knots_}, paras_] :=
  Module[{dim, coeffMat}, dim = Length@paras - 1;
    coeffMat = Function[{u0},
      With[{i = searchSpan[knots, u0]},
        Join[ConstantArray[0, i - deg],
          BSplineBasis[{deg, knots}, #, u0] & /@ Range[i - deg, i],
          ConstantArray[0, dim - i]]
      ]
  ] /@ paras;
 (*generating the LinearSolveFunction[] object*)

LinearSolve[coeffMat]
]
calcParas[pts_, type_] := Which[
  type === Automatic || type === "ChordLength",
  FoldList[Plus, 0, Normalize[(Norm /@ Differences[pts]), Total]] // N,
  type === "Centripetal",
  FoldList[Plus, 0, 
  Normalize[(Norm /@ Differences[pts])^(1/2), Total]] // N,
  type === "EqualSpaced", Range[0, 1, 1/(Length@pts - 1)] // N
]
(*Implementations*)
Options[BSplineSurfaceInterpolate] =
  {Parametrization -> Automatic, InterpolationOrder -> 3,
  OriginPoints -> False, ControlNets -> False};

BSplineSurfaceInterpolate[
  pts : {{{_, _, _} ..} ..},
  opts : OptionsPattern[{BSplineSurfaceInterpolate, Graphics3D}]] /;

ArrayQ[pts, 3, NumericQ] := Module[
{m, n, pz, io, cn, op, deg1, deg2, paras1, paras2,
calcKnots, knots1, knots2, temp, ctrlnets},
{m, n} = Dimensions[pts, 2] - 1;
(*achieve the values of options*)

pz = OptionValue[Parametrization];
io = OptionValue[InterpolationOrder] /. Automatic -> 3;
cn = OptionValue[ControlNets];
op = OptionValue[OriginPoints];
If[Head[io] === List,
{deg1, deg2} = io, deg1 = deg2 = io];
(*reset the value of deg to m if the value of deg1 greater than m*)

If[deg1 > m, deg1 = m]; If[deg2 > n, deg2 = n];
(*calculate the parameters in two directions*)

paras1 = Mean[calcParas[#, pz] & /@ Transpose[pts]];
paras2 = Mean[calcParas[#, pz] & /@ pts];
(*calculate the knots vectors in two diredtions*)
calcKnots =
 Function[{dim, deg, paras},
  Join[ConstantArray[0, deg + 1],
    1/deg (Plus @@ (paras[[# + 1 ;; # + deg]]) & /@ 
       Range[1, dim - deg]),
    ConstantArray[1, deg + 1]] // N
 ];
knots1 = calcKnots[m, deg1, paras1];
knots2 = calcKnots[n, deg2, paras2];
(*compute the control nets of the NURBS surface*)
temp = calcInterpCtrlPoints[{deg1, knots1}, paras1] /@ Transpose[pts];
ctrlnets = calcInterpCtrlPoints[{deg2, knots2}, paras2] /@ Transpose[temp];
(*visualize the ultimate interpolation NURBS surface*)

Graphics3D[
{BSplineSurface[ctrlnets,
 SplineDegree -> {deg1, deg2}, SplineKnots -> {knots1, knots2}],
Sequence @@
 If[cn,
  {Dashed, Green, Line[ctrlnets], Line[Transpose[ctrlnets]],
   Red, PointSize[Medium], Point[#] & /@ ctrlnets}, {}],
Sequence @@
 If[op,
  {PointSize[Medium], Blue, Map[Point, pts], Gray}, {}]},
 Evaluate@FilterRules[{opts}, Options[Graphics3D]]
 ]
]
(*Initialization the interpolation points*)

pts1 = {{{20, 3.01026957638`, 1}, {20, 5.314505776636686`, 2}, {20, 
 10.494223363285943`, 3}, {20, 16.585365853657912`, 4}, {20, 
 22.13093709884531`, 5}, {20, 26.78433889602054`, 6}, {20, 30, 
 7}}, {{40, 2.5609756097561167`, 1}, {40, 3.940949935815186`, 
 2}, {40, 6.103979460847167`, 3}, {40, 8.921694480102463`, 
 4}, {40, 12.50962772785579`, 5}, {40, 17.092426187419257`, 
 6}, {40, 22.13093709884531`, 7}}, {{60, 2.3748395378690628`, 
 1}, {60, 3.177150192554557`, 2}, {60, 4.358151476251605`, 
 3}, {60, 6.059050064184852`, 4}, {60, 8.401797175866495`, 
 5}, {60, 11.206675224646983`, 6}, {60, 14.80744544287548`, 
 7}}, {{80, 2.2657252888318355`, 1}, {80, 2.7856225930680356`, 
 2}, {80, 3.4403080872913994`, 3}, {80, 4.2105263157894735`, 
 4}, {80, 5.7124518613607185`, 5}, {80, 8.318356867779203`, 
 6}, {80, 11.59178433889602`, 7}}, {{100, 2.1822849807445444`, 
 1}, {100, 2.593068035943517`, 2}, {100, 3.0680359435173297`, 
 3}, {100, 3.5879332477535297`, 4}, {100, 4.255455712451861`, 
 5}, {100, 5.423620025673941`, 6}, {100, 8.164313222079588`, 7}}};

Manipulate[
 BSplineSurfaceInterpolate[pts, OriginPoints -> op, ControlNets -> cn,
   Parametrization -> pz, InterpolationOrder -> io,
   PlotRange -> {{20, 100}, {0, 30}, {0, 8}}],
   "data samples",
 {{pts, pts1, ""}, {pts1 -> "sample 1"}},
 Delimiter,
 "interpolation order",
 {{io, 3, ""}, 1, 5, 1, Appearance -> "Labeled", ImageSize -> Tiny},
 Row[{Control[{{op, True, "origin points"}, {True, False}}]}],
 Row[{Control[{{cn, False, "control nets"}, {True, False}}]}],
 "",
 "parametrization",
 {{pz, "ChordLength", ""}, {"ChordLength", "Centripetal", 
 "EqualSpaced"}, ControlType -> PopupMenu},
 SaveDefinitions -> True, ControlPlacement -> Left

]

enter image description here

Improve this answer
$\endgroup$
7
  • $\begingroup$ thank you very much for your idea, that I think it is very useful! I could for example estimate the difference between the guessed polynomial function and the interpolating function to see which is the best approximation. I think your idea is very useful. Thank you very much! :-) $\endgroup$
    – VDF
    Commented Nov 6, 2020 at 22:48
  • $\begingroup$ I have tried to increase the order of the interpolating surface polynomial, but some oscillations appear. Do you have an idea on how can I remove them? $\endgroup$
    – VDF
    Commented Nov 6, 2020 at 23:19
  • $\begingroup$ Choosing another kind of surface to fit. Perhaps a b-spline. $\endgroup$
    – Cesareo
    Commented Nov 6, 2020 at 23:41
  • $\begingroup$ Thank you very much @Cesareo, your suggestions are very useful. $\endgroup$
    – VDF
    Commented Nov 7, 2020 at 8:56
  • 1
    $\begingroup$ Follows a script from a Wolfram repository with minor changes to implement the spline issue. $\endgroup$
    – Cesareo
    Commented Nov 7, 2020 at 11:14
1
$\begingroup$

Try this:

Istim3D["Domain"]
{{200., 1000.}, {2.18228, 26.7843}}

Istim3D["Grid"]
{{200., 3.01027}, {200., 5.31451}, {200., 10.4942}, {200., 
  16.5854}, {400., 2.56098}, {400., 3.94095}, {400., 6.10398}, {400., 
  8.92169}, {400., 12.5096}, {400., 17.0924}, {400., 22.1309}, {600., 
  2.37484}, {600., 3.17715}, {600., 4.35815}, {600., 6.05905}, {600., 
  8.4018}, {600., 11.2067}, {600., 14.8074}, {600., 18.5879}, {600., 
  22.3107}, {600., 26.7843}, {800., 2.26573}, {800., 2.78562}, {800., 
  3.44031}, {800., 4.21053}, {800., 5.71245}, {800., 8.31836}, {800., 
  11.5918}, {800., 14.4416}, {800., 17.7792}, {800., 21.0591}, {800., 
  25.5327}, {1000., 2.18228}, {1000., 2.59307}, {1000., 
  3.06804}, {1000., 3.58793}, {1000., 4.25546}, {1000., 
  5.42362}, {1000., 8.16431}, {1000., 11.0719}, {1000., 
  13.4917}, {1000., 16.0847}, {1000., 19.172}, {1000., 
  22.3556}, {1000., 25.8472}}

Istim3D["ValuesOnGrid"]
{1., 2., 3., 4., 1., 2., 3., 4., 5., 6., 7., 1., 2., 3., 4., 5., 6., \
7., 8., 9., 10., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 1., \
2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13.}

The "Domain" gives the Min/Max of x/y coordinates.

The "Grid" gives the x/y values of the grid points and the "ValuesOnGrid" the belonging z values.

Now I need to guess a bit. I think, because you specified "InterpolationOrder->1", MMA triangulates the given x/y points. Then, given some x,y values, MMA determines the triangle which contains the point and determines the height of the triangle above this point.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you very much @Daniel Huber it is also a good strategy. $\endgroup$
    – VDF
    Commented Nov 6, 2020 at 22:46
  • $\begingroup$ @Cesareo thank you very much for what you have showed. As you have highlighted the points you have chosen are not the one I have. I am sorry, but I am not able to modify what you have proposed me by applying it to my points. Could you please show me how to do it? $\endgroup$
    – VDF
    Commented Nov 8, 2020 at 20:58

Not the answer you're looking for? Browse other questions tagged or ask your own question.