# Finite difference for Lighthill source term

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: • Welcome to SciComp.SE! Could you add more detail (ideally, with some pictures) about what kind of data you have, what you expect to see, and what you are seeing ("wild unconnectivity" is a bit vague...) – Christian Clason Sep 30 '15 at 10:03
• OK, I will try to elaborate something 'visible' (quiver plot 100x200 is not readible in this site). I haven't post anything specific cause I expect a mistake somewhere in process. – Victor Pira Sep 30 '15 at 10:07
• (If it's just applying the product rule, then it's wrong -- you should always get products of derivatives, not derivatives of products, otherwise there'd be no point...) – Christian Clason Sep 30 '15 at 10:44
• That might be the mistake. Could you please be a bit more specific? – Victor Pira Sep 30 '15 at 10:46

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).
$$\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.)