Finite difference for Lighthill source term

Question

I have experimental data containing horizontal and vertical components of speed and I need to evaluate this:

$$ \frac{\partial^2}{\partial x_i \partial x_j}\left(v_iv_j\right) $$

Denoting $U$ as horizontal component and $V$ as vertical, we should obtain:

$$ \frac{\partial^2}{\partial x_i \partial x_j}\left(v_iv_j\right) \rightarrow \frac{\partial^2}{\partial x^2}(U^2)+\frac{\partial^2}{\partial y^2}(V^2) + 2\frac{\partial^2}{\partial x \partial y}(UV) $$

I have used central finite difference of the second order for $\partial^2/\partial x^2$ and finite difference of the first order for mixed derivative. In MATLAB:

%H and WW are height and width respectively
for ii=3:H-2
    for jj=3:WW-2

        %d^2/dx^2
        A(ii,jj)=(-1/12*U(ii,jj-2)^2+4/3*U(ii,jj-1)^2-5/2*U(ii,jj)^2+4/3*U(ii,jj+1)^2-1/12*U(ii,jj+2)^2)/h^2;

        %d^2/dy^2
        B(ii,jj)=(-1/12*V(ii-2,jj)^2+4/3*V(ii-1,jj)^2-5/2*V(ii,jj)^2+4/3*V(ii+1,jj)^2-1/12*V(ii+2,jj)^2)/h^2;

        %Mixed derivative
        C(ii,jj)=(V(ii+1,jj+1)*U(ii+1,jj+1)-V(ii+1,jj-1)*U(ii+1,jj-1)-V(ii-1,jj+1)*U(ii-1,jj+1)+V(ii-1,jj-1)*U(ii-1,jj-1))/(4*h^2);

        %sum of derivatives
        LGHT(ii,jj)=(A(ii,jj)+B(ii,jj)+2*C(ii,jj));

   end
end

Results are not very satisfactory so far: "wild unconnectivity" all around the screen.

• Do I have a mistake somewhere?
• Would you use different kind of scheme?
• What kind of data treatment (spline?) would you recommend for obtaining smoother results?

Example of input data:

Output:

Answer 1

Your reformulation is wrong: The product rule for derivatives is $(uv)' = u'v + uv'$. Applying this in your situation (with $\partial_x$ instead of $\frac{\partial}{\partial x}$ and $\partial_{xy} = \frac{\partial^2}{\partial x\partial y}$ for brevity would yield $$\partial_{xy}(uv) = \partial_y(\partial_x u \,v + u\,\partial_x v) = \partial_{xy}u \,v + \partial_x u \, \partial_y v + \partial_y u\,\partial_x v + u\,\partial_{xy}v $$ (note that you have four different terms, not three).

If you want to have a second-order accurate scheme (mixing schemes of different accuracy is in general a waste of effort), the standard stencils are

$$\partial_x u_{ij} \approx \frac{u_{i+1,j}-u_{i-1,j}}{2h},\qquad\partial_y u_{ij} \approx \frac{u_{i,j+1}-u_{i,j-1}}{2h}$$ and $$\partial_{xy} u_{ij} \approx \frac{u_{i+1,j+1}-u_{i+1,j-1}-u_{i-1,j+1}+u_{i-1,j-1}}{4h^2}.$$ You can either apply the second directly to $(u\,v)_{ij}$, or plug in both stencils on the right-hand side. (Or better yet, both, and compare the results.)