Bivariate Hermite interpolation containing second order derivatives

Question

I need to interpolate a function (say $\Phi(x,y)$) with four known points. I have $\Phi(x_i,y_i)$, $\Phi_x(x_i,y_i)$, $\Phi_y(x_i,y_i)$ $\Phi_{xx} (x_i,y_i)$, $\Phi_{yy} (x_i,y_i)$ and $\Phi_{xy}(x_i,y_i)$ at each four points. I cannot omit any of these twenty four(6 x 4) values. How do I interpolate my function efficiently? Please give some reference.

Answer 1

So we need to find such polynomial $P(x, y)$ containing 24 unknown coefficients to match 24 values of $\Phi, \Phi_x, \Phi_y, \Phi_{xy}, \Phi_{xx}, \Phi_{yy}$.

Let's start with a subproblem — conventional Hermite interpolation: given $\Phi, \Phi_x, \Phi_y, \Phi_{xy}$ find such polynomial $P_H(x,y)$ with $\deg_x P_H(x,y) = \deg_y P_H(x,y) = 3$. This is well-known bicubic interpolation problem.

Now we have extra $24 - 16 = 8$ conditions given by $\Phi_{xx}$ and $\Phi_{yy}$ at every node. To match the second derivatives we need to raise the degree of the polynomial from third to fifth. But buquintic interpolation has $6^2 = 36$ unknown coefficients. We need to eliminate some ($36 - 24 = 12$) coefficients without loosing the interpolation property.

Applying trial and error method I've obtained the following set of $24$ basis monomials for the interpolation: $$ 16 \text{ monomials } x^i y^j, \; i,j \in 0,\dots,3\; \text{ (basis for Hermite's bicubic)}\\ 4 \text{ monomials } x^i y^j, \; i \in 4,5, \; j \in 0,1 \text{ (basis for second derivatives w.r.t $x$)}\\ 4 \text{ monomials } x^i y^j, \; i \in 0,1, \; j \in 4,5 \text{ (basis for second derivatives w.r.t $y$)} $$

So the interpolant is sought in form $$ P(x, y) = \sum_{i=0}^3 \sum_{j=0}^3 c_{i,j} x^i y^j + \sum_{i=4}^5 \sum_{j=0}^1 c_{i,j} x^i y^j + \sum_{i=0}^1 \sum_{j=4}^5 c_{i,j} x^i y^j $$ The $24$ coefficients $c_{i,j}$ are then determined from the linear system $$ P(x_k, y_k) = \Phi(x_k, y_k),\; k = 1,\dots,4\\ P_x(x_k, y_k) = \Phi_x(x_k, y_k),\; k = 1,\dots,4\\ P_y(x_k, y_k) = \Phi_y(x_k, y_k),\; k = 1,\dots,4\\ P_{xy}(x_k, y_k) = \Phi_{xy}(x_k, y_k),\; k = 1,\dots,4\\ P_{xx}(x_k, y_k) = \Phi_{xx}(x_k, y_k),\; k = 1,\dots,4\\ P_{yy}(x_k, y_k) = \Phi_{yy}(x_k, y_k),\; k = 1,\dots,4 $$

Direct check using Wolfram Mathematica shows that the linear system for the $24$ coefficients $c_{i,j}$ is solvable.

l=Flatten@{
  Table[c[i,j] x^i y^j,{i,0,3},{j,0,3}],
  Table[c[i,j] x^i y^j,{i,4,5},{j,0,1}],
  Table[c[i,j] x^i y^j,{i,0,1},{j,4,5}]
};
u[x_,y_]=Total@l;
cs=(l/.{x->1,y->1});
A = Last@CoefficientArrays[Flatten@{
  Table[u[x, y], {x, 0, 1}, {y, 0, 1}],
  Table[Derivative[1, 0][u][x, y], {x, 0, 1}, {y, 0, 1}],
  Table[Derivative[0, 1][u][x, y], {x, 0, 1}, {y, 0, 1}],
  Table[Derivative[1, 1][u][x, y], {x, 0, 1}, {y, 0, 1}],
  Table[Derivative[2, 0][u][x, y], {x, 0, 1}, {y, 0, 1}], 
  Table[Derivative[0, 2][u][x, y], {x, 0, 1}, {y, 0, 1}]
}, cs];
MatrixRank[A]
(* Out = 24 *)

The complete answer (quite lengthy) can be obtained from solving the system (here's the WM code to produce the answer)

Thread[(cs /. c[i_, j_] -> Subscript[c, i, j]) == 
   Simplify@LinearSolve[A, Flatten@{
       Table[F[x, y], {x, 0, 1}, {y, 0, 1}],
       Table[Fx[x, y], {x, 0, 1}, {y, 0, 1}],
       Table[Fy[x, y], {x, 0, 1}, {y, 0, 1}], 
       Table[Fxy[x, y], {x, 0, 1}, {y, 0, 1}], 
       Table[Fxx[x, y], {x, 0, 1}, {y, 0, 1}], 
       Table[Fyy[x, y], {x, 0, 1}, {y, 0, 1}]}]] // Column

Below is an example for $\Phi(x,y) = \sin\left(x + \frac{e^y}{x+1}\right)$. Yellow if $\Phi(x,y)$ and blue is $P(x,y)$.

Here's the difference

Appendix. The Lagrange basis for this problem is $$ \ell _{00}=(x-1) (y-1) \left(\left(4 x^2-2 x+1\right) y^2+x \left(6 \ x^3-9 x^2+x+1\right)+\\+(1-2 x) x y+6 y^4-9 y^3+y+1\right) \\ \ell _{01}=-(x-1) y \left(6 x^4-9 x^3+x^2 \left(4 y^2-6 y+3\right)+x (3-2 y) y+\\+y^2 (3 y (2 y-5)+10)\right) \\ \ell _{10}=-x (y-1) \left(\left(4 x^2-6 x+3\right) y^2+x^2 (3 x (2 \ x-5)+10)+\\+(3-2 x) x y+6 y^4-9 y^3\right) \\ \ell _{11}=x y \left(6 x^4-15 x^3+2 x^2 (y (2 y-3)+6)+x \left(-6 \ y^2+9 y-3\right)+\\+3 (y-1)^2 y (2 y-1)\right) \\ \ell _ {00}^{{x}}=(x-1)^2 x (y-1) (x (3 x-2)+(y-1) (2 y+1)) \\ \ell _ {01}^{{x}}=-(x-1)^2 x y (x (3 x-2)+y (2 y-3)) \\ \ell _ {10}^{{x}}=(x-1) x^2 (y-1) (x (3 x-4)+y (2 y-1)) \\ \ell _ {11}^{{x}}=-(x-1) x^2 y \left(3 x^2-4 x+2 y^2-3 y+1\right) \ \\ \ell _ {00}^{{y}}=(x-1) (y-1)^2 y (x (2 x-1)+(y-1) (3 y+1)) \\ \ell _ {01}^{{y}}=(x-1) (y-1) y^2 (x (2 x-1)+y (3 y-4)) \\ \ell _ {10}^{{y}}=-x (y-1)^2 y (x (2 x-3)+y (3 y-2)) \\ \ell _ {11}^{{y}}=-x (y-1) y^2 \left(2 x^2-3 x+3 y^2-4 y+1\right) \ \\ \ell _ {00}^{{xy}}=(x-1)^2 x (y-1)^2 y \qquad \ell _ {01}^{{xy}}=(x-1)^2 x (y-1) y^2 \\ \ell _ {10}^{{xy}}=(x-1) x^2 (y-1)^2 y \qquad \ell _ {11}^{{xy}}=(x-1) x^2 (y-1) y^2 \\ \ell _ {00}^{{xx}}=\frac{1}{2} (x-1)^3 x^2 (y-1) \qquad \ell _ {01}^{{xx}}=-\frac{1}{2} (x-1)^3 x^2 y \\ \ell _ {10}^{{xx}}=-\frac{1}{2} (x-1)^2 x^3 (y-1) \qquad \ell _ {11}^{{xx}}=\frac{1}{2} (x-1)^2 x^3 y \\ \ell _ {00}^{{yy}}=\frac{1}{2} (x-1) (y-1)^3 y^2 \qquad \ell _ {01}^{{yy}}=-\frac{1}{2} (x-1) (y-1)^2 y^3 \\ \ell _ {10}^{{yy}}=-\frac{1}{2} x (y-1)^3 y^2 \qquad \ell _ {11}^{{yy}}=\frac{1}{2} x (y-1)^2 y^3 $$