0
$\begingroup$

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?

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ 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)? $\endgroup$ Commented Dec 15, 2017 at 16:35
  • $\begingroup$ 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. $\endgroup$ Commented Dec 18, 2017 at 7:50
  • $\begingroup$ Even though an analytic solution is wanted, an explicit numeric example would be useful just to see what exactly is the problem under consideration. $\endgroup$ Commented Jan 14, 2018 at 15:41

1 Answer 1

Reset to default
1
$\begingroup$

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.

Improve this answer
$\endgroup$
8
  • $\begingroup$ 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++. $\endgroup$ Commented Dec 15, 2017 at 11:49
  • $\begingroup$ I see. But you can also implement Newton's method for this specific problem quite easily in C++. $\endgroup$ Commented Dec 15, 2017 at 11:58
  • $\begingroup$ 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. $\endgroup$ Commented Dec 15, 2017 at 12:05
  • $\begingroup$ 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... $\endgroup$ Commented Dec 15, 2017 at 12:39
  • $\begingroup$ And I mean MHz with a big M... $\endgroup$ Commented Dec 15, 2017 at 12:47

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.