How to derive/interpolate a polynomial f(x,y) from several discrete polynomials f(x)?

I have 7 discrete degree 3 polynomials $f_i(x)$. Each of them represents a snapshot of data from an experiment and fits this data quite well.

This graph depicts the 5 largest of the 7 polynomials:

The seemingly simple structural connection between the snapshots is obvious for the observer.

What I actually need is to build a polynomial $f(x,y)$ that performs a reasonable interpolation between the discrete steps $i$. Since the existing polynomials $f_i(x)$ fit the experimental data so well, I want $f(x,y=i)$ to resemble $f_i(x)$ as much as possible, especially the 3 largest ones (blue, yellow and green in the image).

However, all my trials with Fit and FindFit with 2 indeterminates and my experimental data produced quite convoluted surfaces instead of the simple “converging tunnel” that the image suggests, and deviated a lot from the discrete polynomials $f_i(x)$.

Any idea of how to achieve what I need?

I’m new to Mathematica, so please forgive me if I’m overlooking the obvious. TIA for your help!

Edit: Here are the 7 polynomials for the $x$ values:

f1[x_]:=-1.5652680002977166+0.00008966734044623339 x+1.2238004079440088*10^-11 x^2-1.0023031996777336*10^-16 x^3
f2[x_]:=-0.709251613865797+0.00007808877386592759 x-4.318708627177192*10^-11 x^2-2.177549100092401*10^-16 x^3
f3[x_]:=-0.10969451811604027+0.00007655917624397411 x-4.888912463799568*10^-11 x^2-7.008331859901413*10^-15 x^3
f4[x_]:=0.39887872110014777 +0.000017214364803508095 x-1.6581991806280448*10^-10 x^2-6.44037732894833*10^-15 x^3
f5[x_]:=0.2621776923966859 +9.278887731970977*10^-7 x-9.73805099727172*10^-10 x^2-1.3177607015593936*10^-13 x^3
f6[x_]:=0.06690167343301313 +0.000027421328844440038 x-1.4575700564412294*10^-8 x^2-3.974730171304088*10^-11 x^3
f7[x_]:=0.1099576690761461 +0.00046442717100703064 x-2.153183483513223*10^-6 x^2-6.236855134171703*10^-8 x^3


The corresponding $y$ values are:

f1: 1*10^6
f2: 5*10^5
f3: 1*10^5
f4: 5*10^4
f5: 1*10^4
f6: 1*10^3
f7: 1*10^2


In the image above, I have scaled the $y$ values to $7-log_{10}(y)$.

Result values $z < 0$ are to be ignored.

Edit 2: Here is the raw data that I used to construct the polynomials (please note that $y$ / $i$ is the first parameter in each triplet, i.e. the triplets with an identical 1st parameter constitute the data (= 2nd and 3rd parameter) for 1 discrete polynomial $f_i(x)$):

{{1000000, 1000000, 0.10545}, {1000000, 999000, 0.18}, {1000000,
997000, 0.29}, {1000000, 995000, 0.58}, {1000000, 992000,
0.83}, {1000000, 991000, 0.93}, {1000000, 990100, 1}, {1000000,
990000, 7.5}, {1000000, 900000, 12.56}, {1000000, 800000,
26}, {1000000, 700000, 35}, {1000000, 600000, 36}, {1000000, 500000,
32.8}, {1000000, 400000, 30}, {1000000, 200000, 15.55}, {1000000,
100000, 6.79}, {1000000, 50000, 3.75}, {500000, 500000,
0.319}, {500000, 100000, 6.45}, {500000, 50000, 3.06}, {100000,
100000, 0.049}, {100000, 50000, 2.72}, {100000, 10000,
0.644}, {50000, 50000, 0.04}, {50000, 10000, 0.548}, {50000, 5000,
0.48}, {10000, 10000, 0.0423}, {10000, 5000, 0.226}, {10000, 1000,
0.262}, {1000, 1000, 0.04}, {1000, 500, 0.072}, {100, 100,
0.0725}, {100, 50, 0.12}}


There is much more empirical data, but applying Fit to the data above produced optimal discrete polynomials $f_i(x)$ for the empirical data as a whole.

Edit 3: I used logarithmic scaling for the image I posted as I thought that would be easier for the eyes, but in fact it seems to confusingly suggest a logarithmic structure which isn’t actually there. So here are two additional graphics without a log scale. One from above:

This makes it obvious that the “width” of the polynomials is in fact exactly linear; this is simply because in the experiment, $x$ cannot be larger than $y$, and $z$ is always close to $0$ if $x ≈ y$.

The second image is from the side (again, without log scale) and suggests that the connection between the maxima of the discrete functions is also more linear than logarithmic:

Edit 4: A reply suggested a simple linear interpolation, but this does not work, even if the situation were perfectly linear.

Let’s look at an idealized variant of my experimental data: only degree 2 polynomials and a perfectly linear connection between the data samples (as in my experiment, the “width” of the parabolas equals $y$):

testdata1={{0,0},{5,10},{10,0}};
test1F=Fit[testdata1,{1,x,x^2},x]
testdata2={{0,0},{10,20},{20,0}};
test2F=Fit[testdata2,{1,x,x^2},x]
testdata3={{0,0},{20,40},{40,0}};
test3F=Fit[testdata3,{1,x,x^2},x]
ParametricPlot3D[{{test1F,10,x},{test2F,20,x},{test3F,40,x}},{x,0,40},PlotRange->{{0,40},{0,40},{0,40}}]


This gives you the following result:

Out[1]= 4.0x-0.4x^2
Out[2]= 4.0x-0.2x^2
Out[3]= 4.0x-0.1x^2


Now, let’s assume we don’t have testdata2 and test2F and, knowing the connection is perfectly linear, build a linear interpolation between test1F and test3F, as was suggested in one reply:

test[x_, y_] = (40 - y)/30*test1F + (y - 10)/30*test3F


Let’s verify with the known data:

test[5,10]
test[20,40]
Out[4]= 10.0
Out[5]= 40.0


Works fine.

No let’s try testdata2 with $x = 10$, which we know must result in $z=20$:

test[10,20]
Out[6]= 10.0


Fail.

So even an idealized, perfectly linear situation would not work with the suggested form of linear interpolation.

Let alone my not perfectly linear situation. This would be almost impossible to solve “manually” in a reasonable amount of time, so when I asked this question, I was sure Mathematica would offer an algorithmic solution for this, and I just couldn’t find it.

• Why not to add those 7 polynomials to your description? It will allow people to address specific problem. – BlacKow Sep 24 '16 at 16:54
• Welcome! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – user9660 Sep 24 '16 at 17:33
• Would you leave available the data? – LCarvalho Sep 24 '16 at 20:38
• I didn’t list the actual polynomials because I wanted to avoid to be too verbose, and thought my question was generic enough so that the specific polynomials would not matter. I have added them now. – Snoopy Sep 24 '16 at 22:02
• I get this for Fit[] with the raw data. Haven't compared it with your f's, though. – Michael E2 Sep 25 '16 at 16:48

Update: This replaces a generic interpolation of the f's -- see edit history.

data = {{1000000, 1000000, 0.10545}, {1000000, 999000,
0.18}, {1000000, 997000, 0.29}, {1000000, 995000, 0.58}, {1000000,
992000, 0.83}, {1000000, 991000, 0.93}, {1000000, 990100,
1}, {1000000, 990000, 7.5}, {1000000, 900000, 12.56}, {1000000,
800000, 26}, {1000000, 700000, 35}, {1000000, 600000,
36}, {1000000, 500000, 32.8}, {1000000, 400000, 30}, {1000000,
200000, 15.55}, {1000000, 100000, 6.79}, {1000000, 50000,
3.75}, {500000, 500000, 0.319}, {500000, 100000, 6.45}, {500000,
50000, 3.06}, {100000, 100000, 0.049}, {100000, 50000,
2.72}, {100000, 10000, 0.644}, {50000, 50000, 0.04}, {50000,
10000, 0.548}, {50000, 5000, 0.48}, {10000, 10000,
0.0423}, {10000, 5000, 0.226}, {10000, 1000, 0.262}, {1000, 1000,
0.04}, {1000, 500, 0.072}, {100, 100, 0.0725}, {100, 50, 0.12}};
ff = Fit[data, {1, x, y, x^2, x y, y^2, x^3, x^2 y, x y^2, y^3}, {y, x}]
(*
0.0360751 - 9.24559*10^-6 x - 2.92375*10^-10 x^2 -
1.00408*10^-16 x^3 + 0.0000149795 y + 3.34707*10^-10 x y +
3.04925*10^-16 x^2 y - 5.74977*10^-11 y^2 - 2.35943*10^-16 x y^2 +
4.09308*10^-17 y^3
*)

Show[
Plot3D[ff,
{x, -1 - Min@data[[All, 2]], 1 + Max@data[[All, 2]]},
{y, -1 - Min@data[[All, 1]], 1 + Max@data[[All, 1]]},
Mesh -> {0, Union@data[[All, 1]]}, MeshStyle -> Thick],
Graphics3D[{Red, PointSize[Medium], Point[data[[All, {2, 1, 3}]]]}],
PlotRange -> MinMax@data[[All, 3]]
]


The black mesh lines represent the OP's functions f1 etc., but since the lines are the result of a multivariate fitting, they will differ somewhat from them.

• Thanks for your generic solution, but unfortunately, if it doesn’t work so well for plotting, it won’t work so well numerically, either. I mean, it’s very obvious that the 3D result should be a “tunnel”. I applied your solution to the simple test case outlined in Edit 4 of my question above, and indeed, it yields the same incorrect results as @bills suggestion: If I use all 3 test functions, then of course the result is correct for all 3 test points. But again, if I omit test2F and build your suggestion with only test1Fand test3F, I’ll get 10 for if[10, 20] instead of 20. – Snoopy Sep 25 '16 at 16:34
• @Snoopy Yes, that's what I meant about having the actual functions/data/code.... – Michael E2 Sep 25 '16 at 16:40
• But the construction principle of my idealized example is so simple and obvious, an interpolation should be capable of figuring it out, no? Arguably, if I only specify test1F and test3F, it’s too little information for Interpolation to figure it out. But if I specify all 3 test functions, it should be perfectly clear. But it still isn’t: In this case, if I try if[15,30], I’ll get 37.5 instead of 30. – Snoopy Sep 25 '16 at 16:55
• Thanks a lot for your update! This finally looks like a starting point I can build upon, with a “tunnel” as it should be. I tried similar things, but with other parameters for Fit, and had decisively worse results. Maybe this is already good enough for my data, and if it isn’t, I figure I could calculate more data points from my discrete functions that I know work well and then use these additional data points with Fit. Thanks again! – Snoopy Sep 25 '16 at 20:41

I already posted an older version of this answer before, but that contained a stupid error on my part, so I removed it again. This new version is hopefully correct.

OK, this might be 50% answer and 50% question, but I think I should sum up the results I’ve achieved so far without further editing the original question.

Idealized Scenario

I still feel my issue is generic at its core, so the best way to start is with the idealized version of my empirical data which avoids the “empirical noise”.

In a more elaborate form than described in Edit 4 of my question, this idealized version has 7 discrete functions that describe the “data”, just as my actual empirical version does. Each function is simply a symmetric parabola defined by a degree 2 polynomial $f_i(x)$. The width of each parabola (i.e. the distance between the two points where it crosses $0$), its height (i.e. its maximum value) and its position on the $y$ axis are all identical.

So with this very simple, very “linear” setup, we have the following data and corresponding functions ($f_1$ has the position $y=10$, $f_2$ has the position $y=20$ etc.):

data1 = {{0, 0}, {5, 10}, {10, 0}};
data2 = {{0, 0}, {10, 20}, {20, 0}};
data3 = {{0, 0}, {15, 30}, {30, 0}};
data4 = {{0, 0}, {20, 40}, {40, 0}};
data5 = {{0, 0}, {25, 50}, {50, 0}};
data6 = {{0, 0}, {30, 60}, {60, 0}};
data7 = {{0, 0}, {35, 70}, {70, 0}};
f1 = Fit[data1, {1, x, x^2}, x]
f2 = Fit[data2, {1, x, x^2}, x]
f3 = Fit[data3, {1, x, x^2}, x]
f4 = Fit[data4, {1, x, x^2}, x]
f5 = Fit[data5, {1, x, x^2}, x]
f6 = Fit[data6, {1, x, x^2}, x]
f7 = Fit[data7, {1, x, x^2}, x]


Since Fit offers approximations, we’ll get some convoluted floating point numbers in $f_1$ – $f_7$, but eliminating the approximation errors, the polynomials are actually quite simple and can be rewritten as follows:

ClearAll[f1,f2,f3,f4,f5,f6,f7]
f1[x_]:=4x-2/5  x^2
f2[x_]:=4x-2/10 x^2
f3[x_]:=4x-2/15 x^2
f4[x_]:=4x-2/20 x^2
f5[x_]:=4x-2/25 x^2
f6[x_]:=4x-2/30 x^2
f7[x_]:=4x-2/35 x^2
ParametricPlot3D[{{f1[x],10,x},{f2[x],20,x},{f3[x],30,x},{f4[x],40,x},{f5[x],50,x},{f6[x],60,x},{f7[x],70,x}},{x,0,80},PlotRange->{{0,80},{0,80},{0,80}}]


This is how this looks like:

The “tunnel” shape of the data is easy to recognize.

Now, if you look closer at $f_1(x)$ – $f_7(x)$, you’ll realize that the denominator of the coefficient of $x^2$ is always $y/2$. So we can derive the following formula for $f(x,y)$:

f[x_,y_]:=4(x-x^2/y)


Ideally, this would be the kind of solution I was looking for (for the idealized version of my experiment).

A simple 3D plot might look disillusioning at first, as it shows an image that I encountered very often during my experimentation, and which does not resemble a tunnel at all:

Plot3D[f[x,y],{x,10,70},{y,10,70}]


However, this is only because Mathematica uses unsuitable default values for the $z$ axis. With suitable values, we’ll get a perfect tunnel:

Plot3D[f[x,y],{x,0,70},{y,0,70},PlotRange->{{0,80},{0,80},{0,80}},BoxRatios->1]


Again, ideally this is the kind of result I was looking for. Note that the connections between the tops and of both intersections with $z=0$ are all straight lines, as it should be. (Also note that Mathematica inverts the $x$ axis to 80 → 0; I have no idea why nor how to fix that.)

Now, how do the suggested solutions perform in my idealized scenario?

Interpolation

First, let’s try Interpolation:

polys=Through[{f1,f2,f3,f4,f5,f6,f7}[x]];
indcs={10,20,30,40,50,60,70};
if[x_,y_]=Interpolation[Sort@Table[{indcs[[yi]],polys[[yi]]},{yi,Length@polys}]][y];
Plot3D[if[x,y],{x,0,70},{y,-0.1+Min@indcs,0.1+Max@indcs},Mesh->{0,indcs},MeshStyle->Thick,AxesLabel->Automatic,PlotRange->{{0,70},{10,70},{0,80}},BoxRatios->1]


This looks promising indeed (again, if we use the correct PlotRange instead of Mathematica’s default).

However, this result depends on the large amount of data we provided. If we only provide 4 of the 7 functions – still a very clear-cut scenario, the 4 functions evenly spaced –, we’ll get this animal instead:

polys=Through[{f1,f3,f5,f7}[x]];
indcs={10,30,50,70};


It’s very obvious that this result is already quite off, and it would get worse if the provided functions were not evenly spaced on the $y$ axis.

So Interpolation seems to be very demanding when it comes to the amount of data it needs.

Another possibly important disadvantage is that you cannot copy an Interpolation result into your own C program, as you can with Fit results.

But of course, data points from the functions that are provided to Interpolation are 100% exactly reproduced; that’s what interpolation means (blue = interpolation results, green = original discrete parabolas):

g7=Show[Plot[if[y,70],{y,0,70}],Plot[f7[x],{x,0,70},PlotStyle->Green],PlotRange->All,ImageSize->{225,150}];
g6=Show[Plot[if[y,60],{y,0,60}],Plot[f6[x],{x,0,60},PlotStyle->Green],PlotRange->All,ImageSize->{225,150}];
g5=Show[Plot[if[y,50],{y,0,50}],Plot[f5[x],{x,0,50},PlotStyle->Green],PlotRange->All,ImageSize->{225,150}];
g4=Show[Plot[if[y,40],{y,0,40}],Plot[f4[x],{x,0,40},PlotStyle->Green],PlotRange->All,ImageSize->{225,150}];
g3=Show[Plot[if[y,30],{y,0,30}],Plot[f3[x],{x,0,30},PlotStyle->Green],PlotRange->All,ImageSize->{225,150}];
g2=Show[Plot[if[y,20],{y,0,20}],Plot[f2[x],{x,0,20},PlotStyle->Green],PlotRange->All,ImageSize->{225,150}];
g1=Show[Plot[if[y,10],{y,0,10}],Plot[f1[x],{x,0,10},PlotStyle->Green],PlotRange->All,ImageSize->{225,150}];
GraphicsGrid[{{g7,g6,g5},{g4,g3,g2},{g1}}]


Fit

So let’s try Fit:

data={{0,10,0},{5,10,10},{10,10,0},{0,20,0},{10,20,20},{20,20,0},{0,30,0},{15,30,30},{30,30,0},{0,40,0},{20,40,40},{40,40,0},{0,50,0},{25,50,50},{50,50,0},{0,60,0},{30,60,60},{60,60,0},{0,70,0},{35,70,70},{70,70,0}};
ff=Fit[data,{1,x,y,x^2,x y,y^2,x^3,x^2 y,x y^2,y^3},{y,x}]
p=Plot3D[ff,{x,0,70},{y,0,70},PlotRange->{{0,80},{0,80},{0,80}},BoxRatios->1,ImageSize->{400,400}];
GraphicsGrid[{{p,p}}]
Out[3]= -4.09885*10^-14+0.333333 x-0.0135615 x^2+0.000128517 x^3+3.70613*10^-15 y+0.162738 x y-0.00108214 x^2 y-0.162738 y^2+0.000161997 x y^2+0.000920142 y^3


The polynomial looks awfully convoluted, but the 3D surface looks again promising.

However, it’s easy to spot some nonlinearities. So how much does the result (blue) deviate from the original discrete parabolas (green)? Unfortunately, a lot:

g7=Show[Plot[ff/.x->70,{y,0,70}],Plot[f7[x],{x,0,70},PlotStyle->Green],PlotRange->All,ImageSize->{225,150}];
g6=Show[Plot[ff/.x->60,{y,0,60}],Plot[f6[x],{x,0,60},PlotStyle->Green],PlotRange->All,ImageSize->{225,150}];
g5=Show[Plot[ff/.x->50,{y,0,50}],Plot[f5[x],{x,0,50},PlotStyle->Green],PlotRange->All,ImageSize->{225,150}];
g4=Show[Plot[ff/.x->40,{y,0,40}],Plot[f4[x],{x,0,40},PlotStyle->Green],PlotRange->All,ImageSize->{225,150}];
g3=Show[Plot[ff/.x->30,{y,0,30}],Plot[f3[x],{x,0,30},PlotStyle->Green],PlotRange->All,ImageSize->{225,150}];
g2=Show[Plot[ff/.x->20,{y,0,20}],Plot[f2[x],{x,0,20},PlotStyle->Green],PlotRange->All,ImageSize->{225,150}];
g1=Show[Plot[ff/.x->10,{y,0,10}],Plot[f1[x],{x,0,10},PlotStyle->Green],PlotRange->All,ImageSize->{225,150}];
GraphicsGrid[{{g7,g6,g5},{g4,g3,g2},{g1}}]


Since the Fit curves are asymmetric whereas the original parabolas aren’t, I first tried to omit the cubic functions from the list of functions fed to Fit:

ff=Fit[data,{1,x,y,x^2,x y,y^2,x^2 y,x y^2},{y,x}];


This produced better results for some curves, but not for all:

So I figured I should try to improve this result by using much more data points for Fit, derived from my original parabolic functions. I used 11 data points from each parabola:

data=FlattenAt[Table[Sequence @@{{i,j,f[i,j]}},{j,10,70,10},{i,0,j,j/10}],{{1},{2},{3},{4},{5},{6},{7}}];


This results in another improvement, though not by much:

Further improvement is possible by adding more functions to the functions list for Fit:

ff=Fit[data,{1,x,y,x^2,x y,y^2,x^3,x^2 y,x y^2,y^3,x y^3,x^3 y,x^3 y^3,x^2 y^3,x^3 y^2,x^4,y^4,x^4 y,x^4 y^2,x^4 y^3,x^4 y^4,x^3 y^4,x^2 y^4, x y^4},{y,x}]
Out[4]= 5.38885 -0.81743 x+0.0370019 x^2-0.000647046 x^3+3.85345*10^-6 x^4+1.01148 y+0.497555 x y-0.0233103 x^2 y+0.000414295 x^3 y-2.48915*10^-6 x^4 y-0.738136 y^2+0.0259108 x y^2-9.25723*10^-6 x^3 y^2+8.51752*10^-8 x^4 y^2+0.0216425 y^3-0.00138867 x y^3+0.0000271303 x^2 y^3-1.52213*10^-7 x^3 y^3-2.17344*10^-10 x^4 y^3-0.000286201 y^4+0.0000210475 x y^4-5.24626*10^-7 x^2 y^4+5.27399*10^-9 x^3 y^4-1.77042*10^-11 x^4 y^4


However, here is where it stops. Because if we use even more functions for Fit (i.e. above degree 4), the 3D surface begins to degenerate (just when the 7 discrete parabolas are being approximated well):

ff=Fit[data,{1,x,y,x^2,x y,y^2,x^3,x^2 y,x y^2,y^3,x y^3,x^3 y,x^3 y^3,x^2 y^3,x^3 y^2,x^4,y^4,x^4 y,x^4 y^2,x^4 y^3,x^4 y^4,x^3 y^4,x^2 y^4, x y^4,x^5,y^5,x^5 y,x^5 y^2,x^5 y^3,x^5 y^4,x^5 y^5,x^4 y^5,x^3 y^5,x^2 y^5,x y^5},{y,x}]


Using even more data points calculated from the original parabolas doesn’t help, either.

What about using less parabolas, 4 instead of 7, that so clearly damaged the results of Interpolation?

It’s even worse with Fit:

data=FlattenAt[Table[Sequence @@{{i,j,f[i,j]}},{j,10,70,20},{i,0,j,j/10}],{{1},{2},{3},{4}}]


Recap

So here we’re stuck.

With

• an opaque Interpolation algorithm we cannot copy to C code, or else

5.38885 -0.81743 x+0.0370019 x^2-0.000647046 x^3+3.85345*10^-6 x^4+1.01148 y+0.497555 x y-0.0233103 x^2 y+0.000414295 x^3 y-2.48915*10^-6 x^4 y-0.738136 y^2+0.0259108 x y^2-9.25723*10^-6 x^3 y^2+8.51752*10^-8 x^4 y^2+0.0216425 y^3-0.00138867 x y^3+0.0000271303 x^2 y^3-1.52213*10^-7 x^3 y^3-2.17344*10^-10 x^4 y^3-0.000286201 y^4+0.0000210475 x y^4-5.24626*10^-7 x^2 y^4+5.27399*10^-9 x^3 y^4-1.77042*10^-11 x^4 y^4

• some remaining, notable deviations from the original parabolas in the case of Fit,
• and stark deviations for both Interpolation and Fit in the case of less data to start from,

4(x-x^2/y)