# Where do two bicubic-interpolated surfaces intersect a plane?

I have two surfaces, $p_a$ and $p_b$, described by two matrices. I want to find out if both surfaces intersect with the planes $q_a$ and $q_b$ respectively for any pair of x and y coordinates. To do this I set up an equation for the solution within each grid cell, and I traverse both matrices. For bilinear interpolation, Mathematica can find an analytic solution, but for bicubic interpolation I had to solve it numerically. The equation system for bicubic interpolation I want to solve analytically for $x$ and $y$ looks as follows:

$p_i(x,y)=\left(\begin{array}{cccc} 1 & x & x^2 & x^3 \end{array} \right).\left(\begin{array}{cccc} aa_i & ab_i & ac_i & ad_i \\ ba_i & bb_i & bc_i & bd_i \\ ca_i & cb_i & cc_i & cd_i \\ da_i & db_i & dc_i & dd_i \\ \end{array}\right).\left( \begin{array}{c} 1 \\ y \\ y^2 \\ y^3 \\ \end{array} \right)$

$p_a(x,y)=q_a$

$p_b(x,y)=q_b$

All symbols within the matrix are constant within each grid cell and are calculated from the points around the grid cell. I try to solve it in Mathematica by running the following code:

Subscript[p, i_][x_, y_] := {{1, x, x^2, x^3}} .
{{Subscript[aa, i], Subscript[ab, i], Subscript[ac, i],Subscript[ad, i]},
{Subscript[ba, i], Subscript[bb, i], Subscript[bc, i], Subscript[bd, i]},
{Subscript[ca, i], Subscript[cb, i], Subscript[cc, i], Subscript[cd, i]},
{Subscript[da, i], Subscript[db, i], Subscript[dc, i], Subscript[dd, i]}} .
{{1}, {y}, {y^2}, {y^3}}
Solve[Subscript[p, a][x, y] == Subscript[q, a] &&
Subscript[p, b][x, y] == Subscript[q, b],
{x, y}]


This has been running for hours without any result. Is it possible to solve this analytically? Do I have the wrong approach?

• Do you only want to know if the surface intersects the z=const plane or also where? In general i would expect this to be a curve, so which intersection point would you want to compute then, an arbitrary one or all of them? Also why two surfaces? Does it make a difference to your problem instead of solving the two surface and planes separately? Could you give a bit more context to the problem by telling where you want to apply it? What do you know about the input data (the matrix coefficients and qa/qb)? – Thies Heidecke Dec 15 '17 at 16:35
• The intersection between the planes and the surfaces are curves, but the intersections between the curves are points. I need to find those points. I know that the intersection between the planes and the surface are close to a line. I know that the input data normally only gives one intersection point between the curves. Why two surfaces? With only one surface, the intersection between the plane and the surface would give a curve and not a point, and I need a point. – Johan Dec 18 '17 at 7:50
• Even though an analytic solution is wanted, an explicit numeric example would be useful just to see what exactly is the problem under consideration. – Daniel Lichtblau Jan 14 '18 at 15:41

I doubt that this is solvable analytically. You can try it with FindRoot, instead:

SeedRandom[123];
a0 = RandomReal[{-1, 1}, {2}];
A = RandomReal[{-1, 1}, {2, 4, 4}];
Quiet[Block[{x, xx},
xx = {x[[1]], x[[2]]};
f = x \[Function] (A.{1., x[[1]], x[[1]]^2, x[[1]]^3}).{1., x[[2]], x[[2]]^2, x[[2]]^3} - a0;
Df = x \[Function] Evaluate[D[f[xx], {xx, 1}]];
]];
x0 = {0., 0.};
x = FindRoot[f, {x0}, Jacobian -> Df][[1]]
f[x]


{-1.76746, -0.310105}

{4.44089*10^-16, 6.66134*10^-16}

Of course, success is not guaranteed and the outcome of the method depends on the starting point x0.

• I have to implement the solution in C++, and I cannot solve it numerically in Mathematica. I wanted to use Mathematica to find an analytic solution that I can implement in C++. – Johan Dec 15 '17 at 11:49
• I see. But you can also implement Newton's method for this specific problem quite easily in C++. – Henrik Schumacher Dec 15 '17 at 11:58
• I have already implemented a numeric method in C++, but it is too slow for the embedded target. It needs to run at 3 kHz. The bilinear analytic solution is running at 8 kHz, whereas the bicubic numeric solution is running at 32 Hz. No doubt I can optimize it a lot more, but I'm not sure I could reach 3 kHz. – Johan Dec 15 '17 at 12:05
• If you mean that you want to solve at least 3000 of such 2 x 2 root finding problems per second: A basal Newton implementation in C for these kind of problems runs with a rate of about 2 MHz on my machine... – Henrik Schumacher Dec 15 '17 at 12:39
• And I mean MHz with a big M... – Henrik Schumacher Dec 15 '17 at 12:47